diff --git a/src/one-body-sample.org b/src/one-body-sample.org
index 37a4454097fa73e085c4ec17d0a35d91e1765140..8453610219a3104407a16d75e4804c13b8f74f05 100644
--- a/src/one-body-sample.org
+++ b/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