move operators and assembling routines to DOFMatrixBase and DOFVectorBase

src/amdis/Assembler.hpp

 #pragma once
 
-#include <memory>
-#include <tuple>
-
-#include <amdis/DirichletBC.hpp>
-#include <amdis/LocalAssemblerList.hpp>
-#include <amdis/common/Mpl.hpp>
-
 namespace AMDiS
 {
-  template <class Traits>
-  class Assembler
+  template <class GridView, class Element, class Operators, class ElementAssembler>
+  void assembleOperators(
+      GridView const& gridView,
+      Element const& element,
+      Operators& operators,
+      ElementAssembler const& localAssembler)
   {
-    using GlobalBasis = typename Traits::GlobalBasis;
-
-    /// The grid view the global FE basis lives on
-    using GridView = typename GlobalBasis::GridView;
-
-  public:
-
-    /// Constructor, stores a shared-pointer to the feSpaces
-    Assembler(GlobalBasis& globalBasis,
-              MatrixOperators<GlobalBasis>& matrixOperators,
-              VectorOperators<GlobalBasis>& rhsOperators,
-              Constraints<GlobalBasis>& constraints)
-      : globalBasis_(globalBasis)
-      , matrixOperators_(matrixOperators)
-      , rhsOperators_(rhsOperators)
-      , constraints_(constraints)
-    {}
-
-    /// Assemble the linear system
-    template <class SystemMatrixType, class SystemVectorType>
-    void assemble(
-        SystemMatrixType& matrix,
-        SystemVectorType& solution,
-        SystemVectorType& rhs,
-        bool asmMatrix, bool asmVector);
-
-
-  private:
-
-    /// Sets the system to zero and initializes all operators and boundary conditions
-    template <class SystemMatrixType, class SystemVectorType>
-    void initMatrixVector(
-        SystemMatrixType& matrix,
-        SystemVectorType& solution,
-        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,
-        Operators& operators,
-        ElementAssembler const& elementAssembler) const;
-
-    /// Finish insertion into the matrix and assembles boundary conditions
-    /// Return the number of nonzeros assembled into the matrix
-    template <class SystemMatrixType, class SystemVectorType>
-    std::size_t finishMatrixVector(
-        SystemMatrixType& matrix,
-        SystemVectorType& solution,
-        SystemVectorType& rhs,
-        bool asmMatrix, bool asmVector) const;
-
-
-  private:
-
-    GlobalBasis& globalBasis_;
-    MatrixOperators<GlobalBasis>& matrixOperators_;
-    VectorOperators<GlobalBasis>& rhsOperators_;
-    Constraints<GlobalBasis>& constraints_;
-
-    //TODO: add caching of localBases
-  };
+    // assemble element operators
+    localAssembler(element, operators.element);
+
+    // assemble intersection operators
+    if (!operators.intersection.empty()
+        || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
+    {
+      for (auto const& intersection : intersections(gridView, element)) {
+        if (intersection.boundary())
+          localAssembler(intersection, operators.boundary);
+        else
+          localAssembler(intersection, operators.intersection);
+      }
+    }
+  }
 
 } // end namespace AMDiS
-
-#include "Assembler.inc.hpp"

src/amdis/Assembler.inc.hpp

-#pragma once
-
-#include <dune/common/dynmatrix.hh>
-#include <dune/common/dynvector.hh>
-#include <dune/functions/functionspacebases/subspacebasis.hh>
-#include <amdis/utility/TreePath.hpp>
-#include <amdis/utility/Visitor.hpp>
-
-#include <amdis/common/Math.hpp>
-
-namespace AMDiS {
-
-template <class Traits>
-  template <class SystemMatrixType, class SystemVectorType>
-void Assembler<Traits>::assemble(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector)
-{
-  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
-  initMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
-
-  auto localView = globalBasis_.localView();
-
-  // 2. create a local matrix and vector
-  std::size_t localSize = localView.maxSize();
-
-  Dune::DynamicMatrix<double> elementMatrix(localSize, localSize);
-  Dune::DynamicVector<double> elementVector(localSize);
-
-  // 3. traverse grid and assemble operators on the elements
-  for (auto const& element : elements(globalBasis_.gridView()))
-  {
-    elementMatrix = 0;
-    elementVector = 0;
-
-    localView.bind(element);
-    auto geometry = element.geometry();
-
-    // traverse type-tree of global-basis
-    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-    {
-      auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-      auto rowLocalView = rowBasis.localView();
-      rowLocalView.bind(element);
-
-      auto& rhsOp = rhsOperators_[rowNode];
-      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, rowLocalView.tree(), elementVector);
-        };
-
-        this->assembleElementOperators(element, rhsOp, vecAssembler);
-      }
-
-      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-      {
-        auto& matOp = matrixOperators_[rowNode][colNode];
-        if (matOp.doAssemble(asmMatrix) && !matOp.empty()) {
-          auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-          auto colLocalView = colBasis.localView();
-          colLocalView.bind(element);
-
-          matOp.bind(element, geometry);
-
-          auto matAssembler = [&](auto const& context, auto& operator_list) {
-            for (auto scaled : operator_list)
-              scaled.op->assemble(context, rowLocalView.tree(), colLocalView.tree(), elementMatrix);
-          };
-
-          this->assembleElementOperators(element, matOp, matAssembler);
-        }
-      });
-    });
-
-    // add element-matrix to system-matrix and element-vector to rhs
-    matrix.insert(localView, localView, elementMatrix);
-    rhs.insert(localView, elementVector);
-
-    // unbind all operators
-    forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
-      rhsOperators_[rowNode].unbind();
-      forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
-        matrixOperators_[rowNode][colNode].unbind();
-      });
-    });
-
-    localView.unbind();
-  }
-
-  // 4. finish matrix insertion and apply dirichlet boundary conditions
-  std::size_t nnz = finishMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
-
-  msg("fill-in of assembled matrix: {}", nnz);
-}
-
-
-template <class Traits>
-  template <class Element, class Operators, class ElementAssembler>
-void Assembler<Traits>::assembleElementOperators(
-    Element const& element,
-    Operators& operators,
-    ElementAssembler const& localAssembler) const
-{
-  // assemble element operators
-  localAssembler(element, operators.element);
-
-  // assemble intersection operators
-  if (!operators.intersection.empty()
-      || (!operators.boundary.empty() && element.hasBoundaryIntersections()))
-  {
-    for (auto const& intersection : intersections(globalBasis_.gridView(), element)) {
-      if (intersection.boundary())
-        localAssembler(intersection, operators.boundary);
-      else
-        localAssembler(intersection, operators.intersection);
-    }
-  }
-}
-
-
-template <class Traits>
-  template <class SystemMatrixType, class SystemVectorType>
-void Assembler<Traits>::initMatrixVector(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector) const
-{
-  matrix.init(asmMatrix);
-  solution.compress();
-  rhs.compress();
-  if (asmVector)
-    rhs = 0;
-
-  auto localView = globalBasis_.localView();
-  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-  {
-#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
-    if (rowNode.isLeaf)
-      msg("{} DOFs for Basis[{}]",globalBasis_.dimension(rowTreePath), to_string(rowTreePath));
-#endif
-
-    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-
-    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-    {
-      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-
-      for (auto bc : constraints_[rowNode][colNode])
-        bc->init(matrix, solution, rhs, rowBasis, colBasis);
-    });
-  });
-  msg("{} total DOFs", globalBasis_.dimension());
-}
-
-
-template <class Traits>
-  template <class SystemMatrixType, class SystemVectorType>
-std::size_t Assembler<Traits>::finishMatrixVector(
-    SystemMatrixType& matrix,
-    SystemVectorType& solution,
-    SystemVectorType& rhs,
-    bool asmMatrix, bool asmVector) const
-{
-  matrix.finish();
-
-  auto localView = globalBasis_.localView();
-  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
-  {
-    auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
-    auto& rhsOp = rhsOperators_[rowNode];
-    if (rhsOp.doAssemble(asmVector))
-      rhsOp.assembled = true;
-
-    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
-    {
-      auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
-      auto& matOp = matrixOperators_[rowNode][colNode];
-      if (matOp.doAssemble(asmMatrix))
-        matOp.assembled = true;
-
-      // finish boundary condition
-      for (auto bc : constraints_[rowNode][colNode])
-        bc->finish(matrix, solution, rhs, rowBasis, colBasis);
-    });
-  });
-
-  return matrix.nnz();
-}
-
-} // end namespace AMDiS

src/amdis/Boundary.hpp

+#pragma once
+
+namespace AMDiS
+{
+  struct BoundaryType { int b; };
+
+} // end namespace AMDiS

src/amdis/CMakeLists.txt

@@ -20,7 +20,7 @@ install(FILES
     AdaptStationary.hpp
     AMDiS.hpp
     Assembler.hpp
-    Assembler.inc.hpp
+    Boundary.hpp
     ContextGeometry.hpp
     CreatorInterface.hpp
     CreatorMap.hpp

src/amdis/DirichletBC.hpp

@@ -8,6 +8,7 @@
 #include <dune/common/hybridutilities.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
 
+#include <amdis/Boundary.hpp>
 #include <amdis/Output.hpp>
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/ValueCategory.hpp>
@@ -16,8 +17,6 @@
 
 namespace AMDiS
 {
-  struct BoundaryType { int b; };
-
   /// Implements a boundary condition of Dirichlet-type.
   /**
    * By calling the methods \ref init() and \ref finish before and after
@@ -119,17 +118,17 @@ namespace AMDiS
   };
 
 
-  template <class GlobalBasis>
+  template <class Basis>
   struct DirichletData
   {
-    using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
+    using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
 
     template <class RowNode, class ColNode>
     using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RowNode>>> >;
   };
 
-  template <class GlobalBasis>
-  using Constraints = MatrixData<GlobalBasis, DirichletData<GlobalBasis>::template type>;
+  template <class RowBasis, class ColBasis>
+  using Constraints = MatrixData<RowBasis, ColBasis, DirichletData<RowBasis>::template type>;
 
 } // end namespace AMDiS
 

src/amdis/LocalAssemblerList.hpp

@@ -3,11 +3,19 @@
 #include <list>
 #include <memory>
 
+#include <amdis/Boundary.hpp>
 #include <amdis/LocalAssemblerBase.hpp>
 #include <amdis/utility/TreeData.hpp>
 
 namespace AMDiS
 {
+  namespace tag
+  {
+    template <class E> struct element_operator { using type = E; };
+    template <class I> struct intersection_operator { using type = I; };
+    template <class I> struct boundary_operator { using type = I; BoundaryType id; };
+  }
+
   template <class GridView>
   class OperatorLists
   {
@@ -18,8 +26,6 @@ namespace AMDiS
     struct DataElement
     {
       std::shared_ptr<OperatorType> op;
-      double* factor = nullptr;
-      double* estFactor = nullptr;
       BoundaryType b = {0};
     };
 
@@ -39,9 +45,9 @@ namespace AMDiS
       }
 
       /// Test whether to assemble on the node
-      bool doAssemble(bool flag) const
+      bool doAssemble() const
       {
-        return flag && (!assembled || changing);
+        return !assembled || changing;
       }
 
       /// Bind all operators to the grid element and geometry
@@ -61,6 +67,26 @@ namespace AMDiS
         for (auto& scaled : intersection) scaled.op->unbind();
       }
 
+
+      template <class... Args>
+      void push(tag::element_operator<Element>, Args&&... args)
+      {
+        element.push_back({std::forward<Args>(args)...});
+      }
+
+      template <class... Args>
+      void push(tag::boundary_operator<Intersection> b, Args&&... args)
+      {
+        boundary.push_back({std::forward<Args>(args)..., b.id});
+      }
+
+      template <class... Args>
+      void push(tag::intersection_operator<Intersection>, Args&&... args)
+      {
+        intersection.push_back({std::forward<Args>(args)...});
+      }
+
+
       /// List of operators to be assembled on grid elements
       std::list<DataElement<ElementOperator>> element;
       /// List of operators to be assembled on boundary intersections
@@ -87,8 +113,8 @@ namespace AMDiS
   };
 
 
-  template <class GlobalBasis>
-  using MatrixOperators = MatrixData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template MatData>;
+  template <class RowBasis, class ColBasis>
+  using MatrixOperators = MatrixData<RowBasis, ColBasis, OperatorLists<typename RowBasis::GridView>::template MatData>;
 
   template <class GlobalBasis>
   using VectorOperators = VectorData<GlobalBasis, OperatorLists<typename GlobalBasis::GridView>::template VecData>;

src/amdis/ProblemStat.hpp

@@ -13,7 +13,6 @@
 #include <dune/grid/common/grid.hh>
 
 #include <amdis/AdaptInfo.hpp>
-#include <amdis/Assembler.hpp>
 #include <amdis/CreatorInterface.hpp>
 #include <amdis/CreatorMap.hpp>
 #include <amdis/DirichletBC.hpp>
@@ -116,22 +115,34 @@ namespace AMDiS
     /// Adds an operator to \ref A.
     /** @{ */
     template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
-    void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
+    void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {})
+    {
+      systemMatrix_->addOperator(tag::element_operator<Element>{}, op, row, col);
+    }
 
     // operator evaluated on the boundary
     template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
-    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
+    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {})
+    {
+      systemMatrix_->addOperator(tag::boundary_operator<typename GridView::Intersection>{b}, op, row, col);
+    }
     /** @} */
 
 
     /// Adds an operator to \ref rhs.
     /** @{ */
     template <class Operator, class TreePath = RootTreePath>
-    void addVectorOperator(Operator const& op, TreePath = {});
+    void addVectorOperator(Operator const& op, TreePath path = {})
+    {
+      rhs_->addOperator(tag::element_operator<Element>{}, op, path);
+    }
 
     // operator evaluated on the boundary
     template <class Operator, class TreePath = RootTreePath>
-    void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
+    void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {})
+    {
+      rhs_->addOperator(tag::boundary_operator<typename GridView::Intersection>{b}, op, path);
+    }
     /** @} */
 
 
@@ -272,9 +283,7 @@ namespace AMDiS
     void initGlobalBasis(GlobalBasis const& globalBasis)
     {
       localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis_->localView());
-      matrixOperators_.init(localView_->tree(), tag::store{});
-      rhsOperators_.init(localView_->tree(), tag::store{});
-      constraints_.init(localView_->tree(), tag::store{});
+      constraints_.init(localView_->tree(), localView_->tree(), tag::store{});
     }
 
     void createMatricesAndVectors()
@@ -377,9 +386,7 @@ namespace AMDiS
 
   private: // some internal data-structures
 
-    MatrixOperators<GlobalBasis> matrixOperators_;
-    VectorOperators<GlobalBasis> rhsOperators_;
-    Constraints<GlobalBasis> constraints_;
+    Constraints<GlobalBasis, GlobalBasis> constraints_;
   };
 
 #if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION

src/amdis/ProblemStat.inc.hpp

@@ -8,7 +8,6 @@
 #include <dune/typetree/childextraction.hh>
 
 #include <amdis/AdaptInfo.hpp>
-#include <amdis/Assembler.hpp>
 #include <amdis/FileWriter.hpp>
 #include <amdis/LocalAssembler.hpp>
 #include <amdis/GridFunctionOperator.hpp>
@@ -158,94 +157,6 @@ void ProblemStat<Traits>::createFileWriter()
 }
 
 
-// add matrix/vector operator terms
-
-template <class Traits>
-  template <class Operator, class RowTreePath, class ColTreePath>
-void ProblemStat<Traits>::addMatrixOperator(
-    Operator const& preOp,
-    RowTreePath row, ColTreePath col)
-{
-  static_assert( Concepts::PreTreePath<RowTreePath>,
-      "row must be a valid treepath, or an integer/index-constant");
-  static_assert( Concepts::PreTreePath<ColTreePath>,
-      "col must be a valid treepath, or an integer/index-constant");
-
-  auto i = child(localView_->tree(), makeTreePath(row));
-  auto j = child(localView_->tree(), makeTreePath(col));
-
-  auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
-  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);