Commit 36d6a365 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/local_operator' into 'develop'

Feature/local operator

See merge request !24
parents fd552d52 15764e96
Pipeline #1111 passed with stage
in 3 minutes and 24 seconds
......@@ -32,6 +32,8 @@ int main(int argc, char** argv)
auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
prob.addVectorOperator(opForce, _0);
auto opForce2 = makeOperator(tag::test{}, [](auto const& x) { return -2.0; }, 0);
prob.addVectorOperator(BoundaryType{0}, opForce2, _0);
// set boundary condition
auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
......
......@@ -8,7 +8,7 @@
#include <amdis/DirichletBC.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/LocalAssemblerBase.hpp>
#include <amdis/LocalAssemblerList.hpp>
#include <amdis/common/Mpl.hpp>
#include <amdis/common/TypeDefs.hpp>
......@@ -54,6 +54,8 @@ namespace AMDiS
SystemVectorType& rhs,
bool asmMatrix, bool asmVector) const;
/// Assemble operators on an element, by passing the element/intersection to
/// `elementAssembler` functor.
template <class Element, class Operators, class ElementAssembler>
void assembleElementOperators(
Element const& element,
......@@ -70,14 +72,14 @@ namespace AMDiS
bool asmMatrix, bool asmVector) const;
/// Return whether the matrix-block needs to be assembled
/// Return the element the LocalViews are bound to
template <class LocalView0, class... LovalViews>
auto const& getElement(LocalView0 const& localView, LovalViews const&...) const
{
return localView.element();
}
/// Return whether the matrix-block needs to be assembled
/// Return the gridView the localViews are bound to
template <class LocalView0, class... LovalViews>
auto const& getGridView(LocalView0 const& localView, LovalViews const&...) const
{
......
......@@ -45,12 +45,12 @@ void Assembler<Traits>::assemble(
rowLocalView.bind(element);
auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.assemble(asmVector) && !rhsOp.empty()) {
if (rhsOp.doAssemble(asmVector) && !rhsOp.empty()) {
rhsOp.bind(element, geometry);
auto vecAssembler = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list)
scaled.op->assemble(context, elementVector, rowLocalView.tree());
scaled.op->assemble(context, rowLocalView.tree(), elementVector);
};
this->assembleElementOperators(element, rhsOp, vecAssembler);
......@@ -59,7 +59,7 @@ void Assembler<Traits>::assemble(
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto& matOp = matrixOperators_[rowNode][colNode];
if (matOp.assemble(asmMatrix) && !matOp.empty()) {
if (matOp.doAssemble(asmMatrix) && !matOp.empty()) {
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
auto colLocalView = colBasis.localView();
colLocalView.bind(element);
......@@ -68,7 +68,7 @@ void Assembler<Traits>::assemble(
auto matAssembler = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list)
scaled.op->assemble(context, elementMatrix, rowLocalView.tree(), colLocalView.tree());
scaled.op->assemble(context, rowLocalView.tree(), colLocalView.tree(), elementMatrix);
};
this->assembleElementOperators(element, matOp, matAssembler);
......@@ -121,10 +121,10 @@ template <class Traits>
void Assembler<Traits>::assembleElementOperators(
Element const& element,
Operators& operators,
ElementAssembler const& elementAssembler) const
ElementAssembler const& localAssembler) const
{
// assemble element operators
elementAssembler(element, operators.element);
localAssembler(element, operators.element);
// assemble intersection operators
if (!operators.intersection.empty()
......@@ -132,9 +132,9 @@ void Assembler<Traits>::assembleElementOperators(
{
for (auto const& intersection : intersections(globalBasis_.gridView(), element)) {
if (intersection.boundary())
elementAssembler(intersection, operators.boundary);
localAssembler(intersection, operators.boundary);
else
elementAssembler(intersection, operators.intersection);
localAssembler(intersection, operators.intersection);
}
}
}
......@@ -191,14 +191,14 @@ std::size_t Assembler<Traits>::finishMatrixVector(
{
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.assemble(asmVector))
if (rhsOp.doAssemble(asmVector))
rhsOp.assembled = true;
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
auto& matOp = matrixOperators_[rowNode][colNode];
if (matOp.assemble(asmMatrix))
if (matOp.doAssemble(asmMatrix))
matOp.assembled = true;
// finish boundary condition
......
#pragma once
#include <type_traits>
#include <dune/common/typetraits.hh>
#include <dune/common/std/optional.hh>
#include <dune/geometry/type.hh>
namespace AMDiS
{
template <class LocalContext, class Geometry, class LocalGeometry>
namespace Impl
{
template <class E, class = Dune::void_t<>>
struct ContextTypes
{
using Entity = E;
using LocalGeometry = typename E::Geometry;
};
// specialization for intersections
template <class I>
struct ContextTypes<I, Dune::void_t<decltype(std::declval<I>().inside())>>
{
using Entity = typename I::Entity;
using LocalGeometry = typename I::LocalGeometry;
};
} // end namespace Impl
/// \brief Wrapper class for element and geometry
/**
* A LocalContext can be either a grid entity of codim 0 (called an element)
* or an intersection of elements. The element and its geometry may be stored
* externally and can be passed along with the localContext object.
* Since an intersection has a geometry (and localGeometry) different from the
* geometry (and localGeometry) of the entity it belongs to, these objects
* are provided as well.
**/
template <class LocalContextType>
struct ContextGeometry
{
public:
using LocalContext = LocalContextType;
using Element = typename Impl::ContextTypes<LocalContext>::Entity;
using Geometry = typename Element::Geometry;
using LocalGeometry = typename Impl::ContextTypes<LocalContext>::LocalGeometry;
using IsEntity = std::is_same<Element, LocalContext>;
enum {
dim = Geometry::mydimension, //< the dimension of the grid element
dow = Geometry::coorddimension //< the dimension of the world
};
LocalContext const& localContext;
Geometry const& geometry;
LocalGeometry const& localGeometry;
ContextGeometry(LocalContext const& localContext,
Geometry const& geometry,
LocalGeometry const& localGeometry)
: localContext(localContext)
, geometry(geometry)
, localGeometry(localGeometry)
/// Constructor. Stores pointer to localContext, element, and geometry.
ContextGeometry(LocalContext const& localContext, Element const& element, Geometry const& geometry)
: localContext_(&localContext)
, element_(&element)
, geometry_(&geometry)
{}
/// Coordinate `p` given in `localGeometry`, transformed to coordinate in `geometry`.
template <class Coordinate>
decltype(auto) position(Coordinate const& p) const
public:
Element const& element() const
{
return position(p, std::is_same<Geometry, LocalGeometry>{});
return *element_;
}
/// The integration element from the `localGeometry`, the quadrature points are
/// defined in.
LocalContext const& localContext() const
{
return *localContext_;
}
Geometry const& geometry() const
{
return *geometry_;
}
LocalGeometry const& localGeometry() const
{
return localGeometryImpl(IsEntity{});
}
public:
/// Coordinate `p` given in `localGeometry`, transformed to coordinate in `geometry`.
template <class Coordinate>
auto integrationElement(Coordinate const& p) const
decltype(auto) local(Coordinate const& p) const
{
return localGeometry.integrationElement(p);
return localImpl(p, IsEntity{});
}
/// Transformation of coordinate `p` given in `localGeometry` to world space coordinates.
template <class Coordinate>
decltype(auto) global(Coordinate const& p) const
{
return geometry.global(p);
return geometry_->global(p);
}
/// Return the geometry-type of the localContext
Dune::GeometryType type() const
{
return localContext_->type();
}
/// The integration element from the `localGeometry`, the quadrature points are
/// defined in.
template <class Coordinate>
auto integrationElement(Coordinate const& p) const
{
return localGeometry().integrationElement(p);
}
private: // implementation detail
// position for elements
template <class Coordinate>
Coordinate const& position(Coordinate const& p, std::true_type) const
Coordinate const& localImpl(Coordinate const& p, std::true_type) const
{
return p;
}
// position for intersection
template <class Coordinate>
auto position(Coordinate const& p, std::false_type) const
auto localImpl(Coordinate const& p, std::false_type) const
{
return localGeometry().global(p);
}
// local-geometry is the same as geometry
Geometry const& localGeometryImpl(std::true_type) const
{
return *geometry_;
}
// local-geometry of intersection in inside element
LocalGeometry const& localGeometryImpl(std::false_type) const
{
return localGeometry.global(p);
if (!localGeometry_)
localGeometry_.emplace(localContext_->geometryInInside());
return *localGeometry_;
}
private:
LocalContext const* localContext_;
Element const* element_;
Geometry const* geometry_;
// The localGeometry may be constructed only if needed
mutable Dune::Std::optional<LocalGeometry> localGeometry_;
};
} // end namespace AMDiS
#pragma once
#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/IndexSeq.hpp>
#include <amdis/common/Loops.hpp>
namespace AMDiS
{
template <class FeSpace>
struct LocalView
{
using type = typename FeSpace::LocalView;
};
template <class FeSpace>
struct LocalIndexSet
{
using type = typename FeSpace::LocalIndexSet;
};
template <class FeSpaces>
class FiniteElementSpaces
{
template <std::size_t I>
using FeSpace = std::tuple_element_t<I, FeSpaces>;
static constexpr int nComponents = std::tuple_size<FeSpaces>::value;
static_assert( nComponents > 0, "" );
using LocalViews = MapTuple_t<LocalView, FeSpaces>;
using LocalIndexSets = MapTuple_t<LocalIndexSet, FeSpaces>;
/// The grid view the global FE basis lives on
using GridView = typename FeSpace<0>::GridView;
/// Type of the grid element we are bound to
using Element = typename GridView::template Codim<0>::Entity;
public:
explicit FiniteElementSpaces(std::shared_ptr<FeSpaces> const& feSpaces)
: feSpaces_(feSpaces)
, localViews_(mapTuple([](auto const& basis) { return basis.localView(); }, *feSpaces))
, localIndexSets_(mapTuple([](auto const& basis) { return basis.localIndexSet(); }, *feSpaces))
{}
/// Update all global bases. This will update the indexing information of the global basis.
/// NOTE: It must be called if the grid has changed.
void update(GridView const& gv)
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(*feSpaces_).update(gv);
});
}
/// Bind the \ref localViews and \ref localIndexSets to a grid element
void bind(Element const& element)
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
auto& localView = std::get<I>(localViews_);
localView.bind(element);
auto& localIndexSet = std::get<I>(localIndexSets_);
localIndexSet.bind(localView);
});
// NOTE: maybe create element-geometry here
bound_ = true;
}
/// Unbind from the current element
void unbind()
{
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(localIndexSets_).unbind();
std::get<I>(localViews_).unbind();
});
bound_ = false;
}
template <std::size_t I>
auto const& feSpace(const index_t<I> _i = {}) const
{
return std::get<I>(*feSpaces_);
}
template <std::size_t I>
auto const& localView(const index_t<I> _i = {}) const
{
assert( bound_ && "localViews must be bound to an element." );
return std::get<I>(localViews_);
}
template <std::size_t I>
auto const& localIndexSet(const index_t<I> _i = {}) const
{
assert( bound_ && "localIndexSets must be bound to a localView." );
return std::get<I>(localIndexSets_);
}
auto const& element() const
{
assert( bound_ && "localViews must be bound to an element." );
return std::get<0>(localViews_).element();
}
auto const& gridView() const
{
return std::get<0>(*feSpaces_).gridView();
}
private:
/// Tuple of global functionspace bases
std::shared_ptr<FeSpaces> feSpaces_;
/// Tuple of localView objects, obtained from the tuple of global bases
LocalViews localViews_;
/// Tuple of localIndexSet objects, obtained from the tuple of global bases
LocalIndexSets localIndexSets_;
/// True, if localViews and localIndexSets are bound to an element
bool bound_ = false;
};
} // end namespace AMDiS
......@@ -4,6 +4,8 @@
#include <type_traits>
#include <amdis/GridFunctions.hpp>
#include <amdis/LocalOperator.hpp>
#include <amdis/QuadratureFactory.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/utility/FiniteElementType.hpp>
......@@ -14,7 +16,8 @@ namespace AMDiS
* @{
**/
/// \brief The main implementation of an operator to be used in a \ref LocalAssembler.
/// \brief The main implementation of an CRTP-base class for operators using a grid-function
/// coefficient to be used in a \ref LocalAssembler.
/**
* An Operator that takes a GridFunction as coefficient.
* Provides quadrature rules and handles the evaluation of the GridFunction at
......@@ -22,54 +25,58 @@ namespace AMDiS
*
* The class is specialized, by deriving from it, in \ref GridFunctionOperator.
*
* \tparam Derived The class derived from GridFunctionOperatorBase
* \tparam LocalContext The Element or Intersection type
* \tparam GridFunction The GridFunction, a LocalFunction is created from, and
* that is evaluated at quadrature points.
* \tparam QuadratureCreator A functor that provides a \ref Dune::QuadratureRule.
*
* **Requirements:**
* - `LocalContext` models either Entity (of codim 0) or Intersection
* - `GridFunction` models the \ref Concepts::GridFunction
* - `QuadratureCreator` models \ref Concepts::Callable<Dune::GeometryType, LocalFunction, F>
* where `F` is a functor of the signature `int(int)` that calculates the
* degree of the (bi)linear-form. The argument passed to `F` is the polynomial
* order of the GridFunction.
**/
template <class GridFunction, class QuadratureCreator>
template <class Derived, class LocalContext, class GridFunction>
class GridFunctionOperatorBase
: public LocalOperator<Derived, LocalContext>
{
/// The type of the localFunction associated with the GridFunction
using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
using Element = typename GridFunction::EntitySet::Element;
/// The Codim=0 entity of the grid, the localFunction can be bound to
using Element = typename Impl::ContextTypes<LocalContext>::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 GridFunction::EntitySet::LocalCoordinate;
/// A factory for QuadratureRules that incooperate the order of the LocalFunction
using QuadFactory = QuadratureFactory<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>;
public:
/// \brief Constructor. Stores a copy of `expr`.
/// \brief Constructor. Stores a copy of `gridFct`.
/**
* An expression operator takes an expression, following the interface of
* \ref ExpressionBase, and stores a copy. Additionally, it gets the
* differentiation order, to calculate the quadrature degree in \ref getDegree.
* A GridFunctionOperator takes a gridFunction and the
* differentiation order of the operator, to calculate the
* quadrature degree in \ref getDegree.
**/
GridFunctionOperatorBase(GridFunction const& gridFct, QuadratureCreator creator, int order)
GridFunctionOperatorBase(GridFunction const& gridFct, int termOrder)
: gridFct_(gridFct)
, localFct_(localFunction(gridFct_))
, quadCreator_(creator)
, order_(order)
, termOrder_(termOrder)
{}
// Copy constructor
GridFunctionOperatorBase(GridFunctionOperatorBase const& that)
: gridFct_(that.gridFct_)
, localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
, quadCreator_(that.quadCreator_)
, order_(that.order_)
, termOrder_(that.termOrder_)
{}
// Move constructor
GridFunctionOperatorBase(GridFunctionOperatorBase&& that)
: gridFct_(std::move(that.gridFct_))
, localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
, quadCreator_(std::move(that.quadCreator_))
, order_(std::move(that.order_))
, termOrder_(std::move(that.termOrder_))
{}
/// \brief Binds operator to `element` and `geometry`.
......@@ -78,328 +85,197 @@ namespace AMDiS
* 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`.
* By default, it binds the \ref localFct_ to the `element` and the Quadrature
* factory to the localFunction.
**/
void bind(Element const& element, Geometry const& geometry)
void bind_impl(Element const& element, Geometry const& geometry)
{
if (!bound_) {
localFct_.bind(element);
isSimplex_ = geometry.type().isSimplex();
isAffine_ = geometry.affine();
bound_ = true;
}
assert( bool(quadFactory_) );
localFct_.bind(element);
quadFactory_->bind(localFct_);
}
/// Unbinds operator from element.
void unbind()
void unbind_impl()
{
bound_ = false;
localFct_.unbind();
}
/// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
template <class PreQuadFactory>
void setQuadFactory(PreQuadFactory const& pre)
{
quadFactory_ = makeQuadratureFactoryPtr<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>(pre);
}
protected:
/// Return expression value at LocalCoordinates
auto coefficient(LocalCoordinate const& local) const
{
assert( bound_ );
assert( this->bound_ );
return localFct_(local);
}
/// Create a quadrature rule using the \ref QuadratureCreator by calculating the
/// Create a quadrature rule using the \ref QuadratureFactory by calculating the
/// quadrature order necessary to integrate the (bi)linear-form accurately.
template <class... Nodes>
decltype(auto) getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
{
auto degreeFunctor = [&,this](int order) -> int { return this->getDegree(order, nodes...); };
return quadCreator_(type, localFct_, degreeFunctor);
}
protected:
/// \brief Return the quadrature degree for a vector operator.
/**
* The quadrature degree that is necessary, to integrate the expression
* multiplied with (possibly derivatives of) shape functions. Thus it needs
* the order of derivatives, this operator implements.
**/
template <class Node>
int getDegree(int coeffDegree, Node const& node) const
{
assert( bound_ );
int psiDegree = polynomialDegree(node);
int degree = psiDegree + coeffDegree;
if (isSimplex_)
degree -= order_;
if (isAffine_)
degree += 1; // oder += (order+1)