Skip to content
GitLab
Commit f9a757d2 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

started to rewrite Assembler using dune-functions globalBasis

parent beb6bd46
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ f9a757d2
......@@ -7,6 +7,8 @@ if(NOT (dune-common_DIR OR dune-common_ROOT OR
${PROJECT_BINARY_DIR})
endif()
set(CXX_MAX_STANDARD 14)
#find dune-common and set the module path
find_package(dune-common REQUIRED)
list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
......
dune/amdis/Assembler.hpp
View file @ f9a757d2
......@@ -6,139 +6,89 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/amdis/FiniteElementSpaces.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
namespace AMDiS
{
template <class FeSpaces>
template <class GlobalBasis>
class Assembler
{
template <std::size_t I>
using FeSpace = std::tuple_element_t<I, FeSpaces>;
/// Number of problem components
static constexpr int nComponents = std::tuple_size<FeSpaces>::value;
/// The grid view the global FE basis lives on
using GridView = typename FeSpace<0>::GridView;
using SystemMatrixType = SystemMatrix<FeSpaces>;
using SystemVectorType = SystemVector<FeSpaces>;
template <class T>
using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
using GridView = typename GlobalBasis::GridView;
template <class T>
using VectorEntries = Dune::FieldVector<T, nComponents>;
using SystemMatrixType = SystemMatrix<GlobalBasis>;
using SystemVectorType = SystemVector<GlobalBasis>;
using ElementVector = Impl::ElementVector;
using ElementMatrix = Impl::ElementMatrix;
using TypeTree = typename GlobalBasis::LocalView::Tree;
using
public:
/// Constructor, stores a shared-pointer to the feSpaces
Assembler(std::shared_ptr<FeSpaces> const& feSpaces, bool asmMatrix, bool asmVector)
: globalBases_(feSpaces)
Assembler(std::shared_ptr<GlobalBasis> const& globalBasis, bool asmMatrix, bool asmVector)
: globalBasis_(globalBasis)
, asmMatrix_(asmMatrix)
, asmVector_(asmVector)
{}
/// Assemble the linear system
template <class Operators>
template <class MatrixOperators, class VectorOperators>
void assemble(
GridView const& gv,
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators);
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators);
private:
/// Sets the system to zero and initializes all operators and boundary conditions
template <class Operators>
template <class MatrixOperators, class VectorOperators>
void initMatrixVector(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators) const;
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const;
/// Assemble a block-matrix of element-matrices and return a matrix of flags, whether
/// a given block has received some entries.
template <class Operators>
MatrixEntries<bool> assembleElementMatrices(
MatrixEntries<Operators>& operators,
MatrixEntries<ElementMatrix>& elementMatrix) const;
/// Assemble a block-vector of element-vectors and return a vector of flags, whether
/// a given block has received some entries.
template <class Operators>
VectorEntries<bool> assembleElementVectors(
VectorEntries<Operators>& operators,
VectorEntries<ElementVector>& elementVector) const;
template <class MatrixOperators, class VectorOperators>
void assembleElement(
SystemMatrixType& matrix,
SystemVectorType& rhs,
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const;
/// Assemble one block of the block-element-matrix
// The MatrixData argument stores all matrix-operators
template <class Operators, class RowView, class ColView>
bool assembleElementMatrix(
template <class ElementContainer, class Container, class Operators, class... Bases>
void assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
ElementMatrix& elementMatrix,
RowView const& rowLocalView,
ColView const& colLocalView) const;
/// Assemble one block of the block-element-vector
// The VectorData argument stores all vector-operators
template <class Operators, class RowView>
bool assembleElementVector(
Operators& operators,
ElementVector& elementVector,
RowView const& rowLocalView) const;
/// Add the block-element-matrix to the system-matrix
void addElementMatrices(
SystemMatrixType& dofmatrix,
MatrixEntries<bool> const& addMat,
MatrixEntries<ElementMatrix> const& elementMatrix) const;
Bases const&... subBases) const;
/// Add the block-element-vector to the system-vector
void addElementVectors(
SystemVectorType& dofvector,
VectorEntries<bool> const& addVec,
VectorEntries<ElementVector> const& elementVector) const;
/// Finish insertion into the matrix and assembles boundary conditions
/// Return the number of nonzeros assembled into the matrix
template <class Operators>
template <class MatrixOperators, class VectorOperators>
std::size_t finishMatrixVector(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators) const;
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const;
/// Return whether the matrix-block needs to be assembled
template <std::size_t R, std::size_t C, class Operators>
bool assembleMatrix(const index_t<R>, const index_t<C>,
MatrixEntries<Operators> const& matrix_operators) const
{
return asmMatrix_ && (!matrix_operators[R][C].assembled || matrix_operators[R][C].changing);
}
/// Return whether the vector-block needs to be assembled
template <std::size_t R, class Operators>
bool assembleVector(const index_t<R>, VectorEntries<Operators> const& rhs_operators) const
template <class Basis0, class... Bases>
auto getElement(Basis0 const& basis0, Bases const&... bases) const
{
return asmVector_ && (!rhs_operators[R].assembled || rhs_operators[R].changing);
return basis0.localView().element();
}
private:
FiniteElementSpaces<FeSpaces> globalBases_;
std::shared_ptr<GlobalBasis> globalBasis_;
bool asmMatrix_;
bool asmVector_;
};
......
dune/amdis/Assembler.inc.hpp
View file @ f9a757d2
......@@ -2,224 +2,131 @@
namespace AMDiS {
template <class FeSpaces>
template <class Operators>
void Assembler<FeSpaces>::assemble(
template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators>
void Assembler<GlobalBasis>::assemble(
GridView const& gv,
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators)
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators)
{
// 1. Update global bases
globalBases_.update(gv);
globalBasis_.update(gv);
// 2. init matrix and rhs vector and initialize dirichlet boundary conditions
initMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators);
initMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators);
// 3. traverse grid and assemble operators on the elements
for (auto const& element : elements(globalBases_.gridView()))
for (auto const& element : elements(globalBasis_.gridView()))
{
globalBases_.bind(element);
MatrixEntries<ElementMatrix> elementMatrix;
VectorEntries<ElementVector> elementVector;
MatrixEntries<bool> addMat = assembleElementMatrices(matrix_operators, elementMatrix);
VectorEntries<bool> addVec = assembleElementVectors(rhs_operators, elementVector);
addElementMatrices(matrix, addMat, elementMatrix);
addElementVectors(rhs, addVec, elementVector);
globalBases_.unbind();
globalBasis_.bind(element);
assembleElement(matrix, rhs, matrixOperators, rhsOperators);
globalBasis_.unbind();
}
// 4. finish matrix insertion and apply dirichlet boundary conditions
std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators);
std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrixOperators, rhsOperators);
msg("fillin of assembled matrix: ", nnz);
}
template <class FeSpaces>
template <class Operators>
void Assembler<FeSpaces>::initMatrixVector(
template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators>
void Assembler<GlobalBasis>::initMatrixVector(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators) const
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const
{
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
static const int R = decltype(_r)::value;
auto const& rowFeSpace = globalBases_.basis(_r);
msg(rowFeSpace.size(), " DOFs for FeSpace[", R, "]");
if (this->assembleVector(_r, rhs_operators)) {
rhs.compress(_r);
rhs.getVector(_r) = 0.0;
// init vector operators
for (auto& scaled : rhs_operators[R].element)
scaled.op->init(rowFeSpace);
for (auto& scaled : rhs_operators[R].boundary)
scaled.op->init(rowFeSpace);
for (auto& scaled : rhs_operators[R].intersection)
scaled.op->init(rowFeSpace);
}
matrix.init(globalBasis_);
solution.init(globalBasis_);
rhs.init(globalBasis_);
forEach(range_<0, nComponents>, [&,this](auto const _c)
{
static const int C = decltype(_c)::value;
auto const& colFeSpace = globalBases_.basis(_c);
bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators);
matrix(_r, _c).init(asmMatrix);
if (asmMatrix) {
// init matrix operators
for (auto& scaled : matrix_operators[R][C].element)
scaled.op->init(rowFeSpace, colFeSpace);
for (auto& scaled : matrix_operators[R][C].boundary)
scaled.op->init(rowFeSpace, colFeSpace);
for (auto& scaled : matrix_operators[R][C].intersection)
scaled.op->init(rowFeSpace, colFeSpace);
// init boundary condition
for (int c = 0; c < nComponents; ++c)
for (auto bc : matrix_operators[R][c].dirichlet)
bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
}
});
});
}
auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
if (rowNode.isLeaf)
msg(globalBasis_.size(rowTreePath), " DOFs for Basis[", Dune::TypeTree::treePathIndex(rowTreePath), "]");
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
if (rhsOperators[rowNode].assemble(asmVector_))
rhsOperators[rowNode].init(rowBasis);
template <class FeSpaces>
template <class Operators>
typename Assembler<FeSpaces>::template MatrixEntries<bool>
Assembler<FeSpaces>::assembleElementMatrices(
MatrixEntries<Operators>& operators,
MatrixEntries<ElementMatrix>& elementMatrix) const
{
MatrixEntries<bool> addMat;
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
static const std::size_t R = decltype(_r)::value;
forEach(range_<0, nComponents>, [&,this](auto const _c)
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
static const std::size_t C = decltype(_c)::value;
// assemble block of element matrix
addMat[R][C] = this->assembleMatrix(_r, _c, operators)
? this->assembleElementMatrix(operators[R][C], elementMatrix[R][C],
globalBases_.localView(_r), globalBases_.localView(_c))
: false;
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
if (matrixOperators[rowNode][colNode].assemble(asmMatrix_))
matrixOperators[rowNode][colNode].init(rowBasis, colBasis);
// init boundary condition
// for (int c = 0; c < nComponents; ++c)
// for (auto bc : matrixOperators[R][c].dirichlet)
// bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
});
});
return addMat;
}
template <class FeSpaces>
template <class Operators>
typename Assembler<FeSpaces>::template VectorEntries<bool>
Assembler<FeSpaces>::assembleElementVectors(
VectorEntries<Operators>& operators,
VectorEntries<ElementVector>& elementVector) const
template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators>
void Assembler<GlobalBasis>::assembleElement(
SystemMatrixType& matrix,
SystemVectorType& rhs,
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const
{
VectorEntries<bool> addVec;
forEach(range_<0, nComponents>, [&,this](auto const _r)
auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
static const std::size_t R = decltype(_r)::value;
// assemble block of element vector
addVec[R] = this->assembleVector(_r, operators)
? this->assembleElementVector(operators[R], elementVector[R], globalBases_.localView(_r))
: false;
});
return addVec;
}
template <class FeSpaces>
template <class Operators, class RowView, class ColView>
bool Assembler<FeSpaces>::assembleElementMatrix(
Operators& operators,
ElementMatrix& elementMatrix,
RowView const& rowLocalView,
ColView const& colLocalView) const
{
if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
return false; // nothing to do
auto const nRows = rowLocalView.tree().finiteElement().size();
auto const nCols = colLocalView.tree().finiteElement().size();
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto const& element = rowLocalView.element();
auto const& gridView = rowLocalView.globalBasis().gridView();
auto& rhsOp = rhsOperators[rowNode];
if (rhsOp.assemble(asmVector_) && !rhsOp.empty()) {
auto elementVector = makeElementVector(rowNode);
set_to_zero(elementVector);
// fills the entire matrix with zeroes
elementMatrix.change_dim(nRows, nCols);
set_to_zero(elementMatrix);
bool add = false;
auto assemble_operators = [&](auto& e, auto& operator_list) {
for (auto scaled : operator_list) {
bool add_op = scaled.op->getElementMatrix(e, rowLocalView, colLocalView, elementMatrix, scaled.factor);
add = add || add_op;
assembleElementOperators(elementVector, rhs, rhsOp, rowBasis);
}
};
// assemble element operators
assemble_operators(element, operators.element);
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);
// assemble intersection operators
if (!operators.intersection.empty()
|| (!operators.boundary.empty() && element.hasBoundaryIntersections()))
{
for (auto const& intersection : intersections(gridView, element)) {
if (intersection.boundary())
assemble_operators(intersection, operators.boundary);
else
assemble_operators(intersection, operators.intersection);
}
}
auto elementMatrix = makeElementMatrix(rowNode, colNode);
set_to_zero(elementMatrix);
return add;
assembleElementOperators(elementMatrix, matrix, matOp, rowBasis, colBasis);
}
});
});
}
template <class FeSpaces>
template <class Operators, class RowView>
bool Assembler<FeSpaces>::assembleElementVector(
template <class GlobalBasis>
template <class ElementContainer, class Container, class Operators, class... Bases>
void Assembler<GlobalBasis>::assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
ElementVector& elementVector,
RowView const& rowLocalView) const
Bases const&... subBases) const
{
if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
if (operators.empty())
return false;
auto const nRows = rowLocalView.tree().finiteElement().size();
auto const& element = rowLocalView.element();
auto const& gridView = rowLocalView.globalBasis().gridView();
// Set all entries to zero
elementVector.change_dim(nRows);
set_to_zero(elementVector);
auto const& element = getElement(subBases...);
auto const& gridView = subBasis.gridView();
bool add = false;
auto assemble_operators = [&](auto& e, auto& operator_list) {
auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) {
bool add_op = scaled.op->getElementVector(e, rowLocalView, elementVector, scaled.factor);
bool add_op = scaled.op->assemble(gridView, context, subBases.localView()..., elementContainer, scaled.factor);
add = add || add_op;
}
};
......@@ -239,109 +146,50 @@ bool Assembler<FeSpaces>::assembleElementVector(
}
}
return add;
}
if (!add)
return;
template <class FeSpaces>
void Assembler<FeSpaces>::addElementMatrices(
SystemMatrixType& dofmatrix,
MatrixEntries<bool> const& addMat,
MatrixEntries<ElementMatrix> const& elementMatrix) const
{
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
static const std::size_t R = decltype(_r)::value;
forEach(range_<0, nComponents>, [&,this](auto const _c)
{
static const std::size_t C = decltype(_c)::value;
if (!addMat[R][C])
return;
auto const& rowIndexSet = globalBases_.localIndexSet(_r);
auto const& colIndexSet = globalBases_.localIndexSet(_c);
// NOTE: current implementation does not utilize the multi-indices that we get from localIndexSet.
for (std::size_t i = 0; i < num_rows(elementMatrix[R][C]); ++i) {
// The global index of the i−th vertex of the element
auto const row = rowIndexSet.index(i);
for (std::size_t j = 0; j < num_cols(elementMatrix[R][C]); ++j) {
// The global index of the j−th vertex of the element
auto const col = colIndexSet.index(j);
dofmatrix(_r,_c)(row,col) += elementMatrix[R][C](i,j);
}
}
});
});
}
template <class FeSpaces>
void Assembler<FeSpaces>::addElementVectors(
SystemVectorType& dofvector,
VectorEntries<bool> const& addVec,
VectorEntries<ElementVector> const& elementVector) const
{
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
static const std::size_t R = decltype(_r)::value;
if (!addVec[R])
return;
auto const& localIndexSet = globalBases_.localIndexSet(_r);
for (std::size_t i = 0; i < size(elementVector[R]); ++i) {
// The global index of the i-th vertex
auto const row = localIndexSet.index(i);
dofvector[_r][row] += elementVector[R][i];
}
});
elementContainer.apply(subBases.localIndexSet()..., container);
}
template <class FeSpaces>
template <class Operators>
std::size_t Assembler<FeSpaces>::finishMatrixVector(
template <class GlobalBasis>
template <class MatrixOperators, class VectorOperators>
std::size_t Assembler<GlobalBasis>::finishMatrixVector(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators) const
MatrixOperators& matrixOperators,
VectorOperators& rhsOperators) const
{
std::size_t nnz = 0;
forEach(range_<0, nComponents>, [&,this](auto const _r)
matrix.finish();
auto localView = globalBasis_.localView();
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
static const int R = decltype(_r)::value;
if (this->assembleVector(_r, rhs_operators))
rhs_operators[R].assembled = true;
auto& rhsOp = rhsOperators[rowNode];
if (rhsOp.assemble(asmVector_))
rhsOp.assembled = true;
forEach(range_<0, nComponents>, [&,this](auto const _c)
forEachNode(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
static const int C = decltype(_c)::value;
bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators);
matrix(_r, _c).finish();
if (asmMatrix) {
matrix_operators[R][C].assembled = true;
// finish boundary condition
for (int c = 0; c < nComponents; ++c) {