Be smarter about what parts of the coarse grid matrices to reassemble

In the truncated monotone multigrid method, the coarse grid matrices change from
one iteration to the next.  Therefore, they need to be reassembled before every
iteration.  That is pretty expensive, therefore the MonotoneMG class and
corresponding transfer operators have always included a mechanism to only reassemble
parts of the coarse grid matrices.  That saves quite a bit of time.

Previously, a matrix entry would only be recomputed if the hasObstacle field
was set for its row or column.  That works nicely for contact problems,
where the hasObstacle field has only very few bits set.

However, for problems with many obstacles the coarse matrix reassembly becomes
expensive again.  Therefore, the code is now a bit smarter about what to reassemble.
Matrix entries really only need reassemble if the corresponding dofs changed
from non-critical to critical (or vice versa), or if the next-finer matrix
changed in the corresponding entries.

For my Cosserat shell problem, which has 6x6 blocks and all hasObstacle entries
set, I measure a speedup of about 30% for each MMG iteration.

Note that after this patch, code should yield precisely the same output as before.
If you get different results compare to before, then that's a bug.

dune/solvers/iterationsteps/mmgstep.cc

@@ -31,8 +31,19 @@ preprocess()
                 (*hasObstacleHierarchy_[i])[j][0]=false;
     }
 
+    recompute_.resize(this->mgTransfer_.size());
     for (size_t i=0; i<this->mgTransfer_.size(); i++)
-        dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[i])->setRecomputeBitField(hasObstacleHierarchy_[i]);
+    {
+        recompute_[i].resize(hasObstacleHierarchy_[i]->size(), false);
+        dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
+    }
+
+    oldCritical_.resize(this->numLevels());
+    for (size_t i=0; i<this->numLevels(); i++)
+    {
+      oldCritical_[i].resize(hasObstacleHierarchy_[i]->size());
+      oldCritical_[i].unsetAll();
+    }
     
     // /////////////////////////////////////////////
     //   Set up base solver
@@ -183,6 +194,27 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
             }
         }
 
+        ///////////////////////////////////////////////////////////////////////////////
+        // Compute the part of the coarse grid matrix that needs to be recomputed.
+        // There are two reasons why a matrix entry may need to be recomputed:
+        // 1) A corresponding fine grid vertex switched from critical to non-critical or vice versa
+        // 2) A corresponding fine grid matrix entry got recomputed
+        ///////////////////////////////////////////////////////////////////////////////
+
+        Dune::BitSetVector<dim> changed(critical.size());
+        for (size_t i=0; i<changed.size(); i++)
+          for (int j=0; j<dim; j++)
+            changed[i][j] = (critical[i][j] != oldCritical_[level][i][j]);
+
+        if (level < this->numLevels()-1 )
+          for (size_t i=0; i<changed.size(); i++)
+            for (int j=0; j<dim; j++)
+              changed[i][j] = (changed[i][j] || recompute_[level][i][j]);
+
+        dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[level-1])->restrict(changed, recompute_[level-1]);
+
+        oldCritical_[level] = critical;
+
         // Set bitfield of nodes that will be truncated
         dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[level-1])->setCriticalBitField(&critical);
      

dune/solvers/iterationsteps/mmgstep.hh

@@ -64,6 +64,12 @@ protected:
     //! Bitfield hierarchy containing the coarse obstacle nodes
     std::vector<Dune::BitSetVector<dim>* > hasObstacleHierarchy_;
 
+    //! For each level specify the matrix lines that need reassembly
+    std::vector<Dune::BitSetVector<dim> > recompute_;
+
+    //! For each level specify dofs that changed status
+    std::vector<Dune::BitSetVector<dim> > oldCritical_;
+
     //! Hierarchy containing the obstacle values of the coarse obstacles
     std::vector<ObstacleVectorType*>  obstacleHierarchy_;
 };