From f5f0bd788384ee0bbafecc90023aa682370a83ef Mon Sep 17 00:00:00 2001
From: Patrick <patrick.jaap@tu-dresden.de>
Date: Fri, 8 May 2020 17:58:10 +0200
Subject: [PATCH] Fix: TrustRegionSolver: use blocked matrices for fufem scalar
 basis assembler

valgrind complained here and it is right: in the internal loop of the
assembler we loop over localMatrix.N() and M(). Since dune-functions
uses flat matrices for the localMatrix the number of dof's is dim-times
too large for a powerBasis. In order to spare a double implementation of
the local assembler only for powerBasis, we can simply convert it in the
lambda. valgrind is now happy here.
---
 dune/elasticity/common/trustregionsolver.cc | 34 +++++++++++++++------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 5c7f540..ff75423 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/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(),
-- 
GitLab