[Cleanup] Add solveLinearProblem(), move declarations

src/one-body-sample.cc

@@ -247,9 +247,6 @@ int main(int argc, char *argv[]) {
     Vector gravityFunctional;
     myAssembler.assembleBodyForce(gravity, density, zenith, gravityFunctional);
 
-    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-        frictionalBoundary, frictionData);
-
     // Problem formulation: right-hand side
     auto const createRHS = [&](double _relativeTime, Vector &_ell) {
       myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
@@ -263,6 +260,9 @@ int main(int argc, char *argv[]) {
     ScalarVector alpha_initial(fineVertexCount);
     alpha_initial =
         std::log(parset.get<double>("boundary.friction.initialState"));
+    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
+        frictionalBoundary, frictionData);
+    myGlobalNonlinearity->updateLogState(alpha_initial);
 
     using LinearFactory = SolverFactory<
         dims, BlockNonlinearTNNMGProblem<ConvexProblem<
@@ -270,31 +270,35 @@ int main(int argc, char *argv[]) {
         Grid>;
     ZeroNonlinearity<SmallVector, SmallMatrix> zeroNonlinearity;
 
-    // Solve the stationary problem
-    Vector u_initial(fineVertexCount);
-    u_initial = 0.0;
-    {
-      LinearFactory displacementFactory(parset.sub("solver.tnnmg"), // FIXME
-                                        refinements, *grid,
-                                        displacementDirichletNodes);
-      auto multigridStep = displacementFactory.getSolver();
+    // TODO: clean up once generic lambdas arrive
+    auto const solveLinearProblem = [&](
+        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
+        Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> _norm,
+        Dune::ParameterTree const &_localParset) {
+      LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
+                            refinements, *grid, _dirichletNodes);
 
       typename LinearFactory::ConvexProblem convexProblem(
-          1.0, A, zeroNonlinearity, ell, u_initial);
-      typename LinearFactory::BlockProblem initialDisplacementProblem(
-          parset, convexProblem);
-
-      auto const &localParset = parset.sub("u0.solver");
-      multigridStep->setProblem(u_initial, initialDisplacementProblem);
-      LoopSolver<Vector> initialDisplacementProblemSolver(
-          multigridStep, localParset.get<size_t>("maximumIterations"),
-          localParset.get<double>("tolerance"), &ANorm,
-          localParset.get<Solver::VerbosityMode>("verbosity"),
+          1.0, _matrix, zeroNonlinearity, _rhs, _x);
+      typename LinearFactory::BlockProblem problem(parset, convexProblem);
+
+      auto multigridStep = factory.getSolver();
+      multigridStep->setProblem(_x, problem);
+      LoopSolver<Vector> solver(
+          multigridStep, _localParset.get<size_t>("maximumIterations"),
+          _localParset.get<double>("tolerance"), &_norm,
+          _localParset.get<Solver::VerbosityMode>("verbosity"),
           false); // absolute error
 
-      initialDisplacementProblemSolver.preprocess();
-      initialDisplacementProblemSolver.solve();
-    }
+      solver.preprocess();
+      solver.solve();
+    };
+
+    // Solve the stationary problem
+    Vector u_initial(fineVertexCount);
+    u_initial = 0.0;
+    solveLinearProblem(displacementDirichletNodes, A, ell, u_initial, ANorm,
+                       parset.sub("u0.solver"));
     Vector v_initial(fineVertexCount);
     v_initial = 0.0;
     {
@@ -317,32 +321,13 @@ int main(int argc, char *argv[]) {
         accelerationRHS = 0.0;
         Arithmetic::addProduct(accelerationRHS, A, u_initial);
         Arithmetic::addProduct(accelerationRHS, C, v_initial);
-        myGlobalNonlinearity->updateLogState(alpha_initial);
         // NOTE: We assume differentiability of Psi at v0 here!
         myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
         accelerationRHS -= ell;
         accelerationRHS *= -1.0;
       }
-      LinearFactory accelerationFactory(parset.sub("solver.tnnmg"), // FIXME
-                                        refinements, *grid,
-                                        accelerationDirichletNodes);
-      auto multigridStep = accelerationFactory.getSolver();
-
-      typename LinearFactory::ConvexProblem convexProblem(
-          1.0, M, zeroNonlinearity, accelerationRHS, a_initial);
-      typename LinearFactory::BlockProblem initialAccelerationProblem(
-          parset, convexProblem);
-
-      auto const &localParset = parset.sub("a0.solver");
-      multigridStep->setProblem(a_initial, initialAccelerationProblem);
-      LoopSolver<Vector> initialAccelerationProblemSolver(
-          multigridStep, localParset.get<size_t>("maximumIterations"),
-          localParset.get<double>("tolerance"), &MNorm,
-          localParset.get<Solver::VerbosityMode>("verbosity"),
-          false); // absolute error
-
-      initialAccelerationProblemSolver.preprocess();
-      initialAccelerationProblemSolver.solve();
+      solveLinearProblem(accelerationDirichletNodes, M, accelerationRHS,
+                         a_initial, MNorm, parset.sub("a0.solver"));
     }
     // }}}
     FrictionWriter<Dune::BitSetVector<1>> writer(frictionData, frictionalNodes);