diff --git a/src/.gitignore b/src/.gitignore
index c52dff94094dc5e1803403a968c0b9d144941306..aa8174205fc7e75281254e0edd1c8b65231201a9 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -1,3 +1,4 @@
timestepping.cc
timestepping.hh
timestepping_tmpl.cc
+one-body-sample.cc
\ No newline at end of file
diff --git a/src/Makefile.am b/src/Makefile.am
index 6c8d03b969ff6b75b5a4649f0809aa6103c102bd..1ea490f5c9a8c62c5edae19fc7467473a36e9575 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -111,3 +111,8 @@ $(srcdir)/timestepping.hh: $(srcdir)/timestepping.org
emacs -Q --batch --eval \
"(let (vc-handled-backends) \
(org-babel-tangle-file \"$<\" nil 'c++))"
+
+$(srcdir)/one-body-sample.cc: $(srcdir)/one-body-sample.org
+ emacs -Q --batch --eval \
+ "(let (vc-handled-backends) \
+ (org-babel-tangle-file \"$<\" nil 'c++))"
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
deleted file mode 100644
index aca2b01ff42281bd1ec60ffbcee5755489e47990..0000000000000000000000000000000000000000
--- a/src/one-body-sample.cc
+++ /dev/null
@@ -1,484 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#ifndef srcdir
-#error srcdir unset
-#endif
-
-#ifndef DIM
-#error DIM unset
-#endif
-
-#ifndef HAVE_PYTHON
-#error Python is required
-#endif
-
-#if !HAVE_ALUGRID
-#error ALUGRID is required
-#endif
-
-#include <cmath>
-#include <exception>
-#include <iostream>
-
-#include <boost/format.hpp>
-
-#include <Python.h>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-#include <dune/common/shared_ptr.hh>
-#include <dune/common/timer.hh>
-#include <dune/grid/alugrid.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/fufem/assemblers/functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
-#include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/dunepython.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/constantfunction.hh>
-#include <dune/fufem/functionspacebases/p0basis.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/sharedpointermap.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh> // Solver::FULL
-
-#include <dune/tectonic/myblockproblem.hh>
-#include <dune/tectonic/myconvexproblem.hh>
-
-#include "assemblers.hh"
-#include "mysolver.hh"
-#include "vtk.hh"
-
-#include "enums.hh"
-#include "enum_parser.cc"
-#include "enum_state_model.cc"
-#include "enum_scheme.cc"
-
-#include "timestepping.hh"
-
-#include "ruinastateupdater.hh"
-#include "dieterichstateupdater.hh"
-
-int const dim = DIM;
-
-template <class GridView, class GridCorner>
-void setup_boundary(GridView const &gridView,
- Dune::BitSetVector<dim> &ignoreNodes,
- Dune::BitSetVector<1> &neumannNodes,
- Dune::BitSetVector<1> &frictionalNodes,
- GridCorner const &lowerLeft, GridCorner const &upperRight,
- size_t &specialNode) {
- typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
-
- Dune::MultipleCodimMultipleGeomTypeMapper<
- GridView, Dune::MCMGVertexLayout> const myVertexMapper(gridView);
-
- for (auto it = gridView.template begin<dim>();
- it != gridView.template end<dim>(); ++it) {
- assert(it->geometry().corners() == 1);
- Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
- size_t const id = myVertexMapper.map(*it);
-
- // Find the center of the lower face
- switch (dim) {
- case 3:
- if (coordinates[2] != lowerLeft[2])
- break;
- // fallthrough
- case 2:
- if (coordinates[1] == lowerLeft[1] &&
- std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
- specialNode = id;
- break;
- default:
- assert(false);
- }
-
- // lower face
- if (coordinates[1] == lowerLeft[1]) {
- frictionalNodes[id] = true;
- ignoreNodes[id][1] = true;
- }
-
- // upper face
- else if (coordinates[1] == upperRight[1])
- ignoreNodes[id] = true;
-
- // right face (except for both corners)
- else if (coordinates[0] == upperRight[0])
- ;
-
- // left face (except for both corners)
- else if (coordinates[0] == lowerLeft[0])
- ;
- }
-}
-
-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);
-}
-
-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");
-
- { // Assemble mass matrix
- MassAssembler<GridType, P1Basis::LocalFiniteElement,
- P1Basis::LocalFiniteElement> const localMass;
- OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
- .assemble(localMass, massMatrix);
- massMatrix *= density;
- }
-
- { // Compute normal stress
- double volume = 1.0;
- for (int i = 0; i < dim; ++i)
- volume *= (upperRight[i] - lowerLeft[i]);
-
- double area = 1.0;
- for (int i = 0; i < dim; ++i)
- if (i != 1)
- area *= (upperRight[i] - lowerLeft[i]);
-
- // volume * gravity * density / area = normal stress
- // V * g * rho / A = sigma_n
- // m^d * N/kg * kg/m^d / m^(d-1) = N/m^(d-1)
- normalStress = volume * gravity * density / area;
- }
-
- { // Compute gravitational body force
- SmallVector weightedGravitationalDirection(0);
- weightedGravitationalDirection[1] = -density * gravity;
- ConstantFunction<SmallVector, SmallVector> const gravityFunction(
- weightedGravitationalDirection);
- L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler(
- gravityFunction);
- FunctionalAssembler<P1Basis>(p1Basis)
- .assemble(gravityFunctionalAssembler, gravityFunctional, true);
- }
- }
-
- // Assemble elastic force on the body
- MatrixType stiffnessMatrix;
- {
- StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
- P1Basis::LocalFiniteElement> const
- localStiffness(E, nu);
- OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
- .assemble(localStiffness, stiffnessMatrix);
- }
- 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);
- setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
- lowerLeft, upperRight, specialNode);
- assert(specialNode != finestSize);
- assert(frictionalNodes[specialNode][0]);
-
- 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.state.initial");
- 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);
-
- std::fstream state_writer("state", std::fstream::out);
- std::fstream displacement_writer("displacement", std::fstream::out);
- std::fstream velocity_writer("velocity", std::fstream::out);
- std::fstream coefficient_writer("coefficient", std::fstream::out);
-
- auto const L = parset.get<double>("boundary.friction.ruina.L");
- auto const a = parset.get<double>("boundary.friction.ruina.a");
- auto const b = parset.get<double>("boundary.friction.ruina.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) {
- double const time = tau * run;
- {
- stateUpdater->nextTimeStep();
-
- VectorType ell(finestSize);
- assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
- leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
- ell += gravityFunctional;
-
- // Copy everything for now
- 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) {
- auto myGlobalNonlinearity =
- assemble_nonlinearity<MatrixType, VectorType>(
- parset.sub("boundary.friction"), *nodalIntegrals, alpha,
- surfaceNormalStress);
-
- MyConvexProblemType const myConvexProblem(
- problem_A, *myGlobalNonlinearity, problem_rhs);
- MyBlockProblemType myBlockProblem(parset, myConvexProblem);
- auto multigridStep = mySolver.getSolver();
- multigridStep->setProblem(problem_iterate, myBlockProblem);
-
- LoopSolver<VectorType> overallSolver(
- multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
- solver_tolerance, &energyNorm, verbosity,
- false); // absolute error
- overallSolver.solve();
-
- timeSteppingScheme->postProcess(problem_iterate);
- timeSteppingScheme->extractDisplacement(u);
- timeSteppingScheme->extractVelocity(ud);
-
- stateUpdater->solve(ud);
- stateUpdater->extractState(alpha);
-
- if (parset.get<bool>("printProgress")) {
- std::cerr << '.';
- std::cerr.flush();
- }
- 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;
- }
- if (parset.get<bool>("printProgress"))
- std::cerr << std::endl;
- }
-
- { // Write all kinds of data
- state_writer << alpha[specialNode][0] << " " << std::endl;
- displacement_writer << u[specialNode][0] << " " << std::endl;
- velocity_writer << ud[specialNode][0] << " " << std::endl;
- coefficient_writer << mu +
- a *std::log(ud[specialNode].two_norm() / V0) +
- b * (alpha[specialNode] + std::log(V0 / L))
- << " " << std::endl;
- }
-
- // Compute von Mises stress and write everything to a file
- if (parset.get<bool>("writeVTK")) {
- VonMisesStressAssembler<GridType> localStressAssembler(
- E, nu,
- Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>(
- p1Basis, u));
- FunctionalAssembler<P0Basis>(p0Basis)
- .assemble(localStressAssembler, vonMisesStress, true);
-
- writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>(
- p1Basis, u, alpha, p0Basis, vonMisesStress, leafView,
- (boost::format("obs%d") % run).str());
- }
- }
- if (parset.get<bool>("enableTimer"))
- std::cerr << std::endl << "Making " << timesteps << " time steps took "
- << timer.elapsed() << "s" << std::endl;
-
- state_writer.close();
- displacement_writer.close();
- velocity_writer.close();
- coefficient_writer.close();
-
- 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;
- }
-}
diff --git a/src/one-body-sample.org b/src/one-body-sample.org
new file mode 100644
index 0000000000000000000000000000000000000000..325a2829ccdc9b06db56b069f533a99b8f4c0001
--- /dev/null
+++ b/src/one-body-sample.org
@@ -0,0 +1,514 @@
+#+name: assembleMassMatrix
+#+begin_src c++
+ {
+ MassAssembler
+ <GridType, P1Basis::LocalFiniteElement,
+ P1Basis::LocalFiniteElement> const localMass;
+ OperatorAssembler<P1Basis,P1Basis>(p1Basis, p1Basis)
+ .assemble(localMass, massMatrix);
+ massMatrix *= density;
+ }
+#+end_src
+
+#+name: computeNormalStress
+#+begin_src c++
+ {
+ double volume = 1.0;
+ for (int i=0; i<dim; ++i)
+ volume *= (upperRight[i] - lowerLeft[i]);
+
+ double area = 1.0;
+ for (int i=0; i<dim; ++i)
+ if (i != 1)
+ area *= (upperRight[i] - lowerLeft[i]);
+
+ // volume * gravity * density / area = normal stress
+ // V * g * rho / A = sigma_n
+ // m^d * N/kg * kg/m^d / m^(d-1) = N/m^(d-1)
+ normalStress = volume * gravity * density / area;
+ }
+#+end_src
+
+#+name: computeGravitationalBodyForce
+#+begin_src c++
+ {
+ SmallVector weightedGravitationalDirection(0);
+ weightedGravitationalDirection[1] = - density * gravity;
+ ConstantFunction<SmallVector,SmallVector> const gravityFunction
+ (weightedGravitationalDirection);
+ L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler
+ (gravityFunction);
+ FunctionalAssembler<P1Basis>(p1Basis)
+ .assemble(gravityFunctionalAssembler, gravityFunctional, true);
+ }
+#+end_src
+
+#+name: assembleStiffnessMatrix
+#+begin_src c++
+ {
+ StVenantKirchhoffAssembler
+ <GridType, P1Basis::LocalFiniteElement,
+ P1Basis::LocalFiniteElement> const localStiffness(E, nu);
+ OperatorAssembler<P1Basis,P1Basis>(p1Basis, p1Basis)
+ .assemble(localStiffness, stiffnessMatrix);
+ }
+#+end_src c++
+
+#+begin_src c++ :tangle one-body-sample.cc :noweb yes
+ #ifdef HAVE_CONFIG_H
+ #include "config.h"
+ #endif
+
+ #ifdef HAVE_IPOPT
+ #undef HAVE_IPOPT
+ #endif
+
+ #ifndef srcdir
+ #error srcdir unset
+ #endif
+
+ #ifndef DIM
+ #error DIM unset
+ #endif
+
+ #ifndef HAVE_PYTHON
+ #error Python is required
+ #endif
+
+ #if ! HAVE_ALUGRID
+ #error ALUGRID is required
+ #endif
+
+ #include <cmath>
+ #include <exception>
+ #include <iostream>
+
+ #include <boost/format.hpp>
+
+ #include <Python.h>
+
+ #include <dune/common/bitsetvector.hh>
+ #include <dune/common/exceptions.hh>
+ #include <dune/common/fmatrix.hh>
+ #include <dune/common/function.hh>
+ #include <dune/common/fvector.hh>
+ #include <dune/common/parametertree.hh>
+ #include <dune/common/parametertreeparser.hh>
+ #include <dune/common/shared_ptr.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/grid/alugrid.hh>
+ #include <dune/grid/common/mcmgmapper.hh>
+ #include <dune/grid/utility/structuredgridfactory.hh>
+ #include <dune/istl/bcrsmatrix.hh>
+ #include <dune/istl/bvector.hh>
+
+ #include <dune/fufem/assemblers/functionalassembler.hh>
+ #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
+ #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+ #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
+ #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
+ #include <dune/fufem/assemblers/operatorassembler.hh>
+ #include <dune/fufem/dunepython.hh>
+ #include <dune/fufem/functions/basisgridfunction.hh>
+ #include <dune/fufem/functions/constantfunction.hh>
+ #include <dune/fufem/functionspacebases/p0basis.hh>
+ #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+ #include <dune/fufem/sharedpointermap.hh>
+ #include <dune/solvers/norms/energynorm.hh>
+ #include <dune/solvers/solvers/loopsolver.hh>
+ #include <dune/solvers/solvers/solver.hh> // Solver::FULL
+
+ #include <dune/tectonic/myblockproblem.hh>
+ #include <dune/tectonic/myconvexproblem.hh>
+
+ #include "assemblers.hh"
+ #include "mysolver.hh"
+ #include "vtk.hh"
+
+ #include "enums.hh"
+ #include "enum_parser.cc"
+ #include "enum_state_model.cc"
+ #include "enum_scheme.cc"
+
+ #include "timestepping.hh"
+
+ #include "ruinastateupdater.hh"
+ #include "dieterichstateupdater.hh"
+
+ int const dim = DIM;
+
+ template <class GridView, class GridCorner>
+ void
+ setup_boundary(GridView const &gridView,
+ Dune::BitSetVector<dim> &ignoreNodes,
+ Dune::BitSetVector<1> &neumannNodes,
+ Dune::BitSetVector<1> &frictionalNodes,
+ GridCorner const &lowerLeft,
+ GridCorner const &upperRight,
+ size_t &specialNode)
+ {
+ typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
+
+ Dune::MultipleCodimMultipleGeomTypeMapper
+ <GridView, Dune::MCMGVertexLayout> const myVertexMapper(gridView);
+
+ for (auto it = gridView.template begin<dim>();
+ it != gridView.template end<dim>();
+ ++it) {
+ assert(it->geometry().corners() == 1);
+ Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
+ size_t const id = myVertexMapper.map(*it);
+
+ // Find the center of the lower face
+ switch (dim) {
+ case 3:
+ if (coordinates[2] != lowerLeft[2])
+ break;
+ // fallthrough
+ case 2:
+ if (coordinates[1] == lowerLeft[1]
+ && std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
+ specialNode = id;
+ break;
+ default:
+ assert(false);
+ }
+
+ // lower face
+ if (coordinates[1] == lowerLeft[1]) {
+ frictionalNodes[id] = true;
+ ignoreNodes[id][1] = true;
+ }
+
+ // upper face
+ else if (coordinates[1] == upperRight[1])
+ ignoreNodes[id] = true;
+
+ // right face (except for both corners)
+ else if (coordinates[0] == upperRight[0])
+ ;
+
+ // left face (except for both corners)
+ else if (coordinates[0] == lowerLeft[0])
+ ;
+ }
+ }
+
+ 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);
+ }
+
+ 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);
+ setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
+ lowerLeft, upperRight, specialNode);
+ assert(specialNode != finestSize);
+ assert(frictionalNodes[specialNode][0]);
+
+ 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.state.initial");
+ 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);
+
+ std::fstream state_writer("state", std::fstream::out);
+ std::fstream displacement_writer("displacement", std::fstream::out);
+ std::fstream velocity_writer("velocity", std::fstream::out);
+ std::fstream coefficient_writer("coefficient", std::fstream::out);
+
+ auto const L = parset.get<double>("boundary.friction.ruina.L");
+ auto const a = parset.get<double>("boundary.friction.ruina.a");
+ auto const b = parset.get<double>("boundary.friction.ruina.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) {
+ double const time = tau * run;
+ {
+ stateUpdater->nextTimeStep();
+
+ VectorType ell(finestSize);
+ assemble_neumann<GridType, GridView, SmallVector, P1Basis>
+ (leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
+ ell += gravityFunctional;
+
+ // Copy everything for now
+ 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) {
+ auto myGlobalNonlinearity
+ = assemble_nonlinearity<MatrixType, VectorType>
+ (parset.sub("boundary.friction"), *nodalIntegrals, alpha,
+ surfaceNormalStress);
+
+ MyConvexProblemType const myConvexProblem
+ (problem_A, *myGlobalNonlinearity, problem_rhs);
+ MyBlockProblemType myBlockProblem(parset, myConvexProblem);
+ auto multigridStep = mySolver.getSolver();
+ multigridStep->setProblem(problem_iterate, myBlockProblem);
+
+ LoopSolver<VectorType> overallSolver
+ (multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+ solver_tolerance, &energyNorm, verbosity, false); // absolute error
+ overallSolver.solve();
+
+ timeSteppingScheme->postProcess(problem_iterate);
+ timeSteppingScheme->extractDisplacement(u);
+ timeSteppingScheme->extractVelocity(ud);
+
+ stateUpdater->solve(ud);
+ stateUpdater->extractState(alpha);
+
+ if (parset.get<bool>("printProgress")) {
+ std::cerr << '.';
+ std::cerr.flush();
+ }
+ 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;
+ }
+ if (parset.get<bool>("printProgress"))
+ std::cerr << std::endl;
+ }
+
+ { // Write all kinds of data
+ state_writer << alpha[specialNode][0] << " " << std::endl;
+ displacement_writer << u[specialNode][0] << " " << std::endl;
+ velocity_writer << ud[specialNode][0] << " " << std::endl;
+ coefficient_writer << mu
+ + a * std::log(ud[specialNode].two_norm() / V0)
+ + b * (alpha[specialNode] + std::log(V0 / L))
+ << " " << std::endl;
+ }
+
+ // Compute von Mises stress and write everything to a file
+ if (parset.get<bool>("writeVTK")) {
+ VonMisesStressAssembler<GridType> localStressAssembler
+ (E, nu,
+ Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>
+ (p1Basis, u));
+ FunctionalAssembler<P0Basis>(p0Basis)
+ .assemble(localStressAssembler, vonMisesStress, true);
+
+ writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>
+ (p1Basis, u, alpha, p0Basis, vonMisesStress, leafView,
+ (boost::format("obs%d") % run).str());
+ }
+ }
+ if (parset.get<bool>("enableTimer"))
+ std::cerr << std::endl
+ << "Making " << timesteps
+ << " time steps took " << timer.elapsed()
+ << "s" << std::endl;
+
+ state_writer.close();
+ displacement_writer.close();
+ velocity_writer.close();
+ coefficient_writer.close();
+
+ 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