Commit 898cd2af authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Implement Operator and LocalOperator with TypeErasure

simplify ProblemStatTest by just assembling and not solving the linear system

Exptacted GridFunctionOperatorTransposed from GridFunctionOperator
parent bb5cd90d
......@@ -145,11 +145,12 @@ namespace AMDiS
}
/// Assemble the matrix operators on the bound element.
template <class RowLocalView, class ColLocalView,
template <class RowLocalView, class ColLocalView, class LocalOperators,
REQUIRES(Concepts::LocalView<RowLocalView>),
REQUIRES(Concepts::LocalView<ColLocalView>)>
void assemble(RowLocalView const& rowLocalView,
ColLocalView const& colLocalView);
ColLocalView const& colLocalView,
LocalOperators& localOperators);
/// Assemble all matrix operators, TODO: incorporate boundary conditions
void assemble();
......@@ -179,8 +180,19 @@ namespace AMDiS
/// Set the flag that forces an update of the pattern since the underlying
/// basis that defines the indexset has been changed
void updateImpl(event::adapt e, index_t<0> i) override { updatePattern_ = true; }
void updateImpl(event::adapt e, index_t<1> i) override { updatePattern_ = true; }
void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(*rowBasis_); }
void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(*colBasis_); }
auto& operators() { return operators_; }
private:
template <class Basis>
void updateImpl2(Basis const& basis)
{
if (!updatePattern_)
Recursive::forEach(operators_, [&](auto& op) { op.update(basis.gridView()); });
updatePattern_ = true;
}
protected:
/// The symmetry property if the bilinear form
......
......@@ -2,10 +2,10 @@
#include <utility>
#include <amdis/Assembler.hpp>
#include <amdis/ContextGeometry.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalOperator.hpp>
#include <amdis/typetree/Traversal.hpp>
#include <amdis/utility/AssembleOperators.hpp>
namespace AMDiS {
......@@ -21,42 +21,35 @@ addOperator(ContextTag contextTag, Expr const& expr,
"col must be a valid treepath, or an integer/index-constant");
auto i = makeTreePath(row);
auto node_i = child(this->rowBasis().localView().treeCache(), i);
auto j = makeTreePath(col);
auto node_j = child(this->colBasis().localView().treeCache(), j);
using LocalContext = typename ContextTag::type;
using Tr = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
auto op = makeLocalOperator<LocalContext>(expr, this->rowBasis().gridView());
auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), node_i, node_j));
operators_[i][j].push(contextTag, std::move(localAssembler));
auto op = makeOperator<LocalContext>(expr, this->rowBasis().gridView());
operators_[i][j].push(contextTag, std::move(op));
updatePattern_ = true;
}
template <class RB, class CB, class T, class Traits>
template <class RowLocalView, class ColLocalView,
template <class RowLocalView, class ColLocalView, class LocalOperators,
REQUIRES_(Concepts::LocalView<RowLocalView>),
REQUIRES_(Concepts::LocalView<ColLocalView>)>
void BiLinearForm<RB,CB,T,Traits>::
assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
LocalOperators& localOperators)
{
elementMatrix_.resize(rowLocalView.size(), colLocalView.size());
elementMatrix_ = 0;
auto const& gv = this->rowBasis().gridView();
auto const& element = rowLocalView.element();
auto geometry = element.geometry();
GlobalContext context{gv, rowLocalView.element(), rowLocalView.element().geometry()};
Traversal::forEachNode(rowLocalView.treeCache(), [&](auto const& rowNode, auto rowTp) {
Traversal::forEachNode(colLocalView.treeCache(), [&](auto const& colNode, auto colTp) {
auto& matOp = operators_[rowTp][colTp];
if (matOp) {
matOp.bind(element, geometry);
assembleOperators(gv, element, matOp, makeMatrixAssembler(rowNode, colNode, elementMatrix_));
matOp.unbind();
}
auto& matOp = localOperators[rowTp][colTp];
matOp.bind(context.element());
matOp.assemble(context, rowNode, colNode, elementMatrix_);
matOp.unbind();
});
});
......@@ -72,13 +65,14 @@ assemble()
auto colLocalView = this->colBasis().localView();
this->init();
auto localOperators = AMDiS::localOperators(operators_);
for (auto const& element : elements(this->rowBasis().gridView(), typename Traits::PartitionSet{})) {
rowLocalView.bind(element);
if (this->rowBasis_ == this->colBasis_)
this->assemble(rowLocalView, rowLocalView);
this->assemble(rowLocalView, rowLocalView, localOperators);
else {
colLocalView.bind(element);
this->assemble(rowLocalView, colLocalView);
this->assemble(rowLocalView, colLocalView, localOperators);
colLocalView.unbind();
}
rowLocalView.unbind();
......
......@@ -21,8 +21,6 @@ install(FILES
AdaptiveGrid.hpp
AdaptStationary.hpp
AMDiS.hpp
Assembler.hpp
AssemblerInterface.hpp
BackupRestore.hpp
BiLinearForm.hpp
BiLinearForm.inc.hpp
......@@ -56,6 +54,7 @@ install(FILES
MeshCreator.hpp
Observer.hpp
Operations.hpp
Operator.hpp
OperatorList.hpp
Output.hpp
PeriodicBC.hpp
......
......@@ -160,4 +160,51 @@ namespace AMDiS
mutable std::optional<LocalGeometry> localGeometry_;
};
template <class GV>
class GlobalContext
{
public:
using GridView = GV;
using Element = typename GV::template Codim<0>::Entity;
using Geometry = typename Element::Geometry;
enum {
dim = GridView::dimension, //< the dimension of the grid element
dow = GridView::dimensionworld //< the dimension of the world
};
/// Constructor. Stores a copy of gridView and a pointer to element and geometry.
GlobalContext(GridView const& gridView, Element const& element,
Geometry const& geometry)
: gridView_(gridView)
, element_(&element)
, geometry_(&geometry)
{}
public:
/// Return the GridView this context is bound to
GridView const& gridView() const
{
return gridView_;
}
/// Return the bound element (entity of codim 0)
Element const& element() const
{
return *element_;
}
/// Return the geometry of the \ref Element
Geometry const& geometry() const
{
return *geometry_;
}
private:
GridView gridView_;
Element const* element_;
Geometry const* geometry_;
};
} // end namespace AMDiS
......@@ -3,11 +3,9 @@
#include <cassert>
#include <type_traits>
#include <amdis/GridFunctions.hpp>
#include <amdis/LocalOperator.hpp>
#include <amdis/common/Transposed.hpp>
#include <amdis/GridFunctionOperatorTransposed.hpp>
#include <amdis/common/Order.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
#include <amdis/utility/QuadratureFactory.hpp>
namespace AMDiS
......@@ -17,256 +15,246 @@ namespace AMDiS
* @{
**/
/// \brief The main implementation of an CRTP-base class for operators using a grid-function
/// coefficient to be used in an \ref Assembler.
template <class LF, class Imp>
class GridFunctionLocalOperator;
/// \brief The main implementation of an operator depending on a grid-function
/**
* An Operator that takes a GridFunction as coefficient.
* Provides quadrature rules and handles the evaluation of the GridFunction at
* local coordinates.
* An Operator that takes a grid-function as coefficient.
* Generates a \ref GridFunctionLocalOperator on \ref localOperator()
*
* The class is specialized, by deriving from it, in \ref GridFunctionOperator.
* The class implements the interface of an \ref Operator.
*
* \tparam Derived The class derived from GridFunctionOperatorBase
* \tparam LC The Element or Intersection type
* \tparam GF The GridFunction, a LocalFunction is created from, and
* that is evaluated at quadrature points.
* \tparam GF The class type of the grid-function
* \tparam Imp Class providing the local assembling method, forwarded to
* GridFunctionLocalOperator class
*
* **Requirements:**
* - `LC` models either Entity (of codim 0) or Intersection
* - `GF` models the \ref Concepts::GridFunction
**/
template <class Derived, class LC, class GF>
class GridFunctionOperatorBase
: public LocalOperator<Derived, LC>
template <class GF, class Imp>
class GridFunctionOperator
{
template <class, class>
friend class LocalOperator;
using ContextType = Impl::ContextTypes<LC>;
using Super = LocalOperator<Derived, LC>;
private:
public:
using GridFunction = GF;
using Implementation = Imp;
/// The type of the localFunction associated with the GridFunction
using LocalFunction = decltype(localFunction(std::declval<GF>()));
/// The Codim=0 entity of the grid, the localFunction can be bound to
using Element = typename ContextType::Entity;
/// The geometry-type of the grid element
using Geometry = typename Element::Geometry;
/// The type of the local coordinates in the \ref Element
using LocalCoordinate = typename GF::EntitySet::LocalCoordinate;
/// A factory for QuadratureRules that incooperate the order of the LocalFunction
using QuadFactory = QuadratureFactory<typename Geometry::ctype, LC::mydimension, LocalFunction>;
public:
/// \brief Constructor. Stores a copy of `gridFct`.
/// \brief Constructor. Stores a copy of `gridFct` and `impl`.
/**
* A GridFunctionOperator takes a gridFunction and the
* differentiation order of the operator, to calculate the
* quadrature degree in \ref getDegree.
* A GridFunctionOperator takes a grid-function and
* an implementation class for the assemble method.
**/
template <class GridFct>
GridFunctionOperatorBase(GridFct&& gridFct, int termOrder)
template <class GridFct, class Impl>
GridFunctionOperator(GridFct&& gridFct, Impl&& impl,
int derivDeg, int gridFctOrder)
: gridFct_(FWD(gridFct))
, termOrder_(termOrder)
, impl_(FWD(impl))
, derivDeg_(derivDeg)
, gridFctOrder_(gridFctOrder)
{}
/// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
/// \tparam PQF A PreQuadratureFactory
template <class PQF>
void setQuadFactory(PQF&& pre)
{
using ctype = typename Geometry::ctype;
quadFactory_ = makeUniquePtr(
makeQuadratureFactory<ctype, LC::mydimension, LocalFunction>(FWD(pre)));
}
protected:
/// Return expression value at LocalCoordinates
auto coefficient(LocalCoordinate const& local) const
{
assert( this->bound_ );
return (*localFct_)(local);
}
template <class GridView>
void update(GridView const&) { /* do nothing */ }
/// Create a quadrature rule using the \ref QuadratureFactory by calculating the
/// quadrature order necessary to integrate the (bi)linear-form accurately.
template <class... Nodes>
auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
friend auto localOperator(GridFunctionOperator const& op)
{
assert( bool(quadFactory_) );
int quadDegree = this->getDegree(termOrder_, quadFactory_->order(), nodes...);
return quadFactory_->rule(type, quadDegree);
return GridFunctionLocalOperator{localFunction(op.gridFct_), op.impl_,
op.derivDeg_, op.gridFctOrder_};
}
private:
/// \brief Binds operator to `element` and `geometry`.
/**
* Binding an operator to the currently visited element in grid traversal.
* Since all operators need geometry information, the `Element::Geometry`
* object `geometry` is created once, during grid traversal, and passed in.
*
* By default, it binds the \ref localFct_ to the `element` and the Quadrature
* factory to the localFunction.
**/
void bind_impl(Element const& element, Geometry const& geometry)
{
assert( bool(quadFactory_) );
localFct_.emplace(localFunction(gridFct_));
localFct_->bind(element);
quadFactory_->bind(localFct_.value());
}
/// Unbinds operator from element.
void unbind_impl()
{
localFct_->unbind();
localFct_.reset();
}
private:
/// The gridFunction to be used within the operator
/// The grid-function associated to this operator
GridFunction gridFct_;
/// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
std::optional<LocalFunction> localFct_;
/// An implementation class for the assembling
Implementation impl_;
/// Assign each element type a quadrature rule
std::shared_ptr<QuadFactory> quadFactory_;
/// Maximal degree of derivative this operator represents
int derivDeg_;
/// the derivative order of this operator (in {0, 1, 2})
int termOrder_ = 0;
/// Polynomial degree of the grid-function (or -1)
int gridFctOrder_;
};
template <class GridFct, class Impl>
GridFunctionOperator(GridFct const& gridFct, Impl const& impl, int, int)
-> GridFunctionOperator<GridFct, Impl>;
/// \brief The transposed operator, implemented in term of its transposed by
/// calling \ref getElementMatrix with inverted arguments.
template <class Derived, class Transposed>
class GridFunctionOperatorTransposed
: public LocalOperator<Derived, typename Transposed::LocalContext>
/// \brief The main implementation of a local-operator depending on a local-function
/**
* A LocalOperator that takes a local-function as coefficient.
* Provides quadrature rules and passes the local-function, bound to an element,
* to the assemble method of an implementation class.
*
* The class implements the interface of a \ref LocalOperator.
*
* \tparam LF The class type of the local-function
* \tparam Imp Class providing the local assembling method
*
* **Requirements:**
* - `LF` models the \ref Concepts::LocalFunction
**/
template <class LF, class Imp>
class GridFunctionLocalOperator
{
template <class, class>
friend class LocalOperator;
private:
/// The type of the localFunction
using LocalFunction = LF;
template <class T, class... Args>
using Constructable = decltype( new T(std::declval<Args>()...) );
/// Type of the implementation class
using Implementation = Imp;
public:
template <class... Args,
std::enable_if_t<Dune::Std::is_detected<Constructable, Transposed, Args...>::value, int> = 0>
GridFunctionOperatorTransposed(Args&&... args)
: transposedOp_(FWD(args)...)
/// \brief Constructor. Stores a copy of `localFct` and `impl`.
/**
* A GridFunctionLocalOperator takes a local-function, an implementation class,
* the differentiation order of the operator and the local-function polynomial
* degree, to calculate the quadrature degree of the operator
**/
template <class LocalFct, class Impl>
GridFunctionLocalOperator(LocalFct&& localFct, Impl&& impl,
int derivDeg, int localFctOrder)
: localFct_(FWD(localFct))
, impl_(FWD(impl))
, derivDeg_(derivDeg)
, localFctOrder_(localFctOrder)
{}
/// Redirects the setQuadFactory call top the transposed operator
/// \tparam PQF A PreQuadratureFactory
template <class PQF>
void setQuadFactory(PQF&& pre)
/// \brief Binds operator to `element`.
/**
* By default, it binds the \ref localFct_ to the `element`.
**/
template <class Element>
void bind(Element const& element)
{
transposedOp_.setQuadFactory(FWD(pre));
localFct_.bind(element);
}
private:
/// Redirects the bind call top the transposed operator
template <class Element, class Geometry>
void bind_impl(Element const& element, Geometry const& geometry)
/// Unbinds operator from element.
void unbind()
{
transposedOp_.bind(element, geometry);
localFct_.unbind();
}
/// Redirects the unbind call top the transposed operator
void unbind_impl()
/// Assemble a local element matrix on the element that is bound.
/**
* This function calls the assemble method from the impl_ class and
* additionally passes a quadrature rule and the localFct_ to that method.
**/
template <class CG, class RN, class CN, class Mat>
void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
Mat& elementMatrix) const
{
transposedOp_.unbind();
auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(),
derivDeg_, localFctOrder(), rowNode, colNode);
impl().assemble(contextGeo, rowNode, colNode, quad, localFct_, elementMatrix);
}
/// Apply the assembling to the transposed elementMatrix with interchanged row-/colNode
/// Assemble a local element vector on the element that is bound.
/**
* \tparam CG ContextGeometry
* \tparam RN RowNode
* \tparam CN ColNode
* \tparam Mat ElementMatrix
* This function calls the assemble method from the impl_ class and
* additionally passes a quadrature rule and the localFct_ to that method.
**/
template <class CG, class RN, class CN, class Mat>
void getElementMatrix(CG const& contextGeometry, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
template <class CG, class Node, class Vec>
void assemble(CG const& contextGeo, Node const& node,
Vec& elementVector) const
{
auto elementMatrixTransposed = transposed(elementMatrix);
transposedOp_.getElementMatrix(
contextGeometry, colNode, rowNode, elementMatrixTransposed);
auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(),
derivDeg_, localFctOrder(), node);
impl().assemble(contextGeo, node, quad, localFct_, elementVector);
}
Implementation & impl() { return impl_; }
Implementation const& impl() const { return impl_; }
protected:
// calculated polynomial order of local-function. Fallback to localFctOrder_;
int localFctOrder() const
{
if constexpr (Concepts::Polynomial<LF>)
return order(localFct_);
else
return localFctOrder_;
}
private:
Transposed transposedOp_;
/// The local-function to be used within the operator
LocalFunction localFct_;
/// Implementation details of the assembling
Implementation impl_;
/// Maximal degree of derivative this operator represents
int derivDeg_;
/// Polynomial degree of the local-function (or -1)
int localFctOrder_;
};
// deduction guide
template <class LocalFct, class Impl>
GridFunctionLocalOperator(LocalFct const& localFct, Impl const& impl, int, int)
-> GridFunctionLocalOperator<LocalFct, Impl>;
#ifndef DOXYGEN
template <class Tag, class PreGridFct, class PQF>
struct PreGridFunctionOperator
/// Registry to specify a tag for each implementation type
template <class Tag, class LocalContext>
struct GridFunctionOperatorRegistry {};
template <class Tag, class Expr>
struct OperatorTerm
{
Tag tag;
PreGridFct expr;
PQF quadFactory;
Expr expr;
int gridFctDeg;
OperatorTerm(Tag tag, Expr const& expr, int gridFctDeg = -1)
: tag(tag)
, expr(expr)
, gridFctDeg(gridFctDeg)
{}
};
#endif
/// Store tag and expression into a \ref PreGridFunctionOperator to create a \ref GridFunctionOperator
template <class Tag, class Expr, class... QuadratureArgs>
auto makeOperator(Tag tag, Expr&& expr, QuadratureArgs&&... args)
/// Store tag and expression into a \ref OperatorTerm to create
/// a \ref GridFunctionOperator
template <class Tag, class Expr>
auto makeOperator(Tag const& tag, Expr&& expr, int gridFctDeg = -1)
{
auto pqf = makePreQuadratureFactory(FWD(args)...);
return PreGridFunctionOperator<Tag, TYPEOF(expr), TYPEOF(pqf)>{tag, FWD(expr), std::move(pqf)};
return OperatorTerm{tag, FWD(expr), gridFctDeg};
}
template <class Tag, class Expr>
auto operatorTerm(Tag const& tag, Expr&& expr, int gridFctDeg = -1)
{
return OperatorTerm{tag, FWD(expr), gridFctDeg};
}
/** @} **/
#ifndef DOXYGEN
/// The base-template for GridFunctionOperators
/**
* An operator can specialize this class, by deriving from \ref GridFunctionOperatorBase.
* With the generic function \ref makeLocalOperator, an instance is created. To
* distinguisch different GridFunction operators, a tag can be provided that has no
* other effect.
*
* \tparam Tag An Identifier for this GridFunctionOperator
* \tparam LC An Element or Intersection the operator is evaluated on
* \tparam GridFct A GridFunction evaluated in local coordinates on the bound element
**/
template <class Tag, class LC, class GridFct>
class GridFunctionOperator
: public GridFunctionOperatorBase<GridFunctionOperator<Tag, LC, GridFct>, LC, GridFct>
{};
template <class Context, class Tag, class GF, class QF>
auto makeGridFunctionOperator(Tag tag, GF&& gf, QF&& qf)
{
GridFunctionOperator<Tag, Context, TYPEOF(gf)> gfo{tag, FWD(gf)};
gfo.setQuadFactory(FWD(qf));
return gfo;
}
template <class R, class Tag>
using IsTransposed = decltype(R::transposedTag(std::declval<Tag>()));
/// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
/// Generate an \ref GridFunctionOperator from a OperatorTerm (tag, expr).
/// @{
template <class Context, class... Args, class GridView>