src/hdf5/restart-io.hh

-#ifndef SRC_HDF_RESTART_HDF_HH
-#define SRC_HDF_RESTART_HDF_HH
-
-#include <dune/fufem/hdf5/file.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-template <class ProgramState> class RestartIO {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-
-public:
-  RestartIO(HDF5::Grouplike &group, size_t vertexCount);
-
-  void write(ProgramState const &programState);
-
-  void read(size_t timeStep, ProgramState &programState);
-
-private:
-  HDF5::SequenceIO<2> displacementWriter_;
-  HDF5::SequenceIO<2> velocityWriter_;
-  HDF5::SequenceIO<2> accelerationWriter_;
-  HDF5::SequenceIO<1> stateWriter_;
-  HDF5::SequenceIO<1> weightedNormalStressWriter_;
-  HDF5::SequenceIO<0> relativeTimeWriter_;
-  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
-};
-#endif

src/hdf5/restart-io_tmpl.cc

-#include "../explicitvectors.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class RestartIO<MyProgramState>;

src/hdf5/restrict.hh

-#ifndef SRC_HDF5_RESTRICT_HH
-#define SRC_HDF5_RESTRICT_HH
-
-#include <cassert>
-
-#include <dune/common/bitsetvector.hh>
-
-#include "../tobool.hh"
-
-template <class Vector, class Patch>
-Vector restrictToSurface(Vector const &v1, Patch const &patch) {
-  auto const &vertices = *patch.getVertices();
-  assert(vertices.size() == v1.size());
-
-  Vector ret(vertices.count());
-  auto target = ret.begin();
-  for (size_t i = 0; i < v1.size(); ++i)
-    if (toBool(vertices[i]))
-      *(target++) = v1[i];
-  assert(target == ret.end());
-  return ret;
-}
-#endif

src/hdf5/surface-writer.hh

-#ifndef SRC_HDF_SURFACE_WRITER_HH
-#define SRC_HDF_SURFACE_WRITER_HH
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/fufem/hdf5/sequenceio.hh>
-#include <dune/fufem/hdf5/singletonwriter.hh>
-
-template <class ProgramState, class GridView> class SurfaceWriter {
-  using Vector = typename ProgramState::Vector;
-  using Patch = BoundaryPatch<GridView>;
-
-public:
-  SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                Patch const &surface);
-
-  void write(ProgramState const &programState);
-
-private:
-  HDF5::Group group_;
-
-  Patch const &surface_;
-
-  HDF5::SequenceIO<2> surfaceDisplacementWriter_;
-  HDF5::SequenceIO<2> surfaceVelocityWriter_;
-};
-#endif

src/hdf5/time-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "time-writer.hh"
-
-template <class ProgramState>
-TimeWriter<ProgramState>::TimeWriter(HDF5::Grouplike &file)
-    : file_(file),
-      relativeTimeWriter_(file_, "relativeTime"),
-      relativeTimeIncrementWriter_(file_, "relativeTimeIncrement") {}
-
-template <class ProgramState>
-void TimeWriter<ProgramState>::write(ProgramState const &programState) {
-  addEntry(relativeTimeWriter_, programState.timeStep,
-           programState.relativeTime);
-  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
-           programState.relativeTau);
-}
-
-#include "time-writer_tmpl.cc"

src/hdf5/time-writer.hh

-#ifndef SRC_HDF_TIME_WRITER_HH
-#define SRC_HDF_TIME_WRITER_HH
-
-#include <dune/fufem/hdf5/file.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-template <class ProgramState> class TimeWriter {
-public:
-  TimeWriter(HDF5::Grouplike &file);
-  void write(ProgramState const &programState);
-
-private:
-  HDF5::Grouplike &file_;
-  HDF5::SequenceIO<0> relativeTimeWriter_;
-  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
-};
-#endif

src/hdf5/time-writer_tmpl.cc

-#include "../explicitvectors.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class TimeWriter<MyProgramState>;

src/matrices.hh

-#ifndef SRC_MATRICES_HH
-#define SRC_MATRICES_HH
-
-template <class Matrix> struct Matrices {
-  Matrix elasticity;
-  Matrix damping;
-  Matrix mass;
-};
-#endif

src/multi-body-problem/CMakeLists.txt

+add_custom_target(tectonic_src_multi-body-problem SOURCES
+  multi-body-problem.cfg
+  multi-body-problem-2D.cfg
+  multi-body-problem-3D.cfg
+) 
+
+set(MSW_SOURCE_FILES
+  ../../dune/tectonic/assemblers.cc
+  ../../dune/tectonic/data-structures/body/body.cc
+  ../../dune/tectonic/data-structures/network/levelcontactnetwork.cc
+  ../../dune/tectonic/data-structures/network/contactnetwork.cc
+  ../../dune/tectonic/data-structures/enumparser.cc
+  #../../dune/tectonic/factories/cantorfactory.cc
+  ../../dune/tectonic/factories/threeblocksfactory.cc
+  ../../dune/tectonic/factories/stackedblocksfactory.cc
+  #../../dune/tectonic/io/vtk.cc
+  #../../dune/tectonic/io/hdf5/frictionalboundary-writer.cc
+  #../../dune/tectonic/io/hdf5/iteration-writer.cc
+  #../../dune/tectonic/io/hdf5/patchinfo-writer.cc
+  #../../dune/tectonic/io/hdf5/restart-io.cc
+  #../../dune/tectonic/io/hdf5/surface-writer.cc
+  #../../dune/tectonic/io/hdf5/time-writer.cc
+  ../../dune/tectonic/problem-data/grid/cuboidgeometry.cc
+  ../../dune/tectonic/problem-data/grid/mygrids.cc
+  ../../dune/tectonic/problem-data/grid/simplexmanager.cc
+  ../../dune/tectonic/spatial-solving/solverfactory.cc
+  ../../dune/tectonic/spatial-solving/fixedpointiterator.cc
+  ../../dune/tectonic/time-stepping/coupledtimestepper.cc
+  ../../dune/tectonic/time-stepping/adaptivetimestepper.cc
+  ../../dune/tectonic/time-stepping/rate.cc
+  ../../dune/tectonic/time-stepping/rate/rateupdater.cc
+  ../../dune/tectonic/time-stepping/state.cc
+  multi-body-problem.cc
+)
+
+foreach(_dim 2 3)
+  set(_msw_target multi-body-problem-${_dim}D)
+
+  add_executable(${_msw_target} ${MSW_SOURCE_FILES})
+  
+  add_dune_ug_flags(${_msw_target})
+  add_dune_hdf5_flags(${_msw_target})
+
+  set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
+endforeach()

