Make multiple runs

src/one-body-sample.cc

@@ -128,6 +128,7 @@ int main() {
     typedef Dune::BlockVector<SmallVector> VectorType;
 
     // FIXME: Random values
+    size_t const runs = 3;
     double const E = 1;
     double const nu = 0.3;
     int const refinements = 5;
@@ -179,55 +180,60 @@ int main() {
     u1 = 0;
     VectorType u2 = u1;
 
-    VectorType f(grid.size(grid.maxLevel(), dim));
-    f = 0;
-    VectorType neumannTerm(grid.size(grid.maxLevel(), dim));
-    assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
-        leafView, p1Basis, neumannNodes, neumannTerm);
-    f += neumannTerm;
-
-    // {{{ Assemble terms for the nonlinearity
-    std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
-    assemble_frictional<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
-        leafView, p1Basis, frictionalNodes, nodalIntegrals);
-
-    // TODO: populate on S_F
-    std::vector<double> normalStress;
-    normalStress.resize(grid.size(grid.maxLevel(), dim));
-    std::fill(normalStress.begin(), normalStress.end(), 0.0);
-
-    // TODO: populate on S_F
-    std::vector<double> coefficientOfFriction;
-    coefficientOfFriction.resize(grid.size(grid.maxLevel(), dim));
-    std::fill(coefficientOfFriction.begin(), coefficientOfFriction.end(), 0.0);
-
-    Dune::GlobalNonlinearity<dim, Dune::LinearFunction> myGlobalNonlinearity(
-        coefficientOfFriction, normalStress, nodalIntegrals);
-    // }}}
-
-    {
-      MyConvexProblemType myConvexProblem(stiffnessMatrix, myGlobalNonlinearity,
-                                          f, u1);
-      MyBlockProblemType myBlockProblem(myConvexProblem);
-      nonlinearGSStep.setProblem(u1, myBlockProblem);
-
-      LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
-                                    solver_tolerance, &energyNorm,
-                                    Solver::FULL); // Solver::QUIET);
-      solver.solve();
-    }
-
-    // Use a linear solver for comparison; should return roughly the
-    // same results if phi vanishes (e.g. because the normalstress is zero)
-    {
-      // TODO: Why does blockGSStep even provide a default constructor?
-      BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u2, f);
-      blockGSStep.ignoreNodes_ = &ignoreNodes;
-
-      LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
-                                    solver_tolerance, &energyNorm,
-                                    Solver::FULL); // Solver::QUIET);
-      solver.solve();
+    for (size_t run = 1; run <= runs; ++run) {
+      VectorType f(grid.size(grid.maxLevel(), dim));
+      f = 0;
+      VectorType neumannTerm(grid.size(grid.maxLevel(), dim));
+      assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
+          leafView, p1Basis, neumannNodes, neumannTerm);
+      f += neumannTerm;
+
+      // {{{ Assemble terms for the nonlinearity
+      std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
+      assemble_frictional<GridType, GridType::LeafGridView, SmallVector,
+                          P1Basis>(leafView, p1Basis, frictionalNodes,
+                                   nodalIntegrals);
+
+      // TODO: populate on S_F
+      std::vector<double> normalStress;
+      normalStress.resize(grid.size(grid.maxLevel(), dim));
+      std::fill(normalStress.begin(), normalStress.end(), 0.0);
+
+      // TODO: populate on S_F
+      std::vector<double> coefficientOfFriction;
+      coefficientOfFriction.resize(grid.size(grid.maxLevel(), dim));
+      std::fill(coefficientOfFriction.begin(), coefficientOfFriction.end(),
+                0.0);
+
+      Dune::GlobalNonlinearity<dim, Dune::LinearFunction> myGlobalNonlinearity(
+          coefficientOfFriction, normalStress, nodalIntegrals);
+      // }}}
+
+      {
+        MyConvexProblemType myConvexProblem(stiffnessMatrix,
+                                            myGlobalNonlinearity, f, u1);
+        MyBlockProblemType myBlockProblem(myConvexProblem);
+        nonlinearGSStep.setProblem(u1, myBlockProblem);
+
+        LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
+                                      solver_tolerance, &energyNorm,
+                                      Solver::FULL); // Solver::QUIET);
+        solver.solve();
+      }
+
+      // Use a linear solver for comparison; should return roughly the
+      // same results if phi vanishes (e.g. because the normalstress is zero)
+      {
+        // TODO: Why does blockGSStep even provide a default constructor?
+        BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u2,
+                                                          f);
+        blockGSStep.ignoreNodes_ = &ignoreNodes;
+
+        LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,
+                                      solver_tolerance, &energyNorm,
+                                      Solver::FULL); // Solver::QUIET);
+        solver.solve();
+      }
     }
 
     VectorType diff = u2;