Fix ISTLMatrixBackend

Now, also FieldMatrix<K,1,1> will be handled correctly. Moreover,
the old C-style casts were replaced by static_casts.

dune/fufem/assemblers/istlbackend.hh

@@ -14,6 +14,7 @@
 #include <dune/functions/common/indexaccess.hh>
 #include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
 
+#include <dune/matrix-vector/traits/scalartraits.hh>
 
 
 namespace Dune {
@@ -66,7 +67,6 @@ constexpr std::true_type isSingleRowMatrix(const SingleRowMatrix<T>*)
 } // namespace Imp
 
 
-
 /**
  * \brief Traits class to check if matrix is marked as single column matrix
  */
@@ -147,8 +147,29 @@ class ISTLMatrixBackend
       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);
+    }
+
+    // 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 std::is_convertible<C&, Result>::value, int>::type = 0,
+      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)
     {
@@ -157,13 +178,13 @@ class ISTLMatrixBackend
       auto&& subRowIndex = rowIndex_;
       auto&& subColIndex = Dune::Functions::shiftedMultiIndex(colIndex_);
       auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return (Result)(hybridIndexAccess(c, _0, [&](auto&& ci) -> decltype(auto) {
-        return (Result)(hybridIndexAccess(ci, colIndex_[_0], subIndexResolver));
-      }));
+      return static_cast<Result>((hybridIndexAccess(c, _0, [&](auto&& ci) -> decltype(auto) {
+        return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
+      })));
     }
 
     template<class C,
-      typename std::enable_if<not std::is_convertible<C&, Result>::value, int>::type = 0,
+      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)
     {
@@ -172,13 +193,13 @@ class ISTLMatrixBackend
       auto&& subRowIndex = Dune::Functions::shiftedMultiIndex(rowIndex_);
       auto&& subColIndex = colIndex_;
       auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return (Result)(hybridIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
-        return (Result)(hybridIndexAccess(ci, _0, subIndexResolver));
-      }));
+      return static_cast<Result>((hybridIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
+        return static_cast<Result>((hybridIndexAccess(ci, _0, subIndexResolver)));
+      })));
     }
 
     template<class C,
-      typename std::enable_if<not std::is_convertible<C&, Result>::value, int>::type = 0,
+      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)
@@ -188,16 +209,16 @@ class ISTLMatrixBackend
       auto&& subRowIndex = Dune::Functions::shiftedMultiIndex(rowIndex_);
       auto&& subColIndex = Dune::Functions::shiftedMultiIndex(colIndex_);
       auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return (Result)(hybridIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
-        return (Result)(hybridIndexAccess(ci, colIndex_[_0], subIndexResolver));
-      }));
+      return static_cast<Result>((hybridIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
+        return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
+      })));
     }
 
     template<class C,
-      typename std::enable_if<std::is_convertible<C&, Result>::value, int>::type = 0>
+      typename std::enable_if<isReturnable<C>::value, int>::type = 0>
     Result operator()(C&& c)
     {
-      return (Result)(std::forward<C>(c));
+      return static_cast<Result>(toScalar(std::forward<C>(c)));
     }
 
     const RowIndex& rowIndex_;

dune/fufem/test/CMakeLists.txt

@@ -18,6 +18,7 @@ set(GRID_BASED_TESTS
   gridfunctionadaptortest
   h1functionalassemblertest
   integraloperatorassemblertest
+  istlbackendtest
   laplaceassemblertest
   polynomialtest
   secondorderassemblertest

dune/fufem/test/istlbackendtest.cc

+#include <config.h>
+#include <dune/common/test/testsuite.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/indices.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+
+#include <dune/fufem/assemblers/istlbackend.hh>
+
+// Due to a formerly unnoticed bug, activate bounds checking:
+#define DUNE_CHECK_BOUNDS 1
+
+using namespace Dune;
+
+template<class Entry, class Matrix, class RowIndex, class ColIndex>
+TestSuite testISTLBackend(Matrix& mat, RowIndex&& row, ColIndex&& col) {
+  TestSuite suite;
+
+  auto backend = Fufem::istlMatrixBackend<Entry>(mat);
+
+  suite.check(backend(row, col) == 1.0, "Check matrix access") << "Type" << typeid(mat).name() << " failed";
+
+  return suite;
+}
+
+int main(int argc, char** argv) {
+  MPIHelper::instance(argc, argv);
+
+  TestSuite suite;
+
+  using K = double;
+
+  // Matrix<FieldMatrix<K, 1,1>
+  {
+    using Matrix = Matrix<FieldMatrix<K,1,1>>;
+    Matrix mat(2,2);
+    mat =1.0;
+
+    // check with multiindex with length 1
+    {
+      auto idx = std::array<size_t, 1>{{1}};
+
+      suite.subTest(testISTLBackend<K>(mat, idx, idx));
+
+      // check const:
+      const auto& const_mat = mat;
+      suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
+    }
+    // check with multiindex with length 2
+    {
+      // someone might set a trailing zero to account for the matrix
+      // character of FieldMatrix<K, 1,1>:
+      auto idx = std::array<size_t, 2>{{1, 0}};
+
+      suite.subTest(testISTLBackend<K>(mat, idx, idx));
+
+      // check const:
+      const auto& const_mat = mat;
+      suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
+    }
+  }
+
+  // Matrix<FieldMatrix<K, n,n> n>1
+  {
+    constexpr int n = 2;
+    using Matrix = Matrix<FieldMatrix<K, n, n>>;
+    Matrix mat(2,2);
+    mat = 1.0;
+
+    auto idx = std::array<size_t, 2>{{1, 1}};
+
+    suite.subTest(testISTLBackend<K>(mat, idx, idx));
+
+    // check const:
+    const auto& const_mat = mat;
+    suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
+  }
+
+  return suite.exit();
+}