src/one-body-problem-2D.cfg → src/multi-body-problem/multi-body-problem-2D.cfg

 # -*- mode:conf -*-
 [boundary.friction]
-smallestDiameter= 2e-3  # [m]
+smallestDiameter = 0.05  # 0.05 2e-3 [m]
 
 [timeSteps]
-refinementTolerance = 1e-5
+refinementTolerance = 2e-2 # 1e-5
 
 [u0.solver]
 tolerance         = 1e-8
@@ -17,5 +17,8 @@ tolerance         = 1e-8
 [v.fpi]
 tolerance         = 1e-5
 
-[solver.tnnmg.linear]
+[solver.tnnmg.preconditioner.basesolver]
+tolerance          = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
 tolerance          = 1e-10

src/multi-body-problem/multi-body-problem.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#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/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/matrices.hh>
+#include <dune/tectonic/data-structures/program_state.hh>
+#include <dune/tectonic/data-structures/friction/globalfriction.hh>
+
+#include <dune/tectonic/factories/stackedblocksfactory.hh>
+#include <dune/tectonic/factories/threeblocksfactory.hh>
+
+#include <dune/tectonic/io/io-handler.hh>
+#include <dune/tectonic/io/hdf5-writer.hh>
+#include <dune/tectonic/io/hdf5/restart-io.hh>
+#include <dune/tectonic/io/vtk.hh>
+
+#include <dune/tectonic/problem-data/bc.hh>
+#include <dune/tectonic/problem-data/mybody.hh>
+#include <dune/tectonic/problem-data/grid/mygrids.hh>
+
+#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
+//#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
+#include <dune/tectonic/spatial-solving/solverfactory.hh>
+
+#include <dune/tectonic/time-stepping/adaptivetimestepper.hh>
+#include <dune/tectonic/time-stepping/rate.hh>
+#include <dune/tectonic/time-stepping/state.hh>
+#include <dune/tectonic/time-stepping/stepbase.hh>
+#include <dune/tectonic/time-stepping/updaters.hh>
+
+#include <dune/tectonic/utils/debugutils.hh>
+#include <dune/tectonic/utils/diameter.hh>
+#include <dune/tectonic/utils/geocoordinate.hh>
+
+// for getcwd
+#include <unistd.h>
+
+//#include <tbb/tbb.h> //TODO multi threading preconditioner?
+//#include <pthread.h>
+
+size_t const dims = MY_DIM;
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+std::vector<std::vector<double>> allReductionFactors;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree("/home/joscha/software/dune/dune-tectonic/src/multi-body-problem/multi-body-problem.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString("/home/joscha/software/dune/dune-tectonic/src/multi-body-problem/multi-body-problem-%dD.cfg", dims), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+  return parset;
+}
+
+static std::atomic<bool> terminationRequested(false);
+void handleSignal(int signum) { terminationRequested = true; }
+
+int main(int argc, char *argv[]) {
+  try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
+    std::ofstream out("multi-body-problem.log");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+    auto const parset = getParameters(argc, argv);
+
+    using Assembler = MyAssembler<DefLeafGridView, dims>;
+    using field_type = Matrix::field_type;
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+
+    using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
+    //using BlocksFactory = ThreeBlocksFactory<Grid, Vector>;
+    BlocksFactory blocksFactory(parset);
+
+    using ContactNetwork = typename BlocksFactory::ContactNetwork;
+    blocksFactory.build();
+
+    ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
+
+    const size_t bodyCount = contactNetwork.nBodies();
+
+    for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+        //printDofLocation(contactNetwork.body(i)->gridView());
+
+        //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+        //def = 1;
+        //deformedGridComplex.setDeformation(def, i);
+
+        const auto& level = *contactNetwork.level(i);
+
+        for (size_t j=0; j<level.nBodies(); j++) {
+            writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+        }
+    }
+
+    /*for (size_t i=0; i<bodyCount; i++) {
+        printDofLocation(contactNetwork.body(i)->gridView());
+        writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+    } */
+
+    // ----------------------------
+    // assemble contactNetwork
+    // ----------------------------
+    contactNetwork.assemble();
+
+    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+    // -----------------
+    // init input/output
+    // -----------------
+    std::vector<size_t> nVertices(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+        nVertices[i] = contactNetwork.body(i)->nVertices();
+    }
+
+    using MyProgramState = ProgramState<Vector, ScalarVector>;
+    MyProgramState programState(nVertices);
+
+    IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
+
+    bool restartRead = ioHandler.read(programState);
+    if (!restartRead) {
+      programState.setupInitialConditions(parset, contactNetwork);
+    }
+
+    //print(programState.u, "u:");
+
+    auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+    for (size_t i=0; i<bodyCount; i++) {
+      contactNetwork.body(i)->setDeformation(programState.u[i]);
+    }
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // ------------------------
+    // assemble global friction
+    // ------------------------
+    contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+
+    auto& globalFriction = contactNetwork.globalFriction();
+    globalFriction.updateAlpha(programState.alpha);
+
+
+    IterationRegister iterationCount;
+
+    ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, true);
+
+    //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
+    // -------------------
+    // Set up TNNMG solver
+    // -------------------
+
+    BitVector totalDirichletNodes;
+    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+    for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+        bool val = false;
+        for (size_t d=0; d<dims; d++) {
+            val = val || totalDirichletNodes[i][d];
+        }
+
+        totalDirichletNodes[i] = val;
+        for (size_t d=0; d<dims; d++) {
+            totalDirichletNodes[i][d] = val;
+        }
+    }
+
+    //print(totalDirichletNodes, "totalDirichletNodes:");
+
+    //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
+    using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+    using NonlinearFactory = SolverFactory<Functional, BitVector>;
+
+    using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+                               StateUpdater<ScalarVector, Vector>>;
+
+    BoundaryFunctions velocityDirichletFunctions;
+    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+        for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+        print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+        }
+    }*/
+
+    std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+    contactNetwork.frictionNodes(frictionNodes);
+
+    /*for (size_t i=0; i<frictionNodes.size(); i++) {
+        print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+    }*/
+
+    Updaters current(
+        initRateUpdater(
+            parset.get<Config::scheme>("timeSteps.scheme"),
+            velocityDirichletFunctions,
+            dirichletNodes,
+            contactNetwork.matrices(),
+            programState.u,
+            programState.v,
+            programState.a),
+        initStateUpdater<ScalarVector, Vector>(
+            parset.get<Config::stateModel>("boundary.friction.stateModel"),
+            programState.alpha,
+            nBodyAssembler.getContactCouplings(),
+            contactNetwork.couplings())
+            );
+
+
+    auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
+
+    const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+
+    auto const mustRefine = [&](Updaters &coarseUpdater,
+                                Updaters &fineUpdater) {
+
+        //return false;
+      //std::cout << "Step 1" << std::endl;
+
+      std::vector<ScalarVector> coarseAlpha;
+      coarseAlpha.resize(bodyCount);
+      coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+      //print(coarseAlpha, "coarseAlpha:");
+
+      std::vector<ScalarVector> fineAlpha;
+      fineAlpha.resize(bodyCount);
+      fineUpdater.state_->extractAlpha(fineAlpha);
+
+      //print(fineAlpha, "fineAlpha:");
+
+      //std::cout << "Step 3" << std::endl;
+
+      ScalarVector::field_type energyNorm = 0;
+      for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+          //std::cout << "for " << i << std::endl;
+
+          //std::cout << not stateEnergyNorms[i] << std::endl;
+
+          if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
+              continue;
+
+          energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
+      }
+      std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance <<  std::endl;
+      std::cout << "must refine: " << (energyNorm > refinementTolerance) <<  std::endl;
+      return energyNorm > refinementTolerance;
+    };
+
+
+    std::signal(SIGXCPU, handleSignal);
+    std::signal(SIGINT, handleSignal);
+    std::signal(SIGTERM, handleSignal);
+
+/*
+    // set patch preconditioner for linear correction in TNNMG method
+    using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
+
+    const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+
+    auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+    PatchSolver patchSolver(gsStep,
+                               preconditionerParset.get<size_t>("maximumIterations"),
+                               preconditionerParset.get<double>("tolerance"),
+                               nullptr,
+                               preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
+                               false); // absolute error
+
+    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
+    preconditioner.setPatchSolver(patchSolver);
+    preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+*/
+    // set adaptive time stepper
+    typename ContactNetwork::ExternalForces externalForces;
+    contactNetwork.externalForces(externalForces);
+
+    StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
+                 externalForces, stateEnergyNorms);
+
+    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        adaptiveTimeStepper(stepBase, contactNetwork, current,
+                            programState.relativeTime, programState.relativeTau,
+                            mustRefine);
+
+    size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+
+    while (!adaptiveTimeStepper.reachedEnd()) {
+      programState.timeStep++;
+
+      //preconditioner.build();
+      iterationCount = adaptiveTimeStepper.advance();
+
+      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
+      programState.relativeTau = adaptiveTimeStepper.relativeTau_;
+      current.rate_->extractDisplacement(programState.u);
+      current.rate_->extractVelocity(programState.v);
+      current.rate_->extractAcceleration(programState.a);
+      current.state_->extractAlpha(programState.alpha);
+      globalFriction.updateAlpha(programState.alpha);
+
+      /*print(programState.u, "current u:");
+      print(programState.v, "current v:");
+      print(programState.a, "current a:");
+      print(programState.alpha, "current alpha:");*/
+
+      contactNetwork.setDeformation(programState.u);
+
+      ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, false);
+
+      if (programState.timeStep==timeSteps) {
+        std::cout << "limit of timeSteps reached!" << std::endl;
+        break; // TODO remove after debugging
+      }
+
+      if (terminationRequested) {
+        std::cerr << "Terminating prematurely" << std::endl;
+        break;
+      }
+
+
+    }
+
+    // output reduction factors
+    std::vector<double> avgReductionFactors(allReductionFactors.size());
+    //std::cout << "allReductionFactors.size: " << allReductionFactors.size() << std::endl;
+
+    for (size_t i=0; i<avgReductionFactors.size(); i++) {
+        avgReductionFactors[i] = 1;
+        size_t c = 0;
+
+        //std::cout << "allReductionFactors[i].size: " << allReductionFactors[i].size() << std::endl;
+        //print(allReductionFactors[i], "all reduction factors: ");
+
+        for (size_t j=0; j<allReductionFactors[i].size(); j++) {
+            avgReductionFactors[i] *= allReductionFactors[i][j];
+            c++;
+        }
+
+        avgReductionFactors[i] = std::pow(avgReductionFactors[i], 1.0/((double)c));
+    }
+
+    print(avgReductionFactors, "average reduction factors: ");
+
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

