Merge branch 'fix/blocked-local-matrix' into 'master'

Fix: TrustRegionSolver: use blocked matrices for fufem scalar basis assembler

See merge request agnumpde/dune-elasticity!31

dune/elasticity/common/trustregionsolver.cc

@@ -2,6 +2,7 @@
 #include <dune/common/timer.hh>
 
 #include <dune/istl/io.hh>
+#include <dune/istl/scaledidmatrix.hh>
 
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 
@@ -62,6 +63,7 @@ setup(const typename BasisType::GridView::Grid& grid,
     damping_                  = damping;
 
     int numLevels = grid_->maxLevel()+1;
+    const auto dim = VectorType::value_type::dimension;
 
 #if HAVE_DUNE_PARMG
     Dune::Timer setupTimer;
@@ -143,13 +145,21 @@ setup(const typename BasisType::GridView::Grid& grid,
     {
       using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
       // construct a Fufem Basis Assembler
-      auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement>();
+      auto laplaceStiffness = LaplaceAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
       // transform it to a dune-functions assembler
       auto localAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
-          laplaceStiffness.assemble(element,
-                                    localMatrix,
-                                    trialLocalView.tree().child(0).finiteElement(),
-                                    ansatzLocalView.tree().child(0).finiteElement());
+        // the fufem basis assembler assumes a scalar basis but a blocked element matrix!
+        // we can use ScaledIdentityMatrix here
+        Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
+        laplaceStiffness.assemble(element,
+                                  localBlockedMatrix,
+                                  trialLocalView.tree().child(0).finiteElement(),
+                                  ansatzLocalView.tree().child(0).finiteElement());
+        // convert back to flat matrix (we need only the diagonal entries in the blocks)
+        localMatrix = 0;
+        for(size_t i=0; i<localBlockedMatrix.N(); i++)
+          for(size_t j=0; j<localBlockedMatrix.M(); j++)
+            Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
       };
 
 
@@ -185,13 +195,21 @@ setup(const typename BasisType::GridView::Grid& grid,
     {
       // construct a Fufem Basis Assembler
       using FiniteElement = std::decay_t<decltype(basis.localView().tree().child(0).finiteElement())>;
-      auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement>();
+      auto massStiffness = MassAssembler<GridType, FiniteElement, FiniteElement, Dune::ScaledIdentityMatrix<double,dim>>();
       // transform it to a dune-functions assembler
       auto localMassAssembler = [&](const auto& element, auto& localMatrix, auto&& trialLocalView, auto&& ansatzLocalView){
+        // the fufem basis assembler assumes a scalar basis but a blocked element matrix!
+        // we can use ScaledIdentityMatrix here
+        Dune::Matrix<Dune::ScaledIdentityMatrix<double,dim>> localBlockedMatrix( localMatrix.N()/dim, localMatrix.M()/dim );
         massStiffness.assemble(element,
-                               localMatrix,
+                               localBlockedMatrix,
                                trialLocalView.tree().child(0).finiteElement(),
                                ansatzLocalView.tree().child(0).finiteElement());
+        // convert back to flat matrix (we need only the diagonal entries in the blocks)
+        localMatrix = 0;
+        for(size_t i=0; i<localBlockedMatrix.N(); i++)
+          for(size_t j=0; j<localBlockedMatrix.M(); j++)
+            Dune::MatrixVector::addToDiagonal(localMatrix[i*dim][j*dim], localBlockedMatrix[i][j].scalar());
       };
 
       auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
@@ -256,7 +274,6 @@ setup(const typename BasisType::GridView::Grid& grid,
 
     if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
     {
-      const auto dim = VectorType::value_type::dimension;
       using namespace Dune::Functions::BasisFactory;
       auto p1PowerBasis = makeBasis(
           basis.gridView(),
@@ -282,7 +299,6 @@ setup(const typename BasisType::GridView::Grid& grid,
 
       if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
       {
-        const auto dim = VectorType::value_type::dimension;
         using namespace Dune::Functions::BasisFactory;
         auto p1PowerBasis = makeBasis(
             basis.gridView(),