Implement BDMatrix<double>

dune/istl/bdmatrix.hh

@@ -32,7 +32,7 @@ namespace Dune {
     //===== type definitions and constants
 
     //! export the type representing the field
-    typedef typename B::field_type field_type;
+    using field_type = typename Imp::BlockTraits<B>::field_type;
 
     //! export the type representing the components
     typedef B block_type;
@@ -47,7 +47,7 @@ namespace Dune {
     typedef typename A::size_type size_type;
 
     //! increment block level counter
-    enum {blocklevel = B::blocklevel+1};
+    static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
 
     /** \brief Default constructor */
     BDMatrix() : BCRSMatrix<B,A>() {}
@@ -108,8 +108,15 @@ namespace Dune {
 
     /** \brief Inverts the matrix */
     void invert() {
-      for (size_type i=0; i<this->N(); i++)
-        (*this)[i][i].invert();
+      Hybrid::ifElse(IsNumber<B>(),
+        [&](auto id) {
+          for (size_type i=0; i<this->N(); i++)
+            (*this)[i][i] = 1.0 / id((*this)[i][i]);
+        },
+        [&](auto id) {
+          for (size_type i=0; i<this->N(); i++)
+            id((*this)[i][i]).invert();
+        });
     }
 
   private:

dune/istl/test/matrixtest.cc

@@ -458,8 +458,34 @@ int main()
 
   testSuperMatrix(bcrsMatrix);
 
-  // ////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////////////////////////////////////////////////
   //   Test the BDMatrix class -- a dynamic block-diagonal matrix
+  ///////////////////////////////////////////////////////////////////////////
+
+  {
+    BDMatrix<double> bdMatrix(3);
+    bdMatrix = 4.0;
+
+    BlockVector<double> x(3), y(3);
+    testMatrix(bdMatrix, x, y);
+
+    // Test construction from initializer list
+    BDMatrix<double> bdMatrix2 = {1.0, 2.0, 3.0};
+    testMatrix(bdMatrix2, x, y);
+
+    // test whether resizing works
+    bdMatrix2.setSize(5);
+    bdMatrix2 = 4.0;
+    x.resize(5);
+    y.resize(5);
+    testMatrix(bdMatrix2, x, y);
+
+    // Test whether inversion works
+    bdMatrix2.invert();
+  }
+
+  // ////////////////////////////////////////////////////////////////////////
+  //   Test the BDMatrix class with FieldMatrix entries
   // ////////////////////////////////////////////////////////////////////////
 
   BDMatrix<FieldMatrix<double,4,4> > bdMatrix(2);
@@ -469,6 +495,11 @@ int main()
 
   // Test construction from initializer list
   BDMatrix<FieldMatrix<double,2,2> > bdMatrix2 = { {{1,0},{0,1}}, {{0,1},{-1,0}}};
+
+  // Test whether inversion works
+  bdMatrix2.invert();
+
+  // Run matrix tests on this matrix
   testSuperMatrix(bdMatrix2);
 
   // test whether resizing works