src/one-body-problem.cfg → src/multi-body-problem/multi-body-problem.cfg

 # -*- mode:conf -*-
 gravity         = 9.81  # [m/s^2]
 
-[io]
-data.write      = true
-printProgress   = false
-restarts.first  = 0
-restarts.spacing= 20
-restarts.write  = true
-vtk.write       = false
-
-[problem]
-finalTime       = 1000  # [s]
 
 [body]
-bulkModulus     = 0.5e5 # [Pa]
-poissonRatio    = 0.3   # [1]
+bulkModulus     = 1.5e5 # [Pa]
+poissonRatio    = 0.11   # [1]
 [body.elastic]
-density         = 900   # [kg/m^3]
-shearViscosity  = 1e3   # [Pas]
-bulkViscosity   = 1e3   # [Pas]
+density         = 1300   # [kg/m^3]
+shearViscosity  = 0   # [Pas]
+bulkViscosity   = 0   # [Pas]
 [body.viscoelastic]
-density         = 1000  # [kg/m^3]
+density         = 1300  # [kg/m^3]
 shearViscosity  = 1e4   # [Pas]
 bulkViscosity   = 1e4   # [Pas]
 
+
 [boundary.friction]
-C               = 10    # [Pa]
-mu0             = 0.7   # [ ]
-V0              = 5e-5  # [m/s]
-L               = 2.25e-5 # [m]
-initialAlpha    = 0     # [ ]
+C               = 6       # [Pa]
+mu0             = 0.48      # [ ]
+V0              = 1e-3     # [m/s]
+L               = 1e-6  # [m]
+initialAlpha    = 0        # [ ]
 stateModel      = AgeingLaw
