Move functions to separate babel blocks

src/one-body-sample.org

@@ -186,7 +186,273 @@
   }
 #+end_src
 * Functions
-#+name: main
+#+name: defineInitTimeStepper
+#+begin_src c++
+  template<class VectorType, class MatrixType, class FunctionType, int dim>
+  Dune::shared_ptr
+  <TimeSteppingScheme
+   <VectorType,MatrixType,FunctionType,dim> >
+  initTimeStepper(Config::scheme scheme,
+                  FunctionType const &dirichletFunction,
+                  Dune::BitSetVector<dim> const &ignoreNodes,
+                  MatrixType const &massMatrix, MatrixType const &stiffnessMatrix,
+                  VectorType const &u_initial,
+                  VectorType const &ud_initial,
+                  VectorType const &udd_initial)
+  {
+    switch (scheme) {
+    case Config::ImplicitTwoStep:
+      return Dune::make_shared<ImplicitTwoStep<VectorType,MatrixType,
+                                               FunctionType,dim> >
+        (stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
+    case Config::ImplicitEuler:
+      return Dune::make_shared<ImplicitEuler<VectorType,MatrixType,
+                                             FunctionType,dim> >
+        (stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
+    case Config::Newmark:
+      return Dune::make_shared <Newmark <VectorType,MatrixType,FunctionType,dim> >
+        (stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
+         ignoreNodes, dirichletFunction);
+    }
+  }
+#+end_src
+#+name: defineInitStateUpdater
+#+begin_src c++
+  template<class SingletonVectorType, class VectorType>
+  Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType> >
+  initStateUpdater(Config::state_model model,
+                   SingletonVectorType const &alpha_initial,
+                   Dune::BitSetVector<1> const &frictionalNodes,
+                   double L)
+  {
+    switch (model) {
+    case Config::Dieterich:
+      return Dune::make_shared
+        <DieterichStateUpdater<SingletonVectorType, VectorType> >
+        (alpha_initial, frictionalNodes, L);
+    case Config::Ruina:
+      return Dune::make_shared
+        <RuinaStateUpdater<SingletonVectorType, VectorType> >
+        (alpha_initial, frictionalNodes, L);
+    }
+  }
+#+end_src
+#+name: defineInitPython
+#+begin_src c++
+  template <class FunctionMap>
+  void
+  initPython(FunctionMap &functions)
+  {
+    Python::start();
+  
+    Python::run("import sys");
+    Python::run("sys.path.append('" srcdir "')");
+  
+    Python::import("one-body-sample")
+      .get("Functions")
+      .toC<typename FunctionMap::Base>(functions);
+  }
+#+end_src
+#+name: mainBody
+#+begin_src c++
+  try {
+    typedef SharedPointerMap
+      <std::string, Dune::VirtualFunction<double, double> > FunctionMap;
+    FunctionMap functions;
+    initPython(functions);
+  
+    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;
+    typedef Dune::BlockVector<SmallVector> VectorType;
+    typedef Dune::BlockVector<Dune::FieldVector<double,1> > SingletonVectorType;
+  
+    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;
+  
+    // {{{ Set up grid
+    typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming>
+      GridType;
+  
+    Dune::FieldVector<typename GridType::ctype, dim> lowerLeft(0);
+    Dune::FieldVector<typename GridType::ctype, dim> upperRight(1);
+    upperRight[0] = parset.get<size_t>("body.width");
+    upperRight[1] = parset.get<size_t>("body.height");
+  
+    Dune::array<unsigned int, dim> elements;
+    std::fill(elements.begin(), elements.end(), 1);
+    elements[0] = parset.get<size_t>("body.width");
+    elements[1] = parset.get<size_t>("body.height");
+  
+    auto grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid
+      (lowerLeft, upperRight, elements);
+  
+    grid->globalRefine(refinements);
+    size_t const finestSize = grid->size(grid->maxLevel(), dim);
+  
+    typedef GridType::LeafGridView GridView;
+    GridView const leafView = grid->leafView();
+    // }}}
+  
+    // Set up bases
+    typedef P0Basis<GridView, double> P0Basis;
+    typedef P1NodalBasis<GridView, double> P1Basis;
+    P0Basis const p0Basis(leafView);
+    P1Basis const p1Basis(leafView);
+  
+    double normalStress;
+    MatrixType massMatrix;
+    VectorType gravityFunctional;
+    {
+      double const gravity = 9.81;
+      double const density = parset.get<double>("body.density");
+  
+      <<assembleMassMatrix>>;
+      <<computeNormalStress>>;
+      <<computeGravitationalBodyForce>>;
+    }
+  
+    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);
+  
+    // {{{ Initialise vectors
+    VectorType u(finestSize); u = 0.0;
+    VectorType ud(finestSize); ud = 0.0;
+  
+    VectorType u_initial(finestSize); u_initial = 0.0;
+    VectorType ud_initial(finestSize); ud_initial = 0.0;
+    VectorType udd_initial(finestSize); udd_initial = 0.0;
+  
+    SingletonVectorType alpha_initial(finestSize);
+    alpha_initial = parset.get<double>("boundary.friction.initial_log_state");
+    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;
+    // }}}
+  
+    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;
+  
+    auto const &dirichletFunction = functions.get("dirichletCondition");
+    auto const &neumannFunction = functions.get("neumannCondition");
+  
+    auto timeSteppingScheme = initTimeStepper
+      (parset.get<Config::scheme>("timeSteppingScheme"), dirichletFunction,
+       ignoreNodes, massMatrix, stiffnessMatrix, u_initial, ud_initial,
+       udd_initial);
+    auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>
+      (parset.get<Config::state_model>("boundary.friction.state_model"),
+       alpha_initial, frictionalNodes, L);
+  
+    auto const state_fpi_max
+      = parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
+    for (size_t run = 1; run <= timesteps; ++run) {
+      stateUpdater->nextTimeStep();
+      timeSteppingScheme->nextTimeStep();
+  
+      double const time = tau * run;
+  
+      VectorType ell(finestSize);
+      assemble_neumann<GridType, GridView, SmallVector, P1Basis>
+        (leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
+      ell += gravityFunctional;
+  
+      MatrixType problem_A;
+      VectorType problem_rhs(finestSize);
+      VectorType problem_iterate(finestSize);
+  
+      stateUpdater->setup(tau);
+      timeSteppingScheme->setup(ell, tau, time,
+                                problem_rhs, problem_iterate, problem_A);
+  
+      VectorType u_saved;
+      for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
+        <<setupAndSolveProblem>>;
+  
+        timeSteppingScheme->postProcess(problem_iterate);
+        timeSteppingScheme->extractDisplacement(u);
+        timeSteppingScheme->extractVelocity(ud);
+  
+        stateUpdater->solve(ud);
+        stateUpdater->extractState(alpha);
+  
+        <<printInnerProgress>>;
+        if (state_fpi > 1 && energyNorm.diff(u_saved, u)
+            < parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
+          break;
+        else
+          u_saved = u;
+  
+        if (state_fpi == state_fpi_max)
+          std::cerr << "[ref = " << refinements
+                    << "]: FPI did not converge after "
+                    << state_fpi_max << " iterations" << std::endl;
+      }
+      <<printOuterProgress>>;
+  
+      <<writeData>>;
+      <<assembleAndWriteStress>>;
+    }
+    <<printBenchmarkData>>;
+    <<closeWriters>>;
+  
+    Python::stop();
+  }
+  catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  }
+  catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+#+end_src
+* Skeleton
 #+begin_src c++ :tangle one-body-sample.cc :noweb yes
   #ifdef HAVE_CONFIG_H
   #include "config.h"
@@ -271,267 +537,13 @@
   
   int const dim = DIM;
   
-  template<class VectorType, class MatrixType, class FunctionType, int dim>
-  Dune::shared_ptr
-  <TimeSteppingScheme
-   <VectorType,MatrixType,FunctionType,dim> >
-  initTimeStepper(Config::scheme scheme,
-                  FunctionType const &dirichletFunction,
-                  Dune::BitSetVector<dim> const &ignoreNodes,
-                  MatrixType const &massMatrix, MatrixType const &stiffnessMatrix,
-                  VectorType const &u_initial,
-                  VectorType const &ud_initial,
-                  VectorType const &udd_initial)
-  {
-    switch (scheme) {
-    case Config::ImplicitTwoStep:
-      return Dune::make_shared<ImplicitTwoStep<VectorType,MatrixType,
-                                               FunctionType,dim> >
-        (stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
-    case Config::ImplicitEuler:
-      return Dune::make_shared<ImplicitEuler<VectorType,MatrixType,
-                                             FunctionType,dim> >
-        (stiffnessMatrix, u_initial, ud_initial, ignoreNodes, dirichletFunction);
-    case Config::Newmark:
-      return Dune::make_shared <Newmark <VectorType,MatrixType,FunctionType,dim> >
-        (stiffnessMatrix, massMatrix, u_initial, ud_initial, udd_initial,
-         ignoreNodes, dirichletFunction);
-    }
-  }
-  
-  template<class SingletonVectorType, class VectorType>
-  Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType> >
-  initStateUpdater(Config::state_model model,
-                   SingletonVectorType const &alpha_initial,
-                   Dune::BitSetVector<1> const &frictionalNodes,
-                   double L)
-  {
-    switch (model) {
-    case Config::Dieterich:
-      return Dune::make_shared
-        <DieterichStateUpdater<SingletonVectorType, VectorType> >
-        (alpha_initial, frictionalNodes, L);
-    case Config::Ruina:
-      return Dune::make_shared
-        <RuinaStateUpdater<SingletonVectorType, VectorType> >
-        (alpha_initial, frictionalNodes, L);
-    }
-  }
-  
-  template <class FunctionMap>
-  void
-  initPython(FunctionMap &functions)
-  {
-    Python::start();
-  
-    Python::run("import sys");
-    Python::run("sys.path.append('" srcdir "')");
-  
-    Python::import("one-body-sample")
-      .get("Functions")
-      .toC<typename FunctionMap::Base>(functions);
-  }
+  <<defineInitTimeStepper>>
+  <<defineInitStateUpdater>>
+  <<defineInitPython>>
   
   int
   main(int argc, char *argv[])
   {
-    try {
-      typedef SharedPointerMap
-        <std::string, Dune::VirtualFunction<double, double> > FunctionMap;
-      FunctionMap functions;
-      initPython(functions);
-  
-      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;
-      typedef Dune::BlockVector<SmallVector> VectorType;
-      typedef Dune::BlockVector<Dune::FieldVector<double,1> > SingletonVectorType;
-  
-      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;
-  
-      // {{{ Set up grid
-      typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming>
-        GridType;
-  
-      Dune::FieldVector<typename GridType::ctype, dim> lowerLeft(0);
-      Dune::FieldVector<typename GridType::ctype, dim> upperRight(1);
-      upperRight[0] = parset.get<size_t>("body.width");
-      upperRight[1] = parset.get<size_t>("body.height");
-  
-      Dune::array<unsigned int, dim> elements;
-      std::fill(elements.begin(), elements.end(), 1);
-      elements[0] = parset.get<size_t>("body.width");
-      elements[1] = parset.get<size_t>("body.height");
-  
-      auto grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid
-        (lowerLeft, upperRight, elements);
-  
-      grid->globalRefine(refinements);
-      size_t const finestSize = grid->size(grid->maxLevel(), dim);
-  
-      typedef GridType::LeafGridView GridView;
-      GridView const leafView = grid->leafView();
-      // }}}
-  
-      // Set up bases
-      typedef P0Basis<GridView, double> P0Basis;
-      typedef P1NodalBasis<GridView, double> P1Basis;
-      P0Basis const p0Basis(leafView);
-      P1Basis const p1Basis(leafView);
-  
-      double normalStress;
-      MatrixType massMatrix;
-      VectorType gravityFunctional;
-      {
-        double const gravity = 9.81;
-        double const density = parset.get<double>("body.density");
-  
-        <<assembleMassMatrix>>;
-        <<computeNormalStress>>;
-        <<computeGravitationalBodyForce>>;
-      }
-  
-      // Assemble elastic force on the body
-      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);
-  
-      // {{{ Initialise vectors
-      VectorType u(finestSize); u = 0.0;
-      VectorType ud(finestSize); ud = 0.0;
-  
-      VectorType u_initial(finestSize); u_initial = 0.0;
-      VectorType ud_initial(finestSize); ud_initial = 0.0;
-      VectorType udd_initial(finestSize); udd_initial = 0.0;
-  
-      SingletonVectorType alpha_initial(finestSize);
-      alpha_initial = parset.get<double>("boundary.friction.initial_log_state");
-      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;
-      // }}}
-  
-      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;
-  
-      auto const &dirichletFunction = functions.get("dirichletCondition");
-      auto const &neumannFunction = functions.get("neumannCondition");
-  
-      auto timeSteppingScheme = initTimeStepper
-        (parset.get<Config::scheme>("timeSteppingScheme"), dirichletFunction,
-         ignoreNodes, massMatrix, stiffnessMatrix, u_initial, ud_initial,
-         udd_initial);
-      auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>
-        (parset.get<Config::state_model>("boundary.friction.state_model"),
-         alpha_initial, frictionalNodes, L);
-  
-      auto const state_fpi_max
-        = parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
-      for (size_t run = 1; run <= timesteps; ++run) {
-        stateUpdater->nextTimeStep();
-        timeSteppingScheme->nextTimeStep();
-  
-        double const time = tau * run;
-  
-        VectorType ell(finestSize);
-        assemble_neumann<GridType, GridView, SmallVector, P1Basis>
-          (leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
-        ell += gravityFunctional;
-  
-        MatrixType problem_A;
-        VectorType problem_rhs(finestSize);
-        VectorType problem_iterate(finestSize);
-  
-        stateUpdater->setup(tau);
-        timeSteppingScheme->setup(ell, tau, time,
-                                  problem_rhs, problem_iterate, problem_A);
-  
-        VectorType u_saved;
-        for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
-          <<setupAndSolveProblem>>;
-  
-          timeSteppingScheme->postProcess(problem_iterate);
-          timeSteppingScheme->extractDisplacement(u);
-          timeSteppingScheme->extractVelocity(ud);
-  
-          stateUpdater->solve(ud);
-          stateUpdater->extractState(alpha);
-  
-          <<printInnerProgress>>;
-          if (state_fpi > 1 && energyNorm.diff(u_saved, u)
-              < parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
-            break;
-          else
-            u_saved = u;
-  
-          if (state_fpi == state_fpi_max)
-            std::cerr << "[ref = " << refinements
-                      << "]: FPI did not converge after "
-                      << state_fpi_max << " iterations" << std::endl;
-        }
-        <<printOuterProgress>>;
-  
-        <<writeData>>;
-        <<assembleAndWriteStress>>;
-      }
-      <<printBenchmarkData>>;
-      <<closeWriters>>;
-  
-      Python::stop();
-    }
-    catch (Dune::Exception &e) {
-      Dune::derr << "Dune reported error: " << e << std::endl;
-    }
-    catch (std::exception &e) {
-      std::cerr << "Standard exception: " << e.what() << std::endl;
-    }
+    <<mainBody>>
   }
-  
 #+end_src