diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index ddaad29b138ebdfd9a486c6f5f910851025bc160..6baf47bfb8cbba1adf2a4ffc10a5cb46f742a29f 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -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
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index c5f6ab4e59cb64e178cdc8ad6cd4bc1e96612d7f..8473ac32278e44489042ffa87e253b582653c7cf 100644
--- a/src/amdis/Assembler.hpp
+++ b/src/amdis/Assembler.hpp
@@ -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
{
diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp
index 1cf4a4eaf39f721d145601816fee68c248095c87..05501f75e98c37487a374f2bf227eab6e15bf3a0 100644
--- a/src/amdis/Assembler.inc.hpp
+++ b/src/amdis/Assembler.inc.hpp
@@ -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
diff --git a/src/amdis/ContextGeometry.hpp b/src/amdis/ContextGeometry.hpp
index d679715fe4ea90c48f418aa509b0ee0adf68143f..6307603e484c7d19668e360b1a659a6e3809cab0 100644
--- a/src/amdis/ContextGeometry.hpp
+++ b/src/amdis/ContextGeometry.hpp
@@ -1,65 +1,157 @@
#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
diff --git a/src/amdis/FiniteElementSpaces.hpp b/src/amdis/FiniteElementSpaces.hpp
deleted file mode 100644
index 10775960cc76104b8a24d2a8a217c52b4597f22e..0000000000000000000000000000000000000000
--- a/src/amdis/FiniteElementSpaces.hpp
+++ /dev/null
@@ -1,133 +0,0 @@
-#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
diff --git a/src/amdis/GridFunctionOperator.hpp b/src/amdis/GridFunctionOperator.hpp
index 9b828ff797e29843d5daa253c6e69464cc616bec..70a12642d81e1fd311b8957c6833bf0e751378f2 100644
--- a/src/amdis/GridFunctionOperator.hpp
+++ b/src/amdis/GridFunctionOperator.hpp
@@ -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)
-
- return degree;
- }
-
-
- /// \brief Return the quadrature degree for a matrix 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 RowNode, class ColNode>
- int getDegree(int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
+ auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
{
- assert( bound_ );
-
- int psiDegree = polynomialDegree(rowNode);
- int phiDegree = polynomialDegree(colNode);
-
- int degree = psiDegree + phiDegree + coeffDegree;
- if (isSimplex_)
- degree -= order_;
- if (isAffine_)
- degree += 1; // oder += (order+1)
-
- return degree;
+ assert( bool(quadFactory_) );
+ int quadDegree = this->getDegree(termOrder_, quadFactory_->order(), nodes...);
+ return quadFactory_->rule(type, quadDegree);
}
private:
+ /// The gridFunction to be used within the operator
GridFunction gridFct_;
+
+ /// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
LocalFunction localFct_;
- QuadratureCreator quadCreator_; //< a creator to provide a quadrature rule
- int order_; //< the derivative order of this operator
+ /// Assign each element type a quadrature rule
+ std::shared_ptr<QuadFactory> quadFactory_ = std::make_shared<QuadFactory>();
- bool isSimplex_ = false; //< the bound element is a simplex
- bool isAffine_ = false; //< the bound geometry is affine
- bool bound_ = false; //< an element is bound to the operator
+ /// the derivative order of this operator (in {0, 1, 2})
+ int termOrder_ = 0;
};
- /// A default implementation of an GridFunctionOperator if no specialization is available.
- /**
- * An operator must implement either \ref calculateElementVector, or
- * \ref calculateElementMatrix, if it is a vector or matrix operator,
- * respectively.
- **/
- template <class Tag, class GridFct, class QuadratureCreator>
- class GridFunctionOperator
- : public GridFunctionOperatorBase<GridFct, QuadratureCreator>
+ /// \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>
{
+ template <class T, class... Args>
+ using Constructable = decltype( new T(std::declval<Args>()...) );
+
public:
- GridFunctionOperator(Tag, GridFct const& gridFct, QuadratureCreator const& quadCreator)
- : GridFunctionOperatorBase<GridFct, QuadratureCreator>(gridFct, 0, quadCreator)
+ template <class... Args,
+ std::enable_if_t<Dune::Std::is_detected<Constructable, Transposed, Args...>::value, int> = 0>
+ GridFunctionOperatorTransposed(Args&&... args)
+ : transposedOp_(std::forward<Args>(args)...)
{}
- /// \brief Assemble a local element vector on the element that is bound.
- /**
- * \param contextGeometry Container for context related data
- * \param quad A quadrature formula
- * \param elementVector The output element vector
- * \param node The node of the test-basis
- **/
- template <class Context, class QuadratureRule,
- class ElementVector, class Node>
- void calculateElementVector(Context const& contextGeometry,
- QuadratureRule const& quad,
- ElementVector& elementVector,
- Node const& node)
+ /// Redirects the bind call top the transposed operator
+ template <class Element, class Geometry>
+ void bind_impl(Element const& element, Geometry const& geometry)
{
- error_exit("Needs to be implemented!");
+ transposedOp_.bind_impl(element, geometry);
}
- /// \brief Assemble a local element matrix on the element that is bound.
- /**
- * \param contextGeometry Container for context related data
- * \param quad A quadrature formula
- * \param elementMatrix The output element matrix
- * \param rowNode The node of the test-basis
- * \param colNode The node of the trial-basis
- * \param flag1 finiteelement is the same for test and trial
- * \param flag2 nodes are the same in the basis-tree
- **/
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& contextGeometry,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE> flag1,
- std::integral_constant<bool, sameNode> flag2)
+ /// Redirects the unbind call top the transposed operator
+ void unbind_impl()
{
- error_exit("Needs to be implemented!");
+ transposedOp_.unbind_impl();
}
- };
-
-#ifndef DOXYGEN
- namespace Concepts
- {
- namespace Definition
+ /// Redirects the setQuadFactory call top the transposed operator
+ template <class PreQuadFactory>
+ void setQuadFactory(PreQuadFactory const& pre)
{
- struct HasGridFunctionOrder
- {
- template <class F, class GridView>
- auto requires_(F&& f, GridView const& gridView)
- -> decltype( order(localFunction(makeGridFunction(f, gridView))) );
- };
+ transposedOp_.setQuadFactory(pre);
}
- template <class F, class GridView = Dune::YaspGrid<2>::LeafGridView>
- constexpr bool HasGridFunctionOrder = models<Definition::HasGridFunctionOrder(F, GridView)>;
-
- } // end namespace Concepts
-
-
- namespace tag
- {
- struct deduce {};
- struct order {};
- struct rule {};
-
- } // end namespace tag
-
-
- template <class Tag, class Expr, class QuadTag, class Rule = void>
- struct ExpressionPreOperator;
-
- template <class Tag, class Expr>
- struct ExpressionPreOperator<Tag, Expr, tag::deduce>
- {
- Tag tag;
- Expr expr;
- };
+ /// Apply the assembling to the transposed elementMatrix with interchanged row-/colNode
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
+ {
+ auto elementMatrixTransposed = trans(elementMatrix);
+ transposedOp_.getElementMatrix(
+ context, colNode, rowNode, elementMatrixTransposed);
+ }
- template <class Tag, class Expr>
- struct ExpressionPreOperator<Tag, Expr, tag::order>
- {
- Tag tag;
- Expr expr;
- int order;
- Dune::QuadratureType::Enum qt;
+ private:
+ Transposed transposedOp_;
};
- template <class Tag, class Expr, class Rule>
- struct ExpressionPreOperator<Tag, Expr, tag::rule, Rule>
+ template <class Tag, class PreGridFct, class PreQuadFactory>
+ struct PreGridFunctionOperator
{
Tag tag;
- Expr expr;
- Rule const& rule;
+ PreGridFct expr;
+ PreQuadFactory quadFactory;
};
-#endif
-
/// Store tag and expression to create a \ref GridFunctionOperator
- template <class Tag, class Expr>
- auto makeOperator(Tag tag, Expr const& expr)
- {
- using RawExpr = Underlying_t<Expr>;
- static_assert(Concepts::HasGridFunctionOrder<RawExpr>,
- "Polynomial degree of expression can not be deduced. You need to provide a polynomial order or a quadrature rule in 'makeOperator()'.");
-
- return Dune::Hybrid::ifElse(bool_<Concepts::HasGridFunctionOrder<RawExpr>>,
- [&](auto) { return ExpressionPreOperator<Tag, Expr, tag::deduce>{tag, expr}; },
- [&](auto) { return ExpressionPreOperator<Tag, double, tag::deduce>{tag, 0.0}; });
- }
-
- /// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator
- template <class Tag, class Expr>
- auto makeOperator(Tag tag, Expr const& expr, int order,
- Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
- {
- return ExpressionPreOperator<Tag, Expr, tag::order>{tag, expr, order, qt};
- }
-
- /// Store tag and expression and a quadrature rule to create a \ref GridFunctionOperator
- template <class Tag, class Expr, class ctype, int dim>
- auto makeOperator(Tag tag, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule)
+ template <class Tag, class PreGridFct, class... QuadratureArgs>
+ auto makeOperator(Tag tag, PreGridFct const& expr, QuadratureArgs&&... args)
{
- return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{tag, expr, rule};
+ auto preQuadFactory = makePreQuadratureFactory(std::forward<QuadratureArgs>(args)...);
+ return PreGridFunctionOperator<Tag, PreGridFct, decltype(preQuadFactory)>{tag, expr, preQuadFactory};
}
/** @} **/
#ifndef DOXYGEN
- namespace Impl
- {
- // Deduce polynomial order of expression automatically. Create standard quadrature rule with degree
- // build up of operator derivative order, and polynomial degrees of trial/test functions.
- template <class Tag, class Expr, class GridView>
- auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::deduce> const& /*op*/, GridView const& /*gridView*/)
- {
- using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
- return [](auto type, auto&& localFct, auto getDegree) -> auto const&
- {
- return QuadratureRules::rule(type, getDegree(order(localFct)));
- };
- }
-
- // Provide polynomial order of expression explicitly
- template <class Tag, class Expr, class GridView>
- auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::order> const& op, GridView const& /*gridView*/)
- {
- using QuadratureRules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
- return [order=op.order,qt=op.qt](auto type, auto&& /*localFct*/, auto getDegree) -> auto const&
- {
- return QuadratureRules::rule(type, getDegree(order), qt);
- };
- }
-
- // Provide quadrature rule explicitly
- template <class Tag, class Expr, class Rule, class GridView>
- auto makeQuadCreator(ExpressionPreOperator<Tag, Expr, tag::rule, Rule> const& op, GridView const& /*gridView*/)
- {
- return [&rule=op.rule](auto&&...) -> auto const&
- {
- return rule;
- };
- }
-
- } // end namespace Impl
+ /// 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 LocalContext 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 LocalContext, class GridFct>
+ class GridFunctionOperator
+ : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LocalContext, GridFct>, LocalContext, GridFct>
+ {};
/// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
/// @{
- template <class Tag, class... Args, class GridView>
- auto makeGridOperator(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
+ template <class LocalContext, class Tag, class... Args, class GridView>
+ auto makeLocalOperator(PreGridFunctionOperator<Tag, Args...> const& op, GridView const& gridView)
{
auto gridFct = makeGridFunction(op.expr, gridView);
- auto quadCreator = Impl::makeQuadCreator(op, gridView);
- using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
- return GridFctOp{op.tag, gridFct, quadCreator};
+ using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
+ GridFctOp localOperator{op.tag, gridFct};
+ localOperator.setQuadFactory(op.quadFactory);
+ return localOperator;
}
- template <class Tag, class... Args, class GridView>
- auto makeGridOperator(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
+ template <class LocalContext, class Tag, class... Args, class GridView>
+ auto makeLocalOperator(std::reference_wrapper<PreGridFunctionOperator<Tag, Args...>> op, GridView const& gridView)
{
- ExpressionPreOperator<Tag, Args...> const& op_ref = op;
+ PreGridFunctionOperator<Tag, Args...> const& op_ref = op;
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
- auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
- using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
- return GridFctOp{op_ref.tag, gridFct, quadCreator};
+ using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
+ GridFctOp localOperator{op_ref.tag, gridFct};
+ localOperator.setQuadFactory(op_ref.quadFactory);
+ return localOperator;
}
/// @}
/// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr).
/// @{
- template <class Tag, class... Args, class GridView>
- auto makeGridOperatorPtr(ExpressionPreOperator<Tag, Args...> const& op, GridView const& gridView)
+ template <class LocalContext, class Tag, class... Args, class GridView>
+ auto makeLocalOperatorPtr(PreGridFunctionOperator<Tag, Args...> const& op, GridView const& gridView)
{
auto gridFct = makeGridFunction(op.expr, gridView);
- auto quadCreator = Impl::makeQuadCreator(op, gridView);
- using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
- return std::make_shared<GridFctOp>(op.tag, gridFct, quadCreator);
+ using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
+ auto localOperator = std::make_shared<GridFctOp>(op.tag, gridFct);
+ localOperator->setQuadFactory(op.quadFactory);
+ return localOperator;
}
- template <class Tag, class... Args, class GridView>
- auto makeGridOperatorPtr(std::reference_wrapper<ExpressionPreOperator<Tag, Args...>> op, GridView const& gridView)
+ template <class LocalContext, class Tag, class... Args, class GridView>
+ auto makeLocalOperatorPtr(std::reference_wrapper<PreGridFunctionOperator<Tag, Args...>> op, GridView const& gridView)
{
- ExpressionPreOperator<Tag, Args...> const& op_ref = op;
+ PreGridFunctionOperator<Tag, Args...> const& op_ref = op;
auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView);
- auto quadCreator = Impl::makeQuadCreator(op_ref, gridView);
- using GridFctOp = GridFunctionOperator<Tag, decltype(gridFct), decltype(quadCreator)>;
- return std::make_shared<GridFctOp>(op_ref.tag, gridFct, quadCreator);
+ using GridFctOp = GridFunctionOperator<Tag, LocalContext, decltype(gridFct)>;
+ auto localOperator = std::make_shared<GridFctOp>(op_ref.tag, gridFct);
+ localOperator->setQuadFactory(op_ref.quadFactory);
+ return localOperator;
}
/// @}
diff --git a/src/amdis/LocalAssembler.hpp b/src/amdis/LocalAssembler.hpp
index 424cc4560b42df6a9a4e98f5dc69757fbb0ea4ef..1d47532dfd273c401c1c5182483ba0e19ecf9ba2 100644
--- a/src/amdis/LocalAssembler.hpp
+++ b/src/amdis/LocalAssembler.hpp
@@ -23,12 +23,8 @@ namespace AMDiS
private:
using Element = typename Super::Element;
- using ElementMatrixVector = typename Super::ElementMatrixVector;
-
using Geometry = typename Super::Geometry;
- using LocalGeometry = typename Super::LocalGeometry;
-
- using Context = ContextGeometry<LocalContext, Geometry, LocalGeometry>;
+ using ElementMatrixVector = typename Super::ElementMatrixVector;
public:
@@ -69,26 +65,16 @@ namespace AMDiS
/// Implementation of \ref LocalAssemblerBase::assemble
/**
- * Provides a quadrature formula with necessary degree to integrate the operator,
- * stores geometry and localGeometry and a \ref ContextGeometry and calls
+ * Stores geometry and localGeometry and calls
* \ref calculateElementVector or \ref calculateElementMatrix on the
* vector or matrix operator, respectively.
- *
- * The call to operator passes two optimization flags:
- * 1. RowNode and ColNode implement the same local finiteelement
- * 2. RowNode and ColNode are the same node in the globalBasis tree.
**/
- virtual bool assemble(
- LocalContext const& localContext,
- ElementMatrixVector& elementMatrixVector,
- Nodes const&... nodes) final
+ virtual void assemble(LocalContext const& localContext,
+ Nodes const&... nodes,
+ ElementMatrixVector& elementMatrixVector) final
{
- auto&& localGeometry = getLocalGeometry(localContext);
- auto&& quad = op_.getQuadratureRule(localGeometry.type(), nodes...);
-
- Context data{localContext, getGeometry(), localGeometry};
- assembleImpl(data, std::forward<decltype(quad)>(quad), elementMatrixVector, nodes...);
- return true;
+ ContextGeometry<LocalContext> context{localContext, getElement(), getGeometry()};
+ assembleImpl(context, nodes..., elementMatrixVector);
}
@@ -97,96 +83,46 @@ namespace AMDiS
protected: // implementation detail
// matrix assembling
- template <class QuadratureRule, class ElementMatrix, class RowNode, class ColNode>
- void assembleImpl(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode, ColNode const& colNode)
- {
- assembleImpl(context, quad, elementMatrix, rowNode, colNode,
- std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{});
- }
-
- // local finite elements are the same for rowNode and colNode
- template <class QuadratureRule, class ElementMatrix, class RowNode, class ColNode>
- void assembleImpl(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode, ColNode const& colNode,
- std::true_type /*sameFE*/)
- {
- if (rowNode.treeIndex() == colNode.treeIndex())
- op_.calculateElementMatrix(
- context, quad, elementMatrix, rowNode, colNode, std::true_type{}, std::true_type{});
- else
- op_.calculateElementMatrix(
- context, quad, elementMatrix, rowNode, colNode, std::true_type{}, std::false_type{});
- }
-
- // local finite elements are different for rowNode and colNode
- template <class QuadratureRule, class ElementMatrix, class RowNode, class ColNode>
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
void assembleImpl(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
RowNode const& rowNode, ColNode const& colNode,
- std::false_type /*sameFE*/)
+ ElementMatrix& elementMatrix)
{
- op_.calculateElementMatrix(
- context, quad, elementMatrix, rowNode, colNode, std::false_type{}, std::false_type{});
+ op_.calculateElementMatrix(context, rowNode, colNode, elementMatrix);
}
-
// vector assembling
- template <class QuadratureRule, class ElementVector, class Node>
+ template <class Context, class Node, class ElementVector>
void assembleImpl(Context const& context,
- QuadratureRule const& quad,
- ElementVector& elementVector,
- Node const& node)
+ Node const& node,
+ ElementVector& elementVector)
{
- op_.calculateElementVector(context, quad, elementVector, node);
+ op_.calculateElementVector(context, node, elementVector);
}
#endif // DOXYGEN
-
public:
+ /// return the bound entity (of codim 0)
Element const& getElement() const
{
assert( element_ );
return *element_;
}
+ /// return the geometry of the bound element
Geometry const& getGeometry() const
{
assert( geometry_ );
return *geometry_;
}
- decltype(auto) getLocalGeometry(LocalContext const& context) const
- {
- return getLocalGeometry(context, std::is_same<LocalContext, Element>{});
- }
-
-
- private: // implementation detail
-
- // localGeometry for Entities
- Geometry const& getLocalGeometry(Element const& element, std::true_type) const
- {
- return *geometry_;
- }
-
- // localGeometry for intersections
- LocalGeometry getLocalGeometry(LocalContext const& intersection, std::false_type) const
- {
- return intersection.geometryInInside();
- }
-
private:
- std::unique_ptr<Operator> storage_; //< the stored operator, implementing \ref GridFunctionOperatorBase
+ /// the stored operator, implementing \ref GridFunctionOperatorBase
+ std::unique_ptr<Operator> storage_;
Operator& op_;
Element const* element_ = nullptr;
diff --git a/src/amdis/LocalAssemblerBase.hpp b/src/amdis/LocalAssemblerBase.hpp
index 8ea0af9b1dc1d9bb082d0cd199d4b7d7a9d628c9..f6320729f7ddfc14d37557fa8635e03908948379 100644
--- a/src/amdis/LocalAssemblerBase.hpp
+++ b/src/amdis/LocalAssemblerBase.hpp
@@ -1,135 +1,46 @@
#pragma once
-#include <list>
-#include <memory>
#include <type_traits>
-#include <amdis/common/ConceptsBase.hpp>
+#include <amdis/ContextGeometry.hpp>
+// #include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/TypeDefs.hpp>
-#include <amdis/utility/TreeData.hpp>
namespace AMDiS
{
- namespace Impl
- {
- template <class E, class = Void_t<>>
- struct ContextImpl
- {
- using Entity = E;
- using Geometry = typename E::Geometry;
- };
-
- // specialization for intersections
- template <class I>
- struct ContextImpl<I, Void_t<decltype(std::declval<I>().inside())>>
- {
- using Entity = typename I::Entity;
- using Geometry = typename I::LocalGeometry;
- };
-
- } // end namespace Impl
-
-
/// Abstract base-class of a \ref LocalAssembler
template <class LocalContext, class... Nodes>
class LocalAssemblerBase
{
public:
-
- using Element = typename Impl::ContextImpl<LocalContext>::Entity;
+ /// The codim=0 grid entity
+ using Element = typename Impl::ContextTypes<LocalContext>::Entity;
+ /// The geometry of the \ref Element
using Geometry = typename Element::Geometry;
- using LocalGeometry = typename Impl::ContextImpl<LocalContext>::Geometry;
static constexpr int numNodes = sizeof...(Nodes);
static_assert( numNodes == 1 || numNodes == 2,
"VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
+ /// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
using ElementMatrixVector = std::conditional_t<
(sizeof...(Nodes)==1), Impl::ElementVector, std::conditional_t<
(sizeof...(Nodes)==2), Impl::ElementMatrix, void>>;
-
public:
-
+ /// Virtual destructor
virtual ~LocalAssemblerBase() {}
+ /// Bind the local-assembler to the grid-element with its corresponding geometry
virtual void bind(Element const& element, Geometry const& geometry) = 0;
+
+ /// Unbind from the element
virtual void unbind() = 0;
/// Assemble an element matrix or element vector on the test- (and trial-) function node(s)
- virtual bool assemble(LocalContext const& localContext,
- ElementMatrixVector& elementMatrixVector,
- Nodes const&... nodes) = 0;
- };
-
-
- template <class GridView>
- class OperatorLists
- {
- using Element = typename GridView::template Codim<0>::Entity;
-
- template <class OperatorType>
- struct Scaled
- {
- std::shared_ptr<OperatorType> op;
- double* factor = nullptr;
- double* estFactor = nullptr;
- BoundaryType b = {0};
- };
-
- template <class... Nodes>
- struct Data
- {
- using ElementOperator = LocalAssemblerBase<Element, Nodes...>;
- using IntersectionOperator = LocalAssemblerBase<typename GridView::Intersection, Nodes...>;
-
- std::list<Scaled<ElementOperator>> element;
- std::list<Scaled<IntersectionOperator>> boundary;
- std::list<Scaled<IntersectionOperator>> intersection;
-
- bool assembled = false; // if false, do reassemble
- bool changing = false; // if true, or assembled false, do reassemble
-
- bool empty() const
- {
- return element.empty() && boundary.empty() && intersection.empty();
- }
-
- bool assemble(bool flag) const
- {
- return flag && (!assembled || changing);
- }
-
- template <class Geo>
- void bind(Element const& elem, Geo const& geo)
- {
- for (auto& scaled : element) scaled.op->bind(elem,geo);
- for (auto& scaled : boundary) scaled.op->bind(elem,geo);
- for (auto& scaled : intersection) scaled.op->bind(elem,geo);
- }
-
- void unbind()
- {
- for (auto& scaled : element) scaled.op->unbind();
- for (auto& scaled : boundary) scaled.op->unbind();
- for (auto& scaled : intersection) scaled.op->unbind();
- }
- };
-
- public:
-
- template <class RowNode, class ColNode>
- using MatData = Data<RowNode, ColNode>;
-
- template <class Node>
- using VecData = Data<Node>;
+ virtual void assemble(LocalContext const& localContext,
+ Nodes const&... nodes,
+ ElementMatrixVector& elementMatrixVector) = 0;
};
-
- template <class GlobalBasis>
- using MatrixOperators = MatrixData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template MatData>;
-
- template <class GlobalBasis>
- using VectorOperators = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
-
} // end namespace AMDiS
diff --git a/src/amdis/LocalAssemblerList.hpp b/src/amdis/LocalAssemblerList.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e4b1dc796a40c0fb90a6aaea7f48160a599de73c
--- /dev/null
+++ b/src/amdis/LocalAssemblerList.hpp
@@ -0,0 +1,96 @@
+#pragma once
+
+#include <list>
+#include <memory>
+
+#include <amdis/LocalAssemblerBase.hpp>
+#include <amdis/utility/TreeData.hpp>
+
+namespace AMDiS
+{
+ template <class GridView>
+ class OperatorLists
+ {
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
+
+ template <class OperatorType>
+ struct DataElement
+ {
+ std::shared_ptr<OperatorType> op;
+ double* factor = nullptr;
+ double* estFactor = nullptr;
+ BoundaryType b = {0};
+ };
+
+ /// Lists of \ref DataElement on the Element, BoundaryIntersction, and InteriorIntersections
+ template <class... Nodes>
+ struct Data
+ {
+ /// The type of local operators associated with grid elements
+ using ElementOperator = LocalAssemblerBase<Element, Nodes...>;
+ /// The type of local operators associated with grid intersections
+ using IntersectionOperator = LocalAssemblerBase<Intersection, Nodes...>;
+
+ /// Return whether any operators are stored on the node
+ bool empty() const
+ {
+ return element.empty() && boundary.empty() && intersection.empty();
+ }
+
+ /// Test whether to assemble on the node
+ bool doAssemble(bool flag) const
+ {
+ return flag && (!assembled || changing);
+ }
+
+ /// Bind all operators to the grid element and geometry
+ template <class Geo>
+ void bind(Element const& elem, Geo const& geo)
+ {
+ for (auto& scaled : element) scaled.op->bind(elem,geo);
+ for (auto& scaled : boundary) scaled.op->bind(elem,geo);
+ for (auto& scaled : intersection) scaled.op->bind(elem,geo);
+ }
+
+ /// Unbind all operators from the element
+ void unbind()
+ {
+ for (auto& scaled : element) scaled.op->unbind();
+ for (auto& scaled : boundary) scaled.op->unbind();
+ for (auto& scaled : intersection) scaled.op->unbind();
+ }
+
+ /// List of operators to be assembled on grid elements
+ std::list<DataElement<ElementOperator>> element;
+ /// List of operators to be assembled on boundary intersections
+ std::list<DataElement<IntersectionOperator>> boundary;
+ /// List of operators to be assembled on interior intersections
+ std::list<DataElement<IntersectionOperator>> intersection;
+
+ /// if false, do reassemble of all operators
+ bool assembled = false;
+
+ /// if true, or \ref assembled false, do reassemble of all operators
+ bool changing = false;
+ };
+
+ public:
+
+ /// List of operators associated with matrix blocks (RowNode, ColNode)
+ template <class RowNode, class ColNode>
+ using MatData = Data<RowNode, ColNode>;
+
+ /// List of operators associated with vector blocks [Node]
+ template <class Node>
+ using VecData = Data<Node>;
+ };
+
+
+ template <class GlobalBasis>
+ using MatrixOperators = MatrixData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template MatData>;
+
+ template <class GlobalBasis>
+ using VectorOperators = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;
+
+} // end namespace AMDiS
diff --git a/src/amdis/LocalOperator.hpp b/src/amdis/LocalOperator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5509510a9ce52afe2312f8765993c6d2d8276600
--- /dev/null
+++ b/src/amdis/LocalOperator.hpp
@@ -0,0 +1,206 @@
+#pragma once
+
+#include <cassert>
+#include <type_traits>
+
+#include <amdis/GridFunctions.hpp>
+#include <amdis/common/Utility.hpp>
+#include <amdis/utility/FiniteElementType.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \addtogroup operators
+ * @{
+ **/
+
+ /// \brief The main implementation of an operator to be used in a \ref LocalAssembler.
+ /**
+ * The CRTP Base class for local operators.
+ *
+ * \tparam Derived The class that derives from this base class
+ * \tparam LocalContextType The type of the element or intersection the operator is evaluated on
+ **/
+ template <class Derived, class LocalContextType>
+ class LocalOperator
+ {
+ public:
+ /// The element or intersection the operator is assembled on
+ using LocalContext = LocalContextType;
+ /// The codim=0 grid entity
+ using Element = typename Impl::ContextTypes<LocalContext>::Entity;
+ /// The geometry of the \ref Element
+ using Geometry = typename Element::Geometry;
+
+ /// Initialize the local operator on the current gridView
+ template <class GridView>
+ void init(GridView const& gridView)
+ {
+ derived().init_impl(gridView);
+ }
+
+ /// \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`.
+ **/
+ void bind(Element const& element, Geometry const& geometry)
+ {
+ if (!bound_) {
+ isSimplex_ = geometry.type().isSimplex();
+ isAffine_ = geometry.affine();
+ bound_ = true;
+ }
+
+ derived().bind_impl(element, geometry);
+ }
+
+ /// Unbinds operator from element.
+ void unbind()
+ {
+ derived().unbind_impl();
+ bound_ = false;
+ }
+
+ /// \brief Assemble a local element matrix on the element that is bound.
+ /**
+ * \param context Container for geometry related data, \ref ContextGeometry
+ * \param rowNode The node of the test-basis
+ * \param colNode The node of the trial-basis
+ * \param elementMatrix The output element matrix
+ **/
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void calculateElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
+ {
+ derived().getElementMatrix(context, rowNode, colNode, elementMatrix);
+ }
+
+ /// \brief Assemble a local element vector on the element that is bound.
+ /**
+ * \param context Container for geometry related data \ref ContextGeometry
+ * \param node The node of the test-basis
+ * \param elementVector The output element vector
+ **/
+ template <class Context, class Node, class ElementVector>
+ void calculateElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
+ {
+ derived().getElementVector(context, node, elementVector);
+ }
+
+
+ Derived & derived() { return static_cast<Derived&>(*this); }
+ Derived const& derived() const { return static_cast<Derived const&>(*this); }
+
+
+ protected:
+
+ // default implementation. Can be overridden in the derived classes
+ template <class GridView>
+ void init_impl(GridView const& /*gridView*/) {}
+
+ // default implementation. Can be overridden in the derived classes
+ template <class Element, class Geometry>
+ void bind_impl(Element const& /*element*/, Geometry const& /*geometry*/) {}
+
+ // default implementation. Can be overridden in the derived classes
+ void unbind_impl() {}
+
+ // default implementation called by \ref calculateElementMatrix
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& /*context*/,
+ RowNode const& /*rowNode*/,
+ ColNode const& /*colNode*/,
+ ElementMatrix& /*elementMatrix*/)
+ {
+ error_exit("Needs to be implemented by derived class!");
+ }
+
+ // default implementation called by \ref calculateElementVector
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& /*context*/,
+ Node const& /*node*/,
+ ElementVector& /*elementVector*/)
+ {
+ error_exit("Needs to be implemented by derived class!");
+ }
+
+ /// \brief Return the quadrature degree for a matrix 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 RowNode, class ColNode>
+ int getDegree(int order, int coeffDegree, RowNode const& rowNode, ColNode const& colNode) const
+ {
+ assert( bound_ );
+ test_warning(coeffDegree >= 0,
+ "polynomial order of coefficient function not determined. Use order 4 by default.");
+
+ int psiDegree = polynomialDegree(rowNode);
+ int phiDegree = polynomialDegree(colNode);
+
+ int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
+ if (isSimplex_)
+ degree -= order;
+ if (isAffine_)
+ degree += 1; // oder += (order+1)
+
+ return degree;
+ }
+
+ /// \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 order, int coeffDegree, Node const& node) const
+ {
+ assert( bound_ );
+ test_warning(coeffDegree >= 0,
+ "polynomial order of coefficient function not determined. Use order 4 by default.");
+
+ int psiDegree = polynomialDegree(node);
+
+ int degree = psiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
+ if (isSimplex_)
+ degree -= order;
+ if (isAffine_)
+ degree += 1; // oder += (order+1)
+
+ return degree;
+ }
+
+ protected:
+
+ bool isSimplex_ = false; //< the bound element is a simplex
+ bool isAffine_ = false; //< the bound geometry is affine
+ bool bound_ = false; //< an element is bound to the operator
+ };
+
+
+ /// Generate a \ref LocalOperator from a PreOperator.
+ template <class Derived, class LocalContext, class GridView>
+ auto makeLocalOperator(LocalOperator<Derived, LocalContext> const& localOp, GridView const& /*gridView*/)
+ {
+ return localOp.derived();
+ }
+
+ /// Generate a shared_ptr to a \ref LocalOperator from a PreOperator.
+ template <class Derived, class LocalContext, class GridView>
+ auto makeLocalOperatorPtr(LocalOperator<Derived, LocalContext> const& localOp, GridView const& /*gridView*/)
+ {
+ return std::make_shared<Derived>(localOp.derived());
+ }
+
+} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 9352089e83f137088cde0b6136cc5d3922baf007..47e5a63684e9b3e17221be5ce6e2f43e1c70caa1 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -19,6 +19,7 @@
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
+#include <amdis/LocalAssemblerList.hpp>
#include <amdis/Marker.hpp>
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 263b5db3d35308019eb2b2eb6146853381c79815..4f08455f00d4d64a853b81d0f33fe518c3cbd767 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -176,7 +176,7 @@ void ProblemStat<Traits>::addMatrixOperator(
auto i = child(localView_->tree(), makeTreePath(row));
auto j = child(localView_->tree(), makeTreePath(col));
- auto op = makeGridOperator(preOp, globalBasis_->gridView());
+ auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
@@ -201,7 +201,7 @@ void ProblemStat<Traits>::addMatrixOperator(
auto j = child(localView_->tree(), makeTreePath(col));
using Intersection = typename GridView::Intersection;
- auto op = makeGridOperator(preOp, globalBasis_->gridView());
+ auto op = makeLocalOperator<Intersection>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
@@ -221,7 +221,7 @@ void ProblemStat<Traits>::addVectorOperator(
auto i = child(localView_->tree(), makeTreePath(path));
- auto op = makeGridOperator(preOp, globalBasis_->gridView());
+ auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
@@ -243,7 +243,7 @@ void ProblemStat<Traits>::addVectorOperator(
auto i = child(localView_->tree(), makeTreePath(path));
using Intersection = typename GridView::Intersection;
- auto op = makeGridOperator(preOp, globalBasis_->gridView());
+ auto op = makeLocalOperator<Intersection>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
diff --git a/src/amdis/QuadratureFactory.hpp b/src/amdis/QuadratureFactory.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2de45d30e36863b6e87d45ad18a62c2d04c856e0
--- /dev/null
+++ b/src/amdis/QuadratureFactory.hpp
@@ -0,0 +1,197 @@
+#pragma once
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/std/type_traits.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+
+#include <amdis/GridFunctions.hpp>
+
+namespace AMDiS
+{
+ /// Base class for quadrature factories for localFunctions
+ template <class ctype, int dimension, class LocalFunction>
+ class QuadratureFactory
+ {
+ protected:
+ using QuadratureRule = Dune::QuadratureRule<ctype, dimension>;
+
+ public:
+ virtual ~QuadratureFactory() {}
+
+ /// Bind the rule to a localFunction
+ virtual void bind(LocalFunction const& localFct)
+ {}
+
+ /// Return the quadrature order of the coefficient relates to the localFunction
+ virtual int order() const
+ {
+ return 0;
+ }
+
+ /// Construct a quadrature rule, based on the local geometry-type and the
+ /// polynomial degree of the integrand function.
+ virtual QuadratureRule const& rule(Dune::GeometryType const& type, int degree) const
+ {
+ using QuadratureRules = Dune::QuadratureRules<ctype, dimension>;
+ return QuadratureRules::rule(type, degree);
+ }
+ };
+
+
+ template <class LF>
+ using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );
+
+
+ /// \brief Factory for quadrature rule, that calculates the coefficient order from
+ /// a localFunction passed to the bind method.
+ template <class ctype, int dimension, class LocalFunction>
+ class QuadFactoryFromLocalFunction
+ : public QuadratureFactory<ctype, dimension, LocalFunction>
+ {
+ using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>;
+ static_assert(Concept::value,
+ "Polynomial order of GridFunction can not be extracted. Provide an explicit order parameter instead.");
+
+ public:
+ virtual void bind(LocalFunction const& localFct) final
+ {
+ order_ = Dune::Hybrid::ifElse(Concept{},
+ [&](auto id) { return AMDiS::order(id(localFct)); },
+ [] (auto) { return -1; });
+ }
+
+ virtual int order() const final { return order_; }
+
+ private:
+ int order_ = -1;
+ };
+
+ struct PreQuadFactoryFromLocalFunction {};
+
+ auto makePreQuadratureFactory()
+ {
+ return PreQuadFactoryFromLocalFunction{};
+ }
+
+ /// Generator for \ref QuadFactoryFromLocalFunction. \relates QuadFactoryFromLocalFunction
+ template <class ctype, int dimension, class LocalFunction>
+ auto makeQuadratureFactory(PreQuadFactoryFromLocalFunction const& /*pre*/)
+ {
+ return QuadFactoryFromLocalFunction<ctype, dimension, LocalFunction>{};
+ }
+
+ template <class ctype, int dimension, class LocalFunction>
+ auto makeQuadratureFactoryPtr(PreQuadFactoryFromLocalFunction const& /*pre*/)
+ {
+ using Factory = QuadFactoryFromLocalFunction<ctype, dimension, LocalFunction>;
+ return std::shared_ptr<Factory>(new Factory{});
+ }
+
+
+
+ /// \brief Factory for quadrature rule, that takes to order of the localFunction
+ /// and a quadrature-type
+ template <class ctype, int dimension, class LocalFunction>
+ class QuadFactoryFromOrder
+ : public QuadratureFactory<ctype, dimension, LocalFunction>
+ {
+ using Super = QuadratureFactory<ctype, dimension, LocalFunction>;
+ using QuadratureRule = typename Super::QuadratureRule;
+
+ public:
+ /// Constructor. Takes the order of the localCoefficient and a quadrature-type
+ QuadFactoryFromOrder(int order, Dune::QuadratureType::Enum qt)
+ : order_(order)
+ , qt_(qt)
+ {}
+
+ virtual int order() const final { return order_; }
+
+ virtual QuadratureRule const& rule(Dune::GeometryType const& type, int degree) const final
+ {
+ using QuadratureRules = Dune::QuadratureRules<ctype, dimension>;
+ return QuadratureRules::rule(type, degree, qt_);
+ }
+
+ private:
+ int order_;
+ Dune::QuadratureType::Enum qt_;
+ };
+
+ struct PreQuadFactoryFromOrder
+ {
+ int order;
+ Dune::QuadratureType::Enum qt;
+ };
+
+ auto makePreQuadratureFactory(int order, Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
+ {
+ return PreQuadFactoryFromOrder{order, qt};
+ }
+
+ /// Generator for \ref QuadFactoryFromOrder. \relates QuadFactoryFromOrder
+ template <class ctype, int dimension, class LocalFunction>
+ auto makeQuadratureFactory(PreQuadFactoryFromOrder const& pre)
+ {
+ return QuadFactoryFromOrder<ctype, dimension, LocalFunction>{pre.order, pre.qt};
+ }
+
+ template <class ctype, int dimension, class LocalFunction>
+ auto makeQuadratureFactoryPtr(PreQuadFactoryFromOrder const& pre)
+ {
+ using Factory = QuadFactoryFromOrder<ctype, dimension, LocalFunction>;
+ return std::make_shared<Factory>(pre.order, pre.qt);
+ }
+
+
+
+ /// \brief Factory for quadrature rule, that is based on an existing rule
+ template <class ctype, int dimension, class LocalFunction>
+ class QuadFactoryFromRule
+ : public QuadratureFactory<ctype, dimension, LocalFunction>
+ {
+ using Super = QuadratureFactory<ctype, dimension, LocalFunction>;
+ using QuadratureRule = typename Super::QuadratureRule;
+
+ public:
+ QuadFactoryFromRule(QuadratureRule const& rule)
+ : rule_(rule)
+ {}
+
+ virtual QuadratureRule const& rule(Dune::GeometryType const& /*type*/, int /*degree*/) const final
+ {
+ return rule_;
+ }
+
+ private:
+ QuadratureRule const& rule_;
+ };
+
+ template <class QuadRule>
+ struct PreQuadFactoryFromRule
+ {
+ QuadRule const& rule;
+ };
+
+ template <class ctype, int dimension>
+ auto makePreQuadratureFactory(Dune::QuadratureRule<ctype, dimension> const& rule)
+ {
+ return PreQuadFactoryFromRule<Dune::QuadratureRule<ctype, dimension>>{rule};
+ }
+
+ /// Generator for \ref QuadFactoryFromRule. \relates QuadFactoryFromRule
+ template <class ctype, int dimension, class LocalFunction, class QuadRule>
+ auto makeQuadratureFactory(PreQuadFactoryFromRule<QuadRule> const& pre)
+ {
+ return QuadFactoryFromRule<ctype, dimension, LocalFunction>{pre.rule};
+ }
+
+ template <class ctype, int dimension, class LocalFunction, class QuadRule>
+ auto makeQuadratureFactoryPtr(PreQuadFactoryFromRule<QuadRule> const& pre)
+ {
+ using Factory = QuadFactoryFromRule<ctype, dimension, LocalFunction>;
+ return std::make_shared<Factory>(pre.rule);
+ }
+
+} // end namespace AMDiS
diff --git a/src/amdis/assembler/ConvectionDiffusionOperator.hpp b/src/amdis/assembler/ConvectionDiffusionOperator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a2f938e7a718a1c353d2052097b8beefb0f0c232
--- /dev/null
+++ b/src/amdis/assembler/ConvectionDiffusionOperator.hpp
@@ -0,0 +1,231 @@
+#pragma once
+
+#include <type_traits>
+
+#include <amdis/LocalOperator.hpp>
+#include <amdis/Output.hpp>
+#include <amdis/common/Utility.hpp>
+#include <amdis/common/ValueCategory.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \addtogroup operators
+ * @{
+ **/
+
+ /// convection-diffusion operator
+ /// <A*grad(u),grad(v)> - <b*u, grad(v)> + <c*u, v> = <f, v> (conserving) or
+ /// <A*grad(u),grad(v)> + <b*grad(u), v> + <c*u, v> = <f, v> (non conserving)
+ template <class LocalContext, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true>
+ class ConvectionDiffusionOperator
+ : public LocalOperator<ConvectionDiffusionOperator<LocalContext, GridFctA, GridFctB, GridFctC, GridFctF, conserving>,
+ LocalContext>
+ {
+ using A_range_type = typename GridFctA::Range;
+ static_assert( Category::Scalar<A_range_type> || Category::Matrix<A_range_type>,
+ "Expression A(x) must be of scalar or matrix type." );
+ using b_range_type = typename GridFctB::Range;
+ static_assert( Category::Scalar<b_range_type> || Category::Vector<b_range_type>,
+ "Expression b(x) must be of scalar or vector type." );
+ using c_range_type = typename GridFctC::Range;
+ static_assert( Category::Scalar<c_range_type>,
+ "Expression c(x) must be of scalar type." );
+ using f_range_type = typename GridFctF::Range;
+ static_assert( Category::Scalar<f_range_type>,
+ "Expression f(x) must be of scalar type." );
+
+ public:
+ ConvectionDiffusionOperator(GridFctA const& gridFctA, GridFctB const& gridFctB,
+ GridFctC const& gridFctC, GridFctF const& gridFctF)
+ : gridFctA_(gridFctA)
+ , gridFctB_(gridFctB)
+ , gridFctC_(gridFctC)
+ , gridFctF_(gridFctF)
+ {}
+
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode, ColNode const& colNode,
+ ElementMatrix& elementMatrix)
+ {
+ static_assert(RowNode::isLeaf && ColNode::isLeaf,
+ "Operator can be applied to Leaf-Nodes only." );
+
+ static_assert(std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{},
+ "Galerkin operator requires equal finite elements for test and trial space." );
+
+ using LocalBasisType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType;
+
+ using RangeType = typename LocalBasisType::Traits::RangeType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename LocalBasisType::Traits::JacobianType;
+
+ auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element());
+ auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element());
+ auto localFctC = localFunction(gridFctC_); localFctC.bind(context.element());
+
+ auto const& localFE = rowNode.finiteElement();
+
+ int quadDeg = std::max({this->getDegree(2,coeffOrder(localFctA),rowNode,colNode),
+ this->getDegree(1,coeffOrder(localFctB),rowNode,colNode),
+ this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)});
+
+ using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>;
+ auto const& quad = QuadratureRules::rule(context.type(), quadDeg);
+
+ for (auto const& qp : quad) {
+ // Position of the current quadrature point in the reference element
+ decltype(auto) local = context.local(qp.position());
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+
+ // The multiplicative factor in the integral transformation formula
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
+
+ // the values of the shape functions on the reference element at the quadrature point
+ thread_local std::vector<RangeType> shapeValues;
+ localFE.localBasis().evaluateFunction(local, shapeValues);
+
+ // The gradients of the shape functions on the reference element
+ thread_local std::vector<JacobianType> shapeGradients;
+ localFE.localBasis().evaluateJacobian(local, shapeGradients);
+
+ // Compute the shape function gradients on the real element
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
+
+ for (std::size_t i = 0; i < gradients.size(); ++i)
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
+
+ const auto A = localFctA(local);
+ const auto b = localFctB(local);
+ const auto c = localFctC(local);
+
+ IF_CONSTEXPR(conserving) {
+ for (std::size_t i = 0; i < localFE.size(); ++i) {
+ const int local_i = rowNode.localIndex(i);
+ const auto gradAi = A * gradients[i];
+ const auto gradBi = b * gradients[i];
+ for (std::size_t j = 0; j < localFE.size(); ++j) {
+ const int local_j = colNode.localIndex(j);
+ elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
+ }
+ }
+ } else {
+ for (std::size_t i = 0; i < localFE.size(); ++i) {
+ const int local_i = rowNode.localIndex(i);
+ const auto grad_i = A * gradients[i];
+ grad_i += b * shapeValues[i];
+ for (std::size_t j = 0; j < localFE.size(); ++j) {
+ const int local_j = colNode.localIndex(j);
+ elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor;
+ }
+ }
+ }
+ }
+
+ localFctA.unbind();
+ localFctB.unbind();
+ localFctC.unbind();
+ }
+
+
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
+ {
+ static_assert(Node::isLeaf
+ "Operator can be applied to Leaf-Nodes only." );
+
+ using RangeType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType::Traits::RangeType;
+
+ auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());
+
+ auto const& localFE = node.finiteElement();
+
+ int quad_order = this->getDegree(0,coeffOrder(localFctF),node);
+
+ using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>;
+ auto const& quad = QuadratureRules::rule(context.type(), quad_order);
+
+ for (auto const& qp : quad) {
+ // Position of the current quadrature point in the reference element
+ decltype(auto) local = context.local(qp.position());
+
+ // The multiplicative factor in the integral transformation formula
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
+
+ // the values of the shape functions on the reference element at the quadrature point
+ std::vector<RangeType> shapeValues;
+ localFE.localBasis().evaluateFunction(local, shapeValues);
+
+ const auto f = localFctF(local);
+
+ for (std::size_t i = 0; i < localFE.size(); ++i) {
+ const int local_i = rowNode.localIndex(i);
+ elementVector[local_i][local_j] += f * shapeValues[i] * factor;
+ }
+ }
+
+ localFctF.unbind();
+ }
+
+ private:
+
+ template <class LF>
+ using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );
+
+ template <class LocalFct>
+ int coeffOrder(LocalFct const& localFct)
+ {
+ using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>;
+ return Dune::Hybrid::ifElse(Concept{},
+ [&](auto id) { return order(id(localFct)); },
+ [] (auto) { return 0; });
+ }
+
+ private:
+ GridFctA gridFctA;
+ GridFctB gridFctB;
+ GridFctC gridFctC;
+ GridFctF gridFctF;
+ };
+
+ template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, bool conserving>
+ struct PreConvectionDiffusionOperator
+ {
+ PreGridFctA gridFctA;
+ PreGridFctB gridFctB;
+ PreGridFctC gridFctC;
+ PreGridFctF gridFctF;
+ };
+
+ template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, bool conserving = true>
+ auto convectionDiffusion(PreGridFctA const& gridFctA, PreGridFctB const& gridFctB,
+ PreGridFctC const& gridFctC, PreGridFctF const& gridFctF,
+ bool_t<conserving> = {})
+ {
+ using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, conserving>;
+ return Pre{gridFctA, gridFctB, gridFctC, gridFctF};
+ }
+
+ template <class LocalContext, class... GrdFcts, bool conserving, class GridView>
+ auto makeLocalOperator(PreConvectionDiffusionOperator<GridFcts..., conserving> const& pre, GridView const& gridView)
+ {
+ auto gridFctA = makeGridFunction(pre.gridFctA, gridView);
+ auto gridFctB = makeGridFunction(pre.gridFctB, gridView);
+ auto gridFctC = makeGridFunction(pre.gridFctC, gridView);
+ auto gridFctF = makeGridFunction(pre.gridFctF, gridView);
+
+ using GridFctOp = ConvectionDiffusionOperator<LocalContext,
+ decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), conserving>;
+
+ GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF};
+ return localOperator;
+ }
+
+ /** @} **/
+
+} // end namespace AMDiS
diff --git a/src/amdis/assembler/FirstOrderDivTestvecTrial.hpp b/src/amdis/assembler/FirstOrderDivTestvecTrial.hpp
index ecdcc7377f3b5094bd67d5fb738e97cc257029cb..dc49905dec8c8242c1a80952be2bb022275675b4 100644
--- a/src/amdis/assembler/FirstOrderDivTestvecTrial.hpp
+++ b/src/amdis/assembler/FirstOrderDivTestvecTrial.hpp
@@ -18,34 +18,19 @@ namespace AMDiS
/// first-order operator \f$ \langle\nabla\cdot\Psi, c\,\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::divtestvec_trial, GridFct, QuadCreator>
- : public GridFunctionOperator<tag::test_divtrialvec, GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::divtestvec_trial, LocalContext, GridFct>
+ : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::divtestvec_trial, LocalContext, GridFct>,
+ GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>>
{
- using Transposed = GridFunctionOperator<tag::test_divtrialvec, GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Transposed = GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>;
+ using Super = GridFunctionOperatorTransposed<Self, Transposed>;
public:
- GridFunctionOperator(tag::divtestvec_trial, GridFct const& expr, QuadCreator const& quadCreator)
- : Transposed(tag::test_divtrialvec{}, expr, quadCreator)
+ GridFunctionOperator(tag::divtestvec_trial, GridFct const& expr)
+ : Super(tag::test_divtrialvec{}, expr)
{}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE> flagSameFE,
- std::integral_constant<bool, sameNode> flagSameNode)
- {
- static_assert(RowNode::isPower && ColNode::isLeaf,
- "ColNode must be Leaf-Node and RowNode must be a Power-Node.");
-
- auto elementMatrixTransposed = trans(elementMatrix);
- Transposed::calculateElementMatrix(
- context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
- }
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderGradTestTrial.hpp b/src/amdis/assembler/FirstOrderGradTestTrial.hpp
index 1ca7e92a1d81e7651ebb98657170df40933e495f..6808103e4b6013135f7c278caed611dea2776dcc 100644
--- a/src/amdis/assembler/FirstOrderGradTestTrial.hpp
+++ b/src/amdis/assembler/FirstOrderGradTestTrial.hpp
@@ -18,31 +18,19 @@ namespace AMDiS
/// first-order operator \f$ \langle\nabla\psi, \mathbf{b}\,\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::gradtest_trial, GridFct, QuadCreator>
- : public GridFunctionOperator<tag::test_gradtrial, GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::gradtest_trial, LocalContext, GridFct>
+ : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::gradtest_trial, LocalContext, GridFct>,
+ GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>>
{
- using Transposed = GridFunctionOperator<tag::test_gradtrial, GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Transposed = GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>;
+ using Super = GridFunctionOperatorTransposed<Self, Transposed>;
public:
- GridFunctionOperator(tag::gradtest_trial, GridFct const& expr, QuadCreator const& quadCreator)
- : Transposed(tag::test_gradtrial{}, expr, quadCreator)
+ GridFunctionOperator(tag::gradtest_trial, GridFct const& expr)
+ : Super(tag::test_gradtrial{}, expr)
{}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE> flagSameFE,
- std::integral_constant<bool, sameNode> flagSameNode)
- {
- auto elementMatrixTransposed = trans(elementMatrix);
- Transposed::calculateElementMatrix(
- context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
- }
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderGradTestTrialvec.hpp b/src/amdis/assembler/FirstOrderGradTestTrialvec.hpp
index 4d270a11bdbdc14dca7fe3bc8774435ae5228ca4..83c6c9e86c8eb9dceef7638ca56d3f612eae309e 100644
--- a/src/amdis/assembler/FirstOrderGradTestTrialvec.hpp
+++ b/src/amdis/assembler/FirstOrderGradTestTrialvec.hpp
@@ -18,31 +18,19 @@ namespace AMDiS
/// first-order operator \f$ \langle\nabla\psi, c\,\Phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::gradtest_trialvec, GridFct, QuadCreator>
- : public GridFunctionOperator<tag::testvec_gradtrial, GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::gradtest_trialvec, LocalContext, GridFct>
+ : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::gradtest_trialvec, LocalContext, GridFct>,
+ GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>>
{
- using Transposed = GridFunctionOperator<tag::testvec_gradtrial, GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Transposed = GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>;
+ using Super = GridFunctionOperatorTransposed<Self, Transposed>;
public:
- GridFunctionOperator(tag::gradtest_trialvec, GridFct const& expr, QuadCreator const& quadCreator)
- : Transposed(tag::testvec_gradtrial{}, expr, quadCreator)
+ GridFunctionOperator(tag::gradtest_trialvec, GridFct const& expr)
+ : Super(tag::testvec_gradtrial{}, expr)
{}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE> flagSameFE,
- std::integral_constant<bool, sameNode> flagSameNode)
- {
- auto elementMatrixTransposed = trans(elementMatrix);
- Transposed::calculateElementMatrix(
- context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
- }
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderPartialTestTrial.hpp b/src/amdis/assembler/FirstOrderPartialTestTrial.hpp
index a37fc0437933a01237f3be7ee5d712965b48e0ae..0c9435f8e40c0d80328ab42a914fd1fe6f84273b 100644
--- a/src/amdis/assembler/FirstOrderPartialTestTrial.hpp
+++ b/src/amdis/assembler/FirstOrderPartialTestTrial.hpp
@@ -21,31 +21,19 @@ namespace AMDiS
/// first-order operator \f$ \langle\partial_i\psi, c\,\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::partialtest_trial, GridFct, QuadCreator>
- : public GridFunctionOperator<tag::test_partialtrial, GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::partialtest_trial, LocalContext, GridFct>
+ : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::partialtest_trial, LocalContext, GridFct>,
+ GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>>
{
- using Transposed = GridFunctionOperator<tag::test_partialtrial, GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Transposed = GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>;
+ using Super = GridFunctionOperatorTransposed<Self, Transposed>;
public:
- GridFunctionOperator(tag::partialtest_trial tag, GridFct const& expr, QuadCreator const& quadCreator)
- : Transposed(tag::test_partialtrial{tag.comp}, expr, quadCreator)
+ GridFunctionOperator(tag::partialtest_trial tag, GridFct const& expr)
+ : Super(tag::test_partialtrial{tag.comp}, expr)
{}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE> flagSameFE,
- std::integral_constant<bool, sameNode> flagSameNode)
- {
- auto elementMatrixTransposed = trans(elementMatrix);
- Transposed::calculateElementMatrix(
- context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
- }
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderTestDivTrialvec.hpp b/src/amdis/assembler/FirstOrderTestDivTrialvec.hpp
index d571babf97040af02165c1f6674f6dd1467b64a8..5cdbeab30c818ca90454eeffdba526943da80ad5 100644
--- a/src/amdis/assembler/FirstOrderTestDivTrialvec.hpp
+++ b/src/amdis/assembler/FirstOrderTestDivTrialvec.hpp
@@ -20,28 +20,26 @@ namespace AMDiS
/// first-order operator \f$ \langle\psi, c\,\nabla\cdot\Phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::test_divtrialvec, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_divtrialvec, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::test_divtrialvec, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 1)
+ GridFunctionOperator(tag::test_divtrialvec, GridFct const& expr)
+ : Super(expr, 1)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isPower,
"RowNode must be Leaf-Node and ColNode must be a Power-Node.");
@@ -52,35 +50,43 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using RangeType = typename RowLocalBasisType::Traits::RangeType;
+ using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ thread_local std::vector<RangeType> shapeValues;
+ rowLocalFE.localBasis().evaluateFunction(local, shapeValues);
// The gradients of the shape functions on the reference element
- colLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,Context::dow> > colGradients(referenceGradients.size());
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > colGradients(shapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], colGradients[i]);
+ jacobian.mv(shapeGradients[i][0], colGradients[i]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
- double value = factor * rowShapeValues_[i];
+ const auto local_i = rowNode.localIndex(i);
+ const auto value = factor * shapeValues[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_kj = colNode.child(k).localIndex(j);
+ const auto local_kj = colNode.child(k).localIndex(j);
elementMatrix[local_i][local_kj] += value * colGradients[j][k];
}
}
@@ -88,9 +94,6 @@ namespace AMDiS
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderTestGradTrial.hpp b/src/amdis/assembler/FirstOrderTestGradTrial.hpp
index a4ae5383cef98820a238a28ea2d6cc659c1ecf5e..bccf890d72ed0d0289b9df7755692603e6a1e0e3 100644
--- a/src/amdis/assembler/FirstOrderTestGradTrial.hpp
+++ b/src/amdis/assembler/FirstOrderTestGradTrial.hpp
@@ -20,28 +20,26 @@ namespace AMDiS
/// first-order operator \f$ \langle\psi, \mathbf{b}\cdot\nabla\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::test_gradtrial, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_gradtrial, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
public:
- GridFunctionOperator(tag::test_gradtrial, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 1)
+ GridFunctionOperator(tag::test_gradtrial, GridFct const& expr)
+ : Super(expr, 1)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
@@ -49,44 +47,48 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using RangeType = typename RowLocalBasisType::Traits::RangeType;
+ using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
const auto b = Super::coefficient(local);
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ thread_local std::vector<RangeType> shapeValues;
+ rowLocalFE.localBasis().evaluateFunction(local, shapeValues);
// The gradients of the shape functions on the reference element
- colLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,Context::dow> > colGradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > colGradients(shapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], colGradients[i]);
+ jacobian.mv(shapeGradients[i][0], colGradients[i]);
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- const int local_j = colNode.localIndex(j);
- double value = factor * (b * colGradients[j]);
+ const auto local_j = colNode.localIndex(j);
+ const auto value = factor * (b * colGradients[j]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
- elementMatrix[local_i][local_j] += value * rowShapeValues_[i];
+ const auto local_i = rowNode.localIndex(i);
+ elementMatrix[local_i][local_j] += value * shapeValues[i];
}
}
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderTestPartialTrial.hpp b/src/amdis/assembler/FirstOrderTestPartialTrial.hpp
index 2d900f5dbb17a585903486e9ca9275a9136c29e8..d69dcb68b094c22ae90b48a6fc41590e93015f21 100644
--- a/src/amdis/assembler/FirstOrderTestPartialTrial.hpp
+++ b/src/amdis/assembler/FirstOrderTestPartialTrial.hpp
@@ -23,29 +23,27 @@ namespace AMDiS
/// first-order operator \f$ \langle\psi, c\,\partial_i\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::test_partialtrial, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_partialtrial, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::test_partialtrial tag, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 1)
+ GridFunctionOperator(tag::test_partialtrial tag, GridFct const& expr)
+ : Super(expr, 1)
, comp_(tag.comp)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
@@ -53,48 +51,54 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
+
+ using RangeType = typename RowLocalBasisType::Traits::RangeType;
+ using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
assert(jacobian.M() == Context::dim);
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
- double c = Super::coefficient(local);
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
+ const auto c = Super::coefficient(local);
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ thread_local std::vector<RangeType> shapeValues;
+ rowLocalFE.localBasis().evaluateFunction(local, shapeValues);
// The gradients of the shape functions on the reference element
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
- colLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- std::vector<double> colPartial(referenceGradients.size());
-
+ thread_local std::vector<RangeFieldType> colPartial(shapeGradients.size());
for (std::size_t i = 0; i < colPartial.size(); ++i) {
- colPartial[i] = jacobian[comp_][0] * referenceGradients[i][0][0];
+ colPartial[i] = jacobian[comp_][0] * shapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian.M(); ++j)
- colPartial[i] += jacobian[comp_][j] * referenceGradients[i][0][j];
+ colPartial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
}
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- const int local_j = colNode.localIndex(j);
- double value = factor * (c * colPartial[j]);
+ const auto local_j = colNode.localIndex(j);
+ const auto value = factor * (c * colPartial[j]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
- elementMatrix[local_i][local_j] += value * rowShapeValues_[i];
+ const auto local_i = rowNode.localIndex(i);
+ elementMatrix[local_i][local_j] += value * shapeValues[i];
}
}
}
-
}
private:
std::size_t comp_;
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/FirstOrderTestvecGradTrial.hpp b/src/amdis/assembler/FirstOrderTestvecGradTrial.hpp
index 2bfd6d52d1bbda560b094803ea272d6325dd4de4..d06a44dbc79a2538e0933246edf81a30c1e80f98 100644
--- a/src/amdis/assembler/FirstOrderTestvecGradTrial.hpp
+++ b/src/amdis/assembler/FirstOrderTestvecGradTrial.hpp
@@ -20,28 +20,26 @@ namespace AMDiS
/// first-order operator \f$ \langle\Psi, c\,\nabla\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::testvec_gradtrial, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec_gradtrial, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::testvec_gradtrial, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 1)
+ GridFunctionOperator(tag::testvec_gradtrial, GridFct const& expr)
+ : Super(expr, 1)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isPower && ColNode::isLeaf,
"ColNode must be Leaf-Node and RowNode must be a Power-Node.");
@@ -52,35 +50,42 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.finiteElement();
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using RangeType = typename RowLocalBasisType::Traits::RangeType;
+ using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ thread_local std::vector<RangeType> shapeValues;
+ rowLocalFE.localBasis().evaluateFunction(local, shapeValues);
// The gradients of the shape functions on the reference element
- colLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,Context::dow> > colGradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > colGradients(shapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], colGradients[i]);
+ jacobian.mv(shapeGradients[i][0], colGradients[i]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- double value = factor * rowShapeValues_[i];
+ const auto value = factor * shapeValues[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- const int local_j = colNode.localIndex(j);
+ const auto local_j = colNode.localIndex(j);
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_ki = rowNode.child(k).localIndex(i);
+ const auto local_ki = rowNode.child(k).localIndex(i);
elementMatrix[local_ki][local_j] += value * colGradients[j][k];
}
}
@@ -88,9 +93,6 @@ namespace AMDiS
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp b/src/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp
index 6a4a62726ae5190576326379eabe0b658732ace1..0e841d53cfef6760ae18741161538680b1922dcf 100644
--- a/src/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp
+++ b/src/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp
@@ -20,70 +20,93 @@ namespace AMDiS
/// second-order operator \f$ \langle\nabla\cdot\Psi, c\,\nabla\cdot\Phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::divtestvec_divtrialvec, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::divtestvec_divtrialvec, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec_divtrialvec, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::divtestvec_divtrialvec, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 2)
+ GridFunctionOperator(tag::divtestvec_divtrialvec, GridFct const& expr)
+ : Super(expr, 2)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isPower && ColNode::isPower,
"Operator can be applied to Power-Nodes only.");
+ const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+ const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ if (sameFE && sameNode)
+ getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix);
+ else
+ getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix);
+ }
+
+ protected:
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixStandard(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
+ {
static const std::size_t CHILDREN = RowNode::CHILDREN;
static_assert( RowNode::CHILDREN == ColNode::CHILDREN, "" );
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > rowReferenceGradients;
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > colReferenceGradients;
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
+
+ using RowFieldType = typename RowLocalBasisType::Traits::RangeFieldType;
+ using ColFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
+ using RowJacobianType = typename RowLocalBasisType::Traits::JacobianType;
+ using ColJacobianType = typename ColLocalBasisType::Traits::JacobianType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
// The gradients of the shape functions on the reference element
- rowLocalFE.localBasis().evaluateJacobian(local, rowReferenceGradients);
- colLocalFE.localBasis().evaluateJacobian(local, colReferenceGradients);
+ thread_local std::vector<RowJacobianType> rowShapeGradients;
+ rowLocalFE.localBasis().evaluateJacobian(local, rowShapeGradients);
- // Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,Context::dow> > rowGradients(rowReferenceGradients.size());
- std::vector<Dune::FieldVector<double,Context::dow> > colGradients(colReferenceGradients.size());
+ thread_local std::vector<ColJacobianType> colShapeGradients;
+ colLocalFE.localBasis().evaluateJacobian(local, colShapeGradients);
+ // Compute the shape function gradients on the real element
+ thread_local std::vector<Dune::FieldVector<RowFieldType,Context::dow> > rowGradients(rowShapeGradients.size());
for (std::size_t i = 0; i < rowGradients.size(); ++i)
- jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
+ jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
+ thread_local std::vector<Dune::FieldVector<ColFieldType,Context::dow> > colGradients(colShapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i)
- jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
+ jacobian.mv(colShapeGradients[i][0], colGradients[i]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_ki = rowNode.child(k).localIndex(i);
- const int local_kj = colNode.child(k).localIndex(j);
+ const auto local_ki = rowNode.child(k).localIndex(i);
+ const auto local_kj = colNode.child(k).localIndex(j);
elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][k];
}
}
@@ -92,52 +115,48 @@ namespace AMDiS
}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& node,
- ColNode const& /*colNode*/,
- std::true_type /*sameFE*/,
- std::true_type /*sameNode*/)
+ template <class Context, class QuadratureRule, class Node, class ElementMatrix>
+ void getElementMatrixOptimized(Context const& context,
+ QuadratureRule const& quad,
+ Node const& node,
+ Node const& /*colNode*/,
+ ElementMatrix& elementMatrix)
{
- static_assert( std::is_same<typename RowNode::NodeTag, Dune::TypeTree::PowerNodeTag>::value &&
- std::is_same<typename ColNode::NodeTag, Dune::TypeTree::PowerNodeTag>::value, "" );
-
- static const std::size_t CHILDREN = RowNode::CHILDREN;
+ static const std::size_t CHILDREN = Node::CHILDREN;
auto const& localFE = node.child(0).finiteElement();
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename LocalBasisType::Traits::JacobianType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
// The gradients of the shape functions on the reference element
- localFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t k = 0; k < CHILDREN; ++k) {
for (std::size_t i = 0; i < localFE.size(); ++i) {
- const int local_ki = node.child(k).localIndex(i);
- double value = factor * gradients[i][k];
+ const auto local_ki = node.child(k).localIndex(i);
+ const auto value = factor * gradients[i][k];
elementMatrix[local_ki][local_ki] += value * gradients[i][k];
for (std::size_t j = i+1; j < localFE.size(); ++j) {
- const int local_kj = node.child(k).localIndex(j);
+ const auto local_kj = node.child(k).localIndex(j);
elementMatrix[local_ki][local_kj] += value * gradients[j][k];
elementMatrix[local_kj][local_ki] += value * gradients[j][k];
}
diff --git a/src/amdis/assembler/SecondOrderGradTestGradTrial.hpp b/src/amdis/assembler/SecondOrderGradTestGradTrial.hpp
index bfb3943d2eedc12cdbb29491b398e3ac190b9dc8..d8f6ea1ddda6a0418390ee8e290e1c8f56953066 100644
--- a/src/amdis/assembler/SecondOrderGradTestGradTrial.hpp
+++ b/src/amdis/assembler/SecondOrderGradTestGradTrial.hpp
@@ -20,128 +20,131 @@ namespace AMDiS
/// second-order operator \f$ \langle\nabla\psi, c\,\nabla\phi\rangle \f$, or \f$ \langle\nabla\psi, A\,\nabla\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::gradtest_gradtrial, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::gradtest_gradtrial, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest_gradtrial, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
using expr_value_type = typename GridFct::Range;
static_assert( Category::Scalar<expr_value_type> || Category::Matrix<expr_value_type>,
"Expression must be of scalar or matrix type." );
public:
- GridFunctionOperator(tag::gradtest_gradtrial, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 2)
+ GridFunctionOperator(tag::gradtest_gradtrial, GridFct const& expr)
+ : Super(expr, 2)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::false_type /*sameNode*/)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
- auto const& rowLocalFE = rowNode.finiteElement();
- auto const& colLocalFE = colNode.finiteElement();
+ const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+ const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
+ using Category = ValueCategory_t<typename GridFct::Range>;
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ if (sameFE && sameNode)
+ getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix, Category{});
+ else if (sameFE)
+ getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix);
+ else
+ error_exit("Not implemented: currently only the implementation for equal fespaces available");
+ }
+
+ protected:
+
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixStandard(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
+ {
+ auto const& localFE = rowNode.finiteElement();
- test_exit( sameFE, "Not implemented: currently only the implementation for equal fespaces available" );
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename LocalBasisType::Traits::JacobianType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
const auto exprValue = Super::coefficient(local);
// The gradients of the shape functions on the reference element
- thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
- rowLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- thread_local std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
- for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
- for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- const int local_j = colNode.localIndex(j);
+ for (std::size_t i = 0; i < localFE.size(); ++i) {
+ const auto local_i = rowNode.localIndex(i);
+ for (std::size_t j = 0; j < localFE.size(); ++j) {
+ const auto local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
}
}
}
}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& node,
- ColNode const& /*colNode*/,
- std::true_type /*sameFE*/,
- std::true_type /*sameNode*/)
- {
- static_assert(RowNode::isLeaf && ColNode::isLeaf,
- "Operator can be applied to Leaf-Nodes only.");
-
- using Category = ValueCategory_t<typename GridFct::Range>;
- calcElementMatrix(context, quad, elementMatrix, node, Category{});
- }
-
- protected:
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode>
- void calcElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& node,
- tag::scalar)
+ template <class Context, class QuadratureRule, class Node, class ElementMatrix>
+ void getElementMatrixOptimized(Context const& context,
+ QuadratureRule const& quad,
+ Node const& node,
+ Node const& /*colNode*/,
+ ElementMatrix& elementMatrix,
+ tag::scalar)
{
auto const& localFE = node.finiteElement();
- thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename LocalBasisType::Traits::JacobianType;
+
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
// The gradients of the shape functions on the reference element
- localFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- thread_local std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- const int local_i = node.localIndex(i);
+ const auto local_i = node.localIndex(i);
elementMatrix[local_i][local_i] += factor * (gradients[i] * gradients[i]);
for (std::size_t j = i+1; j < localFE.size(); ++j) {
- const int local_j = node.localIndex(j);
- double value = factor * (gradients[i] * gradients[j]);
+ const auto local_j = node.localIndex(j);
+ const auto value = factor * (gradients[i] * gradients[j]);
elementMatrix[local_i][local_j] += value;
elementMatrix[local_j][local_i] += value;
@@ -150,42 +153,44 @@ namespace AMDiS
}
}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode>
- void calcElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& node,
- tag::matrix)
+ template <class Context, class QuadratureRule, class Node, class ElementMatrix>
+ void getElementMatrixOptimized(Context const& context,
+ QuadratureRule const& quad,
+ Node const& node,
+ Node const& /*colNode*/,
+ ElementMatrix& elementMatrix,
+ tag::matrix)
{
auto const& localFE = node.finiteElement();
- thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+ using JacobianType = typename LocalBasisType::Traits::JacobianType;
+
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
const auto exprValue = Super::coefficient(local);
// The gradients of the shape functions on the reference element
- localFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<JacobianType> shapeGradients;
+ localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- thread_local std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- const int local_i = node.localIndex(i);
+ const auto local_i = node.localIndex(i);
for (std::size_t j = 0; j < localFE.size(); ++j) {
- const int local_j = node.localIndex(j);
+ const auto local_j = node.localIndex(j);
elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
}
}
@@ -194,20 +199,22 @@ namespace AMDiS
protected:
- template <int dow>
- double eval(Dune::FieldVector<double,1> const& scalar, double factor,
- Dune::FieldVector<double,dow> const& grad_test,
- Dune::FieldVector<double,dow> const& grad_trial) const
+ template <class S, class F, class T, int dow,
+ std::enable_if_t<Category::Scalar<S>,int> = 0>
+ T eval(S const& scalar, F factor,
+ Dune::FieldVector<T,dow> const& grad_test,
+ Dune::FieldVector<T,dow> const& grad_trial) const
{
return (scalar * factor) * (grad_test * grad_trial);
}
- template <int dow>
- double eval(Dune::FieldMatrix<double,dow,dow> const& M, double factor,
- Dune::FieldVector<double,dow> const& grad_test,
- Dune::FieldVector<double,dow> const& grad_trial) const
+ template <class M, class F, class T, int dow,
+ std::enable_if_t<Category::Matrix<M>,int> = 0>
+ T eval(M const& mat, F factor,
+ Dune::FieldVector<T,dow> const& grad_test,
+ Dune::FieldVector<T,dow> const& grad_trial) const
{
- return factor * (grad_test * (M * grad_trial));
+ return factor * (grad_test * (mat * grad_trial));
}
};
diff --git a/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp b/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
index f5c548add90cd01e37592c8e94de651665caa58b..6539acf8ce1ec785c1a2cc689884f5e9fb8a2295 100644
--- a/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
+++ b/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
@@ -24,30 +24,28 @@ namespace AMDiS
/// second-order operator \f$ \langle\partial_i\psi, c\,\partial_j\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::partialtest_partialtrial, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::partialtest_partialtrial, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest_partialtrial, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::partialtest_partialtrial tag, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 2)
+ GridFunctionOperator(tag::partialtest_partialtrial tag, GridFct const& expr)
+ : Super(expr, 2)
, compTest_(tag.comp_test)
, compTrial_(tag.comp_trial)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
@@ -55,46 +53,52 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > rowReferenceGradients;
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > colReferenceGradients;
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using RowFieldType = typename RowLocalBasisType::Traits::RangeFieldType;
+ using ColFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
+ using RowJacobianType = typename RowLocalBasisType::Traits::JacobianType;
+ using ColJacobianType = typename ColLocalBasisType::Traits::JacobianType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
- double c = Super::coefficient(local);
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
+ const auto c = Super::coefficient(local);
// The gradients of the shape functions on the reference element
- rowLocalFE.localBasis().evaluateJacobian(local, rowReferenceGradients);
- colLocalFE.localBasis().evaluateJacobian(local, colReferenceGradients);
+ thread_local std::vector<RowJacobianType> rowShapeGradients;
+ rowLocalFE.localBasis().evaluateJacobian(local, rowShapeGradients);
+ thread_local std::vector<ColJacobianType> colShapeGradients;
+ colLocalFE.localBasis().evaluateJacobian(local, colShapeGradients);
// Compute the shape function gradients on the real element
- std::vector<double> rowPartial(rowReferenceGradients.size());
- std::vector<double> colPartial(colReferenceGradients.size());
-
- // transform gradients to global element
+ thread_local std::vector<RowFieldType> rowPartial(rowShapeGradients.size());
for (std::size_t i = 0; i < rowPartial.size(); ++i) {
- rowPartial[i] = jacobian[compTest_][0] * rowReferenceGradients[i][0][0];
+ rowPartial[i] = jacobian[compTest_][0] * rowShapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian.cols(); ++j)
- rowPartial[i] += jacobian[compTest_][j] * rowReferenceGradients[i][0][j];
+ rowPartial[i] += jacobian[compTest_][j] * rowShapeGradients[i][0][j];
}
+ thread_local std::vector<ColFieldType> colPartial(colShapeGradients.size());
for (std::size_t i = 0; i < colPartial.size(); ++i) {
- colPartial[i] = jacobian[compTrial_][0] * colReferenceGradients[i][0][0];
+ colPartial[i] = jacobian[compTrial_][0] * colShapeGradients[i][0][0];
for (std::size_t j = 1; j < jacobian.cols(); ++j)
- colPartial[i] += jacobian[compTrial_][j] * colReferenceGradients[i][0][j];
+ colPartial[i] += jacobian[compTrial_][j] * colShapeGradients[i][0][j];
}
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- const int local_j = colNode.localIndex(j);
- double value = factor * (c * colPartial[j]);
+ const auto local_j = colNode.localIndex(j);
+ const auto value = factor * (c * colPartial[j]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
+ const auto local_i = rowNode.localIndex(i);
elementMatrix[local_i][local_j] += value * rowPartial[i];
}
}
diff --git a/src/amdis/assembler/StokesOperator.hpp b/src/amdis/assembler/StokesOperator.hpp
index 46c191a66bfd1d42c3d558a71695fe714781cde0..b9dd503273666303b696d2b6a01effdcf4220366 100644
--- a/src/amdis/assembler/StokesOperator.hpp
+++ b/src/amdis/assembler/StokesOperator.hpp
@@ -20,27 +20,26 @@ namespace AMDiS
/// Stokes operator \f$ \langle\nabla\mathbf{v}, \nu \nabla\mathbf{u}\rangle + \langle\nabla\cdot\mathbf{v}, p\rangle + \langle q, \langle\nabla\cdot\mathbf{u}\rangle \f$
- template <class ViscosityExpr, class QuadCreator>
- class GridFunctionOperator<tag::stokes, ViscosityExpr, QuadCreator>
- : public GridFunctionOperatorBase<ViscosityExpr, QuadCreator>
+ template <class LocalContext, class ViscosityExpr>
+ class GridFunctionOperator<tag::stokes, LocalContext, ViscosityExpr>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::stokes, LocalContext, ViscosityExpr>,
+ LocalContext, ViscosityExpr>
{
- using Super = GridFunctionOperatorBase<ViscosityExpr, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, ViscosityExpr>;
static_assert( Category::Scalar<typename ViscosityExpr::Range>, "Viscosity must be of scalar type." );
public:
- GridFunctionOperator(tag::stokes, ViscosityExpr const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 1)
+ GridFunctionOperator(tag::stokes, ViscosityExpr const& expr)
+ : Super(expr, 1)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class Node, class Flag1, class Flag2>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- Node const& tree,
- Node const& /*colNode*/,
- Flag1, Flag2)
+ template <class Context, class Node, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ Node const& tree,
+ Node const& /*colNode*/,
+ ElementMatrix& elementMatrix)
{
using VelocityNode = typename Node::template Child<0>::type;
using PressureNode = typename Node::template Child<1>::type;
@@ -51,44 +50,50 @@ namespace AMDiS
auto const& velocityLocalFE = tree.child(_0,0).finiteElement();
auto const& pressureLocalFE = tree.child(_1).finiteElement();
- std::vector<Dune::FieldVector<double,1>> pressureValues;
- std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+ using VelocityLocalBasisType = typename std::decay_t<decltype(velocityLocalFE)>::Traits::LocalBasisType;
+ using PressureLocalBasisType = typename std::decay_t<decltype(pressureLocalFE)>::Traits::LocalBasisType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using PressureRange = typename PressureLocalBasisType::Traits::RangeType;
+ using RangeFieldType = typename VelocityLocalBasisType::Traits::RangeFieldType;
+ using VelocityJacobian = typename VelocityLocalBasisType::Traits::JacobianType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), tree, tree);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula
- const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
- const double vel_factor = Super::coefficient(local) * factor;
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
+ const auto vel_factor = Super::coefficient(local) * factor;
// The gradients of the shape functions on the reference element
- velocityLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+ thread_local std::vector<VelocityJacobian> shapeGradients;
+ velocityLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
-
+ thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
+ thread_local std::vector<PressureRange> pressureValues;
pressureLocalFE.localBasis().evaluateFunction(local, pressureValues);
// <viscosity * grad(u), grad(v)>
for (std::size_t i = 0; i < velocityLocalFE.size(); ++i) {
- auto value_ii = vel_factor * (gradients[i] * gradients[i]);
+ const auto value_ii = vel_factor * (gradients[i] * gradients[i]);
for (std::size_t k = 0; k < Context::dow; ++k) {
- std::size_t local_ki = tree.child(_0,k).localIndex(i);
+ const auto local_ki = tree.child(_0,k).localIndex(i);
elementMatrix[local_ki][local_ki] += value_ii;
}
for (std::size_t j = i+1; j < velocityLocalFE.size(); ++j) {
- auto value = vel_factor * (gradients[i] * gradients[j]);
+ const auto value = vel_factor * (gradients[i] * gradients[j]);
for (std::size_t k = 0; k < Context::dow; ++k) {
- std::size_t local_ki = tree.child(_0,k).localIndex(i);
- std::size_t local_kj = tree.child(_0,k).localIndex(j);
+ const auto local_ki = tree.child(_0,k).localIndex(i);
+ const auto local_kj = tree.child(_0,k).localIndex(j);
elementMatrix[local_ki][local_kj] += value;
elementMatrix[local_kj][local_ki] += value;
}
@@ -98,10 +103,10 @@ namespace AMDiS
// <p, div(v)> + <div(u), q>
for (std::size_t i = 0; i < velocityLocalFE.size(); ++i) {
for (std::size_t j = 0; j < pressureLocalFE.size(); ++j) {
- auto value = factor * pressureValues[j];
+ const auto value = factor * pressureValues[j];
for (std::size_t k = 0; k < Context::dow; ++k) {
- std::size_t local_v = tree.child(_0,k).localIndex(i);
- std::size_t local_p = tree.child(_1).localIndex(j);
+ const auto local_v = tree.child(_0,k).localIndex(i);
+ const auto local_p = tree.child(_1).localIndex(j);
elementMatrix[local_v][local_p] += gradients[i][k] * value;
elementMatrix[local_p][local_v] += gradients[i][k] * value;
diff --git a/src/amdis/assembler/ZeroOrderTest.hpp b/src/amdis/assembler/ZeroOrderTest.hpp
index fd6f6786d6379430cff6b018b3df8d758c6c2588..5639f0c57c26e05bf7f862c823ba9d9d29bc22ef 100644
--- a/src/amdis/assembler/ZeroOrderTest.hpp
+++ b/src/amdis/assembler/ZeroOrderTest.hpp
@@ -19,48 +19,50 @@ namespace AMDiS
/// zero-order vector-operator \f$ (c\, \psi) \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::test, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::test, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::test, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::test, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 0)
+ GridFunctionOperator(tag::test, GridFct const& expr)
+ : Super(expr, 0)
{}
- template <class Context, class QuadratureRule,
- class ElementVector, class Node>
- void calculateElementVector(Context const& context,
- QuadratureRule const& quad,
- ElementVector& elementVector,
- Node const& node)
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
{
static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only");
auto const& localFE = node.finiteElement();
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeType = typename LocalBasisType::Traits::RangeType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), node);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
- localFE.localBasis().evaluateFunction(local, shapeValues_);
+ thread_local std::vector<RangeType> shapeValues;
+ localFE.localBasis().evaluateFunction(local, shapeValues);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- const int local_i = node.localIndex(i);
- elementVector[local_i] += factor * shapeValues_[i];
+ const auto local_i = node.localIndex(i);
+ elementVector[local_i] += factor * shapeValues[i];
}
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> shapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/ZeroOrderTestTrial.hpp b/src/amdis/assembler/ZeroOrderTestTrial.hpp
index d7df01aacf32518142e2c33126a8a22f17d46bb4..47b0e4d3676e9165604c1c7ef475f07d5f86e52c 100644
--- a/src/amdis/assembler/ZeroOrderTestTrial.hpp
+++ b/src/amdis/assembler/ZeroOrderTestTrial.hpp
@@ -19,28 +19,44 @@ namespace AMDiS
/// zero-order operator \f$ \langle\psi, c\,\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::test_trial, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::test_trial, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_trial, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Scalar<typename GridFct::Range>, "Expression must be of scalar type." );
public:
- GridFunctionOperator(tag::test_trial, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 0)
+ GridFunctionOperator(tag::test_trial, GridFct const& expr)
+ : Super(expr, 0)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::false_type /*sameNode*/)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
+ {
+ const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+ const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ if (sameFE && sameNode)
+ getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix);
+ else
+ getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix);
+ }
+
+ protected:
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixStandard(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
@@ -48,26 +64,35 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
- auto* rowShapeValues = &rowShapeValues_;
- auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_);
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
+ using RowRangeType = typename RowLocalBasisType::Traits::RangeType;
+ using ColRangeType = typename ColLocalBasisType::Traits::RangeType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ thread_local std::vector<RowRangeType> rowShapeValueBuffer;
+ thread_local std::vector<ColRangeType> colShapeValueBuffer;
+
+ const bool sameFE = std::is_same<decltype(rowLocalFE),decltype(colLocalFE)>::value;
+ auto* rowShapeValues = &rowShapeValueBuffer;
+ auto* colShapeValues = (sameFE ? &rowShapeValueBuffer : &colShapeValueBuffer);
+
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ rowLocalFE.localBasis().evaluateFunction(local, rowShapeValueBuffer);
if (!sameFE)
- colLocalFE.localBasis().evaluateFunction(local, colShapeValues_);
+ colLocalFE.localBasis().evaluateFunction(local, colShapeValueBuffer);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
- double value = factor * (*rowShapeValues)[i];
+ const auto local_i = rowNode.localIndex(i);
+ const auto value = factor * (*rowShapeValues)[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- const int local_j = colNode.localIndex(j);
+ const auto local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += value * (*colShapeValues)[j];
}
}
@@ -76,38 +101,39 @@ namespace AMDiS
template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& node,
- ColNode const& /*colNode*/,
- std::true_type /*sameFE*/,
- std::true_type /*sameNode*/)
+ class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixOptimized(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& node,
+ ColNode const& /*colNode*/,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
auto const& localFE = node.finiteElement();
- auto& shapeValues = rowShapeValues_;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeType = typename LocalBasisType::Traits::RangeType;
+
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- const double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
+ thread_local std::vector<RangeType> shapeValues;
localFE.localBasis().evaluateFunction(local, shapeValues);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- const int local_i = node.localIndex(i);
+ const auto local_i = node.localIndex(i);
- double value = factor * shapeValues[i];
+ const auto value = factor * shapeValues[i];
elementMatrix[local_i][local_i] += value * shapeValues[i];
for (std::size_t j = i+1; j < localFE.size(); ++j) {
- const int local_j = node.localIndex(j);
+ const auto local_j = node.localIndex(j);
elementMatrix[local_i][local_j] += value * shapeValues[j];
elementMatrix[local_j][local_i] += value * shapeValues[j];
@@ -115,10 +141,6 @@ namespace AMDiS
}
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
- std::vector<Dune::FieldVector<double,1>> colShapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/ZeroOrderTestTrialvec.hpp b/src/amdis/assembler/ZeroOrderTestTrialvec.hpp
index fbe2f3575f1ff39d7e940c500ca12b5fc22098cb..352cacfc0b58455043df285da8222924caba2127 100644
--- a/src/amdis/assembler/ZeroOrderTestTrialvec.hpp
+++ b/src/amdis/assembler/ZeroOrderTestTrialvec.hpp
@@ -19,28 +19,26 @@ namespace AMDiS
/// zero-order operator \f$ \langle\psi, \mathbf{b}\cdot\Phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::test_trialvec, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
public:
- GridFunctionOperator(tag::test_trialvec, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 0)
+ GridFunctionOperator(tag::test_trialvec, GridFct const& expr)
+ : Super(expr, 0)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE>,
- std::integral_constant<bool, sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isLeaf && ColNode::isPower,
"RowNode must be Leaf-Node and ColNode must be a Power-Node.");
@@ -51,39 +49,45 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
- auto* rowShapeValues = &rowShapeValues_;
- auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_);
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
+ using RowRangeType = typename RowLocalBasisType::Traits::RangeType;
+ using ColRangeType = typename ColLocalBasisType::Traits::RangeType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ thread_local std::vector<RowRangeType> rowShapeValueBuffer;
+ thread_local std::vector<ColRangeType> colShapeValueBuffer;
+
+ const bool sameFE = std::is_same<decltype(rowLocalFE),decltype(colLocalFE)>::value;
+ auto* rowShapeValues = &rowShapeValueBuffer;
+ auto* colShapeValues = (sameFE ? &rowShapeValueBuffer : &colShapeValueBuffer);
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
const auto b = Super::coefficient(local);
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ rowLocalFE.localBasis().evaluateFunction(local, rowShapeValueBuffer);
if (!sameFE)
- colLocalFE.localBasis().evaluateFunction(local, colShapeValues_);
+ colLocalFE.localBasis().evaluateFunction(local, colShapeValueBuffer);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- const int local_i = rowNode.localIndex(i);
+ const auto local_i = rowNode.localIndex(i);
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- auto value = b * (factor * (*rowShapeValues)[i] * (*colShapeValues)[j]);
+ const auto value = b * (factor * (*rowShapeValues)[i] * (*colShapeValues)[j]);
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_kj = colNode.child(k).localIndex(j);
+ const auto local_kj = colNode.child(k).localIndex(j);
elementMatrix[local_i][local_kj] += value[k];
}
}
}
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
- std::vector<Dune::FieldVector<double,1>> colShapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/ZeroOrderTestvec.hpp b/src/amdis/assembler/ZeroOrderTestvec.hpp
index 49a5d40a90568993957cc2de0a16aa413178d96a..0a6c876c5d5f5a8fc6c26654e4e461fe50cb1ddd 100644
--- a/src/amdis/assembler/ZeroOrderTestvec.hpp
+++ b/src/amdis/assembler/ZeroOrderTestvec.hpp
@@ -19,25 +19,25 @@ namespace AMDiS
/// zero-order vector-operator \f$ (\mathbf{b}\cdot\Psi) \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::testvec, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::testvec, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
static_assert( Category::Vector<typename GridFct::Range>, "Expression must be of vector type." );
public:
- GridFunctionOperator(tag::testvec, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 0)
+ GridFunctionOperator(tag::testvec, GridFct const& expr)
+ : Super(expr, 0)
{}
- template <class Context, class QuadratureRule,
- class ElementVector, class Node>
- void calculateElementVector(Context const& context,
- QuadratureRule const& quad,
- ElementVector& elementVector,
- Node const& node)
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
{
static_assert(Node::isPower,
"Operator can be applied to Power-Nodes only.");
@@ -46,28 +46,30 @@ namespace AMDiS
auto const& localFE = node.child(0).finiteElement();
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeType = typename LocalBasisType::Traits::RangeType;
+
+ auto const& quad = this->getQuadratureRule(context.type(), node);
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
const auto exprValue = Super::coefficient(local);
- localFE.localBasis().evaluateFunction(local, shapeValues_);
+ thread_local std::vector<RangeType> shapeValues;
+ localFE.localBasis().evaluateFunction(local, shapeValues);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- auto value = exprValue * (factor * shapeValues_[i]);
+ const auto value = exprValue * (factor * shapeValues[i]);
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_ki = node.child(k).localIndex(i);
+ const auto local_ki = node.child(k).localIndex(i);
elementVector[local_ki] += value[k];
}
}
}
}
-
- private:
- std::vector<Dune::FieldVector<double,1>> shapeValues_;
};
/** @} **/
diff --git a/src/amdis/assembler/ZeroOrderTestvecTrial.hpp b/src/amdis/assembler/ZeroOrderTestvecTrial.hpp
index d6083d1ae911f0751411b9b10e7f7b8c4012dd46..f37c6aa46cc7bd46e7de3e1983c46b0c22c0b1b4 100644
--- a/src/amdis/assembler/ZeroOrderTestvecTrial.hpp
+++ b/src/amdis/assembler/ZeroOrderTestvecTrial.hpp
@@ -18,31 +18,19 @@ namespace AMDiS
/// zero-order operator \f$ \langle\Psi, \mathbf{b}\,\phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::testvec_trial, GridFct, QuadCreator>
- : public GridFunctionOperator<tag::test_trialvec, GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::testvec_trial, LocalContext, GridFct>
+ : public GridFunctionOperatorTransposed<GridFunctionOperator<tag::testvec_trial, LocalContext, GridFct>,
+ GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>>
{
- using Transposed = GridFunctionOperator<tag::test_trialvec, GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Transposed = GridFunctionOperator<tag::test_trialvec, LocalContext, GridFct>;
+ using Super = GridFunctionOperatorTransposed<Self, Transposed>;
public:
- GridFunctionOperator(tag::testvec_trial, GridFct const& expr, QuadCreator const& quadCreator)
- : Transposed(tag::test_trialvec{}, expr, quadCreator)
+ GridFunctionOperator(tag::testvec_trial, GridFct const& expr)
+ : Super(tag::test_trialvec{}, expr)
{}
-
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- std::integral_constant<bool, sameFE> flagSameFE,
- std::integral_constant<bool, sameNode> flagSameNode)
- {
- auto elementMatrixTransposed = trans(elementMatrix);
- Transposed::calculateElementMatrix(
- context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
- }
};
/** @} **/
diff --git a/src/amdis/assembler/ZeroOrderTestvecTrialvec.hpp b/src/amdis/assembler/ZeroOrderTestvecTrialvec.hpp
index 77a478e8685d749c9e69a451988fc83b585c0732..f8028777533a6387a55f8b1ef12cca1183f03b9b 100644
--- a/src/amdis/assembler/ZeroOrderTestvecTrialvec.hpp
+++ b/src/amdis/assembler/ZeroOrderTestvecTrialvec.hpp
@@ -19,30 +19,28 @@ namespace AMDiS
/// zero-order operator \f$ \langle\Psi, c\,\Phi\rangle \f$, or \f$ \langle\Psi, A\,\Phi\rangle \f$
- template <class GridFct, class QuadCreator>
- class GridFunctionOperator<tag::testvec_trialvec, GridFct, QuadCreator>
- : public GridFunctionOperatorBase<GridFct, QuadCreator>
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::testvec_trialvec, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::testvec_trialvec, LocalContext, GridFct>,
+ LocalContext, GridFct>
{
- using Super = GridFunctionOperatorBase<GridFct, QuadCreator>;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
using expr_value_type = typename GridFct::Range;
static_assert( Category::Scalar<expr_value_type> || Category::Matrix<expr_value_type>,
"Expression must be of scalar or matrix type." );
public:
- GridFunctionOperator(tag::testvec_trialvec, GridFct const& expr, QuadCreator const& quadCreator)
- : Super(expr, quadCreator, 0)
+ GridFunctionOperator(tag::testvec_trialvec, GridFct const& expr)
+ : Super(expr, 0)
{}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calculateElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- bool_t<sameFE>,
- bool_t<sameNode>)
+ template <class Context, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrix(Context const& context,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix)
{
static_assert(RowNode::isPower && ColNode::isPower,
"Operator can be applied to Power-Nodes only.");
@@ -50,54 +48,65 @@ namespace AMDiS
static_assert( (RowNode::CHILDREN == ColNode::CHILDREN), "" );
// theoretically the restriction of the same nr of childs would not be necessary
+ const bool sameFE = std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>::value;
+ const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
+
using Category = ValueCategory_t<expr_value_type>;
- calcElementMatrix(context, quad, elementMatrix, rowNode, colNode,
- bool_<sameFE>,
- bool_<sameNode>,
- Category{});
+
+ auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
+ if (sameFE && sameNode && std::is_same<Category,tag::scalar>::value)
+ getElementMatrixOptimized(context, quad, rowNode, colNode, elementMatrix, Category{});
+ else
+ getElementMatrixStandard(context, quad, rowNode, colNode, elementMatrix, Category{});
}
protected: // specialization for different value-categories
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE>
- void calcElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- bool_t<sameFE>,
- std::false_type /*sameNode*/,
- tag::scalar)
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixStandard(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix,
+ tag::scalar)
{
static const std::size_t CHILDREN = RowNode::CHILDREN;
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
- auto* rowShapeValues = &rowShapeValues_;
- auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_);
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
+ using RowRangeType = typename RowLocalBasisType::Traits::RangeType;
+ using ColRangeType = typename ColLocalBasisType::Traits::RangeType;
+
+ thread_local std::vector<RowRangeType> rowShapeValueBuffer;
+ thread_local std::vector<ColRangeType> colShapeValueBuffer;
+
+ const bool sameFE = std::is_same<decltype(rowLocalFE),decltype(colLocalFE)>::value;
+ auto* rowShapeValues = &rowShapeValueBuffer;
+ auto* colShapeValues = (sameFE ? &rowShapeValueBuffer : &colShapeValueBuffer);
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ rowLocalFE.localBasis().evaluateFunction(local, rowShapeValueBuffer);
if (!sameFE)
- colLocalFE.localBasis().evaluateFunction(local, colShapeValues_);
+ colLocalFE.localBasis().evaluateFunction(local, colShapeValueBuffer);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- double value = factor * (*rowShapeValues)[i];
+ const auto value = factor * (*rowShapeValues)[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- double value0 = value * (*colShapeValues)[j];
+ const auto value0 = value * (*colShapeValues)[j];
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_ki = rowNode.child(k).localIndex(i);
- const int local_kj = colNode.child(k).localIndex(j);
+ const auto local_ki = rowNode.child(k).localIndex(i);
+ const auto local_kj = colNode.child(k).localIndex(j);
elementMatrix[local_ki][local_kj] += value0;
}
@@ -107,46 +116,45 @@ namespace AMDiS
}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode>
- void calcElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& node,
- ColNode const& /*colNode*/,
- std::true_type /*sameFE*/,
- std::true_type /*sameNode*/,
- tag::scalar)
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixOptimized(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& node,
+ ColNode const& /*colNode*/,
+ ElementMatrix& elementMatrix,
+ tag::scalar)
{
static const std::size_t CHILDREN = RowNode::CHILDREN;
auto const& localFE = node.child(0).finiteElement();
- auto& shapeValues = rowShapeValues_;
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeType = typename LocalBasisType::Traits::RangeType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
+ thread_local std::vector<RangeType> shapeValues;
localFE.localBasis().evaluateFunction(local, shapeValues);
for (std::size_t i = 0; i < localFE.size(); ++i) {
- double value = factor * shapeValues[i];
+ const auto value = factor * shapeValues[i];
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_ki = node.child(k).localIndex(i);
+ const auto local_ki = node.child(k).localIndex(i);
elementMatrix[local_ki][local_ki] += value * shapeValues[i];
}
for (std::size_t j = i+1; j < localFE.size(); ++j) {
- double value0 = value * shapeValues[j];
+ const auto value0 = value * shapeValues[j];
for (std::size_t k = 0; k < CHILDREN; ++k) {
- const int local_ki = node.child(k).localIndex(i);
- const int local_kj = node.child(k).localIndex(j);
+ const auto local_ki = node.child(k).localIndex(i);
+ const auto local_kj = node.child(k).localIndex(j);
elementMatrix[local_ki][local_kj] += value0;
elementMatrix[local_kj][local_ki] += value0;
@@ -156,16 +164,13 @@ namespace AMDiS
}
}
- template <class Context, class QuadratureRule,
- class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode>
- void calcElementMatrix(Context const& context,
- QuadratureRule const& quad,
- ElementMatrix& elementMatrix,
- RowNode const& rowNode,
- ColNode const& colNode,
- bool_t<sameFE>,
- bool_t<sameNode>,
- tag::matrix)
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixStandard(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& rowNode,
+ ColNode const& colNode,
+ ElementMatrix& elementMatrix,
+ tag::matrix)
{
static const std::size_t CHILDREN = RowNode::CHILDREN;
static_assert( std::is_same<expr_value_type, Dune::FieldMatrix<double, CHILDREN, CHILDREN>>::value, "" );
@@ -173,32 +178,41 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
- auto* rowShapeValues = &rowShapeValues_;
- auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_);
+ using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
+ using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
+ using RowRangeType = typename RowLocalBasisType::Traits::RangeType;
+ using ColRangeType = typename ColLocalBasisType::Traits::RangeType;
- for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+ thread_local std::vector<RowRangeType> rowShapeValueBuffer;
+ thread_local std::vector<ColRangeType> colShapeValueBuffer;
+
+ const bool sameFE = std::is_same<decltype(rowLocalFE),decltype(colLocalFE)>::value;
+ auto* rowShapeValues = &rowShapeValueBuffer;
+ auto* colShapeValues = (sameFE ? &rowShapeValueBuffer : &colShapeValueBuffer);
+
+ for (auto const& qp : quad) {
// Position of the current quadrature point in the reference element
- decltype(auto) local = context.position(quad[iq].position());
+ decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula
- double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto factor = context.integrationElement(qp.position()) * qp.weight();
const Dune::FieldMatrix<double, CHILDREN, CHILDREN> exprValue = Super::coefficient(local);
- rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_);
+ rowLocalFE.localBasis().evaluateFunction(local, rowShapeValueBuffer);
if (!sameFE)
- colLocalFE.localBasis().evaluateFunction(local, colShapeValues_);
+ colLocalFE.localBasis().evaluateFunction(local, colShapeValueBuffer);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
- double value0 = factor * (*rowShapeValues)[i];
+ const auto value0 = factor * (*rowShapeValues)[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
- double value = value0 * (*colShapeValues)[j];
- auto mat = exprValue * value;
+ const auto value = value0 * (*colShapeValues)[j];
+ const auto mat = exprValue * value;
for (std::size_t k0 = 0; k0 < CHILDREN; ++k0) {
- const int local_ki = rowNode.child(k0).localIndex(i);
+ const auto local_ki = rowNode.child(k0).localIndex(i);
for (std::size_t k1 = 0; k1 < CHILDREN; ++k1) {
- const int local_kj = colNode.child(k1).localIndex(j);
+ const auto local_kj = colNode.child(k1).localIndex(j);
elementMatrix[local_ki][local_kj] += mat[k0][k1];
}
@@ -208,9 +222,16 @@ namespace AMDiS
}
}
- private:
- std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
- std::vector<Dune::FieldVector<double,1>> colShapeValues_;
+ template <class Context, class QuadratureRule, class RowNode, class ColNode, class ElementMatrix>
+ void getElementMatrixOptimized(Context const& context,
+ QuadratureRule const& quad,
+ RowNode const& node,
+ ColNode const& /*colNode*/,
+ ElementMatrix& elementMatrix,
+ tag::matrix)
+ {
+ assert(false && "Not implemented.");
+ }
};
/** @} **/
diff --git a/src/amdis/common/FieldMatVec.hpp b/src/amdis/common/FieldMatVec.hpp
index b8fac3891e22d708fb6fba10616bd42cf935897d..3786be2b455a22901b47b05685c7db605397c1dd 100644
--- a/src/amdis/common/FieldMatVec.hpp
+++ b/src/amdis/common/FieldMatVec.hpp
@@ -217,6 +217,12 @@ namespace AMDiS
template <class T, int N, int M>
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec);
+ template <class T, int N, int M>
+ FieldMatrix<T,N,M> operator*(FieldMatrix<T,N,M> mat, FieldVector<T,1> const& scalar);
+
+ template <class T, int N>
+ FieldMatrix<T,N,1> operator*(FieldMatrix<T,N,1> mat, FieldVector<T,1> const& scalar);
+
template <class T, int M, int N, int L>
diff --git a/src/amdis/common/FieldMatVec.inc.hpp b/src/amdis/common/FieldMatVec.inc.hpp
index 8353b618d2725aef034be1640605f2c59a985dff..83254c16d7663a36a2df2d4135f127650eae0ebf 100644
--- a/src/amdis/common/FieldMatVec.inc.hpp
+++ b/src/amdis/common/FieldMatVec.inc.hpp
@@ -404,6 +404,18 @@ FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const
return Dune::FMatrixHelp::mult(mat, vec);
}
+template <class T, int N, int M>
+FieldMatrix<T,N,M> operator*(FieldMatrix<T,N,M> mat, FieldVector<T,1> const& scalar)
+{
+ return mat *= scalar[0];
+}
+
+template <class T, int N>
+FieldMatrix<T,N,1> operator*(FieldMatrix<T,N,1> mat, FieldVector<T,1> const& scalar)
+{
+ return mat *= scalar[0];
+}
+
template <class T, int M, int N, int L>
diff --git a/src/amdis/common/TupleUtility.hpp b/src/amdis/common/TupleUtility.hpp
index 202622fd4704e49da37365d8f22434fa7110a2cb..0a6f1c03b475fef208d86ded0960f5c212d24fd0 100644
--- a/src/amdis/common/TupleUtility.hpp
+++ b/src/amdis/common/TupleUtility.hpp
@@ -31,7 +31,7 @@ namespace AMDiS
{
static_assert(std::tuple_size<Tuple>::value == sizeof...(args),
"Nr. of argument != tuple-size");
- return {std::forward<Args>(args)...};
+ return Tuple{std::forward<Args>(args)...};
}
};
@@ -43,7 +43,7 @@ namespace AMDiS
{
static_assert(std::tuple_size<Tuple>::value == 0,
"Construction of empty tuples with empty argument list only!");
- return {};
+ return Tuple{};
}
};
diff --git a/src/amdis/common/Utility.hpp b/src/amdis/common/Utility.hpp
index dd897b88aeb70d0833899abd705c9c5909490c7a..08f9683e07d22ad3440048b47f44bff570b8d876 100644
--- a/src/amdis/common/Utility.hpp
+++ b/src/amdis/common/Utility.hpp
@@ -4,6 +4,12 @@
#include <type_traits>
#include <vector>
+#if AMDIS_HAS_CXX_CONSTEXPR_IF
+#define IF_CONSTEXPR if constexpr
+#else
+#define IF_CONSTEXPR if
+#endif
+
namespace AMDiS
{
namespace Impl
@@ -46,6 +52,23 @@ namespace AMDiS
// ---------------------------------------------------------------------------
+ // base template for tuple size, must implement a static value member
+ template <class TupleType>
+ struct TupleSize;
+
+ // utility constant returning the value of the tuple size
+ template <class TupleType>
+ constexpr std::size_t TupleSize_v = TupleSize<TupleType>::value;
+
+
+ // base template for retrieving the tuple element type, must implement the type member `type`
+ template <size_t I, class TupleType>
+ struct TupleElement;
+
+ // base template for retrieving the tuple element type
+ template <size_t I, class TupleType>
+ using TupleElement_t = typename TupleElement<I, TupleType>::type;
+
/// A variadic type list
template <class... Ts>
@@ -55,6 +78,23 @@ namespace AMDiS
using Types_t = Types<std::decay_t<Ts>...>;
+ template <class... T>
+ struct TupleSize<Types<T...>>
+ : std::integral_constant<std::size_t, sizeof...(T)> {};
+
+ template <size_t I, class Head, class... Tail>
+ struct TupleElement<I, Types<Head, Tail...>>
+ {
+ using type = typename TupleElement<I-1, Types<Tail...>>::type;
+ };
+
+ template <class Head, class... Tail>
+ struct TupleElement<0, Types<Head, Tail...>>
+ {
+ using type = Head;
+ };
+
+
/// Alias that indicates ownership of resources
template <class T>
using owner = T;
diff --git a/src/amdis/gridfunctions/AnalyticGridFunction.hpp b/src/amdis/gridfunctions/AnalyticGridFunction.hpp
index 995a07c552c7f687524bd69e0d6f086edabc5536..c40e5af1b748b6f2821d84458814c8f3a43dc258 100644
--- a/src/amdis/gridfunctions/AnalyticGridFunction.hpp
+++ b/src/amdis/gridfunctions/AnalyticGridFunction.hpp
@@ -25,9 +25,11 @@ namespace AMDiS
using Geometry = typename LocalContext::Geometry;
+ enum { hasDerivative = true };
+
public:
AnalyticLocalFunction(Function const& fct)
- : fct_(fct)
+ : fct_{fct}
{}
void bind(LocalContext const& element)
@@ -42,6 +44,7 @@ namespace AMDiS
Range operator()(Domain const& local) const
{
+ assert( bool(geometry_) );
return fct_(geometry_.value().global(local));
}
@@ -99,6 +102,8 @@ namespace AMDiS
using Domain = typename EntitySet::GlobalCoordinate;
using Range = std::decay_t<std::result_of_t<Function(Domain)>>;
+ enum { hasDerivative = true };
+
private:
using Element = typename EntitySet::Element;
using LocalDomain = typename EntitySet::LocalCoordinate;
@@ -107,8 +112,8 @@ namespace AMDiS
public:
/// \brief Constructor. Stores the function `fct` and creates an `EntitySet`.
AnalyticGridFunction(Function const& fct, GridView const& gridView)
- : fct_(fct)
- , entitySet_(gridView)
+ : fct_{fct}
+ , entitySet_{gridView}
{}
/// \brief Return the evaluated functor at global coordinates
@@ -118,12 +123,11 @@ namespace AMDiS
}
/// \brief Return the LocalFunction of the AnalyticGridFunction.
- friend LocalFunction localFunction(AnalyticGridFunction const& gf)
+ LocalFunction localFunction() const
{
- return LocalFunction{gf.fct_};
+ return {fct_};
}
-
EntitySet const& entitySet() const
{
return entitySet_;
@@ -138,6 +142,12 @@ namespace AMDiS
};
+ template <class F, class GV>
+ auto localFunction(AnalyticGridFunction<F,GV> const& gf)
+ {
+ return gf.localFunction();
+ }
+
/// \brief Return a GridFunction representing the derivative of a functor.
// [expects: Functor f has free function derivative(f)]
template <class F, class GV>
diff --git a/src/amdis/gridfunctions/ConstantGridFunction.hpp b/src/amdis/gridfunctions/ConstantGridFunction.hpp
index d139d0f51a1890d94fb5f012d69bd4dc85fec736..275da96fbd7e114572c8ef65427a7a9e708cb23d 100644
--- a/src/amdis/gridfunctions/ConstantGridFunction.hpp
+++ b/src/amdis/gridfunctions/ConstantGridFunction.hpp
@@ -27,6 +27,8 @@ namespace AMDiS
using Geometry = typename LocalContext::Geometry;
+ enum { hasDerivative = true };
+
public:
ConstantLocalFunction(T const& value)
: value_(value)
@@ -41,8 +43,15 @@ namespace AMDiS
return value_;
}
- /// \relates ConstantLocalFunction
- friend int order(ConstantLocalFunction const& /*lf*/)
+ auto derivative() const
+ {
+ using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
+ using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
+ DerivativeRange diff(0);
+ return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
+ }
+
+ int order() const
{
return 0;
}
@@ -55,12 +64,17 @@ namespace AMDiS
template <class R, class D, class LocalContext, class T>
auto derivative(ConstantLocalFunction<R(D), LocalContext, T> const& lf)
{
- using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
- using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
- DerivativeRange diff(0);
- return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
+ return lf.derivative();
}
+ /// \relates ConstantLocalFunction
+ template <class R, class D, class LocalContext, class T>
+ int order(ConstantLocalFunction<R(D), LocalContext, T> const& lf)
+ {
+ return lf.order();
+ }
+
+
/**
* \addtogroup GridFunctions
@@ -85,6 +99,8 @@ namespace AMDiS
using LocalDomain = typename EntitySet::LocalCoordinate;
using Range = Underlying_t<T>;
+ enum { hasDerivative = false };
+
public:
using LocalFunction = ConstantLocalFunction<Range(LocalDomain), Element, T>;
@@ -106,7 +122,10 @@ namespace AMDiS
return entitySet_;
}
- T const& value() const { return value_; }
+ LocalFunction localFunction() const
+ {
+ return {value_};
+ }
private:
T value_;
@@ -115,9 +134,9 @@ namespace AMDiS
/// Return the LocalFunction of the \ref ConstantGridFunction. \relates ConstantGridFunction
template <class T, class GV>
- typename ConstantGridFunction<T,GV>::LocalFunction localFunction(ConstantGridFunction<T,GV> const& gf)
+ auto localFunction(ConstantGridFunction<T,GV> const& gf)
{
- return {gf.value()};
+ return gf.localFunction();
}
/** @} **/
diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp
index f0a5f9376131a55dc31ed1eae9f937674ab706b4..e71139e6cf220778070202bbcc276c99a885aec1 100644
--- a/src/amdis/gridfunctions/DOFVectorView.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.hpp
@@ -49,6 +49,7 @@ namespace AMDiS
template <class Block>
using Flat = Dune::Functions::FlatVectorBackend<Block>;
+ enum { hasDerivative = false };
public: // a local view on the gradients
@@ -59,6 +60,8 @@ namespace AMDiS
using Domain = LocalDomain;
using Range = DerivativeRange;
+ enum { hasDerivative = false };
+
private:
using LocalBasisView = typename GlobalBasis::LocalView;
using LocalIndexSet = typename GlobalBasis::LocalIndexSet;
@@ -135,6 +138,8 @@ namespace AMDiS
using Domain = typename DOFVectorConstView::LocalDomain;
using Range = typename DOFVectorConstView::Range;
+ enum { hasDerivative = true };
+
private:
using LocalBasisView = typename GlobalBasis::LocalView;
using LocalIndexSet = typename GlobalBasis::LocalIndexSet;
@@ -297,6 +302,13 @@ namespace AMDiS
return *this;
}
+ template <class Expr>
+ DOFVectorMutableView& operator<<(Expr&& expr)
+ {
+ return interpolate(expr);
+ }
+
+
/// Return the mutable DOFVector
DOFVector<Traits>& coefficients() { return *mutableDofVector_; }
diff --git a/src/amdis/gridfunctions/DerivativeGridFunction.hpp b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
index 5a7e621f2cc890ff5959c2259cfcd722329e432b..e4c0adc7f64304342260e00b5f0aa3d5f4bee75c 100644
--- a/src/amdis/gridfunctions/DerivativeGridFunction.hpp
+++ b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
@@ -32,6 +32,8 @@ namespace AMDiS
using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
using LocalFunction = std::decay_t<decltype(derivative(localFunction(std::declval<GridFunction>())))>;
+ enum { hasDerivative = false };
+
public:
/// The Range of the derivative of the GridFunction
using Range = typename DerivativeTraits::Range;
@@ -45,7 +47,7 @@ namespace AMDiS
public:
/// Constructor. Stores a copy of gridFct.
explicit DerivativeGridFunction(GridFunction const& gridFct)
- : gridFct_(gridFct)
+ : gridFct_{gridFct}
{
static_assert(isValidRange<DerivativeTraits>(), "Derivative of GridFunction not defined");
}
@@ -77,7 +79,7 @@ namespace AMDiS
#ifndef DOXYGEN
template <class GridFct,
class LocalFct = decltype(localFunction(std::declval<GridFct>())),
- REQUIRES(not Concepts::HasDerivative<GridFct>)>
+ REQUIRES(not GridFct::hasDerivative)>
auto derivative(GridFct const& gridFct)
{
static_assert(Concepts::HasDerivative<LocalFct>, "derivative(LocalFunction) not defined!");
diff --git a/src/amdis/gridfunctions/FunctorGridFunction.hpp b/src/amdis/gridfunctions/FunctorGridFunction.hpp
index c7905acd403c91df98836a5e006694d51f56b8e9..8b12ab2546c0e7580fc5067164c363560494bcf7 100644
--- a/src/amdis/gridfunctions/FunctorGridFunction.hpp
+++ b/src/amdis/gridfunctions/FunctorGridFunction.hpp
@@ -34,9 +34,11 @@ namespace AMDiS
} // end namespace Impl
+ // forward declaration
template <class Signatur, class Functor, class... LocalFunctions>
class FunctorLocalFunction;
+ // implementation
template <class R, class D, class Functor, class... LocalFunctions>
class FunctorLocalFunction<R(D), Functor, LocalFunctions...>
{
@@ -44,12 +46,28 @@ namespace AMDiS
using Range = R;
using Domain = D;
+ enum { hasDerivative = true };
+
+ template <class LocalFct>
+ struct LocalFunctionWrapper
+ {
+ template <class GridFct>
+ LocalFunctionWrapper(GridFct const& gridFct)
+ : localFct_{localFunction(gridFct)}
+ {}
+
+ LocalFct const& operator*() const { return localFct_; }
+ LocalFct & operator*() { return localFct_; }
+
+ LocalFct localFct_;
+ };
+
public:
/// Constructor. Stores copies of the functor and localFunction(gridfunction)s.
- template <class... LocalFcts>
- FunctorLocalFunction(Functor const& fct, LocalFcts&&... localFcts)
- : fct_(fct)
- , localFcts_(std::forward<LocalFcts>(localFcts)...)
+ template <class... GridFcts>
+ FunctorLocalFunction(Functor const& fct, GridFcts&&... gridFcts)
+ : fct_{fct}
+ , localFcts_{std::forward<GridFcts>(gridFcts)...}
{}
/// Calls \ref bind for all localFunctions
@@ -57,7 +75,7 @@ namespace AMDiS
void bind(Element const& element)
{
forEach(localFcts_, [&](auto& localFct) {
- localFct.bind(element);
+ (*localFct).bind(element);
});
}
@@ -65,22 +83,15 @@ namespace AMDiS
void unbind()
{
forEach(localFcts_, [&](auto& localFct) {
- localFct.unbind();
+ (*localFct).unbind();
});
}
/// Applies the functor \ref fct_ to the evaluated localFunctions
Range operator()(Domain const& x) const
{
- return eval([&x](auto const& localFct) { return localFct(x); },
- MakeSeq_t<sizeof...(LocalFunctions)>{});
- }
-
- private:
- template <class Op, std::size_t... I>
- auto eval(Op op, Indices<I...>) const
- {
- return fct_(op(std::get<I>(localFcts_))...);
+ using Dune::Std::apply;
+ return apply([&](auto&&... localFct) { return fct_((*localFct)(x)...); }, localFcts_);
}
public:
@@ -90,14 +101,14 @@ namespace AMDiS
return fct_;
}
- std::tuple<LocalFunctions...> const& localFcts() const
+ auto const& localFcts() const
{
return localFcts_;
}
private:
Functor fct_;
- std::tuple<LocalFunctions...> localFcts_;
+ std::tuple<LocalFunctionWrapper<LocalFunctions>...> localFcts_;
};
// Derivative of the LocalFunction of a FunctorGridFunction, utilizing
@@ -109,21 +120,24 @@ namespace AMDiS
REQUIRES(Concepts::HasPartial<F>)>
auto derivative(FunctorLocalFunction<Sig,F,LFs...> const& lf)
{
- auto index_seq = std::make_index_sequence<sizeof...(LFs)>{};
+ auto index_seq = std::index_sequence_for<LFs...>{};
// d_i(f)[lgfs...] * lgfs_i
auto term_i = [&](auto const _i)
{
- auto di_f = Dune::Std::apply([&](auto const&... lgfs) {
- return makeFunctorGridFunction(partial(lf.fct(), _i), lgfs...);
+ using Dune::Std::apply;
+ auto di_f = apply([&](auto const&... lgfs) {
+ return makeFunctorGridFunction(partial(lf.fct(), _i), (*lgfs)...);
}, lf.localFcts());
- auto const& lgfs_i = std::get<decltype(_i)::value>(lf.localFcts());
+ using std::get;
+ auto const& lgfs_i = *get<decltype(_i)::value>(lf.localFcts());
return makeFunctorGridFunction(Operation::Multiplies{}, di_f, derivative(lgfs_i));
};
- // sum_i [ d_i(f)[lgfs...] * derivative(lgfs_i) ]
- auto gridFct = Dune::Std::apply([&](auto const... _i)
+ // sum_i [ d_i(f)[lgfs...] * derivative(lgfs_i)
+ using Dune::Std::apply;
+ auto gridFct = apply([&](auto const... _i)
{
return makeFunctorGridFunction(Operation::Plus{}, term_i(_i)...);
}, index_seq);
@@ -139,7 +153,8 @@ namespace AMDiS
&& all_of_v<Concepts::HasOrder<LFs>...>)>
int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
{
- return Dune::Std::apply([&lf](auto const&... lgfs) { return order(lf.fct(), order(lgfs)...); },
+ using Dune::Std::apply;
+ return apply([&lf](auto const&... lgfs) { return order(lf.fct(), order(*lgfs)...); },
lf.localFcts());
}
@@ -173,6 +188,8 @@ namespace AMDiS
/// The set of entities this grid-function binds to
using EntitySet = typename Impl::EntitySetType<GridFunctions...>::type;
+ enum { hasDerivative = false };
+
private:
template <class GridFct>
using LocalFct = std::decay_t<decltype(localFunction(std::declval<GridFct>()))>;
@@ -186,31 +203,28 @@ namespace AMDiS
/// Constructor. Stores copies of the functor and gridfunctions.
template <class... GridFcts>
explicit FunctorGridFunction(Functor const& fct, GridFcts&&... gridFcts)
- : fct_(fct)
- , gridFcts_(std::forward<GridFcts>(gridFcts)...)
+ : fct_{fct}
+ , gridFcts_{std::forward<GridFcts>(gridFcts)...}
{}
/// Applies the functor to the evaluated gridfunctions
Range operator()(Domain const& x) const
{
- return eval(fct_, [&x](auto const& gridFct) { return gridFct(x); },
- MakeSeq_t<sizeof...(GridFunctions)>{});
+ using Dune::Std::apply;
+ return apply([&](auto&&... gridFct) { return fct_(gridFct(x)...); }, gridFcts_);
}
/// Return the stored \ref EntitySet of the first GridFunction
EntitySet const& entitySet() const
{
- return std::get<0>(gridFcts_).entitySet();
+ using std::get;
+ return get<0>(gridFcts_).entitySet();
}
- auto const& fct() const { return fct_; }
- auto const& gridFcts() const { return gridFcts_; }
-
- private:
- template <class Outer, class Inner, std::size_t... I>
- auto eval(Outer outer, Inner inner, Indices<I...>) const
+ LocalFunction localFunction() const
{
- return outer(inner(std::get<I>(gridFcts_))...);
+ using Dune::Std::apply;
+ return apply([&](auto const&... gridFcts) { return LocalFunction{fct_, gridFcts...}; }, gridFcts_);
}
private:
@@ -223,16 +237,10 @@ namespace AMDiS
template <class F, class... GFs>
auto localFunction(FunctorGridFunction<F,GFs...> const& gf)
{
- return Dune::Std::apply([&gf](auto const&... gridFcts)
- {
- using LocalFunction = typename FunctorGridFunction<F,GFs...>::LocalFunction;
- return LocalFunction{gf.fct(), localFunction(gridFcts)...};
- },
- gf.gridFcts());
+ return gf.localFunction();
}
-#ifndef DOXYGEN
// Generator function for FunctorGridFunction expressions
template <class Functor, class... GridFcts>
auto makeFunctorGridFunction(Functor const& f, GridFcts const&... gridFcts)
@@ -245,7 +253,7 @@ namespace AMDiS
}
// PreGridFunction related to FunctorGridFunction.
- template <class Functor, class... GridFunctions>
+ template <class Functor, class... PreGridFunctions>
struct FunctorPreGridFunction
{
using Self = FunctorPreGridFunction;
@@ -258,28 +266,27 @@ namespace AMDiS
return Dune::Std::apply([&](auto const&... pgf) {
return makeFunctorGridFunction(self.fct_,
makeGridFunction(pgf, gridView)...);
- }, self.gridFcts_);
+ }, self.preGridFcts_);
}
};
- template <class... GridFcts>
- explicit FunctorPreGridFunction(Functor const& fct, GridFcts&&... gridFcts)
- : fct_(fct)
- , gridFcts_(std::forward<GridFcts>(gridFcts)...)
+ template <class... PreGridFcts>
+ explicit FunctorPreGridFunction(Functor const& fct, PreGridFcts&&... pgfs)
+ : fct_{fct}
+ , preGridFcts_{std::forward<PreGridFcts>(pgfs)...}
{}
private:
Functor fct_;
- std::tuple<GridFunctions...> gridFcts_;
+ std::tuple<PreGridFunctions...> preGridFcts_;
};
namespace Traits
{
- template <class Functor, class... GridFunctions>
- struct IsPreGridFunction<FunctorPreGridFunction<Functor, GridFunctions...>>
+ template <class Functor, class... PreGridFcts>
+ struct IsPreGridFunction<FunctorPreGridFunction<Functor, PreGridFcts...>>
: std::true_type {};
}
-#endif
/// \brief Generator function for FunctorGridFunction. \relates FunctorGridFunction
@@ -291,10 +298,10 @@ namespace AMDiS
* - `invokeAtQP([](double u, auto const& x) { return u + x[0]; }, 1.0, X());`
* - `invokeAtQP(Operation::Plus{}, X(0), X(1));`
**/
- template <class Functor, class... GridFcts>
- auto invokeAtQP(Functor const& f, GridFcts&&... gridFcts)
+ template <class Functor, class... PreGridFcts>
+ auto invokeAtQP(Functor const& f, PreGridFcts&&... gridFcts)
{
- return FunctorPreGridFunction<Functor, std::decay_t<GridFcts>...>{f, std::forward<GridFcts>(gridFcts)...};
+ return FunctorPreGridFunction<Functor, std::decay_t<PreGridFcts>...>{f, std::forward<PreGridFcts>(gridFcts)...};
}
/** @} **/
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 73701a09ed01b65438f4cf3c7bab0ad55c500955..b5eb3aaba4cf03f3d7803773612b3456642963a3 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -4,7 +4,7 @@ dune_add_test(SOURCES AdaptInfoTest.cpp
dune_add_test(SOURCES ClonablePtrTest.cpp
LINK_LIBRARIES amdis)
-
+
dune_add_test(SOURCES ConceptsTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/TupleUtilityTest.cpp b/test/TupleUtilityTest.cpp
index 4f40c2c374d07aa36432f8b08f376e5037f310ac..9d99c3baf104c2e3e7e2025c1c505677dd2d6eeb 100644
--- a/test/TupleUtilityTest.cpp
+++ b/test/TupleUtilityTest.cpp
@@ -14,13 +14,13 @@ int main()
using TupleInt = std::tuple<int,int>;
using Tuple2 = std::tuple<TupleDouble,TupleInt>;
- Tuple1 u = {1.3, 2};
+ Tuple1 u{1.3, 2};
auto v = constructTuple<Tuple1>(1.5);
auto w = foldTuples<Tuple2>(u,v);
- AMDIS_TEST(u == Tuple1({1.3, 2}));
- AMDIS_TEST(v == Tuple1({1.5, 1}));
- AMDIS_TEST(w == Tuple2( {{1.3, 1.5}, {2, 1}} ));
+ AMDIS_TEST((u == Tuple1{1.3, 2}));
+ AMDIS_TEST((v == Tuple1{1.5, 1}));
+ AMDIS_TEST((w == Tuple2{TupleDouble{1.3, 1.5}, TupleDouble{2, 1}} ));
return report_errors();
}
\ No newline at end of file