Commit 3058a8c5 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Cleanup of boundary conditions

Add documentation to boundary conditions
Correct DirichletBC::fillBC to use column basis for solution
parent 6420faa7
#pragma once
#include <list>
#include <memory>
#include <amdis/Boundary.hpp>
#include <amdis/BoundaryManager.hpp>
#include <amdis/typetree/TreeData.hpp>
namespace AMDiS
{
/// \brief Interface class for boundary conditions
/**
* Stores a boundary subset related to the boundary condition. See \ref BoundarySubset.
* \relates DirichletBoundaryCondition
* \relates PeriodicBoundaryCondition
* \tparam Mat Matrix
* \tparam Sol Vector of solution
* \tparam Rhs Vector of rhs
**/
template <class Mat, class Sol, class Rhs>
class BoundaryCondition
{
public:
/// Default constructor.
BoundaryCondition() = default;
BoundaryCondition(std::shared_ptr<BoundaryManagerBase const> const& boundaryManager, BoundaryType id)
: boundaryManager_(boundaryManager)
, id_(id)
{}
/// Return true if intersection is on boundary with id
template <class Intersection>
bool onBoundary(Intersection const& is) const
{
return is.boundary() && (!boundaryManager_ || boundaryManager_->boundaryId(is) == id_);
}
template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis) { /* do nothing */ }
template <class Matrix, class X, class B, class RN, class RTP, class CN, class CTP>
void fillBoundaryCondition(Matrix& A, X& x, B& b, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath) { /* do nothing */ }
protected:
std::shared_ptr<BoundaryManagerBase const> boundaryManager_{nullptr};
BoundaryType id_{0};
/// \brief Initialize the boundary condition.
/**
* This performs setup before starting the matrix and vector assembly independantly of the
* matrix, rhs or solution.
**/
virtual void init() = 0;
/// \brief Apply the boundary condition to matrix and vector.
/**
* This is called after the matrix `A` and rhs `b` are assembled.
* Implementations may alter values or change the structure of `A`, `x` or `b`.
**/
virtual void apply(Mat& A, Sol& x, Rhs& b) = 0;
};
template <class Mat, class Sol, class Rhs>
struct BCData
{
template <class RN, class CN>
using type = std::list<std::shared_ptr<BoundaryCondition<Mat, Sol, Rhs>>>;
};
/// Container for storing shared pointers to boundary conditions indexed by their position in
/// the basis tree, see \ref MatrixData.
template <class Mat, class Sol, class Rhs, class RB, class CB>
using BoundaryConditions = MatrixData<RB, CB, BCData<Mat, Sol, Rhs>::template type>;
} // end namespace AMDiS
#pragma once
#include <functional>
#include <amdis/Boundary.hpp>
#include <amdis/BoundaryManager.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ConceptsBase.hpp>
namespace AMDiS
{
/// \brief Class defining a subset of a domain boundary.
/**
* Stores a predicate identifying boundary segments. This may be given as
* - a boundary manager and a boundary ID, see \ref BoundaryManager,
* - a functor of the form bool(GlobalCoordinate)
* - no argument in which case the whole boundary is used.
*
* \tparam IS Type of the intersection of the elements with the boundary
**/
template <class IS>
class BoundarySubset
{
using Domain = typename IS::GlobalCoordinate;
public:
using Intersection = IS;
/// Default constructor. Uses a predicate that returns true on the complete boundary.
BoundarySubset()
: predicate_([](Intersection const& is) -> bool { return is.boundary(); })
{}
/// Use a boundary manager and id to determine a subset.
BoundarySubset(BoundaryManagerBase& boundaryManager, BoundaryType id)
: predicate_([&boundaryManager, id](Intersection const& is) -> bool {
return is.boundary() && boundaryManager.boundaryId(is) == id;
})
{}
/// Use a predicate of the form bool(GlobalCoordinate) to determine a subset.
template <class Predicate,
REQUIRES(Concepts::Functor<Predicate, bool(Domain)>)>
BoundarySubset(Predicate&& predicate)
: predicate_([predicate](Intersection const& is) -> bool {
return predicate(is.geometry().center());
})
{}
/// Return true if intersection is on boundary segment
bool operator()(Intersection const& is) const
{
return predicate_(is);
}
protected:
std::shared_ptr<BoundaryManagerBase const> boundaryManager_;
BoundaryType id_{0};
std::function<bool(Intersection const&)> predicate_;
};
} // end namespace AMDiS
......@@ -29,6 +29,7 @@ install(FILES
Boundary.hpp
BoundaryCondition.hpp
BoundaryManager.hpp
BoundarySubset.hpp
ContextGeometry.hpp
CreatorInterface.hpp
CreatorMap.hpp
......
#pragma once
#include <functional>
#include <list>
#include <memory>
#include <utility>
#include <vector>
#include <dune/common/std/optional.hh>
#include <dune/common/concept.hh>
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <amdis/Boundary.hpp>
#include <amdis/BoundaryCondition.hpp>
#include <amdis/BoundarySubset.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/RangeType.hpp>
#include <amdis/typetree/Traversal.hpp>
#include <amdis/typetree/TreeData.hpp>
#include <amdis/typetree/TreePath.hpp>
namespace AMDiS
{
/// Implements a boundary condition of Dirichlet-type.
/**
* By calling the methods \ref init() and \ref finish before and after
* By calling the methods \ref init and \ref finish before and after
* assembling the system-matrix, respectively, dirichlet boundary conditions
* can be applied to the matrix and system vector. Therefore, a predicate
* functions indicates the DOFs where values should be enforced and a second
* functor provided in the constructor is responsible for determining the
* values to be set at the DOFs.
*
* In the \ref finish method the matrix is called with \ref applyDirichletBC
* In the \ref finish method the matrix is called with \ref apply
* to erase the corresponding rows and columns for the DOF indices. This
* application of boundary conditions can be symmetric if the matrix does
* support this symmetric modification. As a result, this method returns a list
* of columns values, that should be subtracted from the rhs.
* support this symmetric modification.
*
* \tparam Mat Matrix
* \tparam Sol Vector of solution
* \tparam Rhs Vector of rhs
* \tparam RB Basis of the row FE space
* \tparam RTP Treepath to the row node this constraint applies to
* \tparam CB Basis of the column FE space
* \tparam CTP Treepath to the column node this constraint applies to
**/
template <class Domain, class Range = double>
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP>
class DirichletBC
: public BoundaryCondition
: public BoundaryCondition<Mat, Sol, Rhs>
{
using Super = BoundaryCondition;
using Super = BoundaryCondition<Mat, Sol, Rhs>;
using RowSubBasis = Dune::Functions::SubspaceBasis<RB, RTP>;
using ColSubBasis = Dune::Functions::SubspaceBasis<CB, CTP>;
// We assume CB's GridView to be the same as RB's
using Domain = typename RB::GridView::template Codim<0>::Geometry::GlobalCoordinate;
using Intersection = typename RB::GridView::Intersection;
using Range = RangeType_t<typename ColSubBasis::LocalView::Tree>;
public:
template <class BM, class Values,
REQUIRES(Concepts::Functor<Values, Range(Domain)>) >
DirichletBC(BM&& boundaryManager, BoundaryType id, Values&& values)
: Super(FWD(boundaryManager), id)
/// Constructor accepting subspacebases.
template <class Values,
REQUIRES(Concepts::Functor<Values, Range(Domain)>)>
DirichletBC(RowSubBasis rowBasis, ColSubBasis colBasis,
BoundarySubset<Intersection> boundarySubset, Values&& values)
: rowBasis_(std::move(rowBasis))
, colBasis_(std::move(colBasis))
, boundarySubset_(std::move(boundarySubset))
, values_(FWD(values))
{}
template <class Predicate, class Values,
REQUIRES(Concepts::Functor<Predicate, bool(Domain)>),
/// Constructor accepting global bases and treepaths.
template <class Values,
REQUIRES(Concepts::Functor<Values, Range(Domain)>)>
DirichletBC(Predicate&& predicate, Values&& values)
: predicate_(FWD(predicate))
, values_(FWD(values))
DirichletBC(RB const& rowBasis, RTP const& rowTreePath,
CB const& colBasis, CTP const& colTreePath,
BoundarySubset<Intersection> boundarySubset, Values&& values)
: DirichletBC(Dune::Functions::subspaceBasis(rowBasis, rowTreePath),
Dune::Functions::subspaceBasis(colBasis, colTreePath),
std::move(boundarySubset), FWD(values))
{}
template <class IS>
bool onBoundary(IS const& intersection) const
{
return predicate_ ? (*predicate_)(intersection.geometry().center())
: Super::onBoundary(intersection);
}
// fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not.
template <class RB, class CB>
void init(RB const& rowBasis, CB const& colBasis);
/**
* \see BoundaryCondition::init
**/
void init() override;
/// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
/// to the matrix and optionally delete the corresponding matrix-column.
/// \brief Apply dirichlet BC to matrix and vector
/**
* \tparam Mat Matrix
* \tparam Sol Vector type of solution
* \tparam Rhs Vector type of rhs
* \tparam RN RowNode
* \tparam CN ColNode
* Add a unit-row to the matrix and optionally delete the corresponding matrix-column.
* Uses a backend-specific implementation.
**/
template <class Mat, class Sol, class Rhs, class RN, class RTP, class CN, class CTP>
void fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath);
void apply(Mat& matrix, Sol& solution, Rhs& rhs) override;
private:
Dune::Std::optional<std::function<bool(Domain)>> predicate_;
RowSubBasis rowBasis_;
ColSubBasis colBasis_;
BoundarySubset<Intersection> boundarySubset_;
std::function<Range(Domain)> values_;
std::vector<bool> dirichletNodes_;
};
template <class Basis>
struct DirichletData
/// Make a DirichletBC from subspacebases
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP, class Values>
auto makeDirichletBC(Dune::Functions::SubspaceBasis<RB, RTP> rowBasis,
Dune::Functions::SubspaceBasis<CB, CTP> colBasis,
BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
Values&& values)
{
using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
template <class RN, class CN>
using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RN>>> >;
};
template <class RowBasis, class ColBasis>
using DirichletBCs = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>;
using BcType = DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>;
return BcType(std::move(rowBasis), std::move(colBasis), std::move(boundarySubset), FWD(values));
}
/// Make a DirichletBC from a row- and colbasis with treepath arguments
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP, class Values,
REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename RB::GridView>, RB>()),
REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename CB::GridView>, CB>())>
auto makeDirichletBC(RB const& rowBasis, RTP const& rowTreePath,
CB const& colBasis, CTP const& colTreePath,
BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
Values&& values)
{
using BcType = DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>;
return BcType(rowBasis, rowTreePath, colBasis, colTreePath,
std::move(boundarySubset), FWD(values));
}
/// Make a DirichletBC from a global row- and colbasis
template <class Mat, class Sol, class Rhs, class RB, class CB, class Values,
REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename RB::GridView>, RB>()),
REQUIRES(Dune::models<Dune::Functions::Concept::GlobalBasis<typename CB::GridView>, CB>())>
auto makeDirichletBC(RB const& rowBasis, CB const& colBasis,
BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
Values&& values)
{
return makeDirichletBC<Mat, Sol, Rhs>(rowBasis, treepath(), colBasis, treepath(),
std::move(boundarySubset), FWD(values));
}
} // end namespace AMDiS
......
......@@ -2,18 +2,17 @@
#include <type_traits>
#include <dune/common/hybridutilities.hh>
#include <dune/functions/functionspacebases/subentitydofs.hh>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/functions/Interpolate.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
#include <amdis/LinearAlgebra.hpp>
namespace AMDiS {
namespace Impl {
/// Traverse all interior boundary DOFs and apply the functor f.
/// The Functor must have the signature `void(int, typename Basis::LocalView, typename Basis::GridView::Intersection)`
/// The Functor must have the signature
/// `void(int, typename Basis::LocalView, typename Basis::GridView::Intersection)`
template <class Basis, class F>
void forEachInteriorBoundaryDOF(Basis const& basis, F&& f)
{
......@@ -34,31 +33,26 @@ void forEachInteriorBoundaryDOF(Basis const& basis, F&& f)
} // end namespace Impl
template <class D, class R>
template <class RB, class CB>
void DirichletBC<D,R>::
init(RB const& rowBasis, CB const& colBasis)
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP>
void DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>::
init()
{
using LV = typename CB::LocalView;
using IS = typename CB::GridView::Intersection;
dirichletNodes_.resize(colBasis.dimension());
dirichletNodes_.assign(dirichletNodes_.size(), false);
Impl::forEachInteriorBoundaryDOF(colBasis, [&](int localIndex, LV const& localView, IS const& intersection) {
dirichletNodes_[localView.index(localIndex)] = dirichletNodes_[localView.index(localIndex)] || onBoundary(intersection);
using LV = typename ColSubBasis::LocalView;
dirichletNodes_.assign(colBasis_.dimension(), false);
Impl::forEachInteriorBoundaryDOF(colBasis_,
[&](int localIndex, LV const& localView, Intersection const& intersection) {
dirichletNodes_[localView.index(localIndex)] = dirichletNodes_[localView.index(localIndex)]
|| boundarySubset_(intersection);
});
}
template <class D, class R>
template <class Mat, class Sol, class Rhs, class RN, class RTP, class CN, class CTP>
void DirichletBC<D,R>::
fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath)
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP>
void DirichletBC<Mat, Sol, Rhs, RB, RTP, CB, CTP>::
apply(Mat& matrix, Sol& solution, Rhs& rhs)
{
Dune::Hybrid::ifElse(std::is_same<RangeType_t<RN>, R>{},
[&](auto id) {
auto&& gf = makeGridFunction(values_, solution.basis()->gridView());
AMDiS::interpolate(*solution.basis(), id(solution), gf, rowTreePath, tag::defaulted{}, dirichletNodes_);
});
auto&& gf = makeGridFunction(values_, colBasis_.gridView());
AMDiS::interpolate(colBasis_, solution, gf, tag::defaulted{}, tag::defaulted{}, dirichletNodes_);
dirichletBC(matrix, solution, rhs, dirichletNodes_);
}
......
......@@ -3,7 +3,7 @@
#include <list>
#include <memory>
#include <amdis/BoundaryCondition.hpp>
#include <amdis/BoundarySubset.hpp>
#include <amdis/AssemblerInterface.hpp>
#include <amdis/typetree/TreeData.hpp>
......@@ -13,10 +13,10 @@ namespace AMDiS
{
template <class E> struct element_operator { using type = E; };
template <class I> struct intersection_operator { using type = I; };
template <class I> struct boundary_operator : public BoundaryCondition
template <class I> struct boundary_operator : public BoundarySubset<I>
{
using type = I;
using BoundaryCondition::BoundaryCondition;
using BoundarySubset<I>::BoundarySubset;
};
}
......@@ -30,7 +30,7 @@ namespace AMDiS
struct DataElement
{
std::shared_ptr<Op> op;
BoundaryCondition bc;
BoundarySubset<Intersection> bs;
};
/// Lists of \ref DataElement on the Element, BoundaryIntersction, and InteriorIntersections
......
#pragma once
#include <functional>
#include <list>
#include <map>
#include <type_traits>
#include <utility>
#include <vector>
#include <amdis/Output.hpp>
#include <dune/common/concept.hh>
#include <dune/common/std/optional.hh>
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <amdis/BoundaryCondition.hpp>
#include <amdis/BoundarySubset.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/TreePath.hpp>
namespace AMDiS
{
/// \brief Affine transformation x := A*x + b
template <class ct, int dow>
struct FaceTransformation
{
using WorldMatrix = FieldMatrix<ct, dow, dow>;
using WorldVector = FieldVector<ct, dow>;
WorldMatrix matrix_; // Note: We assume the matrix to be orthogonal, here
/// orthogonal (rotation) matrix
WorldMatrix matrix_;
/// Shift vector
WorldVector shift_;
WorldVector evaluate(WorldVector const& x) const
......@@ -41,55 +51,73 @@ namespace AMDiS
/// Implements a periodic boundary condition
/**
* \tparam D Domain
* \tparam MI Type of the MultiIndex
* \tparam Mat Matrix
* \tparam Sol Vector of solution
* \tparam Rhs Vector of rhs
* \tparam Basis RowBasis, must be the same as ColBasis
* \tparam TP Treepath to the node this constraint applies to
**/
template <class D, class MI>
template <class Mat, class Sol, class Rhs, class Basis, class TP>
class PeriodicBC
: public BoundaryCondition
: public BoundaryCondition<Mat, Sol, Rhs>
{
using Super = BoundaryCondition;
using Super = BoundaryCondition<Mat, Sol, Rhs>;
using SubspaceBasis = Dune::Functions::SubspaceBasis<Basis, TP>;
public:
using Domain = D;
using MultiIndex = MI;
using FaceTrafo = FaceTransformation<typename D::field_type, D::dimension>;
using Domain = typename SubspaceBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MultiIndex = typename SubspaceBasis::MultiIndex;
using FaceTrafo = FaceTransformation<typename Domain::field_type, Domain::dimension>;
using Intersection = typename SubspaceBasis::GridView::Intersection;
public:
template <class BM>
PeriodicBC(BM&& boundaryManager, BoundaryType id, FaceTrafo faceTrafo)
: Super(FWD(boundaryManager), id)
/// Create the BoundaryCondition and store the face transformation in a local variable.
PeriodicBC(SubspaceBasis basis, BoundarySubset<Intersection> boundarySubset,
FaceTrafo faceTrafo)
: basis_(std::move(basis))
, boundarySubset_(std::move(boundarySubset))
, faceTrafo_(std::move(faceTrafo))
{}
template <class RB, class CB>
void init(RB const& rowBasis, CB const& colBasis);
template <class Mat, class Sol, class Rhs, class RN, class RTP, class CN, class CTP>
void fillBoundaryCondition(Mat& A, Sol& x, Rhs& b, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath);
/// Compute the pairs of associated boundary DOFs
/**
* After leaving the function the vector used by periodicNodes() is initialized and its value is
* true for all DOFs on the periodic boundary. Furthermode associations() is initialized and
* contains the index of the DOF associated to a DOF on the periodic boundary.
**/
void init() override;
/// Move matrix rows (and columns) of dependant DOFs to the corresponding master DOFs
/**
* Add rows (and columns) of dependant DOFs to rows (and columns) master DOFs and set rows (and
* columns) of dependant DOFs to unit-rows[columns] with entries of -1 at master DOF indices.
* Uses a backend-specific implementation.
**/
void apply(Mat& matrix, Sol& solution, Rhs& rhs) override;
/// Return the map of DOF associations
auto const& associations() const
{
return associations_;
}
/// Return the boolean boundary indicator
auto const& periodicNodes() const
{
return periodicNodes_;
}
protected:
template <class Basis>
void initAssociations(Basis const& basis);
template <class Basis>
void initAssociations2(Basis const& basis);
void initAssociations();
void initAssociations2();
// Get the coordinates of the DOFs
template <class Node>
std::vector<Domain> coords(Node const& node, std::vector<std::size_t> const& localIndices) const;
auto coords(Node const& node, std::vector<std::size_t> const& localIndices) const;
private:
SubspaceBasis basis_;
BoundarySubset<Intersection> boundarySubset_;
FaceTrafo faceTrafo_;
std::vector<bool> periodicNodes_;
......@@ -97,18 +125,45 @@ namespace AMDiS
Dune::Std::optional<bool> hasConnectedIntersections_;
};
namespace Impl {
template <class Basis>
using FaceTrafo = FaceTransformation<
typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate::field_type,
Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate::dimension>;
}
template <class Basis>
struct PeriodicData
{
using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
template <class RowNode, class ColNode>
using type = std::list< std::shared_ptr<PeriodicBC<WorldVector, typename Basis::MultiIndex>> >;
};
template <class RowBasis, class ColBasis>
using PeriodicBCs = MatrixData<RowBasis, ColBasis, PeriodicData<RowBasis>::template type>;
/// Make a PeriodicBC from a subspacebasis
template <class Mat, class Sol, class Rhs, class B, class TP>
auto makePeriodicBC(Dune::Functions::SubspaceBasis<B, TP> basis,
BoundarySubset<typename B::GridView::Intersection> boundarySubset,
Impl::FaceTrafo<B> trafo)
{
using BcType = PeriodicBC<Mat, Sol, Rhs, B, TP>;