Implement local mass matrix assembler for arbitrary dune-functions bases

The local assembler can then be used in the
DuneFunctionsOperatorAssembler to create a global mass matrix assembler.

dune/fufem/assemblers/localassemblers/dunefunctionsmassassembler.hh

+#pragma once
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+
+#include <dune/typetree/visitor.hh>
+#include <dune/typetree/traversal.hh>
+
+namespace Dune::Fufem
+{
+
+namespace Impl
+{
+
+/** \brief Computes the inner product of the local basis functions
+ *
+ * It computes at each leaf node the inner product (hence, mass matrix) of the
+ * local basis function. It determines the local index and updates the
+ * corresponding entry of the local matrix.
+ */
+template <class LV, class M>
+class LeafNodeMassAssembler
+{
+public:
+
+  using LocalView = LV;
+  using Matrix = M;
+  using field_type = typename M::field_type;
+
+  /** \brief Constructor with all necessary context information
+   *
+   * \param [in] lv (localView) element information and localBasis
+   * \param [out] m (matrix) resulting local mass matrix, wrapped in the ISTLBackend
+
+   */
+  LeafNodeMassAssembler(const LV& lv, M& m)
+  : localView_(lv)
+  , matrix_(m)
+  {}
+
+  template<class Node, class TreePath>
+  void operator()(Node& node, TreePath treePath)
+  {
+    const auto& element = localView_.element();
+    const auto& geometry = element.geometry();
+    const auto& finiteElement = node.finiteElement();
+    const auto& localBasis = finiteElement.localBasis();
+    constexpr int gridDim = LocalView::GridView::dimension;
+    using ctype = typename LocalView::GridView::ctype;
+
+    std::vector<Dune::FieldVector<field_type,1>> values(localBasis.size());
+
+    QuadratureRuleKey quadKey(finiteElement);
+    const auto& quad = QuadratureRuleCache<double, gridDim>::rule(quadKey);
+
+    for( const auto& qp : quad )
+    {
+      const auto& quadPos = qp.position();
+      const auto integrationElement = geometry.integrationElement(quadPos);
+
+      localBasis.evaluateFunction(quadPos, values);
+
+      // compute matrix entries = inner products
+      auto z = qp.weight() * integrationElement;
+      for(int i=0; i<localBasis.size(); ++i)
+      {
+        double zi = values[i]*z;
+
+        auto rowIndex = node.localIndex(i);
+
+        // start off-diagonal, treat the diagonal extra
+        for (int j=i+1; j<localBasis.size(); ++j)
+        {
+          double zij = values[j] * zi;
+
+          auto colIndex = node.localIndex(j);
+
+          // mass matrix is symmetric
+          matrix_[rowIndex][colIndex] += zij;
+          matrix_[colIndex][rowIndex] += zij;
+        }
+
+        // z * values[i]^2
+        matrix_[rowIndex][rowIndex] += values[i] * zi;
+      }
+    }
+  }
+
+private:
+  const LocalView& localView_;
+  Matrix& matrix_;
+};
+
+} // namespace Impl
+
+/** \brief Local mass assembler for dune-functions basis
+ *
+ * We assume the mass matrix to be symmetric and hence ansatz and trial space to be equal!
+ * The implementation allows an arbitrary dune-functions compatible basis, as long as there
+ * is an ISTLBackend available for the resulting matrix.
+ */
+class DuneFunctionsLocalMassAssembler
+{
+public:
+
+  template<class Element, class LocalMatrix, class LocalView>
+  void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& trialLocalView, const LocalView& ansatzLocalView) const
+  {
+    // the mass matrix operator assumes trial == ansatz
+    if (&trialLocalView != &ansatzLocalView)
+      DUNE_THROW(Dune::Exception,"The mass matrix operator assumes equal ansatz and trial functions");
+
+    // matrix was already resized before
+    localMatrix = 0.;
+
+    // create a tree visitor and compute the inner products of the local functions at the leaf nodes
+    Impl::LeafNodeMassAssembler leadNodeMassAssembler(trialLocalView,localMatrix);
+    Dune::TypeTree::forEachLeafNode(trialLocalView.tree(),leadNodeMassAssembler);
+  }
+
+};
+
+} // namespace Dune::Fufem