Add ISTLMatrixBackend::addToEntry(row,col,value)

This method adds a value to an entry. In contrast to
`matrixBackend(row,col) += value`, the new
`matrixBackend.addToEntry(row,col,value)` can ignore
nonexisting entries to support e.g. `DiagonalMatrix`.

dune/fufem/backends/istlmatrixbackend.hh

@@ -63,6 +63,66 @@ namespace Impl {
     });
   }
 
+  template<class K, int n, class RowIndex, class ColIndex, class F>
+  decltype(auto) visitMatrixEntry(Dune::DiagonalMatrix<K,n>& matrix, const RowIndex& i, const ColIndex& j, F&& f)
+  {
+    if constexpr (std::is_same_v<std::invoke_result_t<F, K>, void>)
+    {
+      if (i==j)
+        f(matrix.diagonal(i));
+    }
+    else
+    {
+      assert(i==j);
+      return f(matrix.diagonal(i));
+    }
+  }
+
+  template<class K, int n, class RowIndex, class ColIndex, class F>
+  decltype(auto) visitMatrixEntry(const Dune::DiagonalMatrix<K,n>& matrix, const RowIndex& i, const ColIndex& j, F&& f)
+  {
+    if constexpr (std::is_same_v<std::invoke_result_t<F, K>, void>)
+    {
+      if (i==j)
+        f(matrix.diagonal(i));
+    }
+    else
+    {
+      assert(i==j);
+      return f(matrix.diagonal(i));
+    }
+  }
+
+  template<class K, int n, class RowIndex, class ColIndex, class F>
+  decltype(auto) visitMatrixEntry(Dune::ScaledIdentityMatrix<K,n>& matrix, const RowIndex& i, const ColIndex& j, F&& f)
+  {
+    if constexpr (std::is_same_v<std::invoke_result_t<F, K>, void>)
+    {
+      if ((i==0) and (j==0))
+        f(matrix.scalar());
+    }
+    else
+    {
+      assert((i==0) and (j==0));
+      return f(matrix.diagonal(i));
+    }
+  }
+
+  template<class K, int n, class RowIndex, class ColIndex, class F>
+  decltype(auto) visitMatrixEntry(const Dune::ScaledIdentityMatrix<K,n>& matrix, const RowIndex& i, const ColIndex& j, F&& f)
+  {
+    if constexpr (std::is_same_v<std::invoke_result_t<F, K>, void>)
+    {
+      if (i==j)
+        f(matrix.scalar());
+    }
+    else
+    {
+      assert(i==j);
+      return f(matrix.scalar());
+    }
+  }
+
 
 
   // Call f(matrix[i0][j0]...[in][jm]), by recursively resolving row- and column-multi-indices.
@@ -157,6 +217,25 @@ class ISTLMatrixBackend
     }
   };
 
+  template<class Result>
+  struct AddToScalar;
+
+  template<class R>
+  struct AddToScalar<R&>
+  {
+    const R& value_;
+    AddToScalar(const R& value) : value_(value) {}
+
+    template<class Matrix,
+      std::enable_if_t<std::is_convertible_v<Matrix&, R&>, int> = 0>
+    void operator()(Matrix& matrix) {
+      matrix += value_;
+    }
+
+    void operator()(Dune::FieldMatrix<R, 1, 1>& matrix) {
+      matrix[0][0] += value_;
+    }
+  };
 
 
 public:
@@ -186,6 +265,12 @@ public:
     return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<Entry&>());
   }
 
+  template<class RowMultiIndex, class ColMultiIndex>
+  void addToEntry(const RowMultiIndex& row, const ColMultiIndex& col, const Entry& value)
+  {
+    Impl::visitMatrixEntryRecursive(*matrix_, row, col, AddToScalar<Entry&>(value));
+  }
+
   /**
    * \brief Const access to wrapped matrix
    */