LocalOperator and QuadratureFactory added, Tuple reimplemented with mapped construct of elements

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

src/amdis/ContextGeometry.hpp

 #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
+    Dune::Std::optional<LocalGeometry> localGeometry_;
   };
 
 } // end namespace AMDiS

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>
 
@@ -33,8 +35,9 @@ namespace AMDiS
    *   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 GridFunction>
   class GridFunctionOperatorBase
+      : public LocalOperator<Derived>
   {
     using LocalFunction = decltype(localFunction(std::declval<GridFunction>()));
     using Element = typename GridFunction::EntitySet::Element;
@@ -42,6 +45,8 @@ namespace AMDiS
 
     using LocalCoordinate = typename GridFunction::EntitySet::LocalCoordinate;
 
+    using QuadFactory = QuadratureFactory<typename Geometry::ctype, Element::dimension, LocalFunction>;
+
   public:
     /// \brief Constructor. Stores a copy of `expr`.
     /**
@@ -49,27 +54,24 @@ namespace AMDiS
      * \ref ExpressionBase, and stores a copy. Additionally, it gets the
      * differentiation order, 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`.
@@ -80,303 +82,158 @@ namespace AMDiS
      *
      * By default, it binds the \ref localFct_ to the `element`.
      **/
-    void bind(Element const& element, Geometry const& geometry)
+    void bind_impl(Element const& element, Geometry const& geometry)
     {
-      if (!bound_) {
+      assert( bool(quadFactory_) );
+      if (!this->bound_) {
         localFct_.bind(element);
-        isSimplex_ = geometry.type().isSimplex();
-        isAffine_ = geometry.affine();
-        bound_ = true;
+        quadFactory_->bind(localFct_);
       }
     }
 
     /// Unbinds operator from element.
-    void unbind()
+    void unbind_impl()
     {
-      bound_ = false;
       localFct_.unbind();
     }
 
+    template <class PreQuadFactory>
+    void setQuadFactory(PreQuadFactory const& pre)
+    {
+      quadFactory_ = makeQuadratureFactoryPtr<typename Geometry::ctype, Element::dimension, 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
     /// 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
+    int termOrder_ = 0;    //< the derivative order of this operator
   };
 
 
-  /// 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>
+  template <class Derived, class Transposed>
+  class GridFunctionOperatorTransposed
+      : public LocalOperator<Derived>
   {
+    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)
+    template <class Element, class Geometry>
+    void bind_impl(Element const& element, Geometry const& geometry)
     {
-      error_exit("Needs to be implemented!");
+      transposedOp_.bind(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)
+    void unbind_impl()
     {
-      error_exit("Needs to be implemented!");
+      transposedOp_.unbind();
     }
-  };
-
 
-#ifndef DOXYGEN
-  namespace Concepts
-  {
-    namespace Definition
+    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 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);
+    }
 
+  private:
+    Transposed transposedOp_;
+  };
 
-  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;
-  };
+  /// A default implementation of an GridFunctionOperator if no specialization is available.
+  /**
+   * An operator must implement either \ref getElementVector, or
+   * \ref getElementMatrix, if it is a vector or matrix operator,
+   * respectively.
+   **/
+  template <class Tag, class GridFct>
+  class GridFunctionOperator;
 
-  template <class Tag, class Expr>
-  struct ExpressionPreOperator<Tag, Expr, tag::order>
-  {
-    Tag tag;
-    Expr expr;
-    int order;
-    Dune::QuadratureType::Enum qt;
-  };
 
-  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)
+  template <class Tag, class PreGridFct, class... QuadratureArgs>
+  auto makeOperator(Tag tag, PreGridFct const& expr, QuadratureArgs&&... args)
   {
-    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