Commit 5c70d545 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

dune-functions globalBasis is now basis of problemStat

parent f9a757d2
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include <dune/amdis/LinearAlgebra.hpp> #include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/common/Mpl.hpp> #include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp> #include <dune/amdis/common/TypeDefs.hpp>
...@@ -18,47 +19,37 @@ namespace AMDiS ...@@ -18,47 +19,37 @@ namespace AMDiS
/// The grid view the global FE basis lives on /// The grid view the global FE basis lives on
using GridView = typename GlobalBasis::GridView; using GridView = typename GlobalBasis::GridView;
using SystemMatrixType = SystemMatrix<GlobalBasis>;
using SystemVectorType = SystemVector<GlobalBasis>;
using TypeTree = typename GlobalBasis::LocalView::Tree;
using
public: public:
/// Constructor, stores a shared-pointer to the feSpaces /// Constructor, stores a shared-pointer to the feSpaces
Assembler(std::shared_ptr<GlobalBasis> const& globalBasis, bool asmMatrix, bool asmVector) Assembler(GlobalBasis& globalBasis,
MatrixOperators<GlobalBasis>& matrixOperators,
VectorOperators<GlobalBasis>& rhsOperators)
: globalBasis_(globalBasis) : globalBasis_(globalBasis)
, asmMatrix_(asmMatrix) , matrixOperators_(matrixOperators)
, asmVector_(asmVector) , rhsOperators_(rhsOperators)
{} {}
void update(GridView const& gv)
{
globalBasis_.update(gv);
}
/// Assemble the linear system /// Assemble the linear system
template <class MatrixOperators, class VectorOperators> template <class SystemMatrixType, class SystemVectorType>
void assemble( void assemble(
GridView const& gv,
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
MatrixOperators& matrixOperators, bool asmMatrix, bool asmVector);
VectorOperators& rhsOperators);
private: private:
/// Sets the system to zero and initializes all operators and boundary conditions /// Sets the system to zero and initializes all operators and boundary conditions
template <class MatrixOperators, class VectorOperators> template <class SystemMatrixType, class SystemVectorType>
void initMatrixVector( void initMatrixVector(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
MatrixOperators& matrixOperators, bool asmMatrix, bool asmVector) const;
VectorOperators& rhsOperators) const;
template <class MatrixOperators, class VectorOperators>
void assembleElement(
SystemMatrixType& matrix,
SystemVectorType& rhs,
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const;
template <class ElementContainer, class Container, class Operators, class... Bases> template <class ElementContainer, class Container, class Operators, class... Bases>
...@@ -71,26 +62,32 @@ namespace AMDiS ...@@ -71,26 +62,32 @@ namespace AMDiS
/// Finish insertion into the matrix and assembles boundary conditions /// Finish insertion into the matrix and assembles boundary conditions
/// Return the number of nonzeros assembled into the matrix /// Return the number of nonzeros assembled into the matrix
template <class MatrixOperators, class VectorOperators> template <class SystemMatrixType, class SystemVectorType>
std::size_t finishMatrixVector( std::size_t finishMatrixVector(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
MatrixOperators& matrixOperators, bool asmMatrix, bool asmVector) const;
VectorOperators& rhsOperators) const;
/// Return whether the matrix-block needs to be assembled
template <class LocalView0, class... LovalViews>
auto const& getElement(LocalView0 const& localView, LovalViews const&...) const
{
return localView.element();
}
/// Return whether the matrix-block needs to be assembled /// Return whether the matrix-block needs to be assembled
template <class Basis0, class... Bases> template <class LocalView0, class... LovalViews>
auto getElement(Basis0 const& basis0, Bases const&... bases) const auto const& getGridView(LocalView0 const& localView, LovalViews const&...) const
{ {
return basis0.localView().element(); return globalBasis_.gridView();
} }
private: private:
std::shared_ptr<GlobalBasis> globalBasis_; GlobalBasis& globalBasis_;
bool asmMatrix_; MatrixOperators<GlobalBasis>& matrixOperators_;
bool asmVector_; VectorOperators<GlobalBasis>& rhsOperators_;
}; };
} // end namespace AMDiS } // end namespace AMDiS
......
#pragma once #pragma once
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/typetree/treepath.hh>
#include <dune/amdis/utility/Visitor.hpp>
namespace AMDiS { namespace AMDiS {
template <class GlobalBasis> template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators> template <class SystemMatrixType, class SystemVectorType>
void Assembler<GlobalBasis>::assemble( void Assembler<GlobalBasis>::assemble(
GridView const& gv,
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
MatrixOperators& matrixOperators, bool asmMatrix, bool asmVector)
VectorOperators& rhsOperators)
{ {
// 1. Update global bases // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
globalBasis_.update(gv); initMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
auto localView = globalBasis_.localView();
auto localIndexSet = globalBasis_.localIndexSet();
// 2. init matrix and rhs vector and initialize dirichlet boundary conditions // 2. create a local matrix and vector
initMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators); std::size_t localSize = localView.maxSize();
mtl::dense2D<double> elementMatrix(localSize, localSize);
set_to_zero(elementMatrix);
mtl::dense_vector<double> elementVector(localSize);
set_to_zero(elementVector);
// 3. traverse grid and assemble operators on the elements // 3. traverse grid and assemble operators on the elements
for (auto const& element : elements(globalBasis_.gridView())) for (auto const& element : elements(globalBasis_.gridView()))
{ {
globalBasis_.bind(element); localView.bind(element);
assembleElement(matrix, rhs, matrixOperators, rhsOperators);
globalBasis_.unbind(); // traverse type-tree of global-basis
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto rowLocalView = rowBasis.localView();
rowLocalView.bind(element); // NOTE: Is this necessary
auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.assemble(asmVector) && !rhsOp.empty())
assembleElementOperators(elementVector, rhs, rhsOp, rowLocalView);
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto& matOp = matrixOperators_[rowNode][colNode];
if (matOp.assemble(asmMatrix) && !matOp.empty()) {
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
auto colLocalView = colBasis.localView();
colLocalView.bind(element); // NOTE: Is this necessary
assembleElementOperators(elementMatrix, matrix, matOp, rowLocalView, colLocalView);
}
});
});
localIndexSet.bind(localView);
// add element-matrix to system-matrix
for (std::size_t i = 0; i < localView.size(); ++i) {
auto const row = localIndexSet.index(i);
for (std::size_t j = 0; j < localView.size(); ++j) {
auto const col = localIndexSet.index(j);
matrix(row,col) += elementMatrix(i,j);
}
}
// add element-vector to system-vector
for (std::size_t i = 0; i < localView.size(); ++i) {
auto const idx = localIndexSet.index(i);
rhs[idx] += elementVector[i];
}
localIndexSet.unbind();
localView.unbind();
} }
// 4. finish matrix insertion and apply dirichlet boundary conditions // 4. finish matrix insertion and apply dirichlet boundary conditions
std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators); std::size_t nnz = finishMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
msg("fillin of assembled matrix: ", nnz); msg("fillin of assembled matrix: ", nnz);
} }
template <class GlobalBasis> template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators> template <class SystemMatrixType, class SystemVectorType>
void Assembler<GlobalBasis>::initMatrixVector( void Assembler<GlobalBasis>::initMatrixVector(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
MatrixOperators& matrixOperators, bool asmMatrix, bool asmVector) const
VectorOperators& rhsOperators) const
{ {
matrix.init(globalBasis_); matrix.init(asmMatrix);
solution.init(globalBasis_); solution.compress();
rhs.init(globalBasis_); rhs.compress();
auto localView = globalBasis_.localView(); auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath) forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{ {
if (rowNode.isLeaf) if (rowNode.isLeaf)
msg(globalBasis_.size(rowTreePath), " DOFs for Basis[", Dune::TypeTree::treePathIndex(rowTreePath), "]"); msg(0, " DOFs for Basis[", 1, "]"); // TODO: add right values
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath); auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
if (rhsOperators[rowNode].assemble(asmVector_)) if (rhsOperators_[rowNode].assemble(asmVector))
rhsOperators[rowNode].init(rowBasis); rhsOperators_[rowNode].init(rowBasis);
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath) forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{ {
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath); auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
if (matrixOperators[rowNode][colNode].assemble(asmMatrix_)) if (matrixOperators_[rowNode][colNode].assemble(asmMatrix))
matrixOperators[rowNode][colNode].init(rowBasis, colBasis); matrixOperators_[rowNode][colNode].init(rowBasis, colBasis);
// init boundary condition // init boundary condition
// for (int c = 0; c < nComponents; ++c) // for (int c = 0; c < nComponents; ++c)
...@@ -72,61 +124,21 @@ void Assembler<GlobalBasis>::initMatrixVector( ...@@ -72,61 +124,21 @@ void Assembler<GlobalBasis>::initMatrixVector(
template <class GlobalBasis> template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators> template <class ElementContainer, class Container, class Operators, class... LocalViews>
void Assembler<GlobalBasis>::assembleElement(
SystemMatrixType& matrix,
SystemVectorType& rhs,
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const
{
auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto& rhsOp = rhsOperators[rowNode];
if (rhsOp.assemble(asmVector_) && !rhsOp.empty()) {
auto elementVector = makeElementVector(rowNode);
set_to_zero(elementVector);
assembleElementOperators(elementVector, rhs, rhsOp, rowBasis);
}
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto& matOp = matrixOperators[rowNode][colNode];
if (matOp.assemble(asmMatrix_) && !matOp.empty()) {
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
auto elementMatrix = makeElementMatrix(rowNode, colNode);
set_to_zero(elementMatrix);
assembleElementOperators(elementMatrix, matrix, matOp, rowBasis, colBasis);
}
});
});
}
template <class GlobalBasis>
template <class ElementContainer, class Container, class Operators, class... Bases>
void Assembler<GlobalBasis>::assembleElementOperators( void Assembler<GlobalBasis>::assembleElementOperators(
ElementContainer& elementContainer, ElementContainer& elementContainer,
Container& container, Container& container,
Operators& operators, Operators& operators,
Bases const&... subBases) const LocalViews const&... localViews) const
{ {
if (operators.empty()) auto const& element = getElement(localViews...);
return false; auto const& gridView = getGridView(localViews...);
auto const& element = getElement(subBases...);
auto const& gridView = subBasis.gridView();
bool add = false; bool add = false;
auto assemble_operators = [&](auto const& context, auto& operator_list) { auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) { for (auto scaled : operator_list) {
bool add_op = scaled.op->assemble(gridView, context, subBases.localView()..., elementContainer, scaled.factor); bool add_op = scaled.op->assemble(context, localViews..., elementContainer, scaled.factor);
add = add || add_op; add = add || add_op;
} }
}; };
...@@ -145,36 +157,30 @@ void Assembler<GlobalBasis>::assembleElementOperators( ...@@ -145,36 +157,30 @@ void Assembler<GlobalBasis>::assembleElementOperators(
assemble_operators(intersection, operators.intersection); assemble_operators(intersection, operators.intersection);
} }
} }
if (!add)
return;
elementContainer.apply(subBases.localIndexSet()..., container);
} }
template <class GlobalBasis> template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators> template <class SystemMatrixType, class SystemVectorType>
std::size_t Assembler<GlobalBasis>::finishMatrixVector( std::size_t Assembler<GlobalBasis>::finishMatrixVector(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
MatrixOperators& matrixOperators, bool asmMatrix, bool asmVector) const
VectorOperators& rhsOperators) const
{ {
matrix.finish(); matrix.finish();
auto localView = globalBasis_.localView(); auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath) forEachNode(localView.tree(), [&,this](auto const& rowNode, auto /*rowTreePath*/)
{ {
auto& rhsOp = rhsOperators[rowNode]; auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.assemble(asmVector_)) if (rhsOp.assemble(asmVector))
rhsOp.assembled = true; rhsOp.assembled = true;
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath) forEachNode(localView.tree(), [&,this](auto const& colNode, auto /*colTreePath*/)
{ {
auto& matOp = matrixOperators[rowNode][colNode]; auto& matOp = matrixOperators_[rowNode][colNode];
if (matOp.assemble(asmMatrix_)) if (matOp.assemble(asmMatrix))
matOp.assembled = true; matOp.assembled = true;
// finish boundary condition // finish boundary condition
......
...@@ -11,8 +11,8 @@ dune_add_library("duneamdis" NO_EXPORT ...@@ -11,8 +11,8 @@ dune_add_library("duneamdis" NO_EXPORT
ProblemInstat.cpp ProblemInstat.cpp
ProblemStat.cpp ProblemStat.cpp
StandardProblemIteration.cpp StandardProblemIteration.cpp
linear_algebra/istl/SystemMatrix.cpp #linear_algebra/istl/SystemMatrix.cpp
linear_algebra/istl/SystemVector.cpp #linear_algebra/istl/SystemVector.cpp
linear_algebra/mtl/SystemMatrix.cpp linear_algebra/mtl/SystemMatrix.cpp
linear_algebra/mtl/SystemVector.cpp linear_algebra/mtl/SystemVector.cpp
utility/Filesystem.cpp utility/Filesystem.cpp
......
...@@ -21,17 +21,16 @@ namespace AMDiS ...@@ -21,17 +21,16 @@ namespace AMDiS
{ {
private: // typedefs and static constants private: // typedefs and static constants
using Mesh = typename Traits::Mesh; using MeshView = typename Traits::GridView;
using MeshView = typename Mesh::LeafGridView;
/// Dimension of the mesh /// Dimension of the mesh
static constexpr int dim = Traits::dim; static constexpr int dim = MeshView::dimension;
/// Dimension of the world /// Dimension of the world
static constexpr int dow = Traits::dow; static constexpr int dow = MeshView::dimensionworld;
/// Number of problem components /// Number of problem components
static constexpr int nComponents = Traits::nComponents; static constexpr int nComponents = 1;
public: public:
...@@ -62,7 +61,7 @@ namespace AMDiS ...@@ -62,7 +61,7 @@ namespace AMDiS
// copy dofvector to vertex data // copy dofvector to vertex data
forEach(range_<0, nComponents>, [this, &solutions](const auto _i) forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
{ {
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors)); this->dofVector2vertexVector(solutions, std::get<_i>(data_vectors));
vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]); vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
}); });
vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/); vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
...@@ -77,7 +76,7 @@ namespace AMDiS ...@@ -77,7 +76,7 @@ namespace AMDiS
// copy dofvector to vertex data // copy dofvector to vertex data
forEach(range_<0, nComponents>, [this, &solutions](const auto _i) forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
{ {
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors)); this->dofVector2vertexVector(solutions, std::get<_i>(data_vectors));
vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]); vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
}); });
vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/); vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/);
...@@ -99,7 +98,7 @@ namespace AMDiS ...@@ -99,7 +98,7 @@ namespace AMDiS
auto localView = feSpace.localView(); auto localView = feSpace.localView();
auto localIndexSet = feSpace.localIndexSet(); auto localIndexSet = feSpace.localIndexSet();
#if 0
// copy data to P1-vector // copy data to P1-vector
for (auto const& element : elements(meshView)) { for (auto const& element : elements(meshView)) {
localView.bind(element); localView.bind(element);
...@@ -125,6 +124,7 @@ namespace AMDiS ...@@ -125,6 +124,7 @@ namespace AMDiS
} }
} }
} }
#endif
} }
private: private:
......
...@@ -2,6 +2,8 @@ ...@@ -2,6 +2,8 @@
#include <vector> #include <vector>
#include <dune/typetree/nodetags.hh>
#include <dune/amdis/LocalAssembler.hpp> #include <dune/amdis/LocalAssembler.hpp>
#include <dune/amdis/common/Mpl.hpp> #include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp> #include <dune/amdis/common/TypeDefs.hpp>
...@@ -35,15 +37,24 @@ namespace AMDiS ...@@ -35,15 +37,24 @@ namespace AMDiS
// tag dispatching for FirstOrderType... // tag dispatching for FirstOrderType...
template <class Operator, class RowView, class ColView> template <class Operator, class RowView, class ColView>
void calculateElementMatrix(Operator& op, RowView const& rowView, ColView const& colView, ElementMatrix& elementMatrix) void calculateElementMatrix(Operator& op,
RowView const& rowView, ColView const& colView,
ElementMatrix& elementMatrix, double fac)
{ {
calculateElementMatrix(op, rowView, colView, elementMatrix, FirstOrderType_<type>); using RowNode = typename RowView::Tree;
using ColNode = typename ColView::Tree;
calculateElementMatrix(op, rowView, colView, elementMatrix, fac,
typename RowNode::NodeTag{}, typename ColNode::NodeTag{}, FirstOrderType_<type>);
} }
template <class Operator, class RowView> template <class Operator, class RowView>
void calculateElementVector(Operator& op, RowView const& rowView, ElementVector& elementVector) void calculateElementVector(Operator& op,
RowView const& rowView,
ElementVector& elementVector, double fac)
{ {
calculateElementVector(op, rowView, elementVector, FirstOrderType_<type>); using RowNode = typename RowView::Tree;
calculateElementVector(op, rowView, elementVector, fac,
typename RowNode::NodeTag{}, FirstOrderType_<type>);
} }