-frictionModel   = Regularised
+frictionModel   = Truncated #Regularised
 [boundary.friction.weakening]
-a               = 0.002 # [ ]
-b               = 0.017 # [ ]
+a               = 0.054    # [ ]
+b               = 0.074    # [ ]
 [boundary.friction.strengthening]
-a               = 0.020 # [ ]
-b               = 0.005 # [ ]
+a               = 0.054    # [ ]
+b               = 0.074    # [ ]
+
+
+[boundary.neumann]
+sigmaN          =  0.0     # 200.0 [Pa]
+
+
+[boundary.dirichlet]
+finalVelocity   = 1e-4      # [m/s]
+
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 2e-4 # 1e-6
+
 
 [timeSteps]
 scheme = newmark
+timeSteps = 5
+
+
+[problem]
+finalTime       = 100  # [s] #1000
+bodyCount       = 4
+
+
+[io]
+data.write      = false
+printProgress   = true
+restarts.first  = 0
+restarts.spacing= 1 #20
+restarts.write  = true #true
+vtk.write       = true
+
 
 [u0.solver]
-maximumIterations = 100000
-verbosity         = quiet
+maximumIterations = 100
+verbosity         = full
 
 [a0.solver]
-maximumIterations = 100000
-verbosity         = quiet
+maximumIterations = 100
+verbosity         = full
 
 [v.solver]
-maximumIterations = 100000
+maximumIterations = 100
 verbosity         = quiet
 
 [v.fpi]
 maximumIterations = 10000
 lambda            = 0.5
 
-[solver.tnnmg.linear]
-maxiumumIterations = 100000
-pre                = 3
-cycle              = 1  # 1 = V, 2 = W, etc.
-post               = 3
+[solver.tnnmg.preconditioner]
+mode         = additive
+patchDepth   = 1
+maximumIterations = 2
+verbosity         = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity         = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity         = quiet
 
 [solver.tnnmg.main]
 pre   = 1
 multi = 5 # number of multigrid steps
 post  = 0
+
+
+

src/one-body-problem-data/bc.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
-
-class VelocityDirichletCondition
-    : public Dune::VirtualFunction<double, double> {
-  void evaluate(double const &relativeTime, double &y) const {
-    // Assumed to vanish at time zero
-    double const finalVelocity = -5e-5;
-    y = (relativeTime <= 0.1)
-            ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
-            : finalVelocity;
-  }
-};
-
-class NeumannCondition : public Dune::VirtualFunction<double, double> {
-  void evaluate(double const &relativeTime, double &y) const { y = 0.0; }
-};
-#endif

src/one-body-problem-data/geometry.tex

-\documentclass[tikz]{minimal}
-
-\usepackage{tikz}
-\usetikzlibrary{calc}
-\usetikzlibrary{decorations.pathreplacing}
-
-\usepackage{siunitx}
-
-\begin{document}
-\pgfmathsetmacro{\rightleg}{0.27}
-\pgfmathsetmacro{\leftleg}{1.00}
-\pgfmathsetmacro{\leftangle}{atan(\rightleg/\leftleg)}
-\begin{tikzpicture}[scale=12, rotate=\leftangle]
-  \pgfmathsetmacro{\mysin}{sin(\leftangle)}
-  \pgfmathsetmacro{\mycos}{cos(\leftangle)}
-  \pgfmathsetmacro{\viscoheight}{0.06}
-  \pgfmathsetmacro{\Zx}{0.35}
-  \pgfmathsetmacro{\weaklen}{0.20}
-
-  \coordinate (A) at (0,0);
-  \node at (A) [left] {A};
-  \coordinate (B) at (\leftleg,-\rightleg);
-  \node at (B) [right] {B};
-  \coordinate (C) at (\leftleg,0);
-  \node at (C) [right] {C};
-
-  \draw (A) -- (B) -- (C) -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
-
-  \coordinate (Z) at (\Zx,0);
-  \node at (Z) [above] {Z};
-  \coordinate (Y) at ($(\Zx,-\Zx/\leftleg * \rightleg)$);
-  \node at (Y) [below] {Y};
-  \coordinate (X) at ($(Y) + (-\weaklen*\mycos,\weaklen*\mysin)$);
-  \node at (X) [below] {X};
-  \path let \p1 = (X) in coordinate (U) at ($(\x1, 0)$);
-  \node at (U) [above] {U};
-
-  \path (A) -- node [above=.25cm, sloped] {$\overline{AZ} = \SI{35}{cm}$} (Z);
-
-  \draw[color=red, thick] (X) -- node [below=.25cm] {$\overline{XY}=\SI{20}{cm}$} (Y);
-  \draw[dashed] (Y) -- (Z);
-  \draw[dashed] (U) -- (X);
-
-  \coordinate (K) at ($(B) + (-\leftleg * \viscoheight / \rightleg,\viscoheight)$);
-  \node at (K) [below] {K};
-  \coordinate (M) at ($(B) + (0, \viscoheight)$);
-  \node at (M) [right] {M};
-  \path (C) -- node [right=.5cm] {$\overline{CM} = \SI{21}{cm}$} (M);
-
-  \path[fill=blue] (K) -- (B) -- node [right=.75cm] {$\overline{BM}=\SI{6}{cm}$} (M) -- cycle;
-
-  \coordinate (G) at ($(A) ! 0.5 ! (X)$);
-  \node at (G) [below] {G};
-  \coordinate (H) at ($(X) ! 0.5 ! (Y)$);
-  \node at (H) [below] {H};
-  \coordinate (J) at ($(Y) ! 0.5 ! (B)$);
-  \node at (J) [below] {J};
-
-  \coordinate (I) at ($(Y) + (G)$);
-  \node at (I) [below] {I};
-
-  \node[align=left] at (0.5,-0.225) {
-    $Z$: coast line\\
-    $\overline{XY}$: velocity weakening zone\\
-    $BKM$: visco-elastic domain};
-\end{tikzpicture}
-
-\end{document}

