Commit beb6bd46 authored by Praetorius, Simon's avatar Praetorius, Simon

FiniteElementSpaces as a simplified CompositeBasis

parent a1c2c4f1
......@@ -16,9 +16,15 @@ namespace AMDiS
template <class FeSpaces>
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>;
......@@ -33,26 +39,21 @@ namespace AMDiS
public:
/// Constructor, stores a shared-pointer to the feSpaces
Assembler(std::shared_ptr<FeSpaces> const& feSpaces_)
: feSpaces(feSpaces_)
{
auto gridView = std::get<0>(*feSpaces).gridView();
forEach(range_<0, nComponents>, [&,this](auto const _r) {
static const std::size_t R = decltype(_r)::value;
std::get<R>(*feSpaces).update(gridView);
});
}
Assembler(std::shared_ptr<FeSpaces> const& feSpaces, bool asmMatrix, bool asmVector)
: globalBases_(feSpaces)
, asmMatrix_(asmMatrix)
, asmVector_(asmVector)
{}
/// Assemble the linear system
template <class Operators>
void assemble(
GridView const& gv,
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
VectorEntries<Operators>& rhs_operators);
private:
/// Sets the system to zero and initializes all operators and boundary conditions
......@@ -62,8 +63,7 @@ namespace AMDiS
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
VectorEntries<Operators>& rhs_operators) const;
/// Assemble a block-matrix of element-matrices and return a matrix of flags, whether
......@@ -71,18 +71,14 @@ namespace AMDiS
template <class Operators>
MatrixEntries<bool> assembleElementMatrices(
MatrixEntries<Operators>& operators,
MatrixEntries<ElementMatrix>& elementMatrix,
FiniteElementSpaces<FeSpaces> const& elementBases,
bool asmMatrix_);
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,
FiniteElementSpaces<FeSpaces> const& elementBases,
bool asmVector_);
VectorEntries<ElementVector>& elementVector) const;
/// Assemble one block of the block-element-matrix
......@@ -92,7 +88,7 @@ namespace AMDiS
Operators& operators,
ElementMatrix& elementMatrix,
RowView const& rowLocalView,
ColView const& colLocalView);
ColView const& colLocalView) const;
/// Assemble one block of the block-element-vector
// The VectorData argument stores all vector-operators
......@@ -100,22 +96,20 @@ namespace AMDiS
bool assembleElementVector(
Operators& operators,
ElementVector& elementVector,
RowView const& rowLocalView);
RowView const& rowLocalView) const;
/// Add the block-element-matrix to the system-matrix
void addElementMatrices(
SystemMatrixType& dofmatrix,
FiniteElementSpaces<FeSpaces> const& elementBases,
MatrixEntries<bool> const& addMat,
MatrixEntries<ElementMatrix> const& elementMatrix);
MatrixEntries<ElementMatrix> const& elementMatrix) const;
/// Add the block-element-vector to the system-vector
void addElementVectors(
SystemVectorType& dofvector,
FiniteElementSpaces<FeSpaces> const& elementBases,
VectorEntries<bool> const& addVec,
VectorEntries<ElementVector> const& elementVector);
VectorEntries<ElementVector> const& elementVector) const;
/// Finish insertion into the matrix and assembles boundary conditions
/// Return the number of nonzeros assembled into the matrix
......@@ -125,27 +119,28 @@ namespace AMDiS
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
VectorEntries<Operators>& rhs_operators) 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, bool asmMatrix_) const
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, bool asmVector_) const
bool assembleVector(const index_t<R>, VectorEntries<Operators> const& rhs_operators) const
{
return asmVector_ && (!rhs_operators[R].assembled || rhs_operators[R].changing);
}
private:
std::shared_ptr<FeSpaces> feSpaces;
FiniteElementSpaces<FeSpaces> globalBases_;
bool asmMatrix_;
bool asmVector_;
};
} // end namespace AMDiS
......
......@@ -5,38 +5,38 @@ namespace AMDiS {
template <class FeSpaces>
template <class Operators>
void Assembler<FeSpaces>::assemble(
GridView const& gv,
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_)
VectorEntries<Operators>& rhs_operators)
{
// 1. init matrix and rhs vector and initialize dirichlet boundary conditions
initMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators, asmMatrix_, asmVector_);
// 1. Update global bases
globalBases_.update(gv);
// 2. create a localView and localIndexSet object for each global basis
FiniteElementSpaces<FeSpaces> elementBases(*feSpaces);
// 2. init matrix and rhs vector and initialize dirichlet boundary conditions
initMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators);
// 3. traverse grid and assemble operators on the elements
for (auto const& element : elements(elementBases.gridView()))
for (auto const& element : elements(globalBases_.gridView()))
{
elementBases.bind(element);
#if 0
globalBases_.bind(element);
MatrixEntries<ElementMatrix> elementMatrix;
VectorEntries<ElementVector> elementVector;
MatrixEntries<bool> addMat = assembleElementMatrices(matrix_operators, elementMatrix, elementBases, asmMatrix_);
VectorEntries<bool> addVec = assembleElementVectors(rhs_operators, elementVector, elementBases, asmVector_);
MatrixEntries<bool> addMat = assembleElementMatrices(matrix_operators, elementMatrix);
VectorEntries<bool> addVec = assembleElementVectors(rhs_operators, elementVector);
addElementMatrices(matrix, addMat, elementMatrix);
addElementVectors(rhs, addVec, elementVector);
addElementMatrices(matrix, elementBases, addMat, elementMatrix);
addElementVectors(rhs, elementBases, addVec, elementVector);
#endif
elementBases.unbind();
globalBases_.unbind();
}
// 4. finish matrix insertion and apply dirichlet boundary conditions
std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators, asmMatrix_, asmVector_);
std::size_t nnz = finishMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators);
msg("fillin of assembled matrix: ", nnz);
}
......@@ -49,17 +49,16 @@ void Assembler<FeSpaces>::initMatrixVector(
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_)
VectorEntries<Operators>& rhs_operators) const
{
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
static const int R = decltype(_r)::value;
auto const& rowFeSpace = std::get<R>(*feSpaces);
auto const& rowFeSpace = globalBases_.basis(_r);
msg(rowFeSpace.size(), " DOFs for FeSpace[", R, "]");
if (this->assembleVector(_r, rhs_operators, asmVector_)) {
if (this->assembleVector(_r, rhs_operators)) {
rhs.compress(_r);
rhs.getVector(_r) = 0.0;
......@@ -75,9 +74,9 @@ void Assembler<FeSpaces>::initMatrixVector(
forEach(range_<0, nComponents>, [&,this](auto const _c)
{
static const int C = decltype(_c)::value;
auto const& colFeSpace = std::get<C>(*feSpaces);
auto const& colFeSpace = globalBases_.basis(_c);
bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators, asmMatrix_);
bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators);
matrix(_r, _c).init(asmMatrix);
if (asmMatrix) {
......@@ -104,9 +103,7 @@ template <class FeSpaces>
typename Assembler<FeSpaces>::template MatrixEntries<bool>
Assembler<FeSpaces>::assembleElementMatrices(
MatrixEntries<Operators>& operators,
MatrixEntries<ElementMatrix>& elementMatrix,
FiniteElementSpaces<FeSpaces> const& elementBases,
bool asmMatrix_)
MatrixEntries<ElementMatrix>& elementMatrix) const
{
MatrixEntries<bool> addMat;
forEach(range_<0, nComponents>, [&,this](auto const _r)
......@@ -117,9 +114,9 @@ Assembler<FeSpaces>::assembleElementMatrices(
static const std::size_t C = decltype(_c)::value;
// assemble block of element matrix
addMat[R][C] = this->assembleMatrix(_r, _c, operators, asmMatrix_)
addMat[R][C] = this->assembleMatrix(_r, _c, operators)
? this->assembleElementMatrix(operators[R][C], elementMatrix[R][C],
elementBases.localView(_r), elementBases.localView(_c))
globalBases_.localView(_r), globalBases_.localView(_c))
: false;
});
});
......@@ -133,9 +130,7 @@ template <class FeSpaces>
typename Assembler<FeSpaces>::template VectorEntries<bool>
Assembler<FeSpaces>::assembleElementVectors(
VectorEntries<Operators>& operators,
VectorEntries<ElementVector>& elementVector,
FiniteElementSpaces<FeSpaces> const& elementBases,
bool asmVector_)
VectorEntries<ElementVector>& elementVector) const
{
VectorEntries<bool> addVec;
forEach(range_<0, nComponents>, [&,this](auto const _r)
......@@ -143,8 +138,8 @@ typename Assembler<FeSpaces>::template VectorEntries<bool>
static const std::size_t R = decltype(_r)::value;
// assemble block of element vector
addVec[R] = this->assembleVector(_r, operators, asmVector_)
? this->assembleElementVector(operators[R], elementVector[R], elementBases.localView(_r))
addVec[R] = this->assembleVector(_r, operators)
? this->assembleElementVector(operators[R], elementVector[R], globalBases_.localView(_r))
: false;
});
......@@ -158,7 +153,7 @@ bool Assembler<FeSpaces>::assembleElementMatrix(
Operators& operators,
ElementMatrix& elementMatrix,
RowView const& rowLocalView,
ColView const& colLocalView)
ColView const& colLocalView) const
{
if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
return false; // nothing to do
......@@ -206,7 +201,7 @@ template <class FeSpaces>
bool Assembler<FeSpaces>::assembleElementVector(
Operators& operators,
ElementVector& elementVector,
RowView const& rowLocalView)
RowView const& rowLocalView) const
{
if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
return false;
......@@ -251,9 +246,8 @@ bool Assembler<FeSpaces>::assembleElementVector(
template <class FeSpaces>
void Assembler<FeSpaces>::addElementMatrices(
SystemMatrixType& dofmatrix,
FiniteElementSpaces<FeSpaces> const& elementBases,
MatrixEntries<bool> const& addMat,
MatrixEntries<ElementMatrix> const& elementMatrix)
MatrixEntries<ElementMatrix> const& elementMatrix) const
{
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
......@@ -264,8 +258,8 @@ void Assembler<FeSpaces>::addElementMatrices(
if (!addMat[R][C])
return;
auto const& rowIndexSet = elementBases.localIndexSet(_r);
auto const& colIndexSet = elementBases.localIndexSet(_c);
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.
......@@ -286,9 +280,8 @@ void Assembler<FeSpaces>::addElementMatrices(
template <class FeSpaces>
void Assembler<FeSpaces>::addElementVectors(
SystemVectorType& dofvector,
FiniteElementSpaces<FeSpaces> const& elementBases,
VectorEntries<bool> const& addVec,
VectorEntries<ElementVector> const& elementVector)
VectorEntries<ElementVector> const& elementVector) const
{
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
......@@ -296,7 +289,7 @@ void Assembler<FeSpaces>::addElementVectors(
if (!addVec[R])
return;
auto const& localIndexSet = elementBases.localIndexSet(_r);
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
......@@ -314,20 +307,19 @@ std::size_t Assembler<FeSpaces>::finishMatrixVector(
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_)
VectorEntries<Operators>& rhs_operators) const
{
std::size_t nnz = 0;
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
static const int R = decltype(_r)::value;
if (this->assembleVector(_r, rhs_operators, asmVector_))
if (this->assembleVector(_r, rhs_operators))
rhs_operators[R].assembled = true;
forEach(range_<0, nComponents>, [&,this](auto const _c)
{
static const int C = decltype(_c)::value;
bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators, asmMatrix_);
bool asmMatrix = this->assembleMatrix(_r, _c, matrix_operators);
matrix(_r, _c).finish();
......
......@@ -38,12 +38,23 @@ namespace AMDiS
using Element = typename GridView::template Codim<0>::Entity;
public:
FiniteElementSpaces(FeSpaces const& feSpaces)
: localViews_(mapTuple([](auto const& basis) { return basis.localView(); }, feSpaces))
, localIndexSets_(mapTuple([](auto const& basis) { return basis.localIndexSet(); }, feSpaces))
, gridView_(std::get<0>(feSpaces).gridView())
FiniteElementSpaces(std::shared_ptr<FeSpaces> const& feSpaces)
: feSpaces_(feSpaces)
, localViews_(mapTuple([](auto const& basis) { return basis.localView(); }, *feSpaces))
, localIndexSets_(mapTuple([](auto const& basis) { return basis.localIndexSet(); }, *feSpaces))
{}
/// Update all global bases. This will update the indexing information of the global basis.
/// NOTE: It must be called if the grid has changed.
void update(GridView const& gv)
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(*feSpaces_).update(gv);
});
}
/// Bind the \ref localViews and \ref localIndexSets to a grid element
void bind(Element const& element)
{
......@@ -67,6 +78,12 @@ namespace AMDiS
});
}
template <std::size_t I>
auto const& basis(const index_t<I> _i = {}) const
{
return std::get<I>(*feSpaces_);
}
template <std::size_t I>
auto const& localView(const index_t<I> _i = {}) const
{
......@@ -81,17 +98,17 @@ namespace AMDiS
auto const& gridView() const
{
return gridView_;
return std::get<0>(*feSpaces_).gridView();
}
private:
std::shared_ptr<FeSpaces> feSpaces_;
/// tuple of localView objects, obtained from the tuple of global bases
LocalViews localViews_;
// tuple of localIndexSet objects, obtained from the tuple of global bases
LocalIndexSets localIndexSets_;
GridView gridView_;
};
} // end namespace AMDiS
......@@ -261,14 +261,16 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
template <class Traits>
void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix_, bool asmVector_)
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix, bool asmVector)
{
Timer t;
Assembler<FeSpaces> assembler(getFeSpaces());
assembler.assemble(*systemMatrix, *solution, *rhs,
matrixOperators, rhsOperators,
asmMatrix_, asmVector_);
Assembler<FeSpaces> assembler(getFeSpaces(), asmMatrix, asmVector);
auto gridView = leafGridView();
assembler.assemble(gridView,
*systemMatrix, *solution, *rhs,
matrixOperators, rhsOperators);
msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment