[cleanup]Simplify ISTLMatrixBackend entry access by multi-indices

This cleans up the entry access by a pair of row- and column-multi-indices as follows:

• Factor out a generic utility function visitMatrixEntry() for accessing the (i,j)-th entry of a matrix with dynamic indices. By combining hybridIndexAccess() and hybridRowIndexAccess() this also works with multi-type matrices.
• Factor out a generic utility function visitMatrixEntryRecursive() for recursively accessing matrix entries with dynamic row- and column-multi-indices. The action to apply to the resolved matrix entry is customizable.
• Use visitMatrixEntryRecursive() to return matrix entries.

This is a significant simplification of the old MultiIndexResolver-code and far more readable. It is also more generic, and forms the foundation of possible extensions to other methods that need multi-index-access, like, e.g., matrixEntryExists(i,j) or addToMatrixEntry(i,j).

dune/fufem/backends/istlmatrixbackend.hh

@@ -23,9 +23,9 @@ namespace Dune::Fufem {
 
 
 
-template<class M, class E=typename M::field_type>
-class ISTLMatrixBackend
-{
+namespace Impl {
+
+
 
   // The following is essentially a copy of Dune::Functions::hybridIndexAccess
   // But since matrices do not provide size(), we have to use N() instead.
@@ -49,92 +49,112 @@ class ISTLMatrixBackend
         });
   }
 
-  template<class Result, class RowIndex, class ColIndex>
-  struct MatrixMultiIndexResolver
+
+
+  // Call f(matrix[i][j]) and return the result.
+  // This works with dynamic indices i,j, even for a multi-type matrix.
+  // However, it requires that f() has a unique return type.
+  template<class Matrix, class RowIndex, class ColIndex, class F>
+  decltype(auto) visitMatrixEntry(Matrix&& matrix, const RowIndex& i, const ColIndex& j, F&& f)
   {
-    MatrixMultiIndexResolver(const RowIndex& rowIndex, const ColIndex& colIndex) :
-      rowIndex_(rowIndex),
-      colIndex_(colIndex)
-    {}
-
-    template<class C>
-    using isReturnable =
-    std::integral_constant<bool,
-      std::is_convertible<C&, Result>::value or
-      Dune::MatrixVector::Traits::ScalarTraits<std::decay_t<C>>::isScalar
-    >;
-
-    template<class C>
-    auto& toScalar(C&& c) {
-      return std::forward<C>(c);
-    }
+    using namespace Dune::Functions;
+    return hybridRowIndexAccess(matrix, i, [&](auto&& matrix_i) -> decltype(auto) {
+      return hybridIndexAccess(matrix_i, j, f);
+    });
+  }
 
-    // Catch the cases where a FieldMatrix<K, 1, 1> was supplied
-    auto& toScalar(FieldMatrix<std::decay_t<Result>, 1, 1>& c) {
-      return c[0][0];
-    }
 
-    auto& toScalar(const FieldMatrix<std::decay_t<Result>, 1, 1>& c) {
-      return c[0][0];
-    }
 
-    template<class C,
-      typename std::enable_if<not isReturnable<C>::value, int>::type = 0,
-      typename std::enable_if<IsSingleRowMatrix<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+  // Call f(matrix[i0][j0]...[in][jm]), by recursively resolving row- and column-multi-indices.
+  // Whenenver a SingleRowMatrix or SingleColumnMatrix is encountered, a zero row- or column-index
+  // is inserted, respectively. The recursion ends, whenever f(matrixEntry) can be called.
+  // This works with dynamic indices i,j, even for a multi-type matrix.
+  // However, it requires that f() has a unique return type.
+  template<class Matrix, class RowIndex, class ColIndex, class F>
+  static decltype(auto) visitMatrixEntryRecursive(Matrix& matrix, const RowIndex& iii, const ColIndex& jjj, F&& f)
+  {
+    using namespace Dune::Indices;
+    using namespace Dune::Functions::Imp;
+    static constexpr bool isSingleRow = IsSingleRowMatrix<Matrix>::value;
+    static constexpr bool isSingleCol = IsSingleColumnMatrix<Matrix>::value;
+    auto splitIndex = [] (auto&& multiIndex) { return std::make_tuple(multiIndex[_0], shiftedDynamicMultiIndex<1>(multiIndex)) ;};
+    if constexpr (std::is_invocable_v<F, Matrix&>)
     {
-      using namespace Dune::Indices;
-      using namespace Dune::Functions;
-      auto&& subRowIndex = rowIndex_;
-      auto&& subColIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(colIndex_);
-      auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return static_cast<Result>((hybridRowIndexAccess(c, _0, [&](auto&& ci) -> decltype(auto) {
-        return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
-      })));
+      return f(matrix);
     }
-
-    template<class C,
-      typename std::enable_if<not isReturnable<C>::value, int>::type = 0,
-      typename std::enable_if<IsSingleColumnMatrix<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+    else if constexpr (isSingleRow)
     {
-      using namespace Dune::Indices;
-      using namespace Dune::Functions;
-      auto&& subRowIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(rowIndex_);
-      auto&& subColIndex = colIndex_;
-      auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return static_cast<Result>((hybridRowIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
-        return static_cast<Result>((hybridIndexAccess(ci, _0, subIndexResolver)));
-      })));
+      assert(jjj.size()>0);
+      auto [j, jj] = splitIndex(jjj);
+      // We have to capture jj explicitly by value, because capturing structured bindings
+      // by reference is not allowed before C++20
+      return visitMatrixEntry(matrix, _0, j, [&, jj=jj](auto&& matrix_0j) -> decltype(auto) {
+        return visitMatrixEntryRecursive(matrix_0j, iii, jj, f);
+      });
     }
-
-    template<class C,
-      typename std::enable_if<not isReturnable<C>::value, int>::type = 0,
-      typename std::enable_if<not IsSingleRowMatrix<C>::value, int>::type = 0,
-      typename std::enable_if<not IsSingleColumnMatrix<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+    else if constexpr (isSingleCol)
     {
-      using namespace Dune::Indices;
-      using namespace Dune::Functions;
-      auto&& subRowIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(rowIndex_);
-      auto&& subColIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(colIndex_);
-      auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return static_cast<Result>((hybridRowIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
-        return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
-      })));
+      assert(iii.size()>0);
+      auto [i, ii] = splitIndex(iii);
+      // We have to capture ii explicitly by value, because capturing structured bindings
+      // by reference is not allowed before C++20
+      return visitMatrixEntry(matrix, i, _0, [&, ii=ii](auto&& matrix_i0) -> decltype(auto) {
+        return visitMatrixEntryRecursive(matrix_i0, ii, jjj, f);
+      });
     }
-
-    template<class C,
-      typename std::enable_if<isReturnable<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+    else
     {
-      return static_cast<Result>(toScalar(std::forward<C>(c)));
+      assert(iii.size()>0);
+      assert(jjj.size()>0);
+      auto [i, ii] = splitIndex(iii);
+      auto [j, jj] = splitIndex(jjj);
+      // We have to capture ii and jj explicitly by value, because capturing structured bindings
+      // by reference is not allowed before C++20
+      return visitMatrixEntry(matrix, i, j, [&, ii=ii, jj=jj](auto&& matrix_ij) -> decltype(auto) {
+        return visitMatrixEntryRecursive(matrix_ij, ii, jj, f);
+      });
+    }
+  }
+
+
+
+} // namespace Impl
+
+
+
+template<class M, class E=typename M::field_type>
+class ISTLMatrixBackend
+{
+
+  template<class Result>
+  struct ToScalar;
+
+  template<class R>
+  struct ToScalar<R&>
+  {
+    template<class Matrix,
+      std::enable_if_t<std::is_convertible_v<Matrix&, R&>, int> = 0>
+    R& operator()(Matrix& matrix) {
+      return matrix;
     }
 
+    R& operator()(Dune::FieldMatrix<R, 1, 1>& matrix) {
+      return matrix[0][0];
+    }
+  };
 
+  template<class R>
+  struct ToScalar<const R&>
+  {
+    template<class Matrix,
+      std::enable_if_t<std::is_convertible_v<Matrix&, const R&>, int> = 0>
+    const R& operator()(Matrix& matrix) {
+      return matrix;
+    }
 
-    const RowIndex& rowIndex_;
-    const ColIndex& colIndex_;
+    const R& operator()(const Dune::FieldMatrix<R, 1, 1>& matrix) {
+      return matrix[0][0];
+    }
   };
 
 
@@ -157,15 +177,13 @@ public:
   template<class RowMultiIndex, class ColMultiIndex>
   const Entry& operator()(const RowMultiIndex& row, const ColMultiIndex& col) const
   {
-    MatrixMultiIndexResolver<const Entry&, RowMultiIndex, ColMultiIndex> i(row, col);
-    return i(*matrix_);
+    return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<const Entry&>());
   }
 
   template<class RowMultiIndex, class ColMultiIndex>
   Entry& operator()(const RowMultiIndex& row, const ColMultiIndex& col)
   {
-    MatrixMultiIndexResolver<Entry&, RowMultiIndex, ColMultiIndex> i(row, col);
-    return i(*matrix_);
+    return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<Entry&>());
   }
 
   /**