Commit 681e70ef authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Constraints improvements

parent 34c64903
#pragma once #pragma once
#include <list> #include <type_traits>
#include <memory>
#include <amdis/typetree/TreeContainer.hpp> #include <dune/functions/common/typeerasure.hh>
namespace AMDiS #include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/TypeTraits.hpp>
namespace AMDiS {
namespace Impl {
template <class Mat, class Sol, class Rhs>
struct BoundaryConditionDefinition
{ {
/// \brief Interface class for boundary conditions // Definition of the interface a BoundaryCondition must fulfill
/** struct Interface
* 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: virtual ~Interface() = default;
/// Default constructor.
BoundaryCondition() = default;
/// \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; virtual void init() = 0;
virtual void apply(Mat&, Sol&, Rhs&) = 0;
};
/// \brief Apply the boundary condition to matrix and vector. // Templatized implementation of the interface
/** template <class Impl>
* This is called after the matrix `A` and rhs `b` are assembled. struct Model : public Impl
* Implementations may alter values or change the structure of `A`, `x` or `b`. {
**/ using Impl::Impl;
virtual void apply(Mat& A, Sol& x, Rhs& b) = 0; void init() final { this->get().init(); }
void apply(Mat& A, Sol& x, Rhs& b) final { this->get().apply(A,x,b); }
}; };
template <class Mat, class Sol, class Rhs> // The Concept definition of a BoundaryCondition
struct BCData struct Concept
{ {
template <class RN, class CN> template <class BC>
using type = std::list<std::shared_ptr<BoundaryCondition<Mat, Sol, Rhs>>>; auto require(BC&& bc) -> decltype
(
bc.init(),
bc.apply(std::declval<Mat&>(), std::declval<Sol&>(), std::declval<Rhs&>())
);
}; };
/// Container for storing shared pointers to boundary conditions indexed by their position in using Base = Dune::Functions::TypeErasureBase<Interface, Model>;
/// the basis tree, see \ref MatrixData. };
template <class Mat, class Sol, class Rhs, class RB, class CB>
using BoundaryConditions } // end namespace Impl
= TreeMatrix<BCData<Mat, Sol, Rhs>::template type,
typename RB::LocalView::Tree, typename CB::LocalView::Tree>;
/// \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 Impl::BoundaryConditionDefinition<Mat,Sol,Rhs>::Base
{
using Definition = Impl::BoundaryConditionDefinition<Mat,Sol,Rhs>;
using Super = typename Definition::Base;
public:
/// \brief Constructor. Pass any type supporting the \ref BoundaryConditionInterface
template <class Impl, Dune::disableCopyMove<BoundaryCondition,Impl> = 0>
BoundaryCondition(Impl&& impl)
: Super{FWD(impl)}
{
static_assert(Concepts::models<typename Definition::Concept(Impl)>,
"Implementation does not model the BoundaryCondition concept.");
}
/// \brief Default Constructor.
BoundaryCondition() = default;
/// \brief Initialize the boundary condition.
/**
* This performs setup before starting the matrix and vector assembly independantly of the
* matrix, rhs or solution.
**/
void init()
{
this->asInterface().init();
}
/// \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`.
**/
void apply(Mat& A, Sol& x, Rhs& b)
{
this->asInterface().apply(A,x,b);
}
};
} // end namespace AMDiS } // end namespace AMDiS
...@@ -32,22 +32,14 @@ namespace AMDiS ...@@ -32,22 +32,14 @@ namespace AMDiS
* application of boundary conditions can be symmetric if the matrix does * application of boundary conditions can be symmetric if the matrix does
* support this symmetric modification. * support this symmetric modification.
* *
* \tparam Mat Matrix * \tparam SubBasis SubspaceBasis of the solution FE space
* \tparam Sol Vector of solution
* \tparam Rhs Vector of rhs
* \tparam RowSubBasis SubspaceBasis of the row FE space
* \tparam ColSubBasis SubspaceBasis of the column FE space
* \tparam ValueGridFct Type of the GridFunction representing the Dirichlet values. * \tparam ValueGridFct Type of the GridFunction representing the Dirichlet values.
**/ **/
template <class Mat, class Sol, class Rhs, class RowSubBasis, class ColSubBasis, template <class SubBasis, class ValueGridFct>
class ValueGridFct>
class DirichletBC class DirichletBC
: public BoundaryCondition<Mat, Sol, Rhs>
{ {
using Super = BoundaryCondition<Mat, Sol, Rhs>;
// We assume CB's GridView to be the same as RB's // We assume CB's GridView to be the same as RB's
using GridView = typename RowSubBasis::GridView; using GridView = typename SubBasis::GridView;
using Intersection = typename GridView::Intersection; using Intersection = typename GridView::Intersection;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate; using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
...@@ -57,78 +49,72 @@ namespace AMDiS ...@@ -57,78 +49,72 @@ namespace AMDiS
/// Constructor accepting subspacebases. /// Constructor accepting subspacebases.
template <class Values, template <class Values,
REQUIRES(Concepts::Functor<Values, Range(Domain)>)> REQUIRES(Concepts::Functor<Values, Range(Domain)>)>
DirichletBC(RowSubBasis rowBasis, ColSubBasis colBasis, DirichletBC(SubBasis basis, BoundarySubset<Intersection> boundarySubset, Values&& values)
BoundarySubset<Intersection> boundarySubset, Values&& values) : basis_(std::move(basis))
: rowBasis_(std::move(rowBasis))
, colBasis_(std::move(colBasis))
, boundarySubset_(std::move(boundarySubset)) , boundarySubset_(std::move(boundarySubset))
, values_(FWD(values)) , values_(FWD(values))
{} {}
/// Make a DirichletBC from a basis with treepath arguments
template <class Basis, class TP, class Values,
REQUIRES(Concepts::GlobalBasis<Basis>)>
DirichletBC(Basis const& basis, TP const& treePath,
BoundarySubset<Intersection> boundarySubset, Values&& values)
: DirichletBC(Dune::Functions::subspaceBasis(basis, treePath),
std::move(boundarySubset), FWD(values))
{}
/// Make a DirichletBC from a global basis
template <class Basis, class Values,
REQUIRES(Concepts::GlobalBasis<Basis>)>
DirichletBC(Basis const& basis, BoundarySubset<Intersection> boundarySubset, Values&& values)
: DirichletBC(Dune::Functions::subspaceBasis(basis, makeTreePath()),
std::move(boundarySubset), FWD(values))
{}
/// Fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the /// Fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the
/// boundary or not. /// boundary or not.
/** /**
* \see BoundaryCondition::init * \see BoundaryCondition::init
**/ **/
void init() override; void init();
/// \brief Apply dirichlet BC to matrix and vector /// \brief Apply dirichlet BC to matrix and vector
/** /**
* Add a unit-row to the matrix and optionally delete the corresponding matrix-column. * Add a unit-row to the matrix and optionally delete the corresponding matrix-column.
* Uses a backend-specific implementation. * Uses a backend-specific implementation.
**/ **/
void apply(Mat& matrix, Sol& solution, Rhs& rhs) override; template <class Mat, class Sol, class Rhs>
void apply(Mat& matrix, Sol& solution, Rhs& rhs);
private: private:
RowSubBasis rowBasis_; SubBasis basis_;
ColSubBasis colBasis_;
BoundarySubset<Intersection> boundarySubset_; BoundarySubset<Intersection> boundarySubset_;
ValueGridFct values_; ValueGridFct values_;
std::vector<bool> dirichletNodes_; std::vector<bool> dirichletNodes_;
}; };
// deduction guides
/// Make a DirichletBC from subspacebases
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP, class Values> // Make a DirichletBC from subspacebases
auto makeDirichletBC(Dune::Functions::SubspaceBasis<RB, RTP> rowBasis, template <class B, class TP, class Values>
Dune::Functions::SubspaceBasis<CB, CTP> colBasis, DirichletBC(Dune::Functions::SubspaceBasis<B, TP>,
BoundarySubset<typename RB::GridView::Intersection> boundarySubset, BoundarySubset<typename B::GridView::Intersection>, Values const&)
Values&& values) -> DirichletBC<Dune::Functions::SubspaceBasis<B,TP>, Underlying_t<Values>>;
{
using RowSubBasis = Dune::Functions::SubspaceBasis<RB, RTP>; // Make a DirichletBC from a basis with treepath arguments
using ColSubBasis = Dune::Functions::SubspaceBasis<CB, CTP>; template <class B, class TP, class Values,
using BcType = DirichletBC<Mat, Sol, Rhs, RowSubBasis, ColSubBasis, TYPEOF(values)>; REQUIRES(Concepts::GlobalBasis<B>)>
return BcType(std::move(rowBasis), std::move(colBasis), std::move(boundarySubset), FWD(values)); DirichletBC(B const&, TP const&, BoundarySubset<typename B::GridView::Intersection>, Values const&)
} -> DirichletBC<Dune::Functions::SubspaceBasis<B, TP>, Underlying_t<Values>>;
/// Make a DirichletBC from a row- and colbasis with treepath arguments // Make a DirichletBC from a global basis
template <class Mat, class Sol, class Rhs, class RB, class RTP, class CB, class CTP, class Values, template <class B, class Values,
REQUIRES(Concepts::GlobalBasis<RB>), REQUIRES(Concepts::GlobalBasis<B>)>
REQUIRES(Concepts::GlobalBasis<CB>)> DirichletBC(B const&, BoundarySubset<typename B::GridView::Intersection>, Values const&)
auto makeDirichletBC(RB const& rowBasis, RTP const& rowTreePath, -> DirichletBC<Dune::Functions::SubspaceBasis<B,TYPEOF(makeTreePath())>, Underlying_t<Values>>;
CB const& colBasis, CTP const& colTreePath,
BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
Values&& values)
{
return makeDirichletBC<Mat,Sol,Rhs>(
Dune::Functions::subspaceBasis(rowBasis, rowTreePath),
Dune::Functions::subspaceBasis(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(Concepts::GlobalBasis<RB>),
REQUIRES(Concepts::GlobalBasis<CB>)>
auto makeDirichletBC(RB const& rowBasis, CB const& colBasis,
BoundarySubset<typename RB::GridView::Intersection> boundarySubset,
Values&& values)
{
return makeDirichletBC<Mat, Sol, Rhs>(
rowBasis, makeTreePath(), colBasis, makeTreePath(),
std::move(boundarySubset), FWD(values));
}
} // end namespace AMDiS } // end namespace AMDiS
......
...@@ -8,50 +8,36 @@ ...@@ -8,50 +8,36 @@
#include <amdis/LinearAlgebra.hpp> #include <amdis/LinearAlgebra.hpp>
namespace AMDiS { namespace AMDiS {
namespace Impl {
/// Traverse all interior boundary DOFs and apply the functor f. template <class B, class V>
/// The Functor must have the signature void DirichletBC<B, V>::
/// `void(int, typename Basis::LocalView, typename Basis::GridView::Intersection)` init()
template <class Basis, class F>
void forEachInteriorBoundaryDOF(Basis const& basis, F&& f)
{ {
auto localView = basis.localView(); dirichletNodes_.assign(basis_.dimension(), false);
auto seDOFs = Dune::Functions::subEntityDOFs(basis);
auto const& gridView = basis.gridView(); auto localView = basis_.localView();
for (auto const& element : elements(gridView, typename BackendTraits<Basis>::PartitionSet{})) { auto seDOFs = Dune::Functions::subEntityDOFs(basis_);
auto const& gridView = basis_.gridView();
for (auto const& element : elements(gridView, typename BackendTraits<B>::PartitionSet{})) {
if (element.hasBoundaryIntersections()) { if (element.hasBoundaryIntersections()) {
localView.bind(element); localView.bind(element);
for(auto const& intersection: intersections(gridView, element)) for(auto const& intersection: intersections(gridView, element))
if (intersection.boundary()) if (intersection.boundary())
for(auto localIndex: seDOFs.bind(localView,intersection)) for(auto localIndex: seDOFs.bind(localView,intersection))
f(localIndex, localView, intersection); dirichletNodes_[localView.index(localIndex)]
= dirichletNodes_[localView.index(localIndex)] || boundarySubset_(intersection);
} }
} }
}
} // end namespace Impl
template <class Mat, class Sol, class Rhs, class RB, class CB, class V>
void DirichletBC<Mat, Sol, Rhs, RB, CB, V>::
init()
{
using LV = typename CB::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 Mat, class Sol, class Rhs, class RB, class CB, class V> template <class B, class V>
void DirichletBC<Mat, Sol, Rhs, RB, CB, V>:: template <class Mat, class Sol, class Rhs>
void DirichletBC<B, V>::
apply(Mat& matrix, Sol& solution, Rhs& rhs) apply(Mat& matrix, Sol& solution, Rhs& rhs)
{ {
AMDiS::interpolate(colBasis_, solution, values_, tag::defaulted{}, tag::defaulted{}, dirichletNodes_); AMDiS::interpolate(basis_, solution, values_, tag::defaulted{}, tag::defaulted{}, dirichletNodes_);
dirichletBC(matrix, solution, rhs, dirichletNodes_); dirichletBC(matrix, solution, rhs, dirichletNodes_);
} }
......
...@@ -57,11 +57,9 @@ namespace AMDiS ...@@ -57,11 +57,9 @@ namespace AMDiS
* \tparam Basis RowBasis, must be the same as ColBasis * \tparam Basis RowBasis, must be the same as ColBasis
* \tparam TP Treepath to the node this constraint applies to * \tparam TP Treepath to the node this constraint applies to
**/ **/
template <class Mat, class Sol, class Rhs, class Basis, class TP> template <class Basis, class TP>
class PeriodicBC class PeriodicBC
: public BoundaryCondition<Mat, Sol, Rhs>
{ {
using Super = BoundaryCondition<Mat, Sol, Rhs>;
using SubspaceBasis = Dune::Functions::SubspaceBasis<Basis, TP>; using SubspaceBasis = Dune::Functions::SubspaceBasis<Basis, TP>;
public: public:
...@@ -85,7 +83,7 @@ namespace AMDiS ...@@ -85,7 +83,7 @@ namespace AMDiS
* true for all DOFs on the periodic boundary. Furthermode associations() is initialized and * 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. * contains the index of the DOF associated to a DOF on the periodic boundary.
**/ **/
void init() override; void init();
/// Move matrix rows (and columns) of dependant DOFs to the corresponding master DOFs /// Move matrix rows (and columns) of dependant DOFs to the corresponding master DOFs
/** /**
...@@ -93,7 +91,8 @@ namespace AMDiS ...@@ -93,7 +91,8 @@ namespace AMDiS
* columns) of dependant DOFs to unit-rows[columns] with entries of -1 at master DOF indices. * columns) of dependant DOFs to unit-rows[columns] with entries of -1 at master DOF indices.
* Uses a backend-specific implementation. * Uses a backend-specific implementation.
**/ **/
void apply(Mat& matrix, Sol& solution, Rhs& rhs) override; template <class Mat, class Sol, class Rhs>
void apply(Mat& matrix, Sol& solution, Rhs& rhs);
/// Return the map of DOF associations /// Return the map of DOF associations
auto const& associations() const auto const& associations() const
...@@ -134,35 +133,35 @@ namespace AMDiS ...@@ -134,35 +133,35 @@ namespace AMDiS
/// Make a PeriodicBC from a subspacebasis /// Make a PeriodicBC from a subspacebasis
template <class Mat, class Sol, class Rhs, class B, class TP> template <class B, class TP>
auto makePeriodicBC(Dune::Functions::SubspaceBasis<B, TP> basis, auto makePeriodicBC(Dune::Functions::SubspaceBasis<B, TP> basis,
BoundarySubset<typename B::GridView::Intersection> boundarySubset, BoundarySubset<typename B::GridView::Intersection> boundarySubset,
Impl::FaceTrafo<B> trafo) Impl::FaceTrafo<B> trafo)
{ {
using BcType = PeriodicBC<Mat, Sol, Rhs, B, TP>; using BcType = PeriodicBC<B, TP>;
return BcType(std::move(basis), std::move(boundarySubset), std::move(trafo)); return BcType(std::move(basis), std::move(boundarySubset), std::move(trafo));
} }
/// Make a PeriodicBC from a global basis and treepath /// Make a PeriodicBC from a global basis and treepath
template <class Mat, class Sol, class Rhs, class B, class TP, template <class B, class TP,
REQUIRES(Concepts::GlobalBasis<B>)> REQUIRES(Concepts::GlobalBasis<B>)>
auto makePeriodicBC(B const& basis, TP const& tp, auto makePeriodicBC(B const& basis, TP const& tp,
BoundarySubset<typename B::GridView::Intersection> boundarySubset, BoundarySubset<typename B::GridView::Intersection> boundarySubset,
Impl::FaceTrafo<B> trafo) Impl::FaceTrafo<B> trafo)
{ {
return makePeriodicBC<Mat, Sol, Rhs>(Dune::Functions::subspaceBasis(basis, tp), return makePeriodicBC(Dune::Functions::subspaceBasis(basis, tp),
std::move(boundarySubset), std::move(trafo)); std::move(boundarySubset), std::move(trafo));
} }
/// Make a PeriodicBC from a global basis /// Make a PeriodicBC from a global basis
template <class Mat, class Sol, class Rhs, class B, template <class B,
REQUIRES(Concepts::GlobalBasis<B>)> REQUIRES(Concepts::GlobalBasis<B>)>
auto makePeriodicBC(B const& basis, auto makePeriodicBC(B const& basis,
BoundarySubset<typename B::GridView::Intersection> boundarySubset, BoundarySubset<typename B::GridView::Intersection> boundarySubset,
Impl::FaceTrafo<B> trafo) Impl::FaceTrafo<B> trafo)
{ {
return makePeriodicBC<Mat, Sol, Rhs>(basis, makeTreePath(), std::move(boundarySubset), return makePeriodicBC(basis, makeTreePath(), std::move(boundarySubset),
std::move(trafo)); std::move(trafo));
} }
} // end namespace AMDiS } // end namespace AMDiS
......
...@@ -14,8 +14,8 @@ ...@@ -14,8 +14,8 @@
namespace AMDiS { namespace AMDiS {
template <class Mat, class Sol, class Rhs, class Basis, class TP> template <class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>:: void PeriodicBC<Basis, TP>::
init() init()
{ {
if (!bool(hasConnectedIntersections_)) { if (!bool(hasConnectedIntersections_)) {
...@@ -41,8 +41,8 @@ init() ...@@ -41,8 +41,8 @@ init()
} }
template <class Mat, class Sol, class Rhs, class Basis, class TP> template <class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>:: void PeriodicBC<Basis, TP>::
initAssociations() initAssociations()
{ {
using std::sqrt; using std::sqrt;
...@@ -124,8 +124,8 @@ namespace Impl ...@@ -124,8 +124,8 @@ namespace Impl
} }
template <class Mat, class Sol, class Rhs, class Basis, class TP> template <class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>:: void PeriodicBC<Basis, TP>::
initAssociations2() initAssociations2()
{ {
using std::sqrt; using std::sqrt;
...@@ -191,9 +191,9 @@ initAssociations2() ...@@ -191,9 +191,9 @@ initAssociations2()
} }
template <class Mat, class Sol, class Rhs, class Basis, class TP> template <class Basis, class TP>
template <class Node> template <class Node>
auto PeriodicBC<Mat, Sol, Rhs, Basis, TP>:: auto PeriodicBC<Basis, TP>::
coords(Node const& tree, std::vector<std::size_t> const& localIndices) const coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
{ {
std::vector<Domain> dofCoords(localIndices.size()); std::vector<Domain> dofCoords(localIndices.size());
...@@ -231,8 +231,9 @@ coords(Node const& tree, std::vector<std::size_t> const& localIndices) const ...@@ -231,8 +231,9 @@ coords(Node const& tree, std::vector<std::size_t> const& localIndices) const
} }
template <class Mat, class Sol, class Rhs, class Basis, class TP> template <class Basis, class TP>
void PeriodicBC<Mat, Sol, Rhs, Basis, TP>:: template <class Mat, class Sol, class Rhs>
void PeriodicBC<Basis, TP>::
apply(Mat& matrix, Sol& solution, Rhs& rhs) apply(Mat& matrix, Sol& solution, Rhs& rhs)
{ {
periodicBC(matrix, solution, rhs, periodicNodes_, associations_); periodicBC(matrix, solution, rhs, periodicNodes_, associations_);
......
...@@ -305,6 +305,10 @@ namespace AMDiS ...@@ -305,6 +305,10 @@ namespace AMDiS
void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b); void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b);
/** @} */ /** @} */
void addConstraint(BoundaryCondition<SystemMatrix, SolutionVector, SystemVector> constraint)
{
constraints_.push_back(constraint);
}