Add a second solution: A nonlinear Gauss-Seidel

With a zero nonlinearity

src/one-body-sample.cc

@@ -19,7 +19,7 @@
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 
-//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
 
 #include <dune/solvers/iterationsteps/blockgsstep.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
@@ -28,9 +28,15 @@
 
 #include <dune/grid/common/mcmgmapper.hh>
 
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+#include <dune/tnnmg/problem-classes/nonlinearity.hh>
+#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
+
 #include <exception>
 #include <iostream>
 
+#include "myblockproblem.hh"
+
 int const dim = 2;
 
 int main() {
@@ -131,6 +137,36 @@ int main() {
       std::cout << "Number of bounding nodes: " << bounding_nodes << std::endl;
       std::cout << "Number of extremal nodes: " << extremal_nodes << std::endl;
     }
+
+    VectorType copy_of_u = u;
+
+    { // experiment with convex problems and the like
+
+      typedef ZeroNonlinearity<SmallVector, SmallMatrix> NonlinearityType;
+      NonlinearityType phi;
+
+      typedef ConvexProblem<NonlinearityType, OperatorType> ConvexProblemType;
+
+      VectorType unused_vector(grid.size(grid.maxLevel(), dim), 0);
+      ConvexProblemType myConvexProblem(1, stiffnessMatrix, 0, unused_vector,
+                                        phi, f, copy_of_u);
+
+      typedef MyBlockProblem<ConvexProblemType> MyBlockProblemType;
+      MyBlockProblemType myBlockProblem(myConvexProblem);
+
+      GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep;
+      nonlinearGSStep.ignoreNodes_ = &ignoreNodes;
+      nonlinearGSStep.setProblem(copy_of_u, myBlockProblem);
+
+      // FIXME: Does this make any sense?
+      EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix);
+      LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
+                                    solver_tolerance, &energyNorm,
+                                    // Solver::QUIET);
+                                    Solver::FULL);
+      solver.solve();
+    }
+
     // TODO: Why does blockGSStep even provide a default constructor?
     BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f);
     blockGSStep.ignoreNodes_ = &ignoreNodes;
@@ -138,10 +174,17 @@ int main() {
     // FIXME: Does this make any sense?
     EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix);
     LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
-                                  solver_tolerance, &energyNorm, Solver::FULL);
+                                  solver_tolerance, &energyNorm,
+                                  // Solver::QUIET);
+                                  Solver::FULL);
 
     solver.solve();
-    return 0;
+
+    VectorType diff = u;
+    diff -= copy_of_u;
+    std::cout << diff << std::endl;
+    std::cout << std::endl;
+    std::cout << diff.infinity_norm() << std::endl;
   }
   catch (Dune::Exception& e) {
     Dune::derr << "Dune reported error: " << e << std::endl;