Merge branch 'feature/operator_localoperator_restructuring' into 'master'

Implement Operator and LocalOperator with TypeErasure

See merge request !204

amdis/BiLinearForm.hpp

@@ -145,11 +145,12 @@ namespace AMDiS
     }
 
     /// Assemble the matrix operators on the bound element.
-    template <class RowLocalView, class ColLocalView,
+    template <class RowLocalView, class ColLocalView, class LocalOperators,
       REQUIRES(Concepts::LocalView<RowLocalView>),
       REQUIRES(Concepts::LocalView<ColLocalView>)>
     void assemble(RowLocalView const& rowLocalView,
-                  ColLocalView const& colLocalView);
+                  ColLocalView const& colLocalView,
+                  LocalOperators& localOperators);
 
     /// Assemble all matrix operators, TODO: incorporate boundary conditions
     void assemble();
@@ -179,8 +180,19 @@ namespace AMDiS
 
     /// Set the flag that forces an update of the pattern since the underlying
     /// basis that defines the indexset has been changed
-    void updateImpl(event::adapt e, index_t<0> i) override { updatePattern_ = true; }
-    void updateImpl(event::adapt e, index_t<1> i) override { updatePattern_ = true; }
+    void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(*rowBasis_); }
+    void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(*colBasis_); }
+
+    auto& operators() { return operators_; }
+
+  private:
+    template <class Basis>
+    void updateImpl2(Basis const& basis)
+    {
+      if (!updatePattern_)
+        Recursive::forEach(operators_, [&](auto& op) { op.update(basis.gridView()); });
+      updatePattern_ = true;
+    }
 
   protected:
     /// The symmetry property if the bilinear form

amdis/BiLinearForm.inc.hpp

@@ -2,10 +2,10 @@
 
 #include <utility>
 
-#include <amdis/Assembler.hpp>
+#include <amdis/ContextGeometry.hpp>
+#include <amdis/GridFunctionOperator.hpp>
 #include <amdis/LocalOperator.hpp>
 #include <amdis/typetree/Traversal.hpp>
-#include <amdis/utility/AssembleOperators.hpp>
 
 namespace AMDiS {
 
@@ -21,42 +21,35 @@ addOperator(ContextTag contextTag, Expr const& expr,
       "col must be a valid treepath, or an integer/index-constant");
 
   auto i = makeTreePath(row);
-  auto node_i = child(this->rowBasis().localView().treeCache(), i);
   auto j = makeTreePath(col);
-  auto node_j = child(this->colBasis().localView().treeCache(), j);
 
   using LocalContext = typename ContextTag::type;
-  using Tr = DefaultAssemblerTraits<LocalContext, ElementMatrix>;
-  auto op = makeLocalOperator<LocalContext>(expr, this->rowBasis().gridView());
-  auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), node_i, node_j));
-
-  operators_[i][j].push(contextTag, std::move(localAssembler));
+  auto op = makeOperator<LocalContext>(expr, this->rowBasis().gridView());
+  operators_[i][j].push(contextTag, std::move(op));
   updatePattern_ = true;
 }
 
 
 template <class RB, class CB, class T, class Traits>
-  template <class RowLocalView, class ColLocalView,
+  template <class RowLocalView, class ColLocalView, class LocalOperators,
     REQUIRES_(Concepts::LocalView<RowLocalView>),
     REQUIRES_(Concepts::LocalView<ColLocalView>)>
 void BiLinearForm<RB,CB,T,Traits>::
-assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
+assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView,
+         LocalOperators& localOperators)
 {
   elementMatrix_.resize(rowLocalView.size(), colLocalView.size());
   elementMatrix_ = 0;
 
   auto const& gv = this->rowBasis().gridView();
-  auto const& element = rowLocalView.element();
-  auto geometry = element.geometry();
+  GlobalContext context{gv, rowLocalView.element(), rowLocalView.element().geometry()};
 
   Traversal::forEachNode(rowLocalView.treeCache(), [&](auto const& rowNode, auto rowTp) {
     Traversal::forEachNode(colLocalView.treeCache(), [&](auto const& colNode, auto colTp) {
-      auto& matOp = operators_[rowTp][colTp];
-      if (matOp) {
-        matOp.bind(element, geometry);
-        assembleOperators(gv, element, matOp, makeMatrixAssembler(rowNode, colNode, elementMatrix_));
-        matOp.unbind();
-      }
+      auto& matOp = localOperators[rowTp][colTp];
+      matOp.bind(context.element());
+      matOp.assemble(context, rowNode, colNode, elementMatrix_);
+      matOp.unbind();
     });
   });
 
@@ -72,13 +65,14 @@ assemble()
   auto colLocalView = this->colBasis().localView();
 
   this->init();
+  auto localOperators = AMDiS::localOperators(operators_);
   for (auto const& element : elements(this->rowBasis().gridView(), typename Traits::PartitionSet{})) {
     rowLocalView.bind(element);
     if (this->rowBasis_ == this->colBasis_)
-      this->assemble(rowLocalView, rowLocalView);
+      this->assemble(rowLocalView, rowLocalView, localOperators);
     else {
       colLocalView.bind(element);
-      this->assemble(rowLocalView, colLocalView);
+      this->assemble(rowLocalView, colLocalView, localOperators);
       colLocalView.unbind();
     }
     rowLocalView.unbind();

amdis/CMakeLists.txt

@@ -21,8 +21,6 @@ install(FILES
     AdaptiveGrid.hpp
     AdaptStationary.hpp
     AMDiS.hpp
-    Assembler.hpp
-    AssemblerInterface.hpp
     BackupRestore.hpp
     BiLinearForm.hpp
     BiLinearForm.inc.hpp
@@ -56,6 +54,7 @@ install(FILES
     MeshCreator.hpp
     Observer.hpp
     Operations.hpp
+    Operator.hpp
     OperatorList.hpp
     Output.hpp
     PeriodicBC.hpp

amdis/ContextGeometry.hpp

@@ -160,4 +160,51 @@ namespace AMDiS
     mutable std::optional<LocalGeometry> localGeometry_;
   };
 
+
+  template <class GV>
+  class GlobalContext
+  {
+  public:
+    using GridView = GV;
+    using Element = typename GV::template Codim<0>::Entity;
+    using Geometry = typename Element::Geometry;
+
+    enum {
+      dim = GridView::dimension,      //< the dimension of the grid element
+      dow = GridView::dimensionworld  //< the dimension of the world
+    };
+
+    /// Constructor. Stores a copy of gridView and a pointer to element and geometry.
+    GlobalContext(GridView const& gridView, Element const& element,
+                  Geometry const& geometry)
+      : gridView_(gridView)
+      , element_(&element)
+      , geometry_(&geometry)
+    {}
+
+  public:
+    /// Return the GridView this context is bound to
+    GridView const& gridView() const
+    {
+      return gridView_;
+    }
+
+    /// Return the bound element (entity of codim 0)
+    Element const& element() const
+    {
+      return *element_;
+    }
+
+    /// Return the geometry of the \ref Element
+    Geometry const& geometry() const
+    {
+      return *geometry_;
+    }
+
+  private:
+    GridView gridView_;
+    Element const* element_;
+    Geometry const* geometry_;
+  };
+
 } // end namespace AMDiS

amdis/GridFunctionOperator.hpp

@@ -3,11 +3,9 @@
 #include <cassert>
 #include <type_traits>
 
-#include <amdis/GridFunctions.hpp>
-#include <amdis/LocalOperator.hpp>
-#include <amdis/common/Transposed.hpp>
+#include <amdis/GridFunctionOperatorTransposed.hpp>
+#include <amdis/common/Order.hpp>
 #include <amdis/common/TypeTraits.hpp>
-#include <amdis/typetree/FiniteElementType.hpp>
 #include <amdis/utility/QuadratureFactory.hpp>
 
 namespace AMDiS
@@ -17,256 +15,246 @@ namespace AMDiS
    * @{
    **/
 
-  /// \brief The main implementation of an CRTP-base class for operators using a grid-function
-  /// coefficient to be used in an \ref Assembler.
+  template <class LF, class Imp>
+  class GridFunctionLocalOperator;
+
+
+  /// \brief The main implementation of an operator depending on a grid-function
   /**
-   * An Operator that takes a GridFunction as coefficient.
-   * Provides quadrature rules and handles the evaluation of the GridFunction at
-   * local coordinates.
+   * An Operator that takes a grid-function as coefficient.
+   * Generates a \ref GridFunctionLocalOperator on \ref localOperator()
    *
-   * The class is specialized, by deriving from it, in \ref GridFunctionOperator.
+   * The class implements the interface of an \ref Operator.
    *
-   * \tparam Derived  The class derived from GridFunctionOperatorBase
-   * \tparam LC       The Element or Intersection type
-   * \tparam GF       The GridFunction, a LocalFunction is created from, and
-   *                  that is evaluated at quadrature points.
+   * \tparam GF   The class type of the grid-function
+   * \tparam Imp  Class providing the local assembling method, forwarded to
+   *              GridFunctionLocalOperator class
    *
    * **Requirements:**
-   * - `LC` models either Entity (of codim 0) or Intersection
    * - `GF` models the \ref Concepts::GridFunction
    **/
-  template <class Derived, class LC, class GF>
-  class GridFunctionOperatorBase
-      : public LocalOperator<Derived, LC>
+  template <class GF, class Imp>
+  class GridFunctionOperator
   {
-    template <class, class>
-    friend class LocalOperator;
-
-    using ContextType = Impl::ContextTypes<LC>;
-    using Super = LocalOperator<Derived, LC>;
-
-  private:
+  public:
     using GridFunction = GF;
+    using Implementation = Imp;
 
-    /// The type of the localFunction associated with the GridFunction
-    using LocalFunction = decltype(localFunction(std::declval<GF>()));
-
-    /// The Codim=0 entity of the grid, the localFunction can be bound to
-    using Element = typename ContextType::Entity;
-
-    /// The geometry-type of the grid element
-    using Geometry = typename Element::Geometry;
-
-    /// The type of the local coordinates in the \ref Element
-    using LocalCoordinate = typename GF::EntitySet::LocalCoordinate;
-
-    /// A factory for QuadratureRules that incooperate the order of the LocalFunction
-    using QuadFactory = QuadratureFactory<typename Geometry::ctype, LC::mydimension, LocalFunction>;
-
-  public:
-    /// \brief Constructor. Stores a copy of `gridFct`.
+    /// \brief Constructor. Stores a copy of `gridFct` and `impl`.
     /**
-     * A GridFunctionOperator takes a gridFunction and the
-     * differentiation order of the operator, to calculate the
-     * quadrature degree in \ref getDegree.
+     * A GridFunctionOperator takes a grid-function and
+     * an implementation class for the assemble method.
      **/
-    template <class GridFct>
-    GridFunctionOperatorBase(GridFct&& gridFct, int termOrder)
+    template <class GridFct, class Impl>
+    GridFunctionOperator(GridFct&& gridFct, Impl&& impl,
+                         int derivDeg, int gridFctOrder)
       : gridFct_(FWD(gridFct))
-      , termOrder_(termOrder)
+      , impl_(FWD(impl))
+      , derivDeg_(derivDeg)
+      , gridFctOrder_(gridFctOrder)
     {}
 
-    /// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
-    /// \tparam  PQF  A PreQuadratureFactory
-    template <class PQF>
-    void setQuadFactory(PQF&& pre)
-    {
-      using ctype = typename Geometry::ctype;
-      quadFactory_ = makeUniquePtr(
-        makeQuadratureFactory<ctype, LC::mydimension, LocalFunction>(FWD(pre)));
-    }
-
-  protected:
-    /// Return expression value at LocalCoordinates
-    auto coefficient(LocalCoordinate const& local) const
-    {
-      assert( this->bound_ );
-      return (*localFct_)(local);
-    }
+    template <class GridView>
+    void update(GridView const&) { /* do nothing */ }
 
-    /// Create a quadrature rule using the \ref QuadratureFactory by calculating the
-    /// quadrature order necessary to integrate the (bi)linear-form accurately.
-    template <class... Nodes>
-    auto const& getQuadratureRule(Dune::GeometryType type, Nodes const&... nodes) const
+    friend auto localOperator(GridFunctionOperator const& op)
     {
-      assert( bool(quadFactory_) );
-      int quadDegree = this->getDegree(termOrder_, quadFactory_->order(), nodes...);
-      return quadFactory_->rule(type, quadDegree);
+      return GridFunctionLocalOperator{localFunction(op.gridFct_), op.impl_,
+        op.derivDeg_, op.gridFctOrder_};
     }
 
   private:
-    /// \brief Binds operator to `element` and `geometry`.
-    /**
-     * Binding an operator to the currently visited element in grid traversal.
-     * Since all operators need geometry information, the `Element::Geometry`
-     * object `geometry` is created once, during grid traversal, and passed in.
-     *
-     * By default, it binds the \ref localFct_ to the `element` and the Quadrature
-     * factory to the localFunction.
-     **/
-    void bind_impl(Element const& element, Geometry const& geometry)
-    {
-      assert( bool(quadFactory_) );
-      localFct_.emplace(localFunction(gridFct_));
-      localFct_->bind(element);
-      quadFactory_->bind(localFct_.value());
-    }
-
-    /// Unbinds operator from element.
-    void unbind_impl()
-    {
-      localFct_->unbind();
-      localFct_.reset();
-    }
-
-  private:
-    /// The gridFunction to be used within the operator
+    /// The grid-function associated to this operator
     GridFunction gridFct_;
 
-    /// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
-    std::optional<LocalFunction> localFct_;
+    /// An implementation class for the assembling
+    Implementation impl_;
 
-    /// Assign each element type a quadrature rule
-    std::shared_ptr<QuadFactory> quadFactory_;
+    /// Maximal degree of derivative this operator represents
+    int derivDeg_;
 
-    /// the derivative order of this operator (in {0, 1, 2})
-    int termOrder_ = 0;
+    /// Polynomial degree of the grid-function (or -1)
+    int gridFctOrder_;
   };
 
+  template <class GridFct, class Impl>
+  GridFunctionOperator(GridFct const& gridFct, Impl const& impl, int, int)
+    -> GridFunctionOperator<GridFct, Impl>;
+
 
-  /// \brief The transposed operator, implemented in term of its transposed by
-  /// calling \ref getElementMatrix with inverted arguments.
-  template <class Derived, class Transposed>
-  class GridFunctionOperatorTransposed
-      : public LocalOperator<Derived, typename Transposed::LocalContext>
+  /// \brief The main implementation of a local-operator depending on a local-function
+  /**
+   * A LocalOperator that takes a local-function as coefficient.
+   * Provides quadrature rules and passes the local-function, bound to an element,
+   * to the assemble method of an implementation class.
+   *
+   * The class implements the interface of a \ref LocalOperator.
+   *
+   * \tparam LF   The class type of the local-function
+   * \tparam Imp  Class providing the local assembling method
+   *
+   * **Requirements:**
+   * - `LF` models the \ref Concepts::LocalFunction
+   **/
+  template <class LF, class Imp>
+  class GridFunctionLocalOperator
   {
-    template <class, class>
-    friend class LocalOperator;
+  private:
+    /// The type of the localFunction
+    using LocalFunction = LF;
 
-    template <class T, class... Args>
-    using Constructable = decltype( new T(std::declval<Args>()...) );
+    /// Type of the implementation class
+    using Implementation = Imp;
 
   public:
-    template <class... Args,
-      std::enable_if_t<Dune::Std::is_detected<Constructable, Transposed, Args...>::value, int> = 0>
-    GridFunctionOperatorTransposed(Args&&... args)
-      : transposedOp_(FWD(args)...)
+    /// \brief Constructor. Stores a copy of `localFct` and `impl`.
+    /**
+     * A GridFunctionLocalOperator takes a local-function, an implementation class,
+     * the differentiation order of the operator and the local-function polynomial
+     * degree, to calculate the quadrature degree of the operator
+     **/
+    template <class LocalFct, class Impl>
+    GridFunctionLocalOperator(LocalFct&& localFct, Impl&& impl,
+                              int derivDeg, int localFctOrder)
+      : localFct_(FWD(localFct))
+      , impl_(FWD(impl))
+      , derivDeg_(derivDeg)
+      , localFctOrder_(localFctOrder)
     {}
 
-    /// Redirects the setQuadFactory call top the transposed operator
-    /// \tparam  PQF  A PreQuadratureFactory
-    template <class PQF>
-    void setQuadFactory(PQF&& pre)
+    /// \brief Binds operator to `element`.
+    /**
+     * By default, it binds the \ref localFct_ to the `element`.
+     **/
+    template <class Element>
+    void bind(Element const& element)
     {
-      transposedOp_.setQuadFactory(FWD(pre));
+      localFct_.bind(element);
     }
 
-  private:
-    /// Redirects the bind call top the transposed operator
-    template <class Element, class Geometry>
-    void bind_impl(Element const& element, Geometry const& geometry)
+    /// Unbinds operator from element.
+    void unbind()
     {
-      transposedOp_.bind(element, geometry);
+      localFct_.unbind();
     }
 
-    /// Redirects the unbind call top the transposed operator
-    void unbind_impl()
+    /// Assemble a local element matrix on the element that is bound.
+    /**
+     * This function calls the assemble method from the impl_ class and
+     * additionally passes a quadrature rule and  the localFct_ to that method.
+     **/
+    template <class CG, class RN, class CN, class Mat>
+    void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
+                  Mat& elementMatrix) const
     {
-      transposedOp_.unbind();
+      auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(),
+        derivDeg_, localFctOrder(), rowNode, colNode);
+      impl().assemble(contextGeo, rowNode, colNode, quad, localFct_, elementMatrix);
     }
 
-    /// Apply the assembling to the transposed elementMatrix with interchanged row-/colNode
+    /// Assemble a local element vector on the element that is bound.
     /**
-     * \tparam CG  ContextGeometry
-     * \tparam RN  RowNode
-     * \tparam CN  ColNode
-     * \tparam Mat ElementMatrix
+     * This function calls the assemble method from the impl_ class and
+     * additionally passes a quadrature rule and the localFct_ to that method.
      **/
-    template <class CG, class RN, class CN, class Mat>
-    void getElementMatrix(CG const& contextGeometry, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
+    template <class CG, class Node, class Vec>
+    void assemble(CG const& contextGeo, Node const& node,
+                  Vec& elementVector) const
     {
-      auto elementMatrixTransposed = transposed(elementMatrix);
-      transposedOp_.getElementMatrix(
-        contextGeometry, colNode, rowNode, elementMatrixTransposed);
+      auto const& quad = getQuadratureRule(contextGeo.localContext().geometry(),
+        derivDeg_, localFctOrder(), node);
+      impl().assemble(contextGeo, node, quad, localFct_, elementVector);
+    }
+
+    Implementation      & impl()       { return impl_; }
+    Implementation const& impl() const { return impl_; }
+
+  protected:
+    // calculated polynomial order of local-function. Fallback to localFctOrder_;
+    int localFctOrder() const
+    {
+      if constexpr (Concepts::Polynomial<LF>)
+        return order(localFct_);
+      else
+        return localFctOrder_;
     }
 
   private:
-    Transposed transposedOp_;
+    /// The local-function to be used within the operator
+    LocalFunction localFct_;
+
+    /// Implementation details of the assembling
+    Implementation impl_;
+
+    /// Maximal degree of derivative this operator represents
+    int derivDeg_;
+
+    /// Polynomial degree of the local-function (or -1)
+    int localFctOrder_;
   };
 
+  // deduction guide
+  template <class LocalFct, class Impl>
+  GridFunctionLocalOperator(LocalFct const& localFct, Impl const& impl, int, int)
+    -> GridFunctionLocalOperator<LocalFct, Impl>;
 
-#ifndef DOXYGEN
-  template <class Tag, class PreGridFct, class PQF>
-  struct PreGridFunctionOperator
+
+  /// Registry to specify a tag for each implementation type
+  template <class Tag, class LocalContext>
+  struct GridFunctionOperatorRegistry {};
+
+  template <class Tag, class Expr>
+  struct OperatorTerm
   {
     Tag tag;
-    PreGridFct expr;
-    PQF quadFactory;
+    Expr expr;
+    int gridFctDeg;
+
+    OperatorTerm(Tag tag, Expr const& expr, int gridFctDeg = -1)
+      : tag(tag)
+      , expr(expr)
+      , gridFctDeg(gridFctDeg)
+    {}
   };
-#endif
 
-  /// Store tag and expression into a \ref PreGridFunctionOperator to create a \ref GridFunctionOperator
-  template <class Tag, class Expr, class... QuadratureArgs>
-  auto makeOperator(Tag tag, Expr&& expr, QuadratureArgs&&... args)
+  /// Store tag and expression into a \ref OperatorTerm to create
+  /// a \ref GridFunctionOperator
+  template <class Tag, class Expr>
+  auto makeOperator(Tag const& tag, Expr&& expr, int gridFctDeg = -1)
   {
-    auto pqf = makePreQuadratureFactory(FWD(args)...);
-    return PreGridFunctionOperator<Tag, TYPEOF(expr), TYPEOF(pqf)>{tag, FWD(expr), std::move(pqf)};
+    return OperatorTerm{tag, FWD(expr), gridFctDeg};
+  }
+
+  template <class Tag, class Expr>
+  auto operatorTerm(Tag const& tag, Expr&& expr, int gridFctDeg = -1)
+  {
+    return OperatorTerm{tag, FWD(expr), gridFctDeg};
   }
 
   /** @} **/
 
 #ifndef DOXYGEN
 
-  /// The base-template for GridFunctionOperators
-  /**
-   * An operator can specialize this class, by deriving from \ref GridFunctionOperatorBase.
-   * With the generic function \ref makeLocalOperator, an instance is created. To
-   * distinguisch different GridFunction operators, a tag can be provided that has no
-   * other effect.
-   *
-   * \tparam Tag  An Identifier for this GridFunctionOperator
-   * \tparam LC  An Element or Intersection the operator is evaluated on
-   * \tparam GridFct  A GridFunction evaluated in local coordinates on the bound element
-   **/
-  template <class Tag, class LC, class GridFct>
-  class GridFunctionOperator
-      : public GridFunctionOperatorBase<GridFunctionOperator<Tag, LC, GridFct>, LC, GridFct>
-  {};
-
-  template <class Context, class Tag, class GF, class QF>
-  auto makeGridFunctionOperator(Tag tag, GF&& gf, QF&& qf)
-  {
-    GridFunctionOperator<Tag, Context, TYPEOF(gf)> gfo{tag, FWD(gf)};
-    gfo.setQuadFactory(FWD(qf));
-    return gfo;
-  }
+  template <class R, class Tag>
+  using IsTransposed = decltype(R::transposedTag(std::declval<Tag>()));
 
-
-  /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr).
+  /// Generate an \ref GridFunctionOperator from a OperatorTerm (tag, expr).
   /// @{
-  template <class Context, class... Args, class GridView>
-  auto makeLocalOperator(PreGridFunctionOperator<Args...> op, GridView const& gridView)
+  template <class Context, class Tag, class Expr, class GridView>
+  auto makeOperator(OperatorTerm<Tag,Expr> const& op, GridView const& gridView)
   {