Commit ddee16ab authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

cleaned up assembling and dirichlet BC

parent 0f94f332
Pipeline #1300 passed with stage
in 21 minutes and 10 seconds
......@@ -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
......
......@@ -25,7 +25,6 @@ install(FILES
CreatorInterface.hpp
CreatorMap.hpp
DirichletBC.hpp
DirichletBC.inc.hpp
FileWriter.hpp
Flag.hpp
GridFunctionOperator.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"
#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
......@@ -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_;
......
......@@ -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);
});
});
......
install(FILES
Common.hpp
DOFMatrixBase.hpp
DOFMatrixBase.inc.hpp
DOFVectorBase.hpp
DOFVectorBase.inc.hpp
DOFVectorInterface.hpp
HierarchicWrapper.hpp
LinearSolver.hpp
......
......@@ -7,7 +7,7 @@ namespace AMDiS
template <class T>
struct Triplet
{
std::size_t row, cols;
std::size_t row, col;
T value;
};
......
......@@ -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"
#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 <