Solve a sample problem using plain Gauss-Seidel

src/one-body-sample.cc

@@ -6,8 +6,11 @@
 
 #include <dune/common/exceptions.hh>
 
+#include <dune/common/bitsetvector.hh>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
 
 #include <dune/grid/yaspgrid.hh>
 
@@ -15,13 +18,23 @@
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 
+//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+
+#include <dune/solvers/iterationsteps/blockgsstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/common/numproc.hh> // Solver::FULL
+
 #include <exception>
 
 int const dim = 2;
 
 int main() {
   try {
-    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, dim, dim>> OperatorType;
+    typedef Dune::FieldVector<double, dim> SmallVector;
+    typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix;
+    typedef Dune::BCRSMatrix<SmallMatrix> OperatorType;
+    typedef Dune::BlockVector<SmallVector> VectorType;
 
     // FIXME: Random values
     double const E = 1e8;
@@ -52,6 +65,32 @@ int main() {
     OperatorType stiffnessMatrix;
     globalAssembler.assemble(localStiffness, stiffnessMatrix);
 
+    VectorType f(grid.size(grid.maxLevel(), dim));
+    // Fill right-hand side with semi-random numbers
+    for (size_t i = 0; i < f.size(); ++i)
+      for (size_t j = 0; j < dim; ++j)
+        f[i][j] = i + j;
+
+    VectorType u(grid.size(grid.maxLevel(), dim));
+    // Fill initial guess with other semi-random numbers
+    for (size_t i = 0; i < f.size(); ++i)
+      for (size_t j = 0; j < dim; ++j)
+        f[i][j] = i * j;
+
+    // TODO: Why does blockGSStep even provide a default constructor?
+    BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f);
+    EnergyNorm<OperatorType, VectorType> energyNorm(blockGSStep);
+
+    // TODO: Why does blockGSStep not warn about this?
+    // Why is this not documented?
+    Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(false);
+    blockGSStep.ignoreNodes_ = &ignoreNodes;
+
+    LoopSolver<VectorType> solver(&blockGSStep, 1000, // maxIterations
+                                  1e-5,               // tolerance
+                                  &energyNorm, Solver::FULL);
+
+    solver.solve();
     return 0;
   }
   catch (Dune::Exception& e) {