Use a TruncatedBlockGaussSeidelStep for comparison

It has the advantage of handling nontrivial restrictions to the DOFs
at a node gracefully. We need this for u_n|S_F

src/one-body-sample.cc

@@ -35,6 +35,7 @@
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/solvers/common/numproc.hh> // Solver::FULL
 #include <dune/solvers/iterationsteps/blockgsstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
@@ -188,10 +189,12 @@ int main() {
     VectorType u1(grid.size(grid.maxLevel(), dim));
     u1 = 0;
     VectorType u2 = u1;
+    VectorType u3 = u1;
     CellVectorType vonMisesStress;
 
     VectorType b1;
     VectorType b2;
+    VectorType b3;
 
     std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
     assemble_frictional<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
@@ -204,11 +207,12 @@ int main() {
       // b = neumann
       assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
           leafView, p1Basis, neumannNodes, b1, run);
-      b2 = b1;
+      b3 = b2 = b1;
 
       // b += linear update
       stiffnessMatrix.umv(u1, b1);
       stiffnessMatrix.umv(u2, b2);
+      stiffnessMatrix.umv(u3, b3);
 
       // {{{ Assemble terms for the nonlinearity
       // TODO: Random value
@@ -272,20 +276,37 @@ int main() {
                                       Solver::QUIET); // Solver::QUIET);
         solver.solve();
       }
+
+      {
+        TruncatedBlockGSStep<OperatorType, VectorType> blockGSStep(
+            stiffnessMatrix, u3, b3);
+        blockGSStep.ignoreNodes_ = &ignoreNodes;
+
+        LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
+                                      solver_tolerance, &energyNorm,
+                                      Solver::QUIET); // Solver::QUIET);
+        solver.solve();
+      }
     }
     std::cout << std::endl;
 
-    VectorType diff = u2;
-    diff -= u1;
-    std::cout << "Infinity norm of the difference of the two solutions: "
-              << diff.infinity_norm() << std::endl;
+    VectorType diff2 = u2;
+    diff2 -= u1;
+    std::cout << "sup |u1 - u2| = " << diff2.infinity_norm() << std::endl;
+    std::cout << "|u1 - u2| = " << diff2.two_norm() << std::endl;
+
+    VectorType diff3 = u3;
+    diff3 -= u1;
+    std::cout << "sup |u1 - u3| = " << diff3.infinity_norm() << std::endl;
+    std::cout << "|u1 - u3| = " << diff3.two_norm() << std::endl;
 
     // Print displacement on frictional boundary
     for (size_t i = 0; i < frictionalNodes.size(); ++i)
       if (frictionalNodes[i][0])
-        std::cout << boost::format("u1[%02d] = %+3e, %|40t|u2[%02d] = %+3e,    "
-                                   " s[%02d] = %+3e") %
-                         i % u1[i] % i % u2[i] % i %
+        std::cout << boost::format("u1[%02d] = %+3e, %|40t|u2[%02d] = %+3e, "
+                                   "%|80t|u3[%02d] = %+3e, %|120t|s[%02d] = "
+                                   "%+3e") %
+                         i % u1[i] % i % u2[i] % i % u3[i] % i %
                          vonMisesStress[i] << std::endl;
   }
   catch (Dune::Exception &e) {