Cleanup

Move code around, add a comment

src/one-body-sample.org

@@ -256,18 +256,13 @@
 #+name: mainBody
 #+begin_src c++
   try {
-    typedef SharedPointerMap
-      <std::string, Dune::VirtualFunction<double, double> > FunctionMap;
-    FunctionMap functions;
-    initPython(functions);
+    Dune::Timer timer;
   
     Dune::ParameterTree parset;
     Dune::ParameterTreeParser::readINITree
       (srcdir "/one-body-sample.parset", parset);
     Dune::ParameterTreeParser::readOptions(argc, argv, parset);
   
-    Dune::Timer timer;
-  
     typedef Dune::FieldVector<double, dim> SmallVector;
     typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix;
     typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
@@ -276,13 +271,11 @@
   
     auto const E = parset.get<double>("body.E");
     auto const nu = parset.get<double>("body.nu");
-    auto const solver_tolerance = parset.get<double>("solver.tolerance");
-  
-    auto const refinements = parset.get<size_t>("grid.refinements");
-  
-    auto const verbose = parset.get<bool>("verbose");
-    Solver::VerbosityMode const verbosity
-      = verbose ? Solver::FULL : Solver::QUIET;
+    auto const mu = parset.get<double>("boundary.friction.mu");
+    auto const a = parset.get<double>("boundary.friction.a");
+    auto const b = parset.get<double>("boundary.friction.b");
+    auto const V0 = parset.get<double>("boundary.friction.V0");
+    auto const L = parset.get<double>("boundary.friction.L");
   
     // {{{ Set up grid
     typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming>
@@ -301,6 +294,7 @@
     auto grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid
       (lowerLeft, upperRight, elements);
   
+    auto const refinements = parset.get<size_t>("grid.refinements");
     grid->globalRefine(refinements);
     size_t const finestSize = grid->size(grid->maxLevel(), dim);
   
@@ -314,6 +308,22 @@
     P0Basis const p0Basis(leafView);
     P1Basis const p1Basis(leafView);
   
+    // Set up the boundary
+    size_t specialNode = finestSize;
+    Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
+    Dune::BitSetVector<1> neumannNodes(finestSize, false);
+    Dune::BitSetVector<1> frictionalNodes(finestSize, false);
+    <<setupBoundary>>;
+  
+    // Set up functions for time-dependent boundary conditions
+    typedef SharedPointerMap
+      <std::string, Dune::VirtualFunction<double, double> > FunctionMap;
+    FunctionMap functions;
+    initPython(functions);
+    auto const &dirichletFunction = functions.get("dirichletCondition");
+    auto const &neumannFunction = functions.get("neumannCondition");
+  
+    // Set up normal stress, mass matrix, and gravity functional
     double normalStress;
     MatrixType massMatrix;
     VectorType gravityFunctional;
@@ -325,18 +335,16 @@
       <<computeNormalStress>>;
       <<computeGravitationalBodyForce>>;
     }
+    SingletonVectorType surfaceNormalStress(finestSize);
+    surfaceNormalStress = 0.0;
+    for (size_t i = 0; i < frictionalNodes.size(); ++i)
+      if (frictionalNodes[i][0])
+        surfaceNormalStress[i] = normalStress;
   
     MatrixType stiffnessMatrix;
     <<assembleStiffnessMatrix>>;
     EnergyNorm<MatrixType, VectorType> energyNorm(stiffnessMatrix);
   
-    // Set up the boundary
-    size_t specialNode = finestSize;
-    Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
-    Dune::BitSetVector<1> neumannNodes(finestSize, false);
-    Dune::BitSetVector<1> frictionalNodes(finestSize, false);
-    <<setupBoundary>>;
-  
     auto const nodalIntegrals
       = assemble_frictional<GridType, GridView, SmallVector, P1Basis>
       (leafView, p1Basis, frictionalNodes);
@@ -354,35 +362,22 @@
     SingletonVectorType alpha(alpha_initial);
   
     SingletonVectorType vonMisesStress;
-  
-    SingletonVectorType surfaceNormalStress(finestSize);
-    surfaceNormalStress = 0.0;
-    for (size_t i = 0; i < frictionalNodes.size(); ++i)
-      if (frictionalNodes[i][0])
-        surfaceNormalStress[i] = normalStress;
     // }}}
   
+    // Set up TNNMG solver
+    auto const solver_tolerance = parset.get<double>("solver.tolerance");
     typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
     typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
-  
-    // Set up TNNMG solver
     MySolver<dim, MatrixType, VectorType, GridType, MyBlockProblemType>
       mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance,
                *grid, ignoreNodes);
   
-    <<createWriters>>;
-  
-    auto const L = parset.get<double>("boundary.friction.L");
-    auto const a = parset.get<double>("boundary.friction.a");
-    auto const b = parset.get<double>("boundary.friction.b");
-    auto const V0 = parset.get<double>("boundary.friction.V0");
-    auto const mu = parset.get<double>("boundary.friction.mu");
-    auto const timesteps = parset.get<size_t>("timeSteps");
-    auto const tau = parset.get<double>("endOfTime") / timesteps;
+    Solver::VerbosityMode const verbosity
+      = parset.get<bool>("verbose") ? Solver::FULL : Solver::QUIET;
   
-    auto const &dirichletFunction = functions.get("dirichletCondition");
-    auto const &neumannFunction = functions.get("neumannCondition");
+    <<createWriters>>;
   
+    // Set up time steppers that
     auto timeSteppingScheme = initTimeStepper
       (parset.get<Config::scheme>("timeSteppingScheme"), dirichletFunction,
        ignoreNodes, massMatrix, stiffnessMatrix, u_initial, ud_initial,
@@ -391,6 +386,9 @@
       (parset.get<Config::state_model>("boundary.friction.state_model"),
        alpha_initial, frictionalNodes, L);
   
+    auto const timesteps = parset.get<size_t>("timeSteps");
+    auto const tau = parset.get<double>("endOfTime") / timesteps;
+  
     auto const state_fpi_max
       = parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
     for (size_t run = 1; run <= timesteps; ++run) {
@@ -412,6 +410,8 @@
       timeSteppingScheme->setup(ell, tau, time,
                                 problem_rhs, problem_iterate, problem_A);
   
+      // Since the velocity explodes in the quasistatic case, use the
+      // displacement as a convergence criterion
       VectorType u_saved;
       for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
         <<setupAndSolveProblem>>;