src/one-body-problem-data/mygeometry.hh

-#ifndef SRC_MYGEOMETRY_HH
-#define SRC_MYGEOMETRY_HH
-
-#include <dune/common/fvector.hh>
-
-#include "midpoint.hh"
-
-namespace MyGeometry {
-namespace {
-  using LocalVector2D = Dune::FieldVector<double, 2>;
-  using LocalMatrix2D = Dune::FieldMatrix<double, 2, 2>;
-
-  using LocalVector = Dune::FieldVector<double, MY_DIM>;
-}
-
-namespace reference {
-  double const s = 1.0; // scaling factor
-
-  double const rightLeg = 0.27 * s;
-  double const leftLeg = 1.00 * s;
-  double const leftAngle = atan(rightLeg / leftLeg);
-  double const viscoHeight = 0.06 * s; // Height of the viscous bottom layer
-  double const weakLen = 0.20 * s;     // Length of the weak zone
-
-  double const zDistance = 0.35;
-
-  LocalVector2D const A = {0, 0};
-  LocalVector2D const B = {leftLeg, -rightLeg};
-  LocalVector2D const C = {leftLeg, 0};
-
-  LocalVector2D const Z = {zDistance * s, 0};
-  LocalVector2D const Y = {zDistance * s, -zDistance *s / leftLeg *rightLeg};
-  LocalVector2D const X = {Y[0] - weakLen * std::cos(leftAngle),
-                           Y[1] + weakLen *std::sin(leftAngle)};
-
-  LocalVector2D const U = {X[0], 0};
-
-  LocalVector2D const K = {B[0] - leftLeg * viscoHeight / rightLeg,
-                           B[1] + viscoHeight};
-  LocalVector2D const M = {B[0], B[1] + viscoHeight};
-
-  LocalVector2D const G = midPoint(A, X);
-  LocalVector2D const H = midPoint(X, Y);
-  LocalVector2D const J = midPoint(Y, B);
-
-  LocalVector2D const I = {Y[0] + G[0], Y[1] + G[1]};
-
-  LocalVector2D const zenith = {0, 1};
-
-  LocalMatrix2D const rotation = {{std::cos(leftAngle), -std::sin(leftAngle)},
-                                  {std::sin(leftAngle), std::cos(leftAngle)}};
-}
-
-namespace {
-  LocalVector rotate(LocalVector2D const &x) {
-    LocalVector2D ret2D;
-    reference::rotation.mv(x, ret2D);
-    LocalVector ret(0);
-    ret[0] = ret2D[0];
-    ret[1] = ret2D[1];
-    return ret;
-  }
-}
-
-double const lengthScale = reference::s;
-
-double const depth = 0.60 * lengthScale;
-
-LocalVector const A = rotate(reference::A);
-LocalVector const B = rotate(reference::B);
-LocalVector const C = rotate(reference::C);
-LocalVector const G = rotate(reference::G);
-LocalVector const H = rotate(reference::H);
-LocalVector const I = rotate(reference::I);
-LocalVector const J = rotate(reference::J);
-LocalVector const K = rotate(reference::K);
-LocalVector const M = rotate(reference::M);
-LocalVector const U = rotate(reference::U);
-LocalVector const X = rotate(reference::X);
-LocalVector const Y = rotate(reference::Y);
-LocalVector const Z = rotate(reference::Z);
-
-LocalVector const zenith = rotate(reference::zenith);
-
-void write();
-
-void render();
-}
-#endif

