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

@@ -62,6 +62,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 +144,19 @@ 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::FieldMatrix<double,dim,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!
+        Dune::Matrix<Dune::FieldMatrix<double,dim,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
+        for(size_t i=0; i<localMatrix.N(); i++)
+          for(size_t j=0; j<localMatrix.M(); j++)
+            localMatrix[i][j] = localBlockedMatrix[i/dim][j/dim][i%dim][j%dim];
       };
 
 
@@ -185,13 +192,19 @@ 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::FieldMatrix<double,dim,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!
+        Dune::Matrix<Dune::FieldMatrix<double,dim,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
+        for(size_t i=0; i<localMatrix.N(); i++)
+          for(size_t j=0; j<localMatrix.M(); j++)
+            localMatrix[i][j] = localBlockedMatrix[i/dim][j/dim][i%dim][j%dim];
       };
 
       auto massMatrixBackend = Dune::Fufem::istlMatrixBackend(localMassMatrix);
@@ -256,7 +269,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 +294,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(),