diff --git a/examples/init/stokes.dat.2d b/examples/init/stokes.dat.2d index 17ba453a4ae4a2205d93d5abbf41459312b409a4..5ec83a051cdf892e33f3965ba0d231dc8bb36ab1 100644 --- a/examples/init/stokes.dat.2d +++ b/examples/init/stokes.dat.2d @@ -6,7 +6,7 @@ stokesMesh->num cells: 4 4 stokes->mesh: stokesMesh -stokes->solver->name: bicgstab_ell +stokes->solver->name: direct stokes->solver->ell: 3 stokes->solver->max iteration: 1000 stokes->solver->relative tolerance: 1e-8 diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt index b0efadeb7ff344592f99b0c6313125fb210ec04b..67c004369b094a17af460416d3ad446c08a65bbf 100644 --- a/src/amdis/CMakeLists.txt +++ b/src/amdis/CMakeLists.txt @@ -25,7 +25,6 @@ install(FILES CreatorInterface.hpp CreatorMap.hpp DirichletBC.hpp - DirichletBC.inc.hpp FileWriter.hpp Flag.hpp GridFunctionOperator.hpp diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp index fa63da2c08aece3ce0f92ed252bc91102ac7baff..82a12b18079d030f929963ac5dde39e2ef97a5ad 100644 --- a/src/amdis/DirichletBC.hpp +++ b/src/amdis/DirichletBC.hpp @@ -6,7 +6,9 @@ #include <vector> #include <dune/common/hybridutilities.hh> +#include <dune/functions/functionspacebases/boundarydofs.hh> #include <dune/functions/functionspacebases/interpolate.hh> +#include <dune/functions/functionspacebases/subspacebasis.hh> #include <amdis/Boundary.hpp> #include <amdis/Output.hpp> @@ -14,6 +16,7 @@ #include <amdis/common/ValueCategory.hpp> #include <amdis/utility/RangeType.hpp> #include <amdis/utility/TreeData.hpp> +#include <amdis/utility/Visitor.hpp> namespace AMDiS { @@ -48,66 +51,61 @@ namespace AMDiS /// Prepare the matrix, rhs and solution block for assembling dirichlet /// boundary conditions, e.g. collect the corresponding indices of boundary /// DOFS for later elimination in \ref finish. - template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis> - void init(Matrix& /*matrix*/, VectorX& /*rhs*/, VectorB& /*solution*/, - RowBasis const& rowBasis, - ColBasis const& /*colBasis*/) + template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode> + void init(Matrix& /*matrix*/, VectorX& /*solution*/, VectorB& rhs, + RowNode const& rowNode, ColNode const&) { if (!initialized_) { - using RowNode = typename RowBasis::LocalView::Tree; - initImpl(rowBasis, typename RowNode::NodeTag{}); + forEachLeafNode_(rowNode, [&](auto const& /*node*/, auto const& tp) { + auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), cat(rowNode.treePath(),tp)); + this->initImpl(subBasis); + }); initialized_ = true; } } - /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row /// to the matrix and optionally delete the corresponding matrix-column. - /** If DBC is set for matrix block {R,C}, then \p apply_row means that - * the block-row R is assembled, and \p apply_col means that the block-col C - * is assembled. The \p matrix, \p rhs, and \p solution correspond to the - * currently visited blocks of system-matrix and system-vector. - **/ - template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis> - void finish(Matrix& matrix, - VectorX& rhs, - VectorB& solution, - RowBasis const& rowBasis, - ColBasis const& colBasis) + template <class Matrix, class VectorX, class VectorB, class RowNode, class ColNode> + void finish(Matrix& matrix, VectorX& solution, VectorB& rhs, + RowNode const& rowNode, ColNode const& colNode) { test_exit_dbg(initialized_, "Boundary condition not initialized!"); auto columns = matrix.applyDirichletBC(dirichletNodes_); - using Dune::Functions::interpolate; - Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{}, + Dune::Hybrid::ifElse(std::is_same<RangeType_t<RowNode>, Range>{}, [&](auto id) { - interpolate(id(rowBasis), rhs, values_, dirichletNodes_); + auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowNode.treePath()); + Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_); }); - Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{}, + Dune::Hybrid::ifElse(std::is_same<RangeType_t<ColNode>, Range>{}, [&](auto id) { - interpolate(id(colBasis), solution, values_, dirichletNodes_); + auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colNode.treePath()); + Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_); }); // subtract columns of dirichlet nodes from rhs - // for (auto const& triplet : columns) - // rhs[triplet.row] -= triplet.value * solution[triplet.col]; + for (auto const& triplet : columns) + rhs[triplet.row] -= triplet.value * solution[triplet.col]; + initialized_ = false; } protected: - template <class RowBasis> - void initImpl(RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag); - - template <class RowBasis> - void initImpl(RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag); - - template <class RowBasis> - void initImpl(RowBasis const& rowBasis, Dune::TypeTree::CompositeNodeTag); - - template <class RB, class RowNodeTag> - void initImpl(RB const&, RowNodeTag) {} + // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not. + template <class SubBasis> + void initImpl(SubBasis const& subBasis) + { + using Dune::Functions::forEachBoundaryDOF; + using LV = typename SubBasis::LocalView; + using I = typename SubBasis::GridView::Intersection; + dirichletNodes_.resize(subBasis.dimension()); + forEachBoundaryDOF(subBasis, [&](int localIndex, LV const& localView, I const& intersection) { + dirichletNodes_[localView.index(localIndex)] = predicate_(intersection.geometry().center()); + }); + } private: std::function<bool(Domain)> predicate_; @@ -131,5 +129,3 @@ namespace AMDiS using Constraints = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>; } // end namespace AMDiS - -#include "DirichletBC.inc.hpp" diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp deleted file mode 100644 index c2536f44d9a2c2b181694a47395f5d1151e6aa88..0000000000000000000000000000000000000000 --- a/src/amdis/DirichletBC.inc.hpp +++ /dev/null @@ -1,53 +0,0 @@ -#pragma once - -#include <dune/functions/functionspacebases/boundarydofs.hh> -#include <dune/functions/functionspacebases/interpolate.hh> -#include <dune/functions/functionspacebases/subspacebasis.hh> - -namespace AMDiS -{ - template <class WorldVector, class Range> - template <class RowBasis> - void DirichletBC<WorldVector, Range>::initImpl( - RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag) - { - using Dune::Functions::forEachBoundaryDOF; - using LV = typename RowBasis::LocalView; - using I = typename RowBasis::GridView::Intersection; - dirichletNodes_.resize(rowBasis.dimension()); - forEachBoundaryDOF(rowBasis, [&](int localIndex, LV const& localView, I const& intersection) { - dirichletNodes_[localView.index(localIndex)] = predicate_(intersection.geometry().center()); - }); - } - - template <class WorldVector, class Range> - template <class RowBasis> - void DirichletBC<WorldVector, Range>::initImpl( - RowBasis const& rowBasis, Dune::TypeTree::PowerNodeTag) - { - using Dune::Functions::subspaceBasis; - using Node = typename RowBasis::LocalView::Tree; - using ChildNode = typename Node::template Child<0>::type; - - auto tp = rowBasis.prefixPath(); - auto const& basis = rowBasis.rootBasis(); - for (std::size_t i = 0; i < degree(rowBasis.localView().tree()); ++i) - initImpl(subspaceBasis(basis, push_back(tp,i)), typename ChildNode::NodeTag{}); - } - - template <class WorldVector, class Range> - template <class RowBasis> - void DirichletBC<WorldVector, Range>::initImpl( - RowBasis const& rowBasis, Dune::TypeTree::CompositeNodeTag) - { - using Dune::Functions::subspaceBasis; - using Node = typename RowBasis::LocalView::Tree; - auto tp = rowBasis.prefixPath(); - auto const& basis = rowBasis.rootBasis(); - forEach(range_<0,Node::CHILDREN>, [&,this](auto const _i) { - using ChildNode = typename Node::template Child<decltype(_i)::value>::type; - this->initImpl(subspaceBasis(basis, push_back(tp,_i)), typename ChildNode::NodeTag{}); - }); - } - -} // end namespace AMDiS diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index 4972a98e3eedac95e15f08d1afd790ce58382400..c0a57d87ffac313acdeb96fb756a7de226dcbeeb 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -282,8 +282,7 @@ namespace AMDiS void initGlobalBasis(GlobalBasis const& globalBasis) { - localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis_->localView()); - constraints_.init(localView_->tree(), localView_->tree(), tag::store{}); + constraints_.init(*globalBasis_, *globalBasis_); } void createMatricesAndVectors() @@ -292,7 +291,8 @@ namespace AMDiS solution_ = std::make_shared<SystemVector>(*globalBasis_); rhs_ = std::make_shared<SystemVector>(*globalBasis_); - AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath) + auto localView = globalBasis_->localView(); + AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath) { std::string i = to_string(treePath); estimates_[i].resize(globalBasis_->gridView().indexSet().size(0)); @@ -358,7 +358,6 @@ namespace AMDiS /// FE spaces of this problem. std::shared_ptr<GlobalBasis> globalBasis_; - std::shared_ptr<typename GlobalBasis::LocalView> localView_; /// A FileWriter object std::list<std::shared_ptr<FileWriterInterface>> filewriter_; diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 87c525af86af2ed73ef8c547fc080fc4e31f0452..cff5b8f227afccf643446a65d74b924e3eca3a15 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -120,7 +120,8 @@ void ProblemStat<Traits>::initialize( template <class Traits> void ProblemStat<Traits>::createMarker() { - AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath) + auto localView = globalBasis_->localView(); + AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath) { std::string componentName = name_ + "->marker[" + to_string(treePath) + "]"; @@ -144,7 +145,8 @@ void ProblemStat<Traits>::createMarker() template <class Traits> void ProblemStat<Traits>::createFileWriter() { - forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath) + auto localView = globalBasis_->localView(); + forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath) { std::string componentName = name_ + "->output[" + to_string(treePath) + "]"; @@ -166,8 +168,9 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val static_assert( Concepts::Functor<Predicate, bool(WorldVector)>, "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept"); - auto i = child(localView_->tree(), makeTreePath(row)); - auto j = child(localView_->tree(), makeTreePath(col)); + auto localView = globalBasis_->localView(); + auto i = child(localView.tree(), makeTreePath(row)); + auto j = child(localView.tree(), makeTreePath(col)); auto valueGridFct = makeGridFunction(values, globalBasis_->gridView()); @@ -261,40 +264,36 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as rhs_->init(asmVector); solution_->resize(*globalBasis_); - forEachNode_(localView_->tree(), [&,this](auto const& rowNode, auto rowTreePath) { - auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTreePath); - forEachNode_(localView_->tree(), [&,this](auto const& colNode, auto colTreePath) { - auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTreePath); + auto localView = globalBasis_->localView(); + forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) { + forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) { for (auto bc : constraints_[rowNode][colNode]) - bc->init(*systemMatrix_, *solution_, *rhs_, rowBasis, colBasis); + bc->init(*systemMatrix_, *solution_, *rhs_, rowNode, colNode); }); }); msg("{} total DOFs", globalBasis_->dimension()); // 2. traverse grid and assemble operators on the elements for (auto const& element : elements(gridView())) { - localView_->bind(element); - auto geometry = element.geometry(); + localView.bind(element); if (asmMatrix) - systemMatrix_->assemble(element, geometry, *localView_, *localView_); + systemMatrix_->assemble(localView, localView); if (asmVector) - rhs_->assemble(element, geometry, *localView_); + rhs_->assemble(localView); - localView_->unbind(); + localView.unbind(); } // 3. finish matrix insertion and apply dirichlet boundary conditions systemMatrix_->finish(asmMatrix); rhs_->finish(asmVector); - forEachNode_(localView_->tree(), [&,this](auto const& rowNode, auto rowTreePath) { - auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTreePath); - forEachNode_(localView_->tree(), [&,this](auto const& colNode, auto colTreePath) { - auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTreePath); + forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) { + forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) { // finish boundary condition for (auto bc : constraints_[rowNode][colNode]) - bc->finish(*systemMatrix_, *solution_, *rhs_, rowBasis, colBasis); + bc->finish(*systemMatrix_, *solution_, *rhs_, rowNode, colNode); }); }); diff --git a/src/amdis/linear_algebra/CMakeLists.txt b/src/amdis/linear_algebra/CMakeLists.txt index 5580cbc926dd01f8e9e63381ad4071253601b903..a436f7713b38c10fe2b334a1bc2dcd8541b5cc43 100644 --- a/src/amdis/linear_algebra/CMakeLists.txt +++ b/src/amdis/linear_algebra/CMakeLists.txt @@ -1,7 +1,9 @@ install(FILES Common.hpp DOFMatrixBase.hpp + DOFMatrixBase.inc.hpp DOFVectorBase.hpp + DOFVectorBase.inc.hpp DOFVectorInterface.hpp HierarchicWrapper.hpp LinearSolver.hpp diff --git a/src/amdis/linear_algebra/Common.hpp b/src/amdis/linear_algebra/Common.hpp index e7bf16d923776bf233d6114d4706986d9929d9b3..197cfe993b8088235ce4e5c5a0b1f03c362371ec 100644 --- a/src/amdis/linear_algebra/Common.hpp +++ b/src/amdis/linear_algebra/Common.hpp @@ -7,7 +7,7 @@ namespace AMDiS template <class T> struct Triplet { - std::size_t row, cols; + std::size_t row, col; T value; }; diff --git a/src/amdis/linear_algebra/DOFMatrixBase.hpp b/src/amdis/linear_algebra/DOFMatrixBase.hpp index c1835c6e79f9dc0f3ec304d22d78d85f1ae78423..eef7e4a7c294ae2658e5d207ef2b0fb5221ad820 100644 --- a/src/amdis/linear_algebra/DOFMatrixBase.hpp +++ b/src/amdis/linear_algebra/DOFMatrixBase.hpp @@ -4,10 +4,7 @@ #include <dune/common/timer.hh> -#include <amdis/Assembler.hpp> -#include <amdis/LocalAssembler.hpp> #include <amdis/LocalAssemblerList.hpp> -#include <amdis/LocalOperator.hpp> #include <amdis/common/Math.hpp> #include <amdis/utility/MultiIndex.hpp> #include <amdis/utility/TreePath.hpp> @@ -40,10 +37,8 @@ namespace AMDiS DOFMatrixBase(RowBasis const& rowBasis, ColBasis const& colBasis) : rowBasis_(&rowBasis) , colBasis_(&colBasis) - , rowLocalView_(rowBasis.localView()) - , colLocalView_(colBasis.localView()) { - operators_.init(rowLocalView_.tree(), colLocalView_.tree(), tag::store{}); + operators_.init(rowBasis, colBasis); } /// Return the row-basis \ref rowBasis of the matrix @@ -84,46 +79,15 @@ namespace AMDiS /// Initialize the matrix for insertion, e.g. allocate the non-zero pattern /// If \p setToZero is true, the matrix is set to 0 - void init(bool asmMatrix) - { - backend_.init(*rowBasis_, *colBasis_, asmMatrix); - - auto const rowSize = rowBasis_->localView().maxSize(); - auto const colSize = colBasis_->localView().maxSize(); - elementMatrix_.resize(rowSize, colSize); - } + void init(bool asmMatrix); /// Finish the matrix insertion, e.g. cleanup or final insertion - void finish(bool asmMatrix) - { - backend_.finish(); - - if (asmMatrix) { - forEachNode_(rowBasis_->localView().tree(), [&,this](auto const& rowNode, auto&&) { - forEachNode_(colBasis_->localView().tree(), [&,this](auto const& colNode, auto&&) { - auto& matOp = operators_[rowNode][colNode]; - if (matOp.doAssemble()) - matOp.assembled = true; - }); - }); - } - } + void finish(bool asmMatrix); /// Insert a block of values into the matrix (add to existing values) void insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView, - ElementMatrix const& elementMatrix) - { - for (size_type i = 0; i < rowLocalView.size(); ++i) { - size_type const row = flatMultiIndex(rowLocalView.index(i)); - for (size_type j = 0; j < colLocalView.size(); ++j) { - if (std::abs(elementMatrix[i][j]) > threshold<double>) { - size_type const col = flatMultiIndex(colLocalView.index(j)); - backend_.insert(row, col, elementMatrix[i][j]); - } - } - } - } + ElementMatrix const& elementMatrix); /// Insert a single value into the matrix (add to existing value) template <class RowIndex, class ColIndex> @@ -138,53 +102,18 @@ namespace AMDiS return operators_[row][col]; } + /// Associate a local operator with this DOFMatrix template <class ContextTag, class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath> void addOperator(ContextTag contextTag, Operator const& preOp, - RowTreePath row = {}, ColTreePath col = {}) - { - static_assert( Concepts::PreTreePath<RowTreePath>, - "row must be a valid treepath, or an integer/index-constant"); - static_assert( Concepts::PreTreePath<ColTreePath>, - "col must be a valid treepath, or an integer/index-constant"); + RowTreePath row = {}, ColTreePath col = {}); - auto i = child(rowBasis_->localView().tree(), makeTreePath(row)); - auto j = child(colBasis_->localView().tree(), makeTreePath(col)); + /// Assemble the matrix operators on the bound element. + void assemble(RowLocalView const& rowLocalView, + ColLocalView const& colLocalView); - using Context = typename ContextTag::type; - auto op = makeLocalOperator<Context>(preOp, rowBasis_->gridView()); - auto localAssembler = makeLocalAssemblerPtr<Context>(std::move(op), i, j); - - operators_[i][j].push(contextTag, localAssembler); - operators_[i][j].changing = true; - } - - // Assemble the operators on the `element`. Assume that `rowLocalView` - // and `colLocalView` are already bound to the element. - void assemble(Element const& element, Geometry const& geometry, - RowLocalView const& rowLocalView, - ColLocalView const& colLocalView) - { - elementMatrix_ = 0; - - forEachNode_(rowLocalView.tree(), [&](auto const& rowNode, auto&&) { - forEachNode_(colLocalView.tree(), [&](auto const& colNode, auto&&) { - auto& matOp = operators_[rowNode][colNode]; - if (matOp.doAssemble() && !matOp.empty()) { - matOp.bind(element, geometry); - auto matAssembler = [&](auto const& context, auto& operators) { - for (auto scaled : operators) - scaled.op->assemble(context, rowNode, colNode, elementMatrix_); - }; - - assembleOperators(rowBasis_->gridView(), element, matOp, matAssembler); - matOp.unbind(); - } - }); - }); - - insert(rowLocalView, colLocalView, elementMatrix_); - } + /// Assemble all matrix operators, TODO: incooperate boundary conditions + void assemble(); /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds /// a one on the diagonal of the matrix. @@ -199,13 +128,9 @@ namespace AMDiS } protected: - /// The finite element space / basis associated with the data matrix - RowBasis const* rowBasis_; - ColBasis const* colBasis_; - - /// The LocalViews aassociated with the row/col basis - RowLocalView rowLocalView_; - ColLocalView colLocalView_; + /// The finite element space / basis associated with the rows/columns + RowBasis const* rowBasis_; + ColBasis const* colBasis_; /// Data backend Backend backend_; @@ -218,3 +143,5 @@ namespace AMDiS }; } // end namespace AMDiS + +#include "DOFMatrixBase.inc.hpp" diff --git a/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp b/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp new file mode 100644 index 0000000000000000000000000000000000000000..c39950bcb963634ad066d4febdf893dcd0d3e7e1 --- /dev/null +++ b/src/amdis/linear_algebra/DOFMatrixBase.inc.hpp @@ -0,0 +1,131 @@ +#pragma once + +#include <amdis/Assembler.hpp> +#include <amdis/LocalAssembler.hpp> +#include <amdis/LocalOperator.hpp> +#include <amdis/utility/Visitor.hpp> + +namespace AMDiS { + +template <class RB, class CB, class B> +void DOFMatrixBase<RB,CB,B>:: +init(bool asmMatrix) +{ + backend_.init(*rowBasis_, *colBasis_, asmMatrix); + + auto const rowSize = rowBasis_->localView().maxSize(); + auto const colSize = colBasis_->localView().maxSize(); + elementMatrix_.resize(rowSize, colSize); +} + + +template <class RB, class CB, class B> +void DOFMatrixBase<RB,CB,B>:: +finish(bool asmMatrix) +{ + backend_.finish(); + + if (asmMatrix) { + forEachNode_(rowBasis_->localView().tree(), [&](auto const& rowNode, auto) { + forEachNode_(colBasis_->localView().tree(), [&](auto const& colNode, auto) { + auto& matOp = operators_[rowNode][colNode]; + if (matOp.doAssemble()) + matOp.assembled = true; + }); + }); + } +} + + +template <class RB, class CB, class B> +void DOFMatrixBase<RB,CB,B>:: +insert(RowLocalView const& rowLocalView, ColLocalView const& colLocalView, + ElementMatrix const& elementMatrix) +{ + for (size_type i = 0; i < rowLocalView.size(); ++i) { + size_type const row = flatMultiIndex(rowLocalView.index(i)); + for (size_type j = 0; j < colLocalView.size(); ++j) { + if (std::abs(elementMatrix[i][j]) > threshold<double>) { + size_type const col = flatMultiIndex(colLocalView.index(j)); + backend_.insert(row, col, elementMatrix[i][j]); + } + } + } +} + + +template <class RB, class CB, class B> + template <class ContextTag, class Operator, class RowTreePath, class ColTreePath> +void DOFMatrixBase<RB,CB,B>:: +addOperator(ContextTag contextTag, Operator const& preOp, + RowTreePath row, ColTreePath col) +{ + static_assert( Concepts::PreTreePath<RowTreePath>, + "row must be a valid treepath, or an integer/index-constant"); + static_assert( Concepts::PreTreePath<ColTreePath>, + "col must be a valid treepath, or an integer/index-constant"); + + auto i = child(rowBasis_->localView().tree(), makeTreePath(row)); + auto j = child(colBasis_->localView().tree(), makeTreePath(col)); + + using Context = typename ContextTag::type; + auto op = makeLocalOperator<Context>(preOp, rowBasis_->gridView()); + auto localAssembler = makeLocalAssemblerPtr<Context>(std::move(op), i, j); + + operators_[i][j].push(contextTag, localAssembler); + operators_[i][j].changing = true; +} + + +template <class RB, class CB, class B> +void DOFMatrixBase<RB,CB,B>:: +assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView) +{ + elementMatrix_ = 0; + auto const& gv = rowBasis_->gridView(); + auto const& element = rowLocalView.element(); + auto geometry = element.geometry(); + + forEachNode_(rowLocalView.tree(), [&](auto const& rowNode, auto) { + forEachNode_(colLocalView.tree(), [&](auto const& colNode, auto) { + auto& matOp = operators_[rowNode][colNode]; + if (matOp.doAssemble() && !matOp.empty()) { + matOp.bind(element, geometry); + auto matAssembler = [&](auto const& context, auto& operators) { + for (auto scaled : operators) + scaled.op->assemble(context, rowNode, colNode, elementMatrix_); + }; + assembleOperators(gv, element, matOp, matAssembler); + matOp.unbind(); + } + }); + }); + + insert(rowLocalView, colLocalView, elementMatrix_); +} + + +template <class RB, class CB, class B> +void DOFMatrixBase<RB,CB,B>:: +assemble() +{ + auto rowLocalView = rowBasis_->localView(); + auto colLocalView = colBasis_->localView(); + + init(true); + for (auto const& element : elements(rowBasis_->gridView())) { + rowLocalView.bind(element); + if (rowBasis_ == colBasis_) + assemble(rowLocalView, rowLocalView); + else { + colLocalView.bind(element); + assemble(rowLocalView, colLocalView); + colLocalView.unbind(); + } + rowLocalView.unbind(element); + } + finish(true); +} + + +} // end namespace AMDiS diff --git a/src/amdis/linear_algebra/DOFVectorBase.hpp b/src/amdis/linear_algebra/DOFVectorBase.hpp index f57f6404776198eb0988b9f9d96f04814c7e0f70..7ed39476b8816be0f6a3ab301a5a72d2caa7db7b 100644 --- a/src/amdis/linear_algebra/DOFVectorBase.hpp +++ b/src/amdis/linear_algebra/DOFVectorBase.hpp @@ -4,10 +4,7 @@ #include <dune/functions/functionspacebases/sizeinfo.hh> -#include <amdis/Assembler.hpp> -#include <amdis/LocalAssembler.hpp> #include <amdis/LocalAssemblerList.hpp> -#include <amdis/LocalOperator.hpp> #include <amdis/common/Math.hpp> #include <amdis/common/ScalarTypes.hpp> #include <amdis/linear_algebra/DOFVectorInterface.hpp> @@ -45,10 +42,9 @@ namespace AMDiS /// Constructor. Constructs new BaseVector. DOFVectorBase(BasisType const& basis) : basis_(&basis) - , localView_(basis.localView()) { compress(); - operators_.init(localView_.tree(), tag::store{}); + operators_.init(basis); } DOFVectorBase(Self const&) = default; @@ -57,7 +53,7 @@ namespace AMDiS /// Copy assignment operator Self& operator=(Self const& that) { - assert(basis_ == that.basis_); + basis_ = that.basis_; backend_.resize(that.size()); backend_ = that.backend_; return *this; @@ -129,37 +125,12 @@ namespace AMDiS return backend_[i]; } - void init(bool asmVector) - { - backend_.resize(size()); - if (asmVector) - backend_.set(0); + void init(bool asmVector); - auto const localSize = basis_->localView().maxSize(); - elementVector_.resize(localSize); - } - - void finish(bool asmVector) - { - if (asmVector) { - forEachNode_(basis_->localView().tree(), [&](auto const& node, auto&&) { - auto& rhsOp = operators_[node]; - if (rhsOp.doAssemble()) - rhsOp.assembled = true; - }); - } - } + void finish(bool asmVector); /// Insert a block of values into the matrix (add to existing values) - void insert(LocalView const& localView, ElementVector const& elementVector) - { - for (size_type i = 0; i < localView.size(); ++i) { - if (std::abs(elementVector[i]) > threshold<double>) { - size_type const idx = flatMultiIndex(localView.index(i)); - backend_[idx] += elementVector[i]; - } - } - } + void insert(LocalView const& localView, ElementVector const& elementVector); template <class Node> auto& operators(Node const& node) @@ -167,56 +138,20 @@ namespace AMDiS return operators_[node]; } - template <class ContextTag, class Operator, - class TreePath = RootTreePath> - void addOperator(ContextTag contextTag, Operator const& preOp, - TreePath path = {}) - { - static_assert( Concepts::PreTreePath<TreePath>, - "path must be a valid treepath, or an integer/index-constant"); - - auto i = child(basis_->localView().tree(), makeTreePath(path)); - - using Context = typename ContextTag::type; - auto op = makeLocalOperator<Context>(preOp, basis_->gridView()); - auto localAssembler = makeLocalAssemblerPtr<Context>(std::move(op), i); - - operators_[i].push(contextTag, localAssembler); - operators_[i].changing = true; - } - - // Assemble the operators on the `element`. Assume that `localView` - // is already bound to the element. - void assemble(Element const& element, Geometry const& geometry, - LocalView const& localView) - { - elementVector_ = 0; + /// Associate a local operator with this DOFVector + template <class ContextTag, class Operator, class TreePath = RootTreePath> + void addOperator(ContextTag contextTag, Operator const& preOp, TreePath path = {}); - forEachNode_(localView.tree(), [&](auto const& node, auto&&) { - auto& rhsOp = operators_[node]; - if (rhsOp.doAssemble() && !rhsOp.empty()) { - rhsOp.bind(element, geometry); + /// Assemble the vector operators on the bound element. + void assemble(LocalView const& localView); - auto vecAssembler = [&](auto const& context, auto& operators) { - for (auto scaled : operators) - scaled.op->assemble(context, node, elementVector_); - }; - - assembleOperators(basis_->gridView(), element, rhsOp, vecAssembler); - rhsOp.unbind(); - } - }); - - insert(localView, elementVector_); - } + /// Assemble all vector operators + void assemble(); private: /// The finite element space / basis associated with the data vector Basis const* basis_; - /// The LocalView aassociated with the basis - LocalView localView_; - /// Data backend Backend backend_; @@ -228,3 +163,5 @@ namespace AMDiS }; } // end namespace AMDiS + +#include "DOFVectorBase.inc.hpp" diff --git a/src/amdis/linear_algebra/DOFVectorBase.inc.hpp b/src/amdis/linear_algebra/DOFVectorBase.inc.hpp new file mode 100644 index 0000000000000000000000000000000000000000..4fc018f0322592750d988a81ba91b3a9abeba490 --- /dev/null +++ b/src/amdis/linear_algebra/DOFVectorBase.inc.hpp @@ -0,0 +1,112 @@ +#pragma once + +#include <amdis/Assembler.hpp> +#include <amdis/LocalAssembler.hpp> +#include <amdis/LocalOperator.hpp> +#include <amdis/utility/Visitor.hpp> + +namespace AMDiS { + +template <class BS, class B> +void DOFVectorBase<BS,B>:: +init(bool asmVector) +{ + backend_.resize(size()); + if (asmVector) + backend_.set(0); + + auto const localSize = basis_->localView().maxSize(); + elementVector_.resize(localSize); +} + + +template <class BS, class B> +void DOFVectorBase<BS,B>:: +finish(bool asmVector) +{ + if (asmVector) { + forEachNode_(basis_->localView().tree(), [&](auto const& node, auto) { + auto& rhsOp = operators_[node]; + if (rhsOp.doAssemble()) + rhsOp.assembled = true; + }); + } +} + + +template <class BS, class B> +void DOFVectorBase<BS,B>:: +insert(LocalView const& localView, ElementVector const& elementVector) +{ + for (size_type i = 0; i < localView.size(); ++i) { + if (std::abs(elementVector[i]) > threshold<double>) { + size_type const idx = flatMultiIndex(localView.index(i)); + backend_[idx] += elementVector[i]; + } + } +} + + +template <class BS, class B> + template <class ContextTag, class Operator, class TreePath> +void DOFVectorBase<BS,B>:: +addOperator(ContextTag contextTag, Operator const& preOp, TreePath path) +{ + static_assert( Concepts::PreTreePath<TreePath>, + "path must be a valid treepath, or an integer/index-constant"); + + auto i = child(basis_->localView().tree(), makeTreePath(path)); + + using Context = typename ContextTag::type; + auto op = makeLocalOperator<Context>(preOp, basis_->gridView()); + auto localAssembler = makeLocalAssemblerPtr<Context>(std::move(op), i); + + operators_[i].push(contextTag, localAssembler); + operators_[i].changing = true; +} + + +template <class BS, class B> +void DOFVectorBase<BS,B>:: +assemble(LocalView const& localView) +{ + elementVector_ = 0; + auto const& gv = basis_->gridView(); + auto const& element = localView.element(); + auto geometry = element.geometry(); + + forEachNode_(localView.tree(), [&](auto const& node, auto) { + auto& rhsOp = operators_[node]; + if (rhsOp.doAssemble() && !rhsOp.empty()) { + rhsOp.bind(element, geometry); + + auto vecAssembler = [&](auto const& context, auto& operators) { + for (auto scaled : operators) + scaled.op->assemble(context, node, elementVector_); + }; + + assembleOperators(gv, element, rhsOp, vecAssembler); + rhsOp.unbind(); + } + }); + + insert(localView, elementVector_); +} + + +template <class BS, class B> +void DOFVectorBase<BS,B>:: +assemble() +{ + auto localView = basis_->localView(); + + init(true); + for (auto const& element : elements(basis_->gridView())) { + localView.bind(element); + assemble(localView); + localView.unbind(element); + } + finish(true); +} + +} // end namespace AMDiS diff --git a/src/amdis/utility/TreeData.hpp b/src/amdis/utility/TreeData.hpp index 83a35f1bce22c7d65ca92ee1981b6eee032ecc7d..d3842885d812cb28309fc7a897ad34819ec96411 100644 --- a/src/amdis/utility/TreeData.hpp +++ b/src/amdis/utility/TreeData.hpp @@ -10,11 +10,6 @@ namespace AMDiS { - namespace tag - { - struct store {}; - } - /** * \brief Container allowing to attach data to each node of a tree * @@ -35,16 +30,16 @@ namespace AMDiS * node type is known. Hence the tree will be traversed on initilization, * copy, assignment, and destruction of a TreeData container. * - * \tparam T Type of the tree + * \tparam Basis Type of the basis with basis-tree * \tparam ND The data stored for a node of type Node will be of type ND<Node> * \tparam LO Set this flag if data should only be attached to leaf nodes. */ - template <class T, template<class> class ND, bool LO> + template <class Basis, template<class> class ND, bool LO> class TreeData { public: //! Type of tree the data is associated with - using Tree = T; + using Tree = typename Basis::LocalView::Tree; //! Type used for indices and size information using size_type = typename Tree::size_type; @@ -59,7 +54,7 @@ namespace AMDiS public: //! Default constructor TreeData() - : tree_(nullptr) + : basis_(nullptr) {} /** @@ -70,19 +65,15 @@ namespace AMDiS * of the tree data. * See also \ref init. */ - TreeData(Tree const& tree, tag::store) - : tree_(&tree) + explicit TreeData(Basis const& basis) + : basis_(&basis) { initData(); } - explicit TreeData(Tree& tree) - : TreeData(tree, tag::store{}) - {} - //! Copy constructor TreeData(TreeData const& other) - : tree_(other.tree_) + : basis_(other.basis_) { initData(); copyData(other); @@ -102,18 +93,13 @@ namespace AMDiS * A reference to the tree is stored because it's needed for destruction * of the tree data. */ - void init(Tree const& tree, tag::store) + void init(Basis const& basis) { destroyData(); - tree_ = &tree; + basis_ = &basis; initData(); } - void init(Tree& tree) - { - init(tree, tag::store{}); - } - //! Copy and Move assignment TreeData& operator=(TreeData other) { @@ -127,12 +113,6 @@ namespace AMDiS destroyData(); } - //! Return the attached tree (maybe nullptr) - Tree const* tree() const - { - return tree_; - } - //! Get mutable reference to data associated to given node template<class Node> NodeData<Node>& operator[](const Node& node) @@ -155,11 +135,16 @@ namespace AMDiS void swap(TreeData& other) { using std::swap; - swap(tree_, other.tree_); + swap(basis_, other.basis_); swap(data_, other.data_); swap(initialized_, other.initialized_); } + Basis const* basis() const + { + return basis_; + } + protected: // For each node allocate memory and store the void* in the \ref data_ void initData() @@ -193,7 +178,7 @@ namespace AMDiS auto* p = (NodeData<Node>*)(data_[node.treeIndex()]); delete p; }); - tree_ = nullptr; + basis_ = nullptr; initialized_ = false; } @@ -202,18 +187,18 @@ namespace AMDiS template <class Func> void apply(Func&& func) { - if (tree_) + if (basis_) applyImpl(func, std::integral_constant<bool, leafOnly>{}); } template <class Func> - void applyImpl(Func&& func, std::true_type) { forEachLeafNode_(*tree_, func); } + void applyImpl(Func&& func, std::true_type) { forEachLeafNode_(basis_->localView().tree(), func); } template <class Func> - void applyImpl(Func&& func, std::false_type) { forEachNode_(*tree_, func); } + void applyImpl(Func&& func, std::false_type) { forEachNode_(basis_->localView().tree(), func); } protected: - const Tree* tree_ = nullptr; + Basis const* basis_ = nullptr; std::vector<void*> data_; bool initialized_ = false; @@ -222,7 +207,7 @@ namespace AMDiS namespace Impl { - template <class Tree, template <class,class> class NodeData> + template <class Basis, template <class,class> class NodeData> struct RowNodeData { template <class RowNode> @@ -233,42 +218,34 @@ namespace AMDiS }; template <class RowNode> - using type = TreeData<Tree, ColNodeData<RowNode>::template type, false>; + using type = TreeData<Basis, ColNodeData<RowNode>::template type, false>; }; } // end namespace Impl template <class RowBasis, class ColBasis, template <class,class> class NodeData, - class RowTree = typename RowBasis::LocalView::Tree, - class ColTree = typename ColBasis::LocalView::Tree, - template <class> class ND = Impl::RowNodeData<ColTree, NodeData>::template type> + template <class> class ND = Impl::RowNodeData<ColBasis, NodeData>::template type> class MatrixData - : public TreeData<RowTree, ND, false> + : public TreeData<RowBasis, ND, false> { - using Super = TreeData<RowTree, ND, false>; + using Super = TreeData<RowBasis, ND, false>; public: - void init(RowTree const& row, ColTree const& col, tag::store) + void init(RowBasis const& rowBasis, ColBasis const& colBasis) { - Super::init(row, tag::store{}); - forEachNode_(row, [&,this](auto const& node, auto&&) + Super::init(rowBasis); + forEachNode_(rowBasis.localView().tree(), [&](auto const& node, auto&&) { - (*this)[node].init(col, tag::store{}); + (*this)[node].init(colBasis); }); } - - void init(RowTree& row, ColTree& col) - { - init(row, col, tag::store{}); - } }; - template <class GlobalBasis, template <class> class NodeData, - class Tree = typename GlobalBasis::LocalView::Tree> + template <class Basis, template <class> class NodeData> class VectorData - : public TreeData<Tree, NodeData, false> + : public TreeData<Basis, NodeData, false> {}; } // end namespace AMDiS diff --git a/src/amdis/utility/TreePath.hpp b/src/amdis/utility/TreePath.hpp index 4949a9281ee4025edd45171809451ffddd1a2603..6a1b695d2158e3ba0a6b4d3774851069bc34bd62 100644 --- a/src/amdis/utility/TreePath.hpp +++ b/src/amdis/utility/TreePath.hpp @@ -223,5 +223,11 @@ namespace AMDiS return Impl::popBackImpl(tp, std::make_index_sequence<sizeof...(T)>{}); } + template <class... S, class... T> + Dune::TypeTree::HybridTreePath<S...,T...> cat(Dune::TypeTree::HybridTreePath<S...> const& tp0, + Dune::TypeTree::HybridTreePath<T...> const& tp1) + { + return Dune::TypeTree::HybridTreePath<S...,T...>(std::tuple_cat(tp0._data,tp1._data)); + } } // end namespace AMDiS diff --git a/test/TreeDataTest.cpp b/test/TreeDataTest.cpp index 466681a2f1922702590ff83867918c6952936304..8c65636cf777de1ad0cf340bf446302d8b8eb0c3 100644 --- a/test/TreeDataTest.cpp +++ b/test/TreeDataTest.cpp @@ -19,26 +19,26 @@ using namespace AMDiS; template <class Node> using NodeData = double; -template <class Tree> -bool operator==(TreeData<Tree,NodeData,false> const& t1, TreeData<Tree,NodeData,false> const& t2) +template <class Basis> +bool operator==(TreeData<Basis,NodeData,false> const& t1, TreeData<Basis,NodeData,false> const& t2) { - AMDIS_TEST(t1.tree() == t2.tree() && t1.tree() != nullptr); + AMDIS_TEST(t1.basis() == t2.basis() && t1.basis() != nullptr); bool same = true; - AMDiS::forEachNode_(*t1.tree(), [&](auto const& node, auto) { + AMDiS::forEachNode_(t1.basis()->localView().tree(), [&](auto const& node, auto) { same = same && (t1[node] == t2[node]); }); return same; } -template <class Tree> -bool operator==(TreeData<Tree,NodeData,true> const& t1, TreeData<Tree,NodeData,true> const& t2) +template <class Basis> +bool operator==(TreeData<Basis,NodeData,true> const& t1, TreeData<Basis,NodeData,true> const& t2) { - AMDIS_TEST(t1.tree() == t2.tree() && t1.tree() != nullptr); + AMDIS_TEST(t1.basis() == t2.basis() && t1.basis() != nullptr); bool same = true; - AMDiS::forEachLeafNode_(*t1.tree(), [&](auto const& node, auto) { + AMDiS::forEachLeafNode_(t1.basis()->localView().tree(), [&](auto const& node, auto) { same = same && (t1[node] == t2[node]); }); @@ -56,39 +56,39 @@ int main () YaspGrid<2> grid(L, s); auto gridView = grid.leafGridView(); - auto taylorHoodBasis = makeBasis( + auto basis = makeBasis( gridView, composite( power<2>(lagrange<2>()), lagrange<1>() )); - auto localView = taylorHoodBasis.localView(); + auto localView = basis.localView(); auto const& tree = localView.tree(); - using Tree = std::remove_reference_t<decltype(tree)>; + using Basis = std::remove_reference_t<decltype(basis)>; // test treeData for all nodes { // call default-constructor - TreeData<Tree, NodeData, false> treeData; - treeData.init(tree); + TreeData<Basis, NodeData, false> treeData; + treeData.init(basis); AMDiS::forEachNode_(tree, [&](auto const& node, auto) { treeData[node] = double(node.treeIndex()); }); // call constructor with tree - TreeData<Tree, NodeData, false> treeData1(tree); + TreeData<Basis, NodeData, false> treeData1(basis); // call init on non-empty treeData - treeData1.init(tree); + treeData1.init(basis); // call copy-constructor - TreeData<Tree, NodeData, false> treeData2(treeData); + TreeData<Basis, NodeData, false> treeData2(treeData); AMDIS_TEST(treeData == treeData2); // call copy-assignment operator on empty treeData - TreeData<Tree, NodeData, false> treeData3; + TreeData<Basis, NodeData, false> treeData3; treeData3 = treeData; AMDIS_TEST(treeData == treeData3); @@ -100,31 +100,31 @@ int main () treeData = std::move(treeData2); // call move-constructor - TreeData<Tree, NodeData, false> treeData4(std::move(treeData3)); + TreeData<Basis, NodeData, false> treeData4(std::move(treeData3)); } // test treeData for leaf only { // call default-constructor - TreeData<Tree, NodeData, true> treeData; - treeData.init(tree); + TreeData<Basis, NodeData, true> treeData; + treeData.init(basis); AMDiS::forEachLeafNode_(tree, [&](auto const& node, auto) { treeData[node] = double(node.treeIndex()); }); // call constructor with tree - TreeData<Tree, NodeData, true> treeData1(tree); + TreeData<Basis, NodeData, true> treeData1(basis); // call init on non-empty treeData - treeData1.init(tree); + treeData1.init(basis); // call copy-constructor - TreeData<Tree, NodeData, true> treeData2(treeData); + TreeData<Basis, NodeData, true> treeData2(treeData); AMDIS_TEST(treeData == treeData2); // call copy-assignment operator on empty treeData - TreeData<Tree, NodeData, true> treeData3; + TreeData<Basis, NodeData, true> treeData3; treeData3 = treeData; AMDIS_TEST(treeData == treeData3); @@ -136,19 +136,19 @@ int main () treeData = std::move(treeData2); // call move-constructor - TreeData<Tree, NodeData, true> treeData4(std::move(treeData3)); + TreeData<Basis, NodeData, true> treeData4(std::move(treeData3)); } // test for operations with uninitialized tree { // call default-constructor without initialization - TreeData<Tree, NodeData, true> treeData; + TreeData<Basis, NodeData, true> treeData; // call copy-constructor - TreeData<Tree, NodeData, true> treeData2(treeData); + TreeData<Basis, NodeData, true> treeData2(treeData); // call move-constructor - TreeData<Tree, NodeData, true> treeData3(std::move(treeData)); + TreeData<Basis, NodeData, true> treeData3(std::move(treeData)); } return report_errors();