Bring back BlockGSStep now that it works

Do not print vonMisesStress per node
That doesn't make sense -- we only have it for elements!

src/one-body-sample.cc

@@ -40,6 +40,7 @@
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #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>
@@ -199,19 +200,23 @@ int main(int argc, char *argv[]) {
 
     VectorType u1(grid.size(grid.maxLevel(), dim));
     u1 = 0;
+    VectorType u2 = u1;
     VectorType u3 = u1;
 
     VectorType u1_diff_old(grid.size(grid.maxLevel(), dim));
     u1_diff_old = 0;
+    VectorType u2_diff_old = u1_diff_old;
     VectorType u3_diff_old = u1_diff_old;
 
     VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
     u1_diff_new = 0;
+    VectorType u2_diff_new = u1_diff_new;
     VectorType u3_diff_new = u1_diff_new;
 
     CellVectorType vonMisesStress;
 
     VectorType b1;
+    VectorType b2;
     VectorType b3;
 
     std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
@@ -225,10 +230,12 @@ int main(int argc, char *argv[]) {
       // b = neumann
       assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
           leafView, p1Basis, neumannNodes, b1, run);
+      b2 = b1;
       b3 = b1;
 
       // b += linear update
       stiffnessMatrix.umv(u1_diff_old, b1);
+      stiffnessMatrix.umv(u2_diff_old, b2);
       stiffnessMatrix.umv(u3_diff_old, b3);
 
       // {{{ Assemble terms for the nonlinearity
@@ -309,6 +316,20 @@ int main(int argc, char *argv[]) {
         writer.write(filename.c_str());
       }
 
+      {
+        BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix,
+                                                          u2_diff_new, b2);
+        blockGSStep.ignoreNodes_ = &ignoreNodes;
+
+        LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
+                                      solver_tolerance, &energyNorm,
+                                      Solver::QUIET); // Solver::QUIET);
+        solver.solve();
+      }
+
+      u2 += u2_diff_new;
+      u2_diff_old = u2_diff_new;
+
       // Use a linear solver for comparison; should return roughly the
       // same results if phi vanishes (e.g. because the normalstress is zero)
       {
@@ -327,18 +348,22 @@ int main(int argc, char *argv[]) {
     }
     std::cout << std::endl;
 
+    VectorType diff2 = u2;
+    diff2 -= u1;
+    std::cout << "sup |u1 - u2| = " << diff2.infinity_norm() << ", "
+              << "|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;
+    std::cout << "sup |u1 - u3| = " << diff3.infinity_norm() << ", "
+              << "|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|u3[%02d] = %+3e "
-                                   "%|80t|s[%02d] = %+3e") %
-                         i % u1[i] % i % u3[i] % i %
-                         vonMisesStress[i] << std::endl;
+        std::cout << boost::format("u1[%02d] = %+3e, %|40t|u2[%02d] = %+3e "
+                                   "%|80t|u3[%02d] = %+3e") %
+                         i % u1[i] % i % u2[i] % i % u3[i] << std::endl;
   }
   catch (Dune::Exception &e) {
     Dune::derr << "Dune reported error: " << e << std::endl;