Assembler and LocalFiniteElement added, to separate assembling from problemStat

dune/amdis/Assembler.hpp

 #pragma once
 
-// std c++ headers
 #include <memory>
+#include <tuple>
 
-// AMDiS includes
-#include <dune/amdis/QuadratureRules.hpp>
-#include <dune/amdis/common/TypeDefs.hpp>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/amdis/LinearAlgebra.hpp>
+#include <dune/amdis/LocalFiniteElemenet.hpp>
+#include <dune/amdis/common/Mpl.hpp>
 
 namespace AMDiS
 {
-  // the LocalContext is either an Codim=0-EntityType or an IntersectionType
-  template <class GridView, class LocalContext>
+  template <class FeSpaces>
   class Assembler
   {
-    static constexpr int dim = GridView::dimension;
-    using ctype = typename GridView::ctype;
+    /// Number of problem components
+    static constexpr int nComponents = std::tuple_size<FeSpaces>::value;
+
+    using SystemMatrixType = SystemMatrix<FeSpaces>;
+    using SystemVectorType = SystemVector<FeSpaces>;
+
+    template <class T>
+    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;
 
-    using LocalQuadratureRules = Dune::QuadratureRules<ctype, LocalContext::mydimension>;
-    using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>;
-    using Geometry = typename Impl::Get<LocalContext>::Geometry;
+    template <class T>
+    using VectorEntries = Dune::FieldVector<T, nComponents>;
 
   public:
-    /// Constructor, initializes the geometry of the element and a
-    /// quadrature-rule wrapper
-    template <class Operator>
-    Assembler(Operator& op, LocalContext const& element, int degree, FirstOrderType type = GRD_PHI)
-      : geometry(get_geometry(element))
+    /// Constructor, stores a shared-pointer to the feSpaces
+    Assembler(std::shared_ptr<FeSpaces> const& feSpaces_)
+      : feSpaces(feSpaces_)
     {
-      int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type);
-      auto const& quad = LocalQuadratureRules::rule(geometry.type(), order);
+      auto gridView = std::get<0>(*feSpaces).gridView();
 
-      quadrature.reset(new QuadratureRule(element, quad));
+      forEach(range_<0, nComponents>, [&,this](auto const _r) {
+        static const std::size_t R = decltype(_r)::value;
+        std::get<R>(*feSpaces).update(gridView);
+      });
     }
 
-    /// \brief Returns reference to the transformed quadrature rule
-    /**
-     * For intersections this corresponds to the points in the intersection
-     * geometry, but coordinates extended to the whole inside-entity. For
-     * elements this is a wrapper around the classical element quadrature rule.
-     * Access the underlying dune-quadrature rule, with `quadrature->getRule()`.
-     **/
-    QuadratureRule const& getQuadrature() const
+    /// Assemble the linear system
+    template <class Operators>
+    void assemble(SystemMatrixType& matrix,
+                  SystemVectorType& solution,
+                  SystemVectorType& rhs,
+                  MatrixEntries<Operators>& matrix_operators,
+                  VectorEntries<Operators>& rhs_operators,
+                  bool asmMatrix_, bool asmVector_);
+
+  private:
+    // Return whether the matrix-block needs to be assembled
+    template <int R, int C, class Operators>
+    bool assembleMatrix(const index_t<R>, const index_t<C>,
+                        MatrixEntries<Operators> const& matrix_operators, bool asmMatrix_) const
     {
-      return *quadrature;
+      return asmMatrix_ && (!matrix_operators[R][C].assembled || matrix_operators[R][C].changing);
     }
 
-    /// Return the geometry of the Object
-    /**
-     * The geometry is either the element geometry, or the geometry of the
-     * inside-element for intersections.
-     **/
-    Geometry const& getGeometry() const
+    // Return whether the vector-block needs to be assembled
+    template <int R, class Operators>
+    bool assembleVector(const index_t<R>, VectorEntries<Operators> const& rhs_operators, bool asmVector_) const
     {
-      return geometry;
+      return asmVector_ && (!rhs_operators[R].assembled || rhs_operators[R].changing);
     }
 
-  protected:
+    // Sets the system to zero and initializes all operators and boundary conditions
+    template <class Operators>
+    void initMatrixVector(SystemMatrixType& matrix,
+                          SystemVectorType& solution,
+                          SystemVectorType& rhs,
+                          MatrixEntries<Operators>& matrix_operators,
+                          VectorEntries<Operators>& rhs_operators,
+                          bool asmMatrix_, bool asmVector_);
 
-    // transform coords from intersection coords to element coords
-    template <class Coordinate>
-    decltype(auto) map(Coordinate const& p) const
-    {
-      return get_position<LocalContext>(geometry, p);
-    }
 
-  private:
-    /// transformed quadrature rule
-    std::shared_ptr<QuadratureRule const> quadrature;
+    /// Assemble a block-matrix of element-matrices and return a matrix of flags, whether
+    /// a given block has received some entries.
+    template <class Operators>
+    MatrixEntries<bool> assembleElementMatrices(MatrixEntries<Operators>& operators,
+                                                MatrixEntries<ElementMatrix>& elementMatrix,
+                                                LocalFiniteElement<FeSpaces> const& localFiniteElem,
+                                                bool asmMatrix_);
+
+    /// Assemble a block-vector of element-vectors and return a vector of flags, whether
+    /// a given block has received some entries.
+    template <class Operators>
+    VectorEntries<bool> assembleElementVectors(VectorEntries<Operators>& operators,
+                                               VectorEntries<ElementVector>& elementVector,
+                                               LocalFiniteElement<FeSpaces> const& localFiniteElem,
+                                               bool asmVector_);
 
-    /// geometry() for entities or geometryInInside() for intersections
-    Geometry geometry;
 
-    // TODO: create geometry just once for each element, e.g. in ProblemStat when traversing the grid
+    /// Assemble one block of the block-element-matrix
+    // The MatrixData argument stores all matrix-operators
+    template <class Operators, class RowView, class ColView>
+    bool assembleElementMatrix(Operators& operators,
+                               ElementMatrix& elementMatrix,
+                               RowView const& rowLocalView, ColView const& colLocalView);
+
+    /// Assemble one block of the block-element-vector
+    // The VectorData argument stores all vector-operators
+    template <class Operators, class RowView>
+    bool assembleElementVector(Operators& operators,
+                               ElementVector& elementVector,
+                               RowView const& rowLocalView);
+
+
+    /// Add the block-element-matrix to the system-matrix
+    void addElementMatrices(SystemMatrixType& dofmatrix,
+                            LocalFiniteElement<FeSpaces> const& localFiniteElem,
+                            MatrixEntries<bool> const& addMat,
+                            MatrixEntries<ElementMatrix> const& elementMatrix);
+
+    /// Add the block-element-vector to the system-vector
+    void addElementVectors(SystemVectorType& dofvector,
+                           LocalFiniteElement<FeSpaces> const& localFiniteElem,
+                           VectorEntries<bool> const& addVec,
+                           VectorEntries<ElementVector> const& elementVector);
+
+    /// Finish insertion into the matrix and assembles boundary conditions
+    /// Return the number of nonzeros assembled into the matrix
+    template <class Operators>
+    std::size_t finishMatrixVector(SystemMatrixType& matrix,
+                                   SystemVectorType& solution,
+                                   SystemVectorType& rhs,
+                                   MatrixEntries<Operators>& matrix_operators,
+                                   VectorEntries<Operators>& rhs_operators,
+                                   bool asmMatrix_, bool asmVector_);
+
+  private:
+    std::shared_ptr<FeSpaces> feSpaces;
   };
 
 } // end namespace AMDiS
+
+#include "Assembler.inc.hpp"

dune/amdis/Assembler.inc.hpp

+#pragma once
+
+namespace AMDiS {
+
+template <class FeSpaces>
+  template <class Operators>
+void Assembler<FeSpaces>::
+assemble(SystemMatrixType& matrix,
+         SystemVectorType& solution,
+         SystemVectorType& rhs,
+         MatrixEntries<Operators>& matrix_operators,
+         VectorEntries<Operators>& rhs_operators,
+         bool asmMatrix_, bool asmVector_)
+{
+  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
+  initMatrixVector(matrix, solution, rhs, matrix_operators, rhs_operators, asmMatrix_, asmVector_);
+
+  // 2. create a localView and localIndexSet object for each global basis
+  LocalFiniteElement<FeSpaces> localFiniteElem(*feSpaces);
+
+  // 3. traverse grid and assemble operators on the elements
+  for (auto const& element : elements(gridView))
+  {
+    localFiniteElem.bind(element);
+
+    MatrixEntries<ElementMatrix> elementMatrix;
+    VectorEntries<ElementVector> elementVector;
+
+    MatrixEntries<bool> addMat = assembleElementMatrices(matrix_operators, elementMatrix, localFiniteElem, asmMatrix_);
+    VectorEntries<bool> addVec = assembleElementVectors(rhs_operators, elementVector, localFiniteElem, asmVector_);
+
+    addElementMatrices(matrix, localFiniteElem, addMat, elementMatrix);
+    addElementVectors(rhs, localFiniteElem, addVec, elementVector);
+  }
+
+  // 4. finish matrix insertion and apply dirichlet boundary conditions
+  std::size_t nnz = finishMatrixVector(asmMatrix_, asmVector_);
+
+  msg("fillin of assembled matrix: ", nnz);
+}
+
+
+template <class FeSpaces>
+  template <class Operators>
+void Assembler<FeSpaces>::
+initMatrixVector(SystemMatrixType& matrix,
+                 SystemVectorType& solution,
+                 SystemVectorType& rhs,
+                 MatrixEntries<Operators>& matrix_operators,
+                 VectorEntries<Operators>& rhs_operators,
+                 bool asmMatrix_, bool asmVector_)
+{
+  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  {
+    static const int R = decltype(_r)::value;
+    msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", R, "]");
+
+    if (assembleVector(_r, rhs_operators, asmVector_)) {
+      rhs.compress(_r);
+      rhs.getVector(_r) = 0.0;
+
+      // init vector operators
+      for (auto& op : rhs_operators[R].element)
+        op.init(this->getFeSpace(_r));
+      for (auto& op : rhs_operators[R].boundary)
+        op.init(this->getFeSpace(_r));
+      for (auto& op : rhs_operators[R].intersection)
+        op.init(this->getFeSpace(_r));
+    }
+
+    forEach(range_<0, nComponents>, [&,this](auto const _c)
+    {
+      static const int C = decltype(_c)::value;
+
+      bool asmMatrix = assembleMatrix(_r, _c, matrix_operators, asmMatrix_);
+      matrix(_r, _c).init(asmMatrix);
+
+      if (asmMatrix) {
+        // init matrix operators
+        for (auto& op : matrix_operators[R][C].element)
+          op.init(this->getFeSpace(_r), this->getFeSpace(_c));
+        for (auto& op : matrix_operators[R][C].boundary)
+          op.init(this->getFeSpace(_r), this->getFeSpace(_c));
+        for (auto& op : matrix_operators[R][C].intersection)
+          op.init(this->getFeSpace(_r), this->getFeSpace(_c));
+
+        // init boundary condition
+        for (int c = 0; c < nComponents; ++c)
+          for (auto bc : dirichletBc[R][c])
+            bc->init(c == C, matrix(_r, _c), solution[_c], rhs[_r]);
+      }
+    });
+  });
+}
+
+
+template <class FeSpaces>
+  template <class Operators>
+MatrixEntries<bool> Assembler<FeSpaces>::
+assembleElementMatrices(MatrixEntries<Operators>& operators,
+                        MatrixEntries<ElementMatrix>& elementMatrix,
+                        LocalFiniteElement<FeSpaces> const& localFiniteElem, bool asmMatrix_)
+{
+  MatrixEntries<bool> addMat;
+  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  {
+    static const std::size_t R = decltype(_r)::value;
+    forEach(range_<0, nComponents>, [&,this](auto const _c)
+    {
+      static const std::size_t C = decltype(_c)::value;
+
+      // assemble block of element matrix
+      addMat[R][C] = assembleMatrix(_r, _c, operators, asmMatrix_)
+        ? assembleElementMatrix(operators[R][C], elementMatrix[R][C],
+                                localFiniteElem.localView(_r), localFiniteElem.localView(_c))
+        : false;
+    });
+  });
+
+  return addMat;
+}
+
+
+template <class FeSpaces>
+  template <class Operators>
+VectorEntries<bool> Assembler<FeSpaces>::
+assembleElementVectors(VectorEntries<Operators>& operators,
+                       VectorEntries<ElementVector>& elementVector,
+                       LocalFiniteElement<FeSpaces> const& localFiniteElem, bool asmVector_)
+{
+  VectorEntries<bool> addVec;
+  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  {
+    static const std::size_t R = decltype(_r)::value;
+
+    // assemble block of element vector
+    addVec[R] = assembleVector(_r, operators, asmVector_)
+      ? assembleElementVector(operators[R], elementVector[R], localFiniteElem.localView(_r))
+      : false;
+  });
+
+  return addVec;
+}
+
+
+template <class FeSpaces>
+  template <class Operators, class RowView, class ColView>
+bool Assembler<FeSpaces>::
+assembleElementMatrix(Operators& operators,
+                      ElementMatrix& elementMatrix,
+                      RowView const& rowLocalView,
+                      ColView const& colLocalView)
+{
+  if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
+    return false; // nothing to do
+
+  auto const nRows = rowLocalView.tree().finiteElement().size();
+  auto const nCols = colLocalView.tree().finiteElement().size();
+
+  auto const& element = rowLocalView.element();
+  auto const& gridView = rowLocalView.globalBasis().gridView();
+
+  // fills the entire matrix with zeroes
+  elementMatrix.change_dim(nRows, nCols);
+  set_to_zero(elementMatrix);
+
+  bool add = false;
+
+  auto assemble_operators = [&](auto& e, auto& operator_list) {
+    for (auto scaled : operator_list) {
+      bool add_op = scaled.op->getElementMatrix(e, rowLocalView, colLocalView, elementMatrix, scaled.factor);
+      add = add || add_op;
+    }
+  };
+
+  // assemble element operators
+  assemble_operators(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())
+        assemble_operators(intersection, operators.boundary);
+      else
+        assemble_operators(intersection, operators.intersection);
+    }
+  }
+
+  return add;
+}
+
+
+template <class FeSpaces>
+  template <class Operators, class RowView>
+bool Assembler<FeSpaces>::
+assembleElementVector(Operators& operators, ElementVector& elementVector, RowView const& rowLocalView)
+{
+
+  if (operators.element.empty() && operators.boundary.empty() && operators.intersection.empty())
+    return false;
+
+  auto const nRows = rowLocalView.tree().finiteElement().size();
+
+  auto const& element = rowLocalView.element();
+  auto const& gridView = rowLocalView.globalBasis().gridView();
+
+  // Set all entries to zero
+  elementVector.change_dim(nRows);
+  set_to_zero(elementVector);
+
+  bool add = false;
+
+  auto assemble_operators = [&](auto& e, auto& operator_list) {
+    for (auto scaled : operator_list) {
+      bool add_op = scaled.op->getElementVector(e, rowLocalView, elementVector, scaled.factor);
+      add = add || add_op;
+    }
+  };
+
+  // assemble element operators
+  assemble_operators(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())
+        assemble_operators(intersection, operators.boundary);
+      else
+        assemble_operators(intersection, operators.intersection);
+    }
+  }
+
+  return add;
+}
+
+
+template <class FeSpaces>
+void Assembler<FeSpaces>::
+addElementMatrices(SystemMatrixType& dofmatrix,
+                   LocalFiniteElement<FeSpaces> const& localFiniteElem,
+                   MatrixEntries<bool> const& addMat,
+                   MatrixEntries<ElementMatrix> const& elementMatrix)
+{
+  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  {
+    static const std::size_t R = decltype(_r)::value;
+    forEach(range_<0, nComponents>, [&,this](auto const _c)
+    {
+      static const std::size_t C = decltype(_c)::value;
+      if (!addMat[R][C])
+        return;
+
+      // NOTE: current implementation does not utilize the multi-indices that we get from localIndexSet.
+
+      for (std::size_t i = 0; i < num_rows(elementMatrix[R][C]); ++i) {
+        // The global index of the i−th vertex of the element
+        auto const row = rowIndexSet.index(i);
+        for (std::size_t j = 0; j < num_cols(elementMatrix[R][C]); ++j) {
+          // The global index of the j−th vertex of the element
+          auto const col = colIndexSet.index(j);
+          dofmatrix[R][C](row,col) += elementMatrix[R][C](i,j);
+        }
+      }
+    });
+  });
+}
+
+
+template <class FeSpaces>
+void Assembler<FeSpaces>::
+addElementVectors(SystemVectorType& dofvector,
+                  LocalFiniteElement<FeSpaces> const& localFiniteElem,
+                  VectorEntries<bool> const& addVec,
+                  VectorEntries<ElementVector> const& elementVector)
+{
+  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  {
+    static const std::size_t R = decltype(_r)::value;
+    if (!addVec[R])
+      return;
+
+    for (std::size_t i = 0; i < size(elementVector[R]); ++i) {
+      // The global index of the i-th vertex
+      auto const row = indexSet.index(i);
+      dofvector[R][row] += elementVector[R][i];
+    }
+  });
+}
+
+
+template <class FeSpaces>
+  template <class Operators>
+std::size_t Assembler<FeSpaces>::
+finishMatrixVector(SystemMatrixType& matrix, SystemVectorType& solution, SystemVectorType& rhs,
+                   MatrixEntries<Operators>& matrix_operators, VectorEntries<Operators>& rhs_operators,
+                   bool asmMatrix_, bool asmVector_)
+{
+  std::size_t nnz = 0;
+  forEach(range_<0, nComponents>, [&,this](auto const _r)
+  {
+    static const int R = decltype(_r)::value;
+    if (assembleVector(asmVector_, _r))
+      rhs_operators[R].assembled = true;
+
+    forEach(range_<0, nComponents>, [&,this](auto const _c)
+    {
+      static const int C = decltype(_c)::value;
+      bool asmMatrix = assembleMatrix(asmMatrix_, _r, _c);
+
+      matrix(_r, _c).finish();
+
+      if (asmMatrix)
+        matrix_operators[R][C].assembled = true;
+
+      if (asmMatrix) {
+        // finish boundary condition
+        for (int c = 0; c < nComponents; ++c) {
+          for (int r = 0; r < nComponents; ++r) {
+            if (r != R && c != C)
+              continue;
+            for (auto bc : matrix_operators[r][c].dirichlet)
+              bc->finish(r == R, c == C, matrix(_r, _c), solution[_c], rhs[_r]);
+          }
+        }
+
+        nnz += matrix(_r, _c).getNnz();
+    });
+  });
+
+  return nnz;
+}
+
+} // end namespace AMDiS

dune/amdis/FirstOrderAssembler.hpp

@@ -2,7 +2,7 @@
 
 #include <vector>
 
-#include <dune/amdis/Assembler.hpp>
+#include <dune/amdis/LocalAssembler.hpp>
 #include <dune/amdis/common/Mpl.hpp>
 #include <dune/amdis/common/TypeDefs.hpp>
 #include <dune/amdis/utility/GetEntity.hpp>
@@ -11,9 +11,9 @@ namespace AMDiS
 {
   template <class GridView, class LocalContext, FirstOrderType type>
   class FirstOrderAssembler
-      : public Assembler<GridView, LocalContext>
+      : public LocalAssembler<GridView, LocalContext>
   {
-    using Super = Assembler<GridView, LocalContext>;
+    using Super = LocalAssembler<GridView, LocalContext>;
 
     static constexpr int dim = GridView::dimension;
     static constexpr int dow = GridView::dimensionworld;

dune/amdis/LocalAssembler.hpp

+#pragma once
+
+// std c++ headers
+#include <memory>
+
+// AMDiS includes
+#include <dune/amdis/QuadratureRules.hpp>
+#include <dune/amdis/common/TypeDefs.hpp>
+
+namespace AMDiS
+{
+  // the LocalContext is either an Codim=0-EntityType or an IntersectionType
+  template <class GridView, class LocalContext>
+  class LocalAssembler
+  {
+    static constexpr int dim = GridView::dimension;
+    using ctype = typename GridView::ctype;
+
+    using LocalQuadratureRules = Dune::QuadratureRules<ctype, LocalContext::mydimension>;
+    using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>;
+    using Geometry = typename Impl::Get<LocalContext>::Geometry;
+
+  public:
+    /// Constructor, initializes the geometry of the element and a
+    /// quadrature-rule wrapper
+    template <class Operator>
+    LocalAssembler(Operator& op, LocalContext const& element, int degree, FirstOrderType type = GRD_PHI)
+      : geometry(get_geometry(element))
+    {
+      int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type);
+      auto const& quad = LocalQuadratureRules::rule(geometry.type(), order);
+
+      quadrature.reset(new QuadratureRule(element, quad));
+    }
+
+    /// \brief Returns reference to the transformed quadrature rule
+    /**
+     * For intersections this corresponds to the points in the intersection
+     * geometry, but coordinates extended to the whole inside-entity. For
+     * elements this is a wrapper around the classical element quadrature rule.
+     * Access the underlying dune-quadrature rule, with `quadrature->getRule()`.
+     **/
+    QuadratureRule const& getQuadrature() const
+    {
+      return *quadrature;
+    }
+
+    /// Return the geometry of the Object
+    /**
+     * The geometry is either the element geometry, or the geometry of the
+     * inside-element for intersections.
+     **/
+    Geometry const& getGeometry() const
+    {
+      return geometry;
+    }
+
+  protected:
+