src/hdf5/frictionalboundary-writer.hh

-#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
-#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-  using Patch = BoundaryPatch<GridView>;
-
-public:
-  FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                           Patch const &frictionalBoundary);
-
-  template <class Friction>
-  void write(ProgramState const &programState, Friction &friction);
-
-private:
-  HDF5::Group group_;
-
-  Patch const &frictionalBoundary_;
-
-  HDF5::SequenceIO<2> frictionalBoundaryDisplacementWriter_;
-  HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
-  HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
-  HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
-};
-#endif

src/hdf5/frictionalboundary-writer_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using MyFriction = GlobalFriction<Matrix, Vector>;
-
-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
-    MyProgramState const &programState, MyFriction &friction);

src/hdf5/iteration-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "iteration-writer.hh"
-
-IterationWriter::IterationWriter(HDF5::Grouplike &file)
-    : group_(file, "iterations"),
-      fpiSubGroup_(group_, "fixedPoint"),
-      mgSubGroup_(group_, "multiGrid"),
-      finalMGIterationWriter_(mgSubGroup_, "final"),
-      finalFPIIterationWriter_(fpiSubGroup_, "final"),
-      totalMGIterationWriter_(mgSubGroup_, "total"),
-      totalFPIIterationWriter_(fpiSubGroup_, "total") {}
-
-void IterationWriter::write(size_t timeStep,
-                            IterationRegister const &iterationCount) {
-  addEntry(finalMGIterationWriter_, timeStep,
-           iterationCount.finalCount.multigridIterations);
-  addEntry(finalFPIIterationWriter_, timeStep,
-           iterationCount.finalCount.iterations);
-  addEntry(totalMGIterationWriter_, timeStep,
-           iterationCount.totalCount.multigridIterations);
-  addEntry(totalFPIIterationWriter_, timeStep,
-           iterationCount.totalCount.iterations);
-}

src/hdf5/iteration-writer.hh

-#ifndef SRC_HDF_ITERATION_WRITER_HH
-#define SRC_HDF_ITERATION_WRITER_HH
-
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-#include "../time-stepping/adaptivetimestepper.hh"
-
-class IterationWriter {
-public:
-  IterationWriter(HDF5::Grouplike &file);
-
-  void write(size_t timeStep, IterationRegister const &iterationCount);
-
-private:
-  HDF5::Group group_;
-  HDF5::Group fpiSubGroup_;
-  HDF5::Group mgSubGroup_;
-
-  HDF5::SequenceIO<0, size_t> finalMGIterationWriter_;
-  HDF5::SequenceIO<0, size_t> finalFPIIterationWriter_;
-  HDF5::SequenceIO<0, size_t> totalMGIterationWriter_;
-  HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_;
-};
-#endif

src/hdf5/patchinfo-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/fufem/grid/hierarchic-approximation.hh>
-#include <dune/fufem/hdf5/singletonwriter.hh>
-
-#include "patchinfo-writer.hh"
-
-template <class LocalVector, class GridView>
-GridEvaluator<LocalVector, GridView>::GridEvaluator(
-    ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
-  double const bufferWidth = 0.05 * MyGeometry::lengthScale;
-  auto const xminmax = std::minmax_element(
-      weakPatch.vertices.begin(), weakPatch.vertices.end(),
-      [](LocalVector const &a, LocalVector const &b) { return a[0] < b[0]; });
-  double const xmin = (*xminmax.first)[0] - bufferWidth;
-  double const xspan = (*xminmax.second)[0] + bufferWidth - xmin;
-
-  double const zmin = -MyGeometry::depth / 2.0;
-  double const zspan = MyGeometry::depth;
-
-  double spatialResolution = 0.01 * MyGeometry::lengthScale;
-  size_t const xsteps = std::round(xspan / spatialResolution);
-  size_t const zsteps = std::round(zspan / spatialResolution);
-
-  xCoordinates.resize(xsteps + 1);
-  for (size_t xi = 0; xi <= xsteps; xi++)
-    xCoordinates[xi] = xmin + xi * xspan / xsteps;
-  zCoordinates.resize(zsteps + 1);
-  for (size_t zi = 0; zi <= zsteps; zi++)
-    zCoordinates[zi] = zmin + zi * zspan / zsteps;
-
-  HierarchicApproximation<typename GridView::Grid, GridView> const
-      hApproximation(gridView.grid(), gridView, 1e-6 * MyGeometry::lengthScale);
-
-  LocalVector global(0);
-  localInfo.resize(xsteps + 1);
-  for (size_t xi = 0; xi < xCoordinates.size(); ++xi) {
-    localInfo[xi].resize(zsteps + 1);
-    for (size_t zi = 0; zi < zCoordinates.size(); ++zi) {
-      global[0] = xCoordinates[xi];
-      global[2] = zCoordinates[zi];
-      auto const element = hApproximation.findEntity(global);
-      localInfo[xi][zi] =
-          std::make_pair(element, element.geometry().local(global));
-    }
-  }
-}
-
-template <class LocalVector, class GridView>
-template <class Function>
-Dune::Matrix<typename Function::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(Function const &f) const {
-  Dune::Matrix<typename Function::RangeType> ret(xCoordinates.size(),
-                                                 zCoordinates.size());
-  for (size_t xi = 0; xi < localInfo.size(); ++xi) {
-    auto const &localInfoX = localInfo[xi];
-    for (size_t zi = 0; zi < localInfoX.size(); ++zi) {
-      auto const &localInfoXZ = localInfoX[zi];
-      f.evaluateLocal(localInfoXZ.first, localInfoXZ.second, ret[xi][zi]);
-    }
-  }
-  return ret;
-}
-
-template <class ProgramState, class VertexBasis, class GridView>
-PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
-    HDF5::Grouplike &file, VertexBasis const &vertexBasis,
-    Patch const &frictionalBoundary,
-    ConvexPolyhedron<LocalVector> const &weakPatch)
-    : group_(file, "weakPatchGrid"),
-      vertexBasis_(vertexBasis),
-      gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
-      weakPatchGridVelocityWriter_(
-          group_, "velocity", gridEvaluator_.xCoordinates.size(),
-          gridEvaluator_.zCoordinates.size(), Vector::block_type::dimension) {
-  HDF5::SingletonWriter<1> weakPatchGridXCoordinateWriter(
-      group_, "xCoordinates", gridEvaluator_.xCoordinates.size());
-  setEntry(weakPatchGridXCoordinateWriter, gridEvaluator_.xCoordinates);
-
-  HDF5::SingletonWriter<1> weakPatchGridZCoordinateWriter(
-      group_, "zCoordinates", gridEvaluator_.zCoordinates.size());
-  setEntry(weakPatchGridZCoordinateWriter, gridEvaluator_.zCoordinates);
-}
-
-template <class ProgramState, class VertexBasis, class GridView>
-void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
-    ProgramState const &programState) {
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
-  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
-  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
-}
-
-#include "patchinfo-writer_tmpl.cc"

src/hdf5/patchinfo-writer.hh

-#ifndef SRC_HDF_PATCHINFO_WRITER_HH
-#define SRC_HDF_PATCHINFO_WRITER_HH
-
-#include <dune/istl/matrix.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/geometry/convexpolyhedron.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-#include "../one-body-problem-data/mygeometry.hh"
-
-template <class LocalVector, class GridView> class GridEvaluator {
-  using Element = typename GridView::Grid::template Codim<0>::Entity;
-
-public:
-  GridEvaluator(ConvexPolyhedron<LocalVector> const &weakPatch,
-                GridView const &gridView);
-
-  template <class Function>
-  Dune::Matrix<typename Function::RangeType> evaluate(Function const &f) const;
-
-  Dune::BlockVector<Dune::FieldVector<double, 1>> xCoordinates;
-  Dune::BlockVector<Dune::FieldVector<double, 1>> zCoordinates;
-
-private:
-  std::vector<std::vector<std::pair<Element, LocalVector>>> localInfo;
-};
-
-template <class ProgramState, class VertexBasis, class GridView>
-class PatchInfoWriter {
-  using Vector = typename ProgramState::Vector;
-  using LocalVector = typename Vector::block_type;
-  using Patch = BoundaryPatch<GridView>;
-
-public:
-  PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
-                  Patch const &frictionalBoundary,
-                  ConvexPolyhedron<LocalVector> const &weakPatch);
-
-  void write(ProgramState const &programState);
-
-private:
-  HDF5::Group group_;
-
-  VertexBasis const &vertexBasis_;
-
-  GridEvaluator<LocalVector, GridView> const gridEvaluator_;
-  HDF5::SequenceIO<3> weakPatchGridVelocityWriter_;
-};
-#endif

src/hdf5/patchinfo-writer_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<GridView, double>;
-using MyFunction = BasisGridFunction<P1Basis, Vector>;
-
-template class GridEvaluator<LocalVector, GridView>;
-template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
-
-template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;

src/hdf5/restart-io.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "restart-io.hh"
-
-template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
-    : displacementWriter_(group, "displacement", vertexCount,
-                          Vector::block_type::dimension),
-      velocityWriter_(group, "velocity", vertexCount,
-                      Vector::block_type::dimension),
-      accelerationWriter_(group, "acceleration", vertexCount,
-                          Vector::block_type::dimension),
-      stateWriter_(group, "state", vertexCount),
-      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
-      relativeTimeWriter_(group, "relativeTime"),
-      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
-
-template <class ProgramState>
-void RestartIO<ProgramState>::write(ProgramState const &programState) {
-  addEntry(displacementWriter_, programState.timeStep, programState.u);
-  addEntry(velocityWriter_, programState.timeStep, programState.v);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha);
-  addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress);
-  addEntry(relativeTimeWriter_, programState.timeStep,
-           programState.relativeTime);
-  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
-           programState.relativeTau);
-}
-
-template <class ProgramState>
-void RestartIO<ProgramState>::read(size_t timeStep,
-                                   ProgramState &programState) {
-  programState.timeStep = timeStep;
-  readEntry(displacementWriter_, timeStep, programState.u);
-  readEntry(velocityWriter_, timeStep, programState.v);
-  readEntry(accelerationWriter_, timeStep, programState.a);
-  readEntry(stateWriter_, timeStep, programState.alpha);
-  readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress);
-  readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
-  readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
-}
-
-#include "restart-io_tmpl.cc"

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;
+  }
+}