src/one-body-problem-data/mygrid.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/fufem/geometry/polyhedrondistance.hh>
-
-#include "mygrid.hh"
-#include "midpoint.hh"
-#include "../diameter.hh"
-
-#if MY_DIM == 3
-SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
-#endif
-
-// back-to-front, front-to-back, front-to-back
-void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
-                                        unsigned int W) {
-#if MY_DIM == 3
-  unsigned int const U2 = U + shift_;
-  unsigned int const V2 = V + shift_;
-  unsigned int const W2 = W + shift_;
-
-  simplices_.push_back({ U, V, W, U2 });
-  simplices_.push_back({ V, V2, W2, U2 });
-  simplices_.push_back({ W, W2, U2, V });
-#else
-  simplices_.push_back({ U, V, W });
-#endif
-}
-
-// back-to-front, back-to-front, front-to-back
-void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
-                                        unsigned int W) {
-#if MY_DIM == 3
-  unsigned int const U2 = U + shift_;
-  unsigned int const V2 = V + shift_;
-  unsigned int const W2 = W + shift_;
-
-  simplices_.push_back({ U, V, W, U2 });
-  simplices_.push_back({ V, V2, W, U2 });
-  simplices_.push_back({ V2, W, U2, W2 });
-#else
-  simplices_.push_back({ U, V, W });
-#endif
-}
-
-auto SimplexManager::getSimplices() -> SimplexList const & {
-  return simplices_;
-}
-
-template <class Grid> GridConstructor<Grid>::GridConstructor() {
-  auto const &A = MyGeometry::A;
-  auto const &B = MyGeometry::B;
-  auto const &C = MyGeometry::C;
-
-  unsigned int const vc = 3;
-
-#if MY_DIM == 3
-  Dune::FieldMatrix<double, 2 * vc, MY_DIM> vertices;
-#else
-  Dune::FieldMatrix<double, vc, MY_DIM> vertices;
-#endif
-  for (size_t i = 0; i < 2; ++i) {
-#if MY_DIM == 3
-    size_t numXYplanes = 2;
-#else
-    size_t numXYplanes = 1;
-#endif
-    size_t k = 0;
-    for (size_t j = 1; j <= numXYplanes; ++j) {
-      vertices[k++][i] = A[i];
-      vertices[k++][i] = B[i];
-      vertices[k++][i] = C[i];
-      assert(k == j * vc);
-    }
-  }
-
-#if MY_DIM == 3
-  for (size_t k = 0; k < vc; ++k) {
-    vertices[k][2] = -MyGeometry::depth / 2.0;
-    vertices[k + vc][2] = MyGeometry::depth / 2.0;
-  }
-#endif
-
-  for (size_t i = 0; i < vertices.N(); ++i)
-    gridFactory.insertVertex(vertices[i]);
-
-  Dune::GeometryType cell;
-#if MY_DIM == 3
-  cell.makeTetrahedron();
-#else
-  cell.makeTriangle();
-#endif
-
-#if MY_DIM == 3
-  SimplexManager sm(vc);
-#else
-  SimplexManager sm;
-#endif
-  sm.addFromVerticesFFB(1, 2, 0);
-  auto const &simplices = sm.getSimplices();
-
-  // sanity-check choices of simplices
-  for (size_t i = 0; i < simplices.size(); ++i) {
-    Dune::FieldMatrix<double, MY_DIM, MY_DIM> check;
-    for (size_t j = 0; j < MY_DIM; ++j)
-      check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
-    assert(check.determinant() > 0);
-    gridFactory.insertElement(cell, simplices[i]);
-  }
-}
-
-template <class Grid> std::shared_ptr<Grid> GridConstructor<Grid>::getGrid() {
-  return std::shared_ptr<Grid>(gridFactory.createGrid());
-}
-
-template <class Grid>
-template <class GridView>
-MyFaces<GridView> GridConstructor<Grid>::constructFaces(
-    GridView const &gridView) {
-  return MyFaces<GridView>(gridView);
-}
-
-template <class GridView>
-template <class Vector>
-bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
-                                    Vector const &c) {
-  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
-}
-
-template <class GridView>
-template <class Vector>
-bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
-                                Vector const &x) {
-  auto const minmax0 = std::minmax(v1[0], v2[0]);
-  auto const minmax1 = std::minmax(v1[1], v2[1]);
-
-  if (minmax0.first - 1e-14 * MyGeometry::lengthScale > x[0] or
-      x[0] > minmax0.second + 1e-14 * MyGeometry::lengthScale)
-    return false;
-  if (minmax1.first - 1e-14 * MyGeometry::lengthScale > x[1] or
-      x[1] > minmax1.second + 1e-14 * MyGeometry::lengthScale)
-    return false;
-
-  return true;
-}
-
-template <class GridView>
-template <class Vector>
-bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
-                                  Vector const &x) {
-  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
-}
-
-template <class GridView>
-MyFaces<GridView>::MyFaces(GridView const &gridView)
-    :
-#if MY_DIM == 3
-      lower(gridView),
-      right(gridView),
-      upper(gridView),
-      front(gridView),
-      back(gridView)
-#else
-      lower(gridView),
-      right(gridView),
-      upper(gridView)
-#endif
-{
-  assert(isClose(MyGeometry::A[1], 0));
-  assert(isClose(MyGeometry::B[1], 0));
-  lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-    return isClose(0, in.geometry().center()[1]);
-  });
-#if MY_DIM == 3
-  front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-    return isClose(MyGeometry::depth / 2.0, in.geometry().center()[2]);
-  });
-
-  back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-    return isClose(-MyGeometry::depth / 2.0, in.geometry().center()[2]);
-  });
-#endif
-
-  upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-    return xyBetween(MyGeometry::A, MyGeometry::C, in.geometry().center());
-  });
-
-  right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-    return xyBetween(MyGeometry::B, MyGeometry::C, in.geometry().center());
-  });
-}
-
-double computeAdmissibleDiameter(double distance, double smallestDiameter) {
-  return (distance / 0.0125 / MyGeometry::lengthScale + 1.0) * smallestDiameter;
-}
-
-template <class Grid, class LocalVector>
-void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
-            double smallestDiameter) {
-  bool needRefine = true;
-  while (true) {
-    needRefine = false;
-    for (auto &&e : elements(grid.leafGridView())) {
-      auto const geometry = e.geometry();
-
-      auto const weakeningRegionDistance =
-          distance(weakPatch, geometry, 1e-6 * MyGeometry::lengthScale);
-      auto const admissibleDiameter =
-          computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter);
-
-      if (diameter(geometry) <= admissibleDiameter)
-        continue;
-
-      needRefine = true;
-      grid.mark(1, e);
-    }
-    if (!needRefine)
-      break;
-
-    grid.preAdapt();
-    grid.adapt();
-    grid.postAdapt();
-  }
-}
-
-#include "mygrid_tmpl.cc"

src/one-body-problem-data/mygrid_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-template class GridConstructor<Grid>;
-
-template struct MyFaces<GridView>;
-
-template MyFaces<GridView> GridConstructor<Grid>::constructFaces(
-    GridView const &gridView);
-
-template void refine<Grid, LocalVector>(
-    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
-    double smallestDiameter);

src/one-body-problem-data/weakpatch.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
-
-template <class LocalVector>
-ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset) {
-  ConvexPolyhedron<LocalVector> weakPatch;
-#if MY_DIM == 3
-  weakPatch.vertices.resize(4);
-  weakPatch.vertices[0] = weakPatch.vertices[2] = MyGeometry::X;
-  weakPatch.vertices[1] = weakPatch.vertices[3] = MyGeometry::Y;
-  for (size_t k = 0; k < 2; ++k) {
-    weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0;
-    weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0;
-  }
-  switch (parset.get<Config::PatchType>("patchType")) {
-    case Config::Rectangular:
-      break;
-    case Config::Trapezoidal:
-      weakPatch.vertices[1][0] += 0.05 * MyGeometry::lengthScale;
-      weakPatch.vertices[3][0] -= 0.05 * MyGeometry::lengthScale;
-      break;
-    default:
-      assert(false);
-  }
-#else
-  weakPatch.vertices.resize(2);
-  weakPatch.vertices[0] = MyGeometry::X;
-  weakPatch.vertices[1] = MyGeometry::Y;
-#endif
-  return weakPatch;
-};
-#endif

src/one-body-problem.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#include <atomic>
-#include <cmath>
-#include <csignal>
-#include <exception>
-#include <fstream>
-#include <iostream>
-#include <iomanip>
-
-#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/parallel/mpihelper.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/formatstring.hh>
-
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/geocoordinate.hh>
-#include <dune/tectonic/myblockproblem.hh>
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/fufem/hdf5/file.hh>
-
-#include "assemblers.hh"
-#include "diameter.hh"
-#include "enumparser.hh"
-#include "enums.hh"
-#include "gridselector.hh"
-#include "hdf5-writer.hh"
-#include "hdf5/restart-io.hh"
-#include "matrices.hh"
-#include "program_state.hh"
-#include "one-body-problem-data/bc.hh"
-#include "one-body-problem-data/mybody.hh"
-#include "one-body-problem-data/mygeometry.hh"
-#include "one-body-problem-data/myglobalfrictiondata.hh"
-#include "one-body-problem-data/mygrid.hh"
-#include "one-body-problem-data/weakpatch.hh"
-#include "spatial-solving/solverfactory.hh"
-#include "time-stepping/adaptivetimestepper.hh"
-#include "time-stepping/rate.hh"
-#include "time-stepping/state.hh"
-#include "time-stepping/updaters.hh"
-#include "vtk.hh"
-
-size_t const dims = MY_DIM;
-
-Dune::ParameterTree getParameters(int argc, char *argv[]) {
-  Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
-  Dune::ParameterTreeParser::readINITree(
-      Dune::Fufem::formatString("one-body-problem-%dD.cfg", dims), parset);
-  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
-  return parset;
-}
-
-static std::atomic<bool> terminationRequested(false);
-void handleSignal(int signum) { terminationRequested = true; }
-
-int main(int argc, char *argv[]) {
-  try {
-    Dune::MPIHelper::instance(argc, argv);
-    auto const parset = getParameters(argc, argv);
-
-    MyGeometry::render();
-    MyGeometry::write();
-
-    using GridView = Grid::LeafGridView;
-    using MyAssembler = MyAssembler<GridView, dims>;
-    using Matrix = MyAssembler::Matrix;
-    using Vector = MyAssembler::Vector;
-    using LocalVector = Vector::block_type;
-    using ScalarMatrix = MyAssembler::ScalarMatrix;
-    using ScalarVector = MyAssembler::ScalarVector;
-
-    auto const weakPatch =
-        getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"));
-
-    // {{{ Set up grid
-    GridConstructor<Grid> gridConstructor;
-    auto grid = gridConstructor.getGrid();
-
-    refine(*grid, weakPatch,
-           parset.get<double>("boundary.friction.smallestDiameter"));
-
-    double minDiameter = std::numeric_limits<double>::infinity();
-    double maxDiameter = 0.0;
-    for (auto &&e : elements(grid->leafGridView())) {
-      auto const geometry = e.geometry();
-      auto const diam = diameter(geometry);
-      minDiameter = std::min(minDiameter, diam);
-      maxDiameter = std::max(maxDiameter, diam);
-    }
-    std::cout << "min diameter: " << minDiameter << std::endl;
-    std::cout << "max diameter: " << maxDiameter << std::endl;
-
-    auto const leafView = grid->leafGridView();
-    auto const leafVertexCount = leafView.size(dims);
-
-    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
-
-    auto myFaces = gridConstructor.constructFaces(leafView);
-
-    BoundaryPatch<GridView> const neumannBoundary(leafView);
-    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    BoundaryPatch<GridView> const &surface = myFaces.upper;
-
-    // Dirichlet Boundary
-    Dune::BitSetVector<dims> noNodes(leafVertexCount);
-    Dune::BitSetVector<dims> dirichletNodes(leafVertexCount);
-    for (size_t i = 0; i < leafVertexCount; ++i) {
-      if (myFaces.right.containsVertex(i))
-        dirichletNodes[i][0] = true;
-
-      if (myFaces.lower.containsVertex(i))
-        dirichletNodes[i][1] = true;
-
-#if MY_DIM == 3
-      if (myFaces.front.containsVertex(i) || myFaces.back.containsVertex(i))
-        dirichletNodes[i][2] = true;
-#endif
-    }
-
-    // Set up functions for time-dependent boundary conditions
-    using Function = Dune::VirtualFunction<double, double>;
-    Function const &velocityDirichletFunction = VelocityDirichletCondition();
-    Function const &neumannFunction = NeumannCondition();
-
-    MyAssembler const myAssembler(leafView);
-
-    MyBody<dims> const body(parset);
-
-    Matrices<Matrix> matrices;
-    myAssembler.assembleElasticity(body.getYoungModulus(),
-                                   body.getPoissonRatio(), matrices.elasticity);
-    myAssembler.assembleViscosity(body.getShearViscosityField(),
-                                  body.getBulkViscosityField(),
-                                  matrices.damping);
-    myAssembler.assembleMass(body.getDensityField(), matrices.mass);
-
-    ScalarMatrix relativeFrictionalBoundaryMass;
-    myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
-                                               relativeFrictionalBoundaryMass);
-    relativeFrictionalBoundaryMass /= frictionalBoundary.area();
-    EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
-        relativeFrictionalBoundaryMass);
-
-    // Assemble forces
-    Vector gravityFunctional;
-    myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
-
-    // Problem formulation: right-hand side
-    std::function<void(double, Vector &)> computeExternalForces =
-        [&](double _relativeTime, Vector &_ell) {
-          myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
-                                      _relativeTime);
-          _ell += gravityFunctional;
-        };
-
-    using MyProgramState = ProgramState<Vector, ScalarVector>;
-    MyProgramState programState(leafVertexCount);
-    auto const firstRestart = parset.get<size_t>("io.restarts.first");
-    auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
-    auto const writeRestarts = parset.get<bool>("io.restarts.write");
-    auto const writeData = parset.get<bool>("io.data.write");
-    bool const handleRestarts = writeRestarts or firstRestart > 0;
-
-    auto dataFile =
-        writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
-    auto restartFile = handleRestarts
-                           ? std::make_unique<HDF5::File>(
-                                 "restarts.h5",
-                                 writeRestarts ? HDF5::Access::READWRITE
-                                               : HDF5::Access::READONLY)
-                           : nullptr;
-    auto restartIO = handleRestarts
-                         ? std::make_unique<RestartIO<MyProgramState>>(
-                               *restartFile, leafVertexCount)
-                         : nullptr;
-
-    if (firstRestart > 0) // automatically adjusts the time and timestep
-      restartIO->read(firstRestart, programState);
-    else
-      programState.setupInitialConditions(parset, computeExternalForces,
-                                          matrices, myAssembler, dirichletNodes,
-                                          noNodes, frictionalBoundary, body);
-
-    MyGlobalFrictionData<LocalVector> frictionInfo(
-        parset.sub("boundary.friction"), weakPatch);
-    auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
-        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
-        frictionalBoundary, frictionInfo, programState.weightedNormalStress);
-    myGlobalFriction->updateAlpha(programState.alpha);
-
-    Vector vertexCoordinates(leafVertexCount);
-    {
-      Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
-      for (auto &&v : vertices(leafView))
-        vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
-    }
-
-    using MyVertexBasis = typename MyAssembler::VertexBasis;
-    auto dataWriter =
-        writeData ? std::make_unique<
-                        HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
-                        *dataFile, vertexCoordinates, myAssembler.vertexBasis,
-                        surface, frictionalBoundary, weakPatch)
-                  : nullptr;
-    MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
-        myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
-
-    IterationRegister iterationCount;
-    auto const report = [&](bool initial = false) {
-      if (writeData) {
-        dataWriter->reportSolution(programState, *myGlobalFriction);
-        if (!initial)
-          dataWriter->reportIterations(programState, iterationCount);
-        dataFile->flush();
-      }
-
-      if (writeRestarts and !initial and
-          programState.timeStep % restartSpacing == 0) {
-        restartIO->write(programState);
-        restartFile->flush();
-      }
-
-      if (parset.get<bool>("io.printProgress"))
-        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
-                  << ", time = " << std::setw(12) << programState.relativeTime
-                  << ", tau = " << std::setw(12) << programState.relativeTau
-                  << std::endl;
-
-      if (parset.get<bool>("io.vtk.write")) {
-        ScalarVector stress;
-        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                           body.getPoissonRatio(),
-                                           programState.u, stress);
-        vtkWriter.write(programState.timeStep, programState.u, programState.v,
-                        programState.alpha, stress);
-      }
-    };
-    report(true);
-
-    // Set up TNNMG solver
-    using NonlinearFactory = SolverFactory<
-        dims,
-        MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-        Grid>;
-    NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
-
-    using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
-                               StateUpdater<ScalarVector, Vector>>;
-    MyUpdater current(
-        initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
-                        velocityDirichletFunction, dirichletNodes, matrices,
-                        programState.u, programState.v, programState.a),
-        initStateUpdater<ScalarVector, Vector>(
-            parset.get<Config::stateModel>("boundary.friction.stateModel"),
-            programState.alpha, *frictionalBoundary.getVertices(),
-            parset.get<double>("boundary.friction.L"),
-            parset.get<double>("boundary.friction.V0")));
-
-    auto const refinementTolerance =
-        parset.get<double>("timeSteps.refinementTolerance");
-    auto const mustRefine = [&](MyUpdater &coarseUpdater,
-                                MyUpdater &fineUpdater) {
-      ScalarVector coarseAlpha;
-      coarseUpdater.state_->extractAlpha(coarseAlpha);
-
-      ScalarVector fineAlpha;
-      fineUpdater.state_->extractAlpha(fineAlpha);
-
-      return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
-    };
-
-    std::signal(SIGXCPU, handleSignal);
-    std::signal(SIGINT, handleSignal);
-    std::signal(SIGTERM, handleSignal);
-
-    AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
-                        EnergyNorm<ScalarMatrix, ScalarVector>>
-        adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
-                            programState.relativeTime, programState.relativeTau,
-                            computeExternalForces, stateEnergyNorm, mustRefine);
-    while (!adaptiveTimeStepper.reachedEnd()) {
-      programState.timeStep++;
-
-      iterationCount = adaptiveTimeStepper.advance();
-
-      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
-      programState.relativeTau = adaptiveTimeStepper.relativeTau_;
-      current.rate_->extractDisplacement(programState.u);
-      current.rate_->extractVelocity(programState.v);
-      current.rate_->extractAcceleration(programState.a);
-      current.state_->extractAlpha(programState.alpha);
-
-      report();
-
-      if (terminationRequested) {
-        std::cerr << "Terminating prematurely" << std::endl;
-        break;
-      }
-    }
-  } catch (Dune::Exception &e) {
-    Dune::derr << "Dune reported error: " << e << std::endl;
-  } catch (std::exception &e) {
-    std::cerr << "Standard exception: " << e.what() << std::endl;
-  }
-}