[Algorit] FIX: Calculate the initial acceleration

src/one-body-sample.cc

@@ -59,6 +59,7 @@
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/sharedpointermap.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/norms/sumnorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
@@ -342,9 +343,33 @@ int main(int argc, char *argv[]) {
     u_initial = 0.0;
     VectorType v_initial(finestSize);
     v_initial = 0.0;
-    // FIXME: This only happens to be correct
+
+    // From here on: calculate the initial acceleration
     VectorType a_initial(finestSize);
     a_initial = 0.0;
+    {
+      // We solve Au + Ma + Psi(v) = ell, thus Ma = - (Au + Psi(v) - ell)
+      VectorType problem_rhs_initial(finestSize);
+      {
+        problem_rhs_initial = 0.0;
+        Arithmetic::addProduct(problem_rhs_initial, stiffnessMatrix, u_initial);
+        myGlobalNonlinearity->updateState(alpha_initial);
+        // NOTE: We assume differentiability of Psi at v0 here!
+        myGlobalNonlinearity->addGradient(v_initial, problem_rhs_initial);
+        createRHS(0.0, ell);
+        problem_rhs_initial -= ell;
+        problem_rhs_initial *= -1.0;
+      }
+      TruncatedBlockGSStep<MatrixType, VectorType> accelerationSolverStep(
+          massMatrix, a_initial, problem_rhs_initial);
+      accelerationSolverStep.ignoreNodes_ = &ignoreNodes;
+      LoopSolver<VectorType> accelerationSolver(
+          &accelerationSolverStep, 100000, // FIXME
+          1e-12,                           // FIXME
+          &massMatrixNorm, Solver::FULL,   // FIXME
+          false);                          // absolute error
+      accelerationSolver.solve();
+    }
     // }}}
 
     // Set up TNNMG solver