diff --git a/data/tools/2d-performance.py b/data/tools/2d-performance.py
new file mode 100644
index 0000000000000000000000000000000000000000..1d2a019d1e3a3b105baa7e9cba168a36d42569ac
--- /dev/null
+++ b/data/tools/2d-performance.py
@@ -0,0 +1,76 @@
+import ConfigParser as cp
+import os
+import numpy as np
+import csv
+import h5py
+
+from support.find_quakes import find_quakes
+
+NBODIES = 2
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# create output directory
+out_dir = os.path.join(out_path)
+if not os.path.exists(out_dir):
+    os.mkdir(out_dir)
+
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+real_time = relative_time * FINAL_TIME
+
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  ## We are interested in an enlarged time range around actual events here,
+  ## (and no other quantities!) hence we pass a very low velocity here.
+  quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
+                       indices = 1:dim(velocityProxy)[1], 1)
+  quakes$beginning <- realTime[quakes$beginningIndex]
+  quakes$ending    <- realTime[quakes$endingIndex]
+  quakes$duration  <- quakes$ending - quakes$beginning
+  numQuakes        <- nrow(quakes)
+
+  relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
+                               quakes[[numQuakes,  'ending']]), f=0.02)
+  threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+  plotMask   <- threeQuakeTimeMask
+  plotIndices<- which(plotMask)
+  plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
+
+  write(relaxedTime[[1]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'min', 'threequakes',
+                                    basedir), 'tex')))
+  write(relaxedTime[[2]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'max', 'threequakes',
+                                    basedir), 'tex')))
+
+  timeWindow = realTime[plotIndices]
+
+  relativeTau          <- h5file['relativeTimeIncrement'][]
+  fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
+  multiGridIterations  <- h5file['/iterations/multiGrid/final'][]
+  write.csv(data.frame(time = timeWindow,
+                       timeIncrement = finalTime * relativeTau[plotIndices],
+                       fixedPointIterations = fixedPointIterations[plotIndices],
+                       multiGridIterations = multiGridIterations[plotIndices]),
+            file.path(directories[['output']],
+                      paste.(pasteColon('2d-performance', basedir), 'csv')),
+            row.names = FALSE, quote = FALSE)
+  h5::h5close(h5file)
diff --git a/dune/tectonic/CMakeLists.txt b/dune/tectonic/CMakeLists.txt
index de041accbf959ae3928646a83ab796ed97e1545a..90dce1b6cb5af6f35ac3201bb070cf85bf68a4d4 100644
--- a/dune/tectonic/CMakeLists.txt
+++ b/dune/tectonic/CMakeLists.txt
@@ -1,6 +1,7 @@
 add_subdirectory("data-structures")
 add_subdirectory("factories")
 add_subdirectory("io")
+add_subdirectory("multi-threading")
 add_subdirectory("problem-data")
 add_subdirectory("spatial-solving")
 add_subdirectory("tests")
diff --git a/dune/tectonic/assemblers.cc b/dune/tectonic/assemblers.cc
index 8cfc399ef56a1ea0a0f06d1cf922e20f669ca44a..37c8cb39aa662cfd2d78d5c0b095d1b7ffa9aabe 100644
--- a/dune/tectonic/assemblers.cc
+++ b/dune/tectonic/assemblers.cc
@@ -4,7 +4,8 @@
 
 #include <dune/istl/scaledidmatrix.hh>
 
-#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/boundaryoperatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
@@ -45,11 +46,9 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
         BoundaryPatch<GridView> const &frictionalBoundary,
         ScalarMatrix &frictionalBoundaryMass) const {
 
-  BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
-                        LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
-      frictionalBoundaryMassAssembler(frictionalBoundary);
-  vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
-                                   frictionalBoundaryMass);
+  MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis> localMassAssembler;
+  const BoundaryOperatorAssembler<VertexBasis, VertexBasis> boundaryAssembler(vertexBasis, vertexBasis, frictionalBoundary);
+  boundaryAssembler.assemble(localMassAssembler, frictionalBoundaryMass);
 }
 
 template <class GridView, int dimension>
diff --git a/dune/tectonic/assemblers.hh b/dune/tectonic/assemblers.hh
index 9ec78078b313d2ba302a2b8bdadf09d924ce449c..e5354141074f9646a6a0407b41d66c978018e707 100644
--- a/dune/tectonic/assemblers.hh
+++ b/dune/tectonic/assemblers.hh
@@ -22,6 +22,7 @@
 template <class GridView, int dimension>
 class MyAssembler {
 public:
+    using GV = GridView;
     using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using Matrix =
       Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
diff --git a/dune/tectonic/io/CMakeLists.txt b/dune/tectonic/io/CMakeLists.txt
index 07045307f72849e735029517fce006ffa6a95150..007f2ad985c8c547faadb507aa50ebbdc044e8e3 100644
--- a/dune/tectonic/io/CMakeLists.txt
+++ b/dune/tectonic/io/CMakeLists.txt
@@ -2,15 +2,18 @@ add_subdirectory("hdf5")
 
 add_custom_target(tectonic_dune_io SOURCES
   hdf5-bodywriter.hh
-  hdf5-writer.hh 
+  hdf5-writer.hh
+  io-handler.hh
   uniform-grid-writer.cc
+  vtk-bodywriter.hh
   vtk.hh
-  vtk.cc
 )
 
 #install headers
 install(FILES
   hdf5-bodywriter.hh
-  hdf5-levelwriter.hh 
+  hdf5-levelwriter.hh
+  io-handler.hh
+  vtk-bodywriter.hh
   vtk.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/io/hdf5-bodywriter.hh b/dune/tectonic/io/hdf5-bodywriter.hh
index 91dd5edeb45db5bf93bddc0f37a8937ba566c08c..cf47eaad468a2be645b08b490ac40a6354e1203b 100644
--- a/dune/tectonic/io/hdf5-bodywriter.hh
+++ b/dune/tectonic/io/hdf5-bodywriter.hh
@@ -18,7 +18,7 @@ class HDF5BodyWriter {
 public:
     using VertexCoordinates = typename ProgramState::Vector;
     using Patch = typename FrictionalBoundaryWriter<GridView>::Patch;
-    using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
+    //using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
 
     using LocalVector = typename VertexCoordinates::block_type;
 
@@ -29,8 +29,8 @@ class HDF5BodyWriter {
             const size_t bodyID,
             const VertexCoordinates& vertexCoordinates,
             const VertexBasis& vertexBasis,
-            const Patch& frictionPatch,
-            const WeakPatches& weakPatches) :
+            const Patch& frictionPatch) :
+           // const WeakPatches& weakPatches) :
         id_(bodyID),
         /*#if MY_DIM == 3
                 patchInfoWriters_(patchCount_),
diff --git a/dune/tectonic/io/hdf5-writer.hh b/dune/tectonic/io/hdf5-writer.hh
index 2a8ecc2fa05b495fd4f9f16f449a14dc051d3534..006201221902a8131c95ffb7b29e5fa4b71cc211 100644
--- a/dune/tectonic/io/hdf5-writer.hh
+++ b/dune/tectonic/io/hdf5-writer.hh
@@ -1,12 +1,14 @@
 #ifndef SRC_HDF5_WRITER_HH
 #define SRC_HDF5_WRITER_HH
 
-#include <dune/fufem/hdf5/file.hh>
+#include <vector>
 
 #include "hdf5/iteration-writer.hh"
 #include "hdf5/time-writer.hh"
 #include "hdf5-bodywriter.hh"
 
+#include <dune/fufem/hdf5/file.hh>
+
 template <class ProgramState, class VertexBasis, class GridView>
 class HDF5Writer {
 public:
@@ -14,15 +16,15 @@ class HDF5Writer {
     using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
     using VertexBases = std::vector<const VertexBasis*>;
     using FrictionPatches = std::vector<const typename HDF5BodyWriter::Patch*>;
-    using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
+    //using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
 
     //friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
 
     HDF5Writer(HDF5::File& file,
              const VertexCoordinates& vertexCoordinates,
              const VertexBases& vertexBases,
-             const FrictionPatches& frictionPatches,
-             const WeakPatches& weakPatches)
+             const FrictionPatches& frictionPatches)
+             //const WeakPatches& weakPatches)
       : file_(file),
         iterationWriter_(file_),
         timeWriter_(file_),
@@ -30,7 +32,7 @@ class HDF5Writer {
 
         for (size_t i=0; i<frictionPatches_.size(); i++) {
             if (frictionPatches_[i]->numVertices() > 0)
-                bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i], weakPatches[i]));
+                bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i])); //, weakPatches[i]));
         }
     }
 
diff --git a/dune/tectonic/io/hdf5/CMakeLists.txt b/dune/tectonic/io/hdf5/CMakeLists.txt
index 511ce68219d2969ae308ba866715f16c6d2fdb68..0483a7049df256b6e507d48e277cc846f5a1468a 100644
--- a/dune/tectonic/io/hdf5/CMakeLists.txt
+++ b/dune/tectonic/io/hdf5/CMakeLists.txt
@@ -1,24 +1,23 @@
 add_custom_target(tectonic_dune_io_hdf5 SOURCES
   frictionalboundary-writer.hh
-  frictionalboundary-writer.cc
+  #frictionalboundary-writer.cc
   iteration-writer.hh 
-  iteration-writer.cc
-  patchinfo-writer.hh 
-  patchinfo-writer.cc
+  #iteration-writer.cc
+  #patchinfo-writer.hh
+  #patchinfo-writer.cc
+  restartbody-io.hh
   restart-io.hh
-  restart-io.cc
   restrict.hh
   surface-writer.hh
-  surface-writer.cc
   time-writer.hh
-  time-writer.cc
 )
 
 #install headers
 install(FILES
   frictionalboundary-writer.hh
   iteration-writer.hh 
-  patchinfo-writer.hh 
+  #patchinfo-writer.hh
+  restartbody-io.hh
   restart-io.hh
   restrict.hh
   surface-writer.hh
diff --git a/dune/tectonic/io/hdf5/frictionalboundary-writer.cc b/dune/tectonic/io/hdf5/frictionalboundary-writer.cc
deleted file mode 100644
index 6bcfa353a4bece1ad77cc148d052afb65875a36a..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer.cc
+++ /dev/null
@@ -1,53 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/fufem/hdf5/singletonwriter.hh>
-
-#include "frictionalboundary-writer.hh"
-#include "restrict.hh"
-
-template <class GridView>
-template <class Vector>
-FrictionalBoundaryWriter<GridView>::FrictionalBoundaryWriter(
-    HDF5::Group& group, const Vector& vertexCoordinates,
-    const Patch& frictionalBoundary)
-    : group_(group),
-      frictionalBoundary_(frictionalBoundary),
-      frictionalBoundaryDisplacementWriter_(group_, "displacement",
-                                            frictionalBoundary.numVertices(),
-                                            Vector::block_type::dimension),
-      frictionalBoundaryVelocityWriter_(group_, "velocity",
-                                        frictionalBoundary.numVertices(),
-                                        Vector::block_type::dimension),
-      frictionalBoundaryStateWriter_(group_, "state",
-                                     frictionalBoundary.numVertices()),
-      frictionalBoundaryCoefficientWriter_(group_, "coefficient",
-                                           frictionalBoundary.numVertices()) {
-  auto const frictionalBoundaryCoordinates =
-      restrictToSurface(vertexCoordinates, frictionalBoundary);
-  HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
-      group_, "coordinates", frictionalBoundaryCoordinates.size(),
-      Vector::block_type::dimension);
-  setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
-}
-
-template <class GridView>
-template <class Vector, class ScalarVector>
-void FrictionalBoundaryWriter<GridView>::write(
-    const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff) {
-
-  auto const frictionalBoundaryDisplacements = restrictToSurface(u, frictionalBoundary_);
-  addEntry(frictionalBoundaryDisplacementWriter_, timeStep, frictionalBoundaryDisplacements);
-
-  auto const frictionalBoundaryVelocities = restrictToSurface(v, frictionalBoundary_);
-  addEntry(frictionalBoundaryVelocityWriter_, timeStep, frictionalBoundaryVelocities);
-
-  auto const frictionalBoundaryStates = restrictToSurface(alpha, frictionalBoundary_);
-  addEntry(frictionalBoundaryStateWriter_, timeStep, frictionalBoundaryStates);
-
-  auto const frictionalBoundaryCoefficient = restrictToSurface(frictionCoeff, frictionalBoundary_);
-  addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
-}
-
-#include "frictionalboundary-writer_tmpl.cc"
diff --git a/dune/tectonic/io/hdf5/frictionalboundary-writer.hh b/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
index a7ca0855da22a302d4130fd9223f7f90050a895d..bc8ca848b0c05e3aec854ae34b946b7b834b2665 100644
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
+++ b/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
@@ -1,9 +1,12 @@
 #ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
 #define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
 
-
 #include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/hdf5/file.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+#include "restrict.hh"
 
 template <class GridView> class FrictionalBoundaryWriter {
 public:
@@ -11,10 +14,42 @@ template <class GridView> class FrictionalBoundaryWriter {
 
   template <class Vector>
   FrictionalBoundaryWriter(HDF5::Group& group, const Vector& vertexCoordinates,
-                           const Patch& frictionalBoundary);
+                           const Patch& frictionalBoundary) : group_(group),
+      frictionalBoundary_(frictionalBoundary),
+      frictionalBoundaryDisplacementWriter_(group_, "displacement",
+                                            frictionalBoundary.numVertices(),
+                                            Vector::block_type::dimension),
+      frictionalBoundaryVelocityWriter_(group_, "velocity",
+                                        frictionalBoundary.numVertices(),
+                                        Vector::block_type::dimension),
+      frictionalBoundaryStateWriter_(group_, "state",
+                                     frictionalBoundary.numVertices()),
+      frictionalBoundaryCoefficientWriter_(group_, "coefficient",
+                                           frictionalBoundary.numVertices()) {
+
+  auto const frictionalBoundaryCoordinates =
+      restrictToSurface(vertexCoordinates, frictionalBoundary);
+  HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
+      group_, "coordinates", frictionalBoundaryCoordinates.size(),
+      Vector::block_type::dimension);
+  setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
+}
 
   template <class Vector, class ScalarVector>
-  void write(const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff);
+  void write(const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff) {
+
+      auto const frictionalBoundaryDisplacements = restrictToSurface(u, frictionalBoundary_);
+      addEntry(frictionalBoundaryDisplacementWriter_, timeStep, frictionalBoundaryDisplacements);
+
+      auto const frictionalBoundaryVelocities = restrictToSurface(v, frictionalBoundary_);
+      addEntry(frictionalBoundaryVelocityWriter_, timeStep, frictionalBoundaryVelocities);
+
+      auto const frictionalBoundaryStates = restrictToSurface(alpha, frictionalBoundary_);
+      addEntry(frictionalBoundaryStateWriter_, timeStep, frictionalBoundaryStates);
+
+      auto const frictionalBoundaryCoefficient = restrictToSurface(frictionCoeff, frictionalBoundary_);
+      addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
+    }
 
 private:
   HDF5::Group& group_;
@@ -26,4 +61,5 @@ template <class GridView> class FrictionalBoundaryWriter {
   HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
   HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
 };
+
 #endif
diff --git a/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc b/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc
deleted file mode 100644
index 51633c69e27c34effb24f3925e162c9c3c5b5d5b..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc
+++ /dev/null
@@ -1,19 +0,0 @@
-#include "../../explicitvectors.hh"
-#include "../../explicitgrid.hh"
-
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-//#include "../../data-structures/program_state.hh"
-
-//using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class FrictionalBoundaryWriter<DefLeafGridView>;
-
-template FrictionalBoundaryWriter<DefLeafGridView>::FrictionalBoundaryWriter(
-    HDF5::Group& group, const Vector& vertexCoordinates, const BoundaryPatch<DefLeafGridView>& frictionalBoundary);
-template void FrictionalBoundaryWriter<DefLeafGridView>::write(
-    const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff);
-
-
diff --git a/dune/tectonic/io/hdf5/iteration-writer.cc b/dune/tectonic/io/hdf5/iteration-writer.cc
deleted file mode 100644
index 568c1cf5056e9fae81d43f594dc7b7704b4740e8..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/iteration-writer.cc
+++ /dev/null
@@ -1,26 +0,0 @@
-#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);
-}
diff --git a/dune/tectonic/io/hdf5/iteration-writer.hh b/dune/tectonic/io/hdf5/iteration-writer.hh
index d0e010f18e5c7723fb7e1704bcd0d91cd11c899d..23a4d81e9a189355425808e41c419824f7c293d9 100644
--- a/dune/tectonic/io/hdf5/iteration-writer.hh
+++ b/dune/tectonic/io/hdf5/iteration-writer.hh
@@ -1,15 +1,32 @@
-#ifndef SRC_HDF_ITERATION_WRITER_HH
-#define SRC_HDF_ITERATION_WRITER_HH
+#ifndef SRC_HDF5_ITERATION_WRITER_HH
+#define SRC_HDF5_ITERATION_WRITER_HH
 
-#include <dune/fufem/hdf5/sequenceio.hh>
 
 #include "../../time-stepping/adaptivetimestepper.hh"
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
 
 class IterationWriter {
 public:
-  IterationWriter(HDF5::Grouplike &file);
+  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 write(size_t timeStep, IterationRegister const &iterationCount);
+  void write(size_t timeStep, const IterationRegister& 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);
+  }
 
 private:
   HDF5::Group group_;
@@ -21,4 +38,5 @@ class IterationWriter {
   HDF5::SequenceIO<0, size_t> totalMGIterationWriter_;
   HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_;
 };
+
 #endif
diff --git a/dune/tectonic/io/hdf5/patchinfo-writer.cc b/dune/tectonic/io/hdf5/patchinfo-writer.cc
deleted file mode 100644
index c7a38dda5f00b73272669717a12692bed56a8b8d..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/patchinfo-writer.cc
+++ /dev/null
@@ -1,98 +0,0 @@
-#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(const CuboidGeometry& cuboidGeometry,
-    ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
-  double const bufferWidth = 0.05 * cuboidGeometry.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 = -cuboidGeometry.depth() / 2.0;
-  double const zspan = cuboidGeometry.depth();
-
-  double spatialResolution = 0.01 * cuboidGeometry.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 * cuboidGeometry.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,
-    size_t const id)
-    : id_(id),
-      group_(file, "weakPatchGrid"+std::to_string(id_)),
-      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) {
-    // TODO
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]);
-  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
-  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
-}
-
-#include "patchinfo-writer_tmpl.cc"
diff --git a/dune/tectonic/io/hdf5/patchinfo-writer.hh b/dune/tectonic/io/hdf5/patchinfo-writer.hh
deleted file mode 100644
index 567c36e987002807c9797d53cd5b8fea81e9b2cc..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/patchinfo-writer.hh
+++ /dev/null
@@ -1,56 +0,0 @@
-#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 "../../problem-data/grid/cuboidgeometry.hh"
-
-template <class LocalVector, class GridView> class GridEvaluator {
-  using Element = typename GridView::Grid::template Codim<0>::Entity;
-  using CuboidGeometry = CuboidGeometry<typename LocalVector::field_type>;
-
-public:
-  GridEvaluator(const CuboidGeometry& cuboidGeometry,
-                const ConvexPolyhedron<LocalVector>& weakPatch,
-                const GridView& 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,
-                  size_t const id);
-
-  void write(ProgramState const &programState);
-
-private:
-  const size_t id_;
-
-  HDF5::Group group_;
-
-  VertexBasis const &vertexBasis_;
-
-  GridEvaluator<LocalVector, GridView> const gridEvaluator_;
-  HDF5::SequenceIO<3> weakPatchGridVelocityWriter_;
-};
-#endif
diff --git a/dune/tectonic/io/hdf5/patchinfo-writer_tmpl.cc b/dune/tectonic/io/hdf5/patchinfo-writer_tmpl.cc
deleted file mode 100644
index 0b281fead6022b1d4616db26be02fc99a08598dd..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/patchinfo-writer_tmpl.cc
+++ /dev/null
@@ -1,15 +0,0 @@
-#include "../../explicitvectors.hh"
-#include "../../explicitgrid.hh"
-
-#include "../../data-structures/program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
-using MyFunction = BasisGridFunction<P1Basis, Vector>;
-
-template class GridEvaluator<LocalVector, DeformedGrid::LeafGridView>;
-
-template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, DeformedGrid::LeafGridView>::evaluate(MyFunction const &f) const;
-
-template class PatchInfoWriter<MyProgramState, P1Basis, DeformedGrid::LeafGridView>;
diff --git a/dune/tectonic/io/hdf5/restart-io.cc b/dune/tectonic/io/hdf5/restart-io.cc
deleted file mode 100644
index d194095926e8728d6e0ced5c4a799f3aeffb29d4..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/restart-io.cc
+++ /dev/null
@@ -1,89 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "restart-io.hh"
-
-template <class ProgramState>
-RestartBodyIO<ProgramState>::RestartBodyIO(HDF5::Grouplike &group, const size_t vertexCount, const size_t id)
-    : id_(id),
-      displacementWriter_(group, "displacement"+std::to_string(id_), vertexCount,
-                          ProgramState::Vector::block_type::dimension),
-      velocityWriter_(group, "velocity"+std::to_string(id_), vertexCount,
-                      ProgramState::Vector::block_type::dimension),
-      accelerationWriter_(group, "acceleration"+std::to_string(id_), vertexCount,
-                          ProgramState::Vector::block_type::dimension),
-      stateWriter_(group, "state"+std::to_string(id_), vertexCount),
-      weightedNormalStressWriter_(group, "weightedNormalStress"+std::to_string(id_), vertexCount) {}
-
-template <class ProgramState>
-void RestartBodyIO<ProgramState>::write(const ProgramState& programState) {
-
-  addEntry(displacementWriter_, programState.timeStep, programState.u[id_]);
-  addEntry(velocityWriter_, programState.timeStep, programState.v[id_]);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a[id_]);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha[id_]);
-  addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress[id_]);
-}
-
-template <class ProgramState>
-void RestartBodyIO<ProgramState>::read(size_t timeStep,
-                                   ProgramState& programState) {
-
-  readEntry(displacementWriter_, timeStep, programState.u[id_]);
-  readEntry(velocityWriter_, timeStep, programState.v[id_]);
-  readEntry(accelerationWriter_, timeStep, programState.a[id_]);
-  readEntry(stateWriter_, timeStep, programState.alpha[id_]);
-  readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress[id_]);
-}
-
-template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts)
-    : bodyIO_(vertexCounts.size()),
-      relativeTimeWriter_(group, "relativeTime"),
-      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {
-
-  for (size_t i=0; i<bodyIO_.size(); i++) {
-    bodyIO_[i] = new RestartBodyIO<ProgramState>(group, vertexCounts[i], i);
-  }
-}
-
-template <class ProgramState>
-RestartIO<ProgramState>::~RestartIO() {
-  for (size_t i=0; i<bodyIO_.size(); i++) {
-    delete bodyIO_[i];
-  }
-}
-
-template <class ProgramState>
-void RestartIO<ProgramState>::write(ProgramState const &programState) {
-  assert(programState.size() == bodyIO_.size());
-
-  for (size_t i=0; i<bodyIO_.size(); i++) {
-    bodyIO_[i]->write(programState);
-  }
-
-  addEntry(relativeTimeWriter_, programState.timeStep,
-           programState.relativeTime);
-  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
-           programState.relativeTau);
-}
-
-template <class ProgramState>
-void RestartIO<ProgramState>::read(size_t timeStep,
-                                   ProgramState& programState) {
-  assert(programState.size() == bodyIO_.size());
-
-  programState.timeStep = timeStep;
-
-  for (size_t i=0; i<bodyIO_.size(); i++) {
-    bodyIO_[i]->read(timeStep, programState);
-  }
-
-  readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
-  readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
-}
-
-#include "restart-io_tmpl.cc"
diff --git a/dune/tectonic/io/hdf5/restart-io.hh b/dune/tectonic/io/hdf5/restart-io.hh
index 6b2286329ac223a738a4547d294acc2618ed84a1..5fcbe09b9d42c30259c6cc5667e4aa70a0f73bd0 100644
--- a/dune/tectonic/io/hdf5/restart-io.hh
+++ b/dune/tectonic/io/hdf5/restart-io.hh
@@ -1,41 +1,62 @@
-#ifndef SRC_HDF_RESTART_HDF_HH
-#define SRC_HDF_RESTART_HDF_HH
+#ifndef SRC_IO_HDF_RESTART_IO_HH
+#define SRC_IO_HDF_RESTART_IO_HH
+
+#include <vector>
 
 #include <dune/fufem/hdf5/file.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>
 
-template <class ProgramState> class RestartBodyIO {
-public:
-  RestartBodyIO(HDF5::Grouplike& group, const size_t vertexCount, const size_t id);
-
-  void write(const ProgramState& programState);
-
-  void read(size_t timeStep, ProgramState& programState);
-
-private:
-  const size_t id_;
-
-  HDF5::SequenceIO<2> displacementWriter_;
-  HDF5::SequenceIO<2> velocityWriter_;
-  HDF5::SequenceIO<2> accelerationWriter_;
-  HDF5::SequenceIO<1> stateWriter_;
-  HDF5::SequenceIO<1> weightedNormalStressWriter_;
-};
+#include "restartbody-io.hh"
 
 template <class ProgramState> class RestartIO {
-
-public:
-  RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts);
-  ~RestartIO();
-
-  void write(const ProgramState& programState);
-
-  void read(size_t timeStep, ProgramState& programState);
-
 private:
   std::vector<RestartBodyIO<ProgramState>* > bodyIO_;
 
   HDF5::SequenceIO<0> relativeTimeWriter_;
   HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
+
+public:
+  RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts)
+      : bodyIO_(vertexCounts.size()),
+        relativeTimeWriter_(group, "relativeTime"),
+        relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {
+
+    for (size_t i=0; i<bodyIO_.size(); i++) {
+      bodyIO_[i] = new RestartBodyIO<ProgramState>(group, vertexCounts[i], i);
+    }
+  }
+
+  ~RestartIO() {
+    for (size_t i=0; i<bodyIO_.size(); i++) {
+      delete bodyIO_[i];
+    }
+  }
+
+  void write(ProgramState const &programState) {
+    assert(programState.size() == bodyIO_.size());
+
+    for (size_t i=0; i<bodyIO_.size(); i++) {
+      bodyIO_[i]->write(programState);
+    }
+
+    addEntry(relativeTimeWriter_, programState.timeStep,
+             programState.relativeTime);
+    addEntry(relativeTimeIncrementWriter_, programState.timeStep,
+             programState.relativeTau);
+  }
+
+  void read(size_t timeStep, ProgramState& programState) {
+    assert(programState.size() == bodyIO_.size());
+
+    programState.timeStep = timeStep;
+
+    for (size_t i=0; i<bodyIO_.size(); i++) {
+      bodyIO_[i]->read(timeStep, programState);
+    }
+
+    readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
+    readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
+  }
 };
+
 #endif
diff --git a/dune/tectonic/io/hdf5/restart-io_tmpl.cc b/dune/tectonic/io/hdf5/restart-io_tmpl.cc
deleted file mode 100644
index dadc25014587373b3ba8d6a0962909958a259b2c..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/restart-io_tmpl.cc
+++ /dev/null
@@ -1,8 +0,0 @@
-#include "../../explicitvectors.hh"
-
-#include "../../data-structures/program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class RestartBodyIO<MyProgramState>;
-template class RestartIO<MyProgramState>;
diff --git a/dune/tectonic/io/hdf5/restartbody-io.hh b/dune/tectonic/io/hdf5/restartbody-io.hh
new file mode 100644
index 0000000000000000000000000000000000000000..76bfae1b9223435e2367fdb1981527a7ee8eb097
--- /dev/null
+++ b/dune/tectonic/io/hdf5/restartbody-io.hh
@@ -0,0 +1,52 @@
+#ifndef SRC_HDF_RESTART_BODY_HDF_HH
+#define SRC_HDF_RESTART_BODY_HDF_HH
+
+#include <vector>
+
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+template <class ProgramState> class RestartBodyIO {
+private:
+  const size_t id_;
+
+  HDF5::SequenceIO<2> displacementWriter_;
+  HDF5::SequenceIO<2> velocityWriter_;
+  HDF5::SequenceIO<2> accelerationWriter_;
+  HDF5::SequenceIO<1> stateWriter_;
+  HDF5::SequenceIO<1> weightedNormalStressWriter_;
+
+public:
+RestartBodyIO(HDF5::Grouplike &group, const size_t vertexCount, const size_t id)
+    : id_(id),
+      displacementWriter_(group, "displacement"+std::to_string(id_), vertexCount,
+                          ProgramState::Vector::block_type::dimension),
+      velocityWriter_(group, "velocity"+std::to_string(id_), vertexCount,
+                      ProgramState::Vector::block_type::dimension),
+      accelerationWriter_(group, "acceleration"+std::to_string(id_), vertexCount,
+                          ProgramState::Vector::block_type::dimension),
+      stateWriter_(group, "state"+std::to_string(id_), vertexCount),
+      weightedNormalStressWriter_(group, "weightedNormalStress"+std::to_string(id_), vertexCount) {}
+
+
+void write(const ProgramState& programState) {
+  addEntry(displacementWriter_, programState.timeStep, programState.u[id_]);
+  addEntry(velocityWriter_, programState.timeStep, programState.v[id_]);
+  addEntry(accelerationWriter_, programState.timeStep, programState.a[id_]);
+  addEntry(stateWriter_, programState.timeStep, programState.alpha[id_]);
+  addEntry(weightedNormalStressWriter_, programState.timeStep,
+           programState.weightedNormalStress[id_]);
+}
+
+
+void read(size_t timeStep, ProgramState& programState) {
+  readEntry(displacementWriter_, timeStep, programState.u[id_]);
+  readEntry(velocityWriter_, timeStep, programState.v[id_]);
+  readEntry(accelerationWriter_, timeStep, programState.a[id_]);
+  readEntry(stateWriter_, timeStep, programState.alpha[id_]);
+  readEntry(weightedNormalStressWriter_, timeStep,
+            programState.weightedNormalStress[id_]);
+}
+
+};
+#endif
diff --git a/dune/tectonic/io/hdf5/surface-writer.cc b/dune/tectonic/io/hdf5/surface-writer.cc
deleted file mode 100644
index 651bb7c264a2450ffc981fd958ecf5f310bb4983..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/surface-writer.cc
+++ /dev/null
@@ -1,37 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "restrict.hh"
-#include "surface-writer.hh"
-
-template <class ProgramState, class GridView>
-SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
-    HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface, size_t const id)
-    : id_(id),
-      group_(file, "surface"+std::to_string(id_)),
-      surface_(surface),
-      surfaceDisplacementWriter_(group_, "displacement", surface.numVertices(),
-                                 Vector::block_type::dimension),
-      surfaceVelocityWriter_(group_, "velocity", surface.numVertices(),
-                             Vector::block_type::dimension) {
-  auto const surfaceCoordinates = restrictToSurface(vertexCoordinates, surface);
-  HDF5::SingletonWriter<2> surfaceCoordinateWriter(group_, "coordinates",
-                                                 surfaceCoordinates.size(),
-                                                 Vector::block_type::dimension);
-  setEntry(surfaceCoordinateWriter, surfaceCoordinates);
-}
-
-template <class ProgramState, class GridView>
-void SurfaceWriter<ProgramState, GridView>::write(
-    ProgramState const &programState) {
-
-  auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_);
-  addEntry(surfaceDisplacementWriter_, programState.timeStep,
-           surfaceDisplacements);
-
-  auto const surfaceVelocities = restrictToSurface(programState.v[id_], surface_);
-  addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
-}
-
-#include "surface-writer_tmpl.cc"
diff --git a/dune/tectonic/io/hdf5/surface-writer.hh b/dune/tectonic/io/hdf5/surface-writer.hh
index b4b47a781c94a246ffc2edda9216769515567e66..15d70d01aa4b62bb09c54894240dbf9129809634 100644
--- a/dune/tectonic/io/hdf5/surface-writer.hh
+++ b/dune/tectonic/io/hdf5/surface-writer.hh
@@ -1,8 +1,10 @@
-#ifndef SRC_HDF_SURFACE_WRITER_HH
-#define SRC_HDF_SURFACE_WRITER_HH
+#ifndef SRC_HDF5_SURFACE_WRITER_HH
+#define SRC_HDF5_SURFACE_WRITER_HH
 
 #include <dune/fufem/boundarypatch.hh>
+#include "restrict.hh"
 
+#include <dune/fufem/hdf5/file.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>
 #include <dune/fufem/hdf5/singletonwriter.hh>
 
@@ -12,9 +14,31 @@ template <class ProgramState, class GridView> class SurfaceWriter {
 
 public:
   SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                Patch const &surface, size_t const id);
-
-  void write(ProgramState const &programState);
+                Patch const &surface, size_t const id) :
+      id_(id),
+      group_(file, "surface"+std::to_string(id_)),
+      surface_(surface),
+      surfaceDisplacementWriter_(group_, "displacement", surface.numVertices(),
+                                 Vector::block_type::dimension),
+      surfaceVelocityWriter_(group_, "velocity", surface.numVertices(),
+                             Vector::block_type::dimension) {
+
+  auto const surfaceCoordinates = restrictToSurface(vertexCoordinates, surface);
+  HDF5::SingletonWriter<2> surfaceCoordinateWriter(group_, "coordinates",
+                                                 surfaceCoordinates.size(),
+                                                 Vector::block_type::dimension);
+  setEntry(surfaceCoordinateWriter, surfaceCoordinates);
+}
+
+  void write(ProgramState const &programState) {
+
+      auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_);
+      addEntry(surfaceDisplacementWriter_, programState.timeStep,
+               surfaceDisplacements);
+
+      auto const surfaceVelocities = restrictToSurface(programState.v[id_], surface_);
+      addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
+    }
 
 private:
   size_t const id_;
diff --git a/dune/tectonic/io/hdf5/surface-writer_tmpl.cc b/dune/tectonic/io/hdf5/surface-writer_tmpl.cc
deleted file mode 100644
index 9d0b3b05262ad4ceaea0904fec4f64830ba8eae1..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/surface-writer_tmpl.cc
+++ /dev/null
@@ -1,8 +0,0 @@
-#include "../../explicitvectors.hh"
-#include "../../explicitgrid.hh"
-
-#include "../../data-structures/program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class SurfaceWriter<MyProgramState, DeformedGrid::LeafGridView>;
diff --git a/dune/tectonic/io/hdf5/time-writer.cc b/dune/tectonic/io/hdf5/time-writer.cc
deleted file mode 100644
index 56e2c9325e89896454b6a60034ffddf11d449faa..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/time-writer.cc
+++ /dev/null
@@ -1,21 +0,0 @@
-#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"
diff --git a/dune/tectonic/io/hdf5/time-writer.hh b/dune/tectonic/io/hdf5/time-writer.hh
index 02c2c73bcd530bb2087bdd86d6756a3bee237cb9..a0959b80e5e1a98a6fd583805d022bae8a617613 100644
--- a/dune/tectonic/io/hdf5/time-writer.hh
+++ b/dune/tectonic/io/hdf5/time-writer.hh
@@ -1,16 +1,26 @@
-#ifndef SRC_HDF_TIME_WRITER_HH
-#define SRC_HDF_TIME_WRITER_HH
+#ifndef DUNE_TECTONIC_HDF5_TIME_WRITER_HH
+#define DUNE_TECTONIC_HDF5_TIME_WRITER_HH
 
 #include <dune/fufem/hdf5/file.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>
 
-template <class ProgramState> class TimeWriter {
+template <class ProgramState>
+class TimeWriter {
 public:
-  TimeWriter(HDF5::Grouplike &file);
-  void write(ProgramState const &programState);
+  TimeWriter(HDF5::Grouplike &file) : file_(file),
+      relativeTimeWriter_(file_, "relativeTime"),
+      relativeTimeIncrementWriter_(file_, "relativeTimeIncrement") {}
+
+  void write(ProgramState const &programState)  {
+      addEntry(relativeTimeWriter_, programState.timeStep,
+               programState.relativeTime);
+      addEntry(relativeTimeIncrementWriter_, programState.timeStep,
+               programState.relativeTau);
+  }
 
 private:
   HDF5::Grouplike &file_;
+
   HDF5::SequenceIO<0> relativeTimeWriter_;
   HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
 };
diff --git a/dune/tectonic/io/hdf5/time-writer_tmpl.cc b/dune/tectonic/io/hdf5/time-writer_tmpl.cc
deleted file mode 100644
index eb4486f8e2311c2fc2ee9284f4cdf6993898cf0f..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/hdf5/time-writer_tmpl.cc
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "../../explicitvectors.hh"
-
-#include "../../data-structures/program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class TimeWriter<MyProgramState>;
diff --git a/dune/tectonic/io/io-handler.hh b/dune/tectonic/io/io-handler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6f38724678fc0faa18daf9884bd7ca6262d22d86
--- /dev/null
+++ b/dune/tectonic/io/io-handler.hh
@@ -0,0 +1,170 @@
+#ifndef SRC_IO_HANDLER_HH
+#define SRC_IO_HANDLER_HH
+
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+//#include <dune/tectonic/explicitgrid.hh>
+
+#include <dune/fufem/hdf5/file.hh>
+
+#include "../time-stepping/adaptivetimestepper.hh"
+
+#include "hdf5/restart-io.hh"
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/time-writer.hh"
+#include "hdf5-bodywriter.hh"
+
+#include "hdf5-writer.hh"
+#include "vtk.hh"
+
+template <class Assembler, class ContactNetwork>
+class IOHandler {
+private:
+    using Vector = typename Assembler::Vector;
+    using LocalVector = typename Vector::block_type;
+
+    using MyVertexBasis = typename Assembler::VertexBasis;
+    using MyCellBasis = typename Assembler::CellBasis;
+
+    const size_t bodyCount_;
+    std::vector<size_t> nVertices_;
+
+    std::vector<Vector> vertexCoordinates_;
+    std::vector<const MyVertexBasis* > vertexBases_;
+    std::vector<const MyCellBasis* > cellBases_;
+    std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches_;
+
+    bool writeVTK_;
+    bool writeData_;
+
+    bool writeRestarts_;
+    size_t restartIdx_;
+    size_t restartSpacing_;
+
+    bool printProgress_;
+
+public:
+    IOHandler(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork)
+        : bodyCount_(contactNetwork.nBodies()),
+          nVertices_(bodyCount_),
+          vertexCoordinates_(bodyCount_),
+          vertexBases_(bodyCount_),
+          cellBases_(bodyCount_),
+          weakPatches_(bodyCount_),
+          writeVTK_(parset.get<bool>("vtk.write")),
+          writeData_(parset.get<bool>("data.write")),
+          writeRestarts_(parset.get<bool>("restarts.write")),
+          restartIdx_(parset.get<size_t>("restarts.first")),
+          restartSpacing_(parset.get<size_t>("restarts.spacing")),
+          printProgress_(parset.get<bool>("printProgress")) {
+
+        for (size_t i=0; i<bodyCount_; i++) {
+            nVertices_[i] = contactNetwork.body(i)->nVertices();
+        }
+
+        for (size_t i=0; i<bodyCount_; i++) {
+          const auto& body = contactNetwork.body(i);
+          vertexBases_[i] = &(body->assembler()->vertexBasis);
+          cellBases_[i] = &(body->assembler()->cellBasis);
+
+          auto& vertexCoords = vertexCoordinates_[i];
+          vertexCoords.resize(nVertices_[i]);
+
+          auto hostLeafView = body->grid()->hostGrid().leafGridView();
+          Dune::MultipleCodimMultipleGeomTypeMapper<
+              decltype(hostLeafView), Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+          for (auto &&v : vertices(hostLeafView))
+            vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+        }
+    }
+
+    template <class ProgramState, class GlobalFriction>
+    void write(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial = false) {
+      if (writeData_) {
+        writeData(programState, contactNetwork, friction, iterationCount, initial);
+      }
+
+      if (writeVTK_) {
+        writeVTK(programState, contactNetwork);
+      }
+
+      if (writeRestarts_) {
+        writeRestarts(programState);
+      }
+
+      if (printProgress_) {
+        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
+                  << ", time = " << std::setw(12) << programState.relativeTime
+                  << ", tau = " << std::setw(12) << programState.relativeTau
+                  << std::endl;
+      }
+    }
+
+    template <class ProgramState>
+    bool read(ProgramState& programState) {
+        using ADS = RestartIO<ProgramState>;
+
+        if (restartIdx_ > 0) {// automatically adjusts the time and timestep
+            auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READONLY);
+            auto restartIO = std::make_unique<ADS>(*restartFile, nVertices_);
+            restartIO->read(restartIdx_, programState);
+
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+private:
+    template <class ProgramState>
+    void writeVTK(const ProgramState& programState, const ContactNetwork& contactNetwork) {
+        using ScalarVector = typename Assembler::ScalarVector;
+        std::vector<ScalarVector> stress(bodyCount_);
+
+        for (size_t i=0; i<bodyCount_; i++) {
+          const auto& body = contactNetwork.body(i);
+          body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+                                           body->data()->getPoissonRatio(),
+                                           programState.u[i], stress[i]);
+
+        }
+
+        const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases_, vertexBases_, "../debug_print/vtk/");
+        vtkWriter.write(programState.timeStep, programState.u, programState.v,
+                        programState.alpha, stress);
+    }
+
+    template <class ProgramState>
+    void writeRestarts(const ProgramState& programState) {
+        if (programState.timeStep % restartSpacing_ == 0) {
+          auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READWRITE);
+          auto restartIO = std::make_unique<RestartIO<ProgramState>>(*restartFile, nVertices_);
+
+          restartIO->write(programState);
+          restartFile->flush();
+        }
+    }
+
+    template <class ProgramState, class GlobalFriction>
+    void writeData(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial) {
+      std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
+      contactNetwork.frictionPatches(frictionPatches);
+
+      auto dataFile = std::make_unique<HDF5::File>("output.h5");
+      auto dataWriter = std::make_unique<
+                            HDF5Writer<ProgramState, MyVertexBasis, typename Assembler::GV>>(
+                            *dataFile, vertexCoordinates_, vertexBases_,
+                            frictionPatches);
+
+      dataWriter->reportSolution(programState, contactNetwork, friction);
+      if (!initial)
+        dataWriter->reportIterations(programState, iterationCount);
+      dataFile->flush();
+    }
+
+};
+#endif
diff --git a/dune/tectonic/io/uniform-grid-writer.cc b/dune/tectonic/io/uniform-grid-writer.cc
index cf7d466d9e182e6be020581de737fd616b6f8385..6ad279cc5c39464d167d0d01bbe703a773ebac3c 100644
--- a/dune/tectonic/io/uniform-grid-writer.cc
+++ b/dune/tectonic/io/uniform-grid-writer.cc
@@ -11,13 +11,13 @@
 // #include <dune/common/parametertreeparser.hh>
 
 #include "../assemblers.hh"
-#include "gridselector.hh"
+#include "../gridselector.hh"
 
-#include "io/vtk.hh"
+#include "vtk.hh"
 
-#include "one-body-problem-data/mygrid.hh"
+#include "../problem-data/grid/gridconstructor.hh"
 
-#include "utils/diameter.hh"
+#include "../utils/diameter.hh"
 
 size_t const dims = MY_DIM;
 size_t const refinements = 5; // FIXME?
diff --git a/dune/tectonic/io/vtk-bodywriter.hh b/dune/tectonic/io/vtk-bodywriter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8ddc3cf014d13debf50533bfe53ff659b37d8e1c
--- /dev/null
+++ b/dune/tectonic/io/vtk-bodywriter.hh
@@ -0,0 +1,58 @@
+#ifndef SRC_VTK_BODYWRITER_HH
+#define SRC_VTK_BODYWRITER_HH
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/fufem/functions/vtkbasisgridfunction.hh>
+#include "../utils/debugutils.hh"
+
+template <class VertexBasis, class CellBasis> class VTKBodyWriter {
+private:
+  CellBasis const &cellBasis_;
+  VertexBasis const &vertexBasis_;
+  std::string const prefix_;
+
+public:
+  VTKBodyWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis, std::string prefix) :
+    cellBasis_(cellBasis), vertexBasis_(vertexBasis), prefix_(prefix)
+  {}
+
+  template <class Vector, class ScalarVector>
+  void write(size_t record, Vector const &u, Vector const &v,
+             ScalarVector const &alpha, ScalarVector const &stress) const {
+      Dune::VTKWriter<typename VertexBasis::GridView> writer(
+          vertexBasis_.getGridView());
+
+      auto const displacementPointer =
+          std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
+              vertexBasis_, u, "displacement");
+      writer.addVertexData(displacementPointer);
+
+      auto const velocityPointer =
+          std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
+              vertexBasis_, v, "velocity");
+      writer.addVertexData(velocityPointer);
+
+      auto const AlphaPointer =
+          std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
+              vertexBasis_, alpha, "Alpha");
+      writer.addVertexData(AlphaPointer);
+
+      auto const stressPointer =
+          std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(
+              cellBasis_, stress, "stress");
+      writer.addCellData(stressPointer);
+
+      std::string const filename = prefix_ + "_" + std::to_string(record);
+      writer.write(filename.c_str(), Dune::VTK::appendedraw);
+  }
+
+  void writeGrid() const {
+      Dune::VTKWriter<typename VertexBasis::GridView> writer(
+          vertexBasis_.getGridView());
+
+      std::string const filename = prefix_ + "_grid";
+      writer.write(filename.c_str(), Dune::VTK::appendedraw);
+  }
+};
+
+#endif
diff --git a/dune/tectonic/io/vtk.cc b/dune/tectonic/io/vtk.cc
deleted file mode 100644
index 5eaa52ba6641257c1a88b359a4cf2ca1e0afd39e..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/vtk.cc
+++ /dev/null
@@ -1,91 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
-
-#include <dune/fufem/functions/vtkbasisgridfunction.hh>
-
-#include "vtk.hh"
-
-#include "../utils/debugutils.hh"
-
-template <class VertexBasis, class CellBasis>
-BodyVTKWriter<VertexBasis, CellBasis>::BodyVTKWriter(
-    CellBasis const & cellBasis, VertexBasis const & vertexBasis,
-    std::string prefix)
-    : cellBasis_(cellBasis), vertexBasis_(vertexBasis), prefix_(prefix) {}
-
-template <class VertexBasis, class CellBasis>
-template <class Vector, class ScalarVector>
-void BodyVTKWriter<VertexBasis, CellBasis>::write(
-    size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha,
-    ScalarVector const &stress) const {
-  Dune::VTKWriter<typename VertexBasis::GridView> writer(
-      vertexBasis_.getGridView());
-
-  auto const displacementPointer =
-      std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
-          vertexBasis_, u, "displacement");
-  writer.addVertexData(displacementPointer);
-
-  auto const velocityPointer =
-      std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
-          vertexBasis_, v, "velocity");
-  writer.addVertexData(velocityPointer);
-
-  auto const AlphaPointer =
-      std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
-          vertexBasis_, alpha, "Alpha");
-  writer.addVertexData(AlphaPointer);
-
-  auto const stressPointer =
-      std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(
-          cellBasis_, stress, "stress");
-  writer.addCellData(stressPointer);
-
-  std::string const filename = prefix_ + "_" + std::to_string(record);
-  writer.write(filename.c_str(), Dune::VTK::appendedraw);
-}
-
-template <class VertexBasis, class CellBasis>
-void BodyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
-  Dune::VTKWriter<typename VertexBasis::GridView> writer(
-      vertexBasis_.getGridView());
-
-  std::string const filename = prefix_ + "_grid";
-  writer.write(filename.c_str(), Dune::VTK::appendedraw);
-}
-
-template <class VertexBasis, class CellBasis>
-MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
-  std::string prefix): bodyVTKWriters_(vertexBases.size()), prefix_(prefix) {
-
-  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
-    bodyVTKWriters_[i] = new BodyVTKWriter<VertexBasis, CellBasis>(*cellBases[i], *vertexBases[i], prefix_+std::to_string(i));
-  }
-}
-
-template <class VertexBasis, class CellBasis>
-MyVTKWriter<VertexBasis, CellBasis>::~MyVTKWriter() {
-  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
-    delete bodyVTKWriters_[i];
-  }
-}
-
-template <class VertexBasis, class CellBasis>
-template <class Vector, class ScalarVector>
-void MyVTKWriter<VertexBasis, CellBasis>::write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
-                        const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const {
-
-  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
-    bodyVTKWriters_[i]->write(record, u[i], v[i], alpha[i], stress[i]);
-  }
-}
-template <class VertexBasis, class CellBasis>
-void MyVTKWriter<VertexBasis, CellBasis>::writeGrids() const {
-  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
-    bodyVTKWriters_[i]->writeGrid();
-  }
-}
-#include "vtk_tmpl.cc"
diff --git a/dune/tectonic/io/vtk.hh b/dune/tectonic/io/vtk.hh
index ee21f8993b324aba1f83e27743b167b5ed7395e2..5634e4ef27c7ca6490dd49d12f660b0d300acca9 100644
--- a/dune/tectonic/io/vtk.hh
+++ b/dune/tectonic/io/vtk.hh
@@ -3,38 +3,47 @@
 
 #include <string>
 
-template <class VertexBasis, class CellBasis> class BodyVTKWriter {
-private:
-  CellBasis const &cellBasis_;
-  VertexBasis const &vertexBasis_;
-  std::string const prefix_;
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
 
-public:
-  BodyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis,
-              std::string prefix);
+#include <dune/fufem/functions/vtkbasisgridfunction.hh>
 
-  template <class Vector, class ScalarVector>
-  void write(size_t record, Vector const &u, Vector const &v,
-             ScalarVector const &alpha, ScalarVector const &stress) const;
-
-  void writeGrid() const;
-};
+#include "../utils/debugutils.hh"
+#include "vtk-bodywriter.hh"
 
 template <class VertexBasis, class CellBasis> class MyVTKWriter {
 private:
-  std::vector<BodyVTKWriter<VertexBasis, CellBasis>* > bodyVTKWriters_;
+  std::vector<VTKBodyWriter<VertexBasis, CellBasis>* > bodyVTKWriters_;
   std::string const prefix_;
 
 public:
   MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
-              std::string prefix);
+              std::string prefix): bodyVTKWriters_(vertexBases.size()), prefix_(prefix) {
+
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        bodyVTKWriters_[i] = new VTKBodyWriter<VertexBasis, CellBasis>(*cellBases[i], *vertexBases[i], prefix_+std::to_string(i));
+      }
+  }
 
-  ~MyVTKWriter();
+  ~MyVTKWriter() {
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        delete bodyVTKWriters_[i];
+      }
+  }
 
   template <class Vector, class ScalarVector>
   void write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
-             const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;
-
-  void writeGrids() const;
+             const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const {
+
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        bodyVTKWriters_[i]->write(record, u[i], v[i], alpha[i], stress[i]);
+      }
+    }
+
+  void writeGrids() const {
+      for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+        bodyVTKWriters_[i]->writeGrid();
+      }
+  }
 };
+
 #endif
diff --git a/dune/tectonic/io/vtk_tmpl.cc b/dune/tectonic/io/vtk_tmpl.cc
deleted file mode 100644
index d1d9f0a77484660bddac03480fc83a8aeb6b7a1c..0000000000000000000000000000000000000000
--- a/dune/tectonic/io/vtk_tmpl.cc
+++ /dev/null
@@ -1,18 +0,0 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-#include <dune/fufem/functionspacebases/p0basis.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
-using MyP0Basis = P0Basis<DeformedGrid::LeafGridView, double>;
-using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
-
-template class MyVTKWriter<P1Basis, MyP0Basis>;
-
-template void MyVTKWriter<P1Basis, MyP0Basis>::write<Vector, ScalarVector>(
-    size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
-        const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;
diff --git a/dune/tectonic/multi-threading/CMakeLists.txt b/dune/tectonic/multi-threading/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4ad03636ec336e1e0b21f21d48449d70c5d764bd
--- /dev/null
+++ b/dune/tectonic/multi-threading/CMakeLists.txt
@@ -0,0 +1,11 @@
+
+add_custom_target(tectonic_dune_multi-threading SOURCES
+  stoppable.hh
+  task.hh 
+)
+
+#install headers
+install(FILES
+  stoppable.hh
+  task.hh 
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/multi-threading/stoppable.hh b/dune/tectonic/multi-threading/stoppable.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/dune/tectonic/multi-threading/task.hh b/dune/tectonic/multi-threading/task.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d90f72761f13c6b569e3c93410e66fba110a7c84
--- /dev/null
+++ b/dune/tectonic/multi-threading/task.hh
@@ -0,0 +1,49 @@
+#include <future>
+#include <utility>
+
+/*
+ * Class that encapsulates promise and future object,
+ * provides API to set exit signal for thread
+ */
+class Task {
+protected:
+    std::promise<void> exitSignal;
+    std::future<void> futureObj;
+
+public:
+    Task() :
+        futureObj(exitSignal.get_future())
+    {}
+
+    Task(Task&& obj) :
+        exitSignal(std::move(obj.exitSignal)), futureObj(std::move(obj.futureObj))
+    {}
+
+    Task& operator=(Task&& obj) {
+        exitSignal = std::move(obj.exitSignal);
+        futureObj = std::move(obj.futureObj);
+        return *this;
+    }
+
+    virtual void run_task() = 0;
+
+    // thread function to be executed by thread
+    void operator()() {
+        run_task();
+    }
+
+    bool stopRequested() {
+        // checks if value in future object is available
+        if (futureObj.wait_for(std::chrono::milliseconds(0)) == std::future_status::timeout)
+            return false;
+        return true;
+    }
+
+    // stop thread by setting value in promise object
+    void stop() {
+        exitSignal.set_value();
+    }
+};
+
+
+
diff --git a/dune/tectonic/problem-data/grid/gridconstructor.hh b/dune/tectonic/problem-data/grid/gridconstructor.hh
index 4faa6fe1254c0745eff55850217d562b1e4ceb7b..e93d3c0131471db3d6fe29cf49c5ecc976ed470d 100644
--- a/dune/tectonic/problem-data/grid/gridconstructor.hh
+++ b/dune/tectonic/problem-data/grid/gridconstructor.hh
@@ -6,7 +6,7 @@
 #include <dune/fufem/geometry/convexpolyhedron.hh>
 #include <dune/fufem/geometry/polyhedrondistance.hh>
 
-#include "../utils/diameter.hh"
+#include "../../utils/diameter.hh"
 
 template <class GridType>
 class GridConstructor {
diff --git a/dune/tectonic/tests/contacttest/CMakeLists.txt b/dune/tectonic/tests/contacttest/CMakeLists.txt
index fc6644b88b6e9c7a81c65fef660fd431dc6129c5..1b36825411e9d918394b7eb263232d53471105e7 100644
--- a/dune/tectonic/tests/contacttest/CMakeLists.txt
+++ b/dune/tectonic/tests/contacttest/CMakeLists.txt
@@ -13,7 +13,7 @@ set(CONTACTTEST_SOURCE_FILES
   ../../factories/twoblocksfactory.cc
   ../../factories/threeblocksfactory.cc
   ../../factories/stackedblocksfactory.cc
-  ../../io/vtk.cc
+  #../../io/vtk.cc
   ../../problem-data/grid/cuboidgeometry.cc
   ../../problem-data/grid/mygrids.cc
   ../../problem-data/grid/simplexmanager.cc
diff --git a/dune/tectonic/tests/hdf5test/CMakeLists.txt b/dune/tectonic/tests/hdf5test/CMakeLists.txt
index 27378cc4d6a359ccc7e42502545c8b99db104732..0676da18247d5ae16f704c957d3c4a481d106063 100644
--- a/dune/tectonic/tests/hdf5test/CMakeLists.txt
+++ b/dune/tectonic/tests/hdf5test/CMakeLists.txt
@@ -10,10 +10,10 @@ set(HDF5TEST_SOURCE_FILES
   ../../data-structures/network/contactnetwork.cc
   ../../data-structures/enumparser.cc
   ../../factories/strikeslipfactory.cc
-  ../../io/vtk.cc
-  ../../io/hdf5/frictionalboundary-writer.cc
-  ../../io/hdf5/iteration-writer.cc
-  ../../io/hdf5/time-writer.cc
+  #../../io/vtk.cc
+  #../../io/hdf5/frictionalboundary-writer.cc
+  #../../io/hdf5/iteration-writer.cc
+  #../../io/hdf5/time-writer.cc
   ../../problem-data/grid/simplexmanager.cc
   ../../problem-data/grid/trianglegeometry.cc
   ../../problem-data/grid/trianglegrids.cc
diff --git a/dune/tectonic/time-stepping/adaptivetimestepper.cc b/dune/tectonic/time-stepping/adaptivetimestepper.cc
index 4c5b9185a232747da2b1d1e02c9ce2053bd903ae..f11a7673bcbf95086d58bd140dbed03486968806 100644
--- a/dune/tectonic/time-stepping/adaptivetimestepper.cc
+++ b/dune/tectonic/time-stepping/adaptivetimestepper.cc
@@ -2,7 +2,10 @@
 #include "config.h"
 #endif
 
+#include <thread>
+
 #include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
 #include <dune/solvers/iterationsteps/cgstep.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
@@ -12,6 +15,8 @@
 
 #include "adaptivetimestepper.hh"
 
+const unsigned int N_THREADS = std::thread::hardware_concurrency();
+
 template<class IterationStepType, class NormType, class ReductionFactorContainer>
 Dune::Solvers::Criterion reductionFactorCriterion(
       IterationStepType& iterationStep,
@@ -221,13 +226,54 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
    linearSolver->addCriterion(reductionFactorCriterion(*linearMultigridStep, norm, reductionFactors));
    //linearSolver->addCriterion(reductionFactorCriterion(*cgStep, norm, reductionFactors));
 
-  if (R1_.updaters == Updaters()) {
-    //setDeformation(current_);
-    R1_ = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, relativeTau_);
-    updateReductionFactors(reductionFactors);
-  }
+  auto refine = [&](UpdatersWithCount& F1, UpdatersWithCount& F2){
+      setDeformation(current_);
+      F1 = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, relativeTau_ / 2.0);
+      updateReductionFactors(reductionFactors);
+      std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
+
+      setDeformation(F1.updaters);
+      auto&& nBodyAssembler = step(currentNBodyAssembler);
+      F2 = step(F1.updaters, nBodyAssembler, linearSolver, relativeTime_ + relativeTau_ / 2.0,
+                relativeTau_ / 2.0);
+      updateReductionFactors(reductionFactors);
+  };
+
+  auto coarsen = [&](UpdatersWithCount& C, UpdatersWithCount& R2){
+      setDeformation(current_);
+      C = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, 2 * relativeTau_).get();
+      updateReductionFactors(reductionFactors);
+      std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
+
+      /*using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
+      std::vector<ScalarVector> cAlpha(contactNetwork_.nBodies());
+      C.updaters.state_->extractAlpha(cAlpha);
+      print(cAlpha, "cAlpha: ");*/
+
+      setDeformation(R1_.updaters);
+      auto&& nBodyAssembler = step(currentNBodyAssembler);
+      R2 = step(R1_.updaters, nBodyAssembler, linearSolver, relativeTime_ + relativeTau_, relativeTau_).get();
+      updateReductionFactors(reductionFactors);
+      std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
+  };
+
+  auto canCoarsen = [&]{
+    std::cout << N_THREADS << " concurrent threads are supported." << std::endl;
 
+    if (R1_.updaters == Updaters()) {
+      //setDeformation(current_);
+      R1_ = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, relativeTau_);
+      updateReductionFactors(reductionFactors);
+    }
+
+    if (N_THREADS<3) {
+
+    } else {
+
+    }
 
+    return ;
+  };
 
   //std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
 
@@ -242,28 +288,12 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
   while (relativeTime_ + relativeTau_ <= 1.0) {
     std::cout << "tau: " << relativeTau_ << std::endl;
 
-    setDeformation(current_);
-    C = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, 2 * relativeTau_);
-    updateReductionFactors(reductionFactors);
-    std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
-
-    /*using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
-    std::vector<ScalarVector> cAlpha(contactNetwork_.nBodies());
-    C.updaters.state_->extractAlpha(cAlpha);
-    print(cAlpha, "cAlpha: ");*/
-
-    setDeformation(R1_.updaters);
-    auto&& nBodyAssembler = step(currentNBodyAssembler);
-    R2 = step(R1_.updaters, nBodyAssembler, linearSolver, relativeTime_ + relativeTau_, relativeTau_);
-    updateReductionFactors(reductionFactors);
-    std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
-
+    coarsen(C, R2);
 
     /*std::vector<ScalarVector> rAlpha(contactNetwork_.nBodies());
     R2.updaters.state_->extractAlpha(rAlpha);
     print(rAlpha, "rAlpha: ");*/
 
-
     if (mustRefine_(C.updaters, R2.updaters))
       break;
 
@@ -278,17 +308,8 @@ IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNo
   UpdatersWithCount F2;
   if (!didCoarsen) {
     while (true) {
-      setDeformation(current_);
-      F1 = step(current_, currentNBodyAssembler, linearSolver, relativeTime_, relativeTau_ / 2.0);
-      updateReductionFactors(reductionFactors);
-      std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
+      refine(F1, F2);
 
-      setDeformation(F1.updaters);
-      auto&& nBodyAssembler = step(currentNBodyAssembler);
-      F2 = step(F1.updaters, nBodyAssembler, linearSolver, relativeTime_ + relativeTau_ / 2.0,
-                relativeTau_ / 2.0);
-      updateReductionFactors(reductionFactors);
-      std::cout << "AdaptiveTimeStepper F2 computed!" << std::endl << std::endl;
       if (!mustRefine_(R1_.updaters, F2.updaters)) {
         std::cout << "Sufficiently refined!" << std::endl;
         break;
@@ -347,21 +368,26 @@ typename AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::Upd
 AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::step(
     const Updaters& oldUpdaters, const NBodyAssembler& nBodyAssembler, std::shared_ptr<LinearSolver>& linearSolver, double rTime, double rTau) {
 
-  UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
+  auto doStep = [&]{
+    UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
 
-  MyCoupledTimeStepper coupledTimeStepper(finalTime_, parset_, nBodyAssembler,
+    MyCoupledTimeStepper coupledTimeStepper(finalTime_, parset_, nBodyAssembler,
                                           ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
                                           newUpdatersAndCount.updaters, errorNorms_, externalForces_);
 
-  newUpdatersAndCount.count = coupledTimeStepper.step(linearSolver, rTime, rTau);
-  iterationRegister_.registerCount(newUpdatersAndCount.count);
+    newUpdatersAndCount.count = coupledTimeStepper.step(linearSolver, rTime, rTau);
+    iterationRegister_.registerCount(newUpdatersAndCount.count);
 
+    return newUpdatersAndCount;
+  };
+
+  return doStep();
   /*using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
   std::vector<ScalarVector> alpha(contactNetwork_.nBodies());
   newUpdatersAndCount.updaters.state_->extractAlpha(alpha);
   print(alpha, "step alpha: ");
   */
-  return newUpdatersAndCount;
+
 }
 
 #include "adaptivetimestepper_tmpl.cc"
diff --git a/dune/tectonic/time-stepping/adaptivetimestepper.hh b/dune/tectonic/time-stepping/adaptivetimestepper.hh
index 391170a0f94f47dafa8a5f97be0967167cd3047a..fa00188d147a590112aa02ab511b3484c0f725bd 100644
--- a/dune/tectonic/time-stepping/adaptivetimestepper.hh
+++ b/dune/tectonic/time-stepping/adaptivetimestepper.hh
@@ -3,6 +3,7 @@
 
 #include <fstream>
 
+//#include "../multi-threading/task.hh"
 #include "coupledtimestepper.hh"
 
 struct IterationRegister {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 080ed3581c0d5814bed561a24fd69c15ac1c473b..414d9ed805a59ea1350083f637d0c87809e0264a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,17 +4,17 @@ add_subdirectory("strikeslip")
 
 set(UGW_SOURCE_FILES
   ../dune/tectonic/assemblers.cc # FIXME
-  ../dune/tectonic/io/uniform-grid-writer.cc
+  #../dune/tectonic/io/uniform-grid-writer.cc
   ../dune/tectonic/io/vtk.cc
   ../dune/tectonic/problem-data/grid/mygrids.cc
 )
 
 foreach(_dim 2 3)
-  set(_ugw_target uniform-grid-writer-${_dim}D)
+  #set(_ugw_target uniform-grid-writer-${_dim}D)
 
-  add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+  #add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
   
-  add_dune_ug_flags(${_ugw_target})
+  #add_dune_ug_flags(${_ugw_target})
 
-  set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()
diff --git a/src/foam/CMakeLists.txt b/src/foam/CMakeLists.txt
index ab24974cb7acf4b09d54775f2ab0f54963a0c01d..5173c53fec345a744f37c7549a79527af0d3fe95 100644
--- a/src/foam/CMakeLists.txt
+++ b/src/foam/CMakeLists.txt
@@ -10,13 +10,13 @@ set(FOAM_SOURCE_FILES
   ../../dune/tectonic/data-structures/network/contactnetwork.cc
   ../../dune/tectonic/data-structures/enumparser.cc
   ../../dune/tectonic/factories/twoblocksfactory.cc
-  ../../dune/tectonic/io/vtk.cc
-  ../../dune/tectonic/io/hdf5/frictionalboundary-writer.cc
-  ../../dune/tectonic/io/hdf5/iteration-writer.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/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
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index d78fc2709deb2b4185253797732f3cf20db11b5c..46e8bc657ba6197618dd6dfb930aeb6b48edbdf5 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -56,7 +56,7 @@
 
 #include <dune/tectonic/factories/twoblocksfactory.hh>
 
-#include <dune/tectonic/io/hdf5-levelwriter.hh>
+#include <dune/tectonic/io/hdf5-writer.hh>
 #include <dune/tectonic/io/hdf5/restart-io.hh>
 #include <dune/tectonic/io/vtk.hh>
 
@@ -249,12 +249,12 @@ int main(int argc, char *argv[]) {
 
     auto dataWriter =
         writeData ? std::make_unique<
-                        HDF5LevelWriter<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
                         *dataFile, vertexCoordinates, vertexBases,
                         frictionPatches, weakPatches)
                   : nullptr;
 
-    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
+    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/home/joscha/software/dune/dune-tectonic/body");
 
     IterationRegister iterationCount;
 
diff --git a/src/multi-body-problem/CMakeLists.txt b/src/multi-body-problem/CMakeLists.txt
index 15a48074443bbe6fd2c7eddb2b8642b430148c77..272477e6da80977fb420868e4791bd02056e0ecc 100644
--- a/src/multi-body-problem/CMakeLists.txt
+++ b/src/multi-body-problem/CMakeLists.txt
@@ -13,13 +13,13 @@ set(MSW_SOURCE_FILES
   #../../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/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/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
diff --git a/src/multi-body-problem/multi-body-problem-2D.cfg b/src/multi-body-problem/multi-body-problem-2D.cfg
index 641078a70876d7286c6691333b63f81740c6898d..d17965a838358b1ed23731805907260c5bf5075e 100644
--- a/src/multi-body-problem/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem/multi-body-problem-2D.cfg
@@ -1,9 +1,9 @@
 # -*- mode:conf -*-
 [boundary.friction]
-smallestDiameter = 0.05  # 2e-3 [m]
+smallestDiameter = 0.05  # 0.05 2e-3 [m]
 
 [timeSteps]
-refinementTolerance = 1e-5 # 1e-5
+refinementTolerance = 2e-2 # 1e-5
 
 [u0.solver]
 tolerance         = 1e-8
diff --git a/src/multi-body-problem/multi-body-problem.cc b/src/multi-body-problem/multi-body-problem.cc
index b94a7bebf15b27cb9c4d52c1e027215aa581b73f..d7a0009ab06ca48e7080e3f11680d27e474a347a 100644
--- a/src/multi-body-problem/multi-body-problem.cc
+++ b/src/multi-body-problem/multi-body-problem.cc
@@ -57,6 +57,7 @@
 #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>
@@ -92,9 +93,9 @@ std::vector<std::vector<double>> allReductionFactors;
 
 Dune::ParameterTree getParameters(int argc, char *argv[]) {
   Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem/multi-body-problem.cfg", 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/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem/multi-body-problem-%dD.cfg", dims), parset);
+      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;
 }
@@ -126,8 +127,8 @@ int main(int argc, char *argv[]) {
     // set up contact network
     // ----------------------
 
-    //using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
-    using BlocksFactory = ThreeBlocksFactory<Grid, Vector>;
+    using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
+    //using BlocksFactory = ThreeBlocksFactory<Grid, Vector>;
     BlocksFactory blocksFactory(parset);
 
     using ContactNetwork = typename BlocksFactory::ContactNetwork;
@@ -144,11 +145,11 @@ int main(int argc, char *argv[]) {
         //def = 1;
         //deformedGridComplex.setDeformation(def, i);
 
-        /*const auto& level = *contactNetwork.level(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++) {
@@ -173,33 +174,13 @@ int main(int argc, char *argv[]) {
 
     using MyProgramState = ProgramState<Vector, ScalarVector>;
     MyProgramState programState(nVertices);
-    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;
 
+    IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
 
-    auto restartIO = handleRestarts
-                         ? std::make_unique<RestartIO<MyProgramState>>(
-                               *restartFile, nVertices)
-                         : nullptr;
-
-    if (firstRestart > 0) // automatically adjusts the time and timestep
-      restartIO->read(firstRestart, programState);
-    else
-     programState.setupInitialConditions(parset, contactNetwork);
+    bool restartRead = ioHandler.read(programState);
+    if (!restartRead) {
+      programState.setupInitialConditions(parset, contactNetwork);
+    }
 
     //print(programState.u, "u:");
 
@@ -218,87 +199,10 @@ int main(int argc, char *argv[]) {
     auto& globalFriction = contactNetwork.globalFriction();
     globalFriction.updateAlpha(programState.alpha);
 
-    using MyVertexBasis = typename Assembler::VertexBasis;
-    using MyCellBasis = typename Assembler::CellBasis;
-    std::vector<Vector> vertexCoordinates(bodyCount);
-    std::vector<const MyVertexBasis* > vertexBases(bodyCount);
-    std::vector<const MyCellBasis* > cellBases(bodyCount);
-
-    auto& wPatches = blocksFactory.weakPatches();
-    std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
-
-
-    for (size_t i=0; i<bodyCount; i++) {
-      const auto& body = contactNetwork.body(i);
-      vertexBases[i] = &(body->assembler()->vertexBasis);
-      cellBases[i] = &(body->assembler()->cellBasis);
-
-      weakPatches[i].resize(1);
-      weakPatches[i][0] = wPatches[i].get();
-
-      auto& vertexCoords = vertexCoordinates[i];
-      vertexCoords.resize(nVertices[i]);
-
-      auto hostLeafView = body->grid()->hostGrid().leafGridView();
-      Dune::MultipleCodimMultipleGeomTypeMapper<
-          LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
-      for (auto &&v : vertices(hostLeafView))
-        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
-    }
-
-    typename ContactNetwork::BoundaryPatches frictionBoundaries;
-    contactNetwork.boundaryPatches("friction", frictionBoundaries);
-
-    std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
-    contactNetwork.frictionPatches(frictionPatches);
-
-    auto dataWriter =
-        writeData ? std::make_unique<
-                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
-                        *dataFile, vertexCoordinates, vertexBases,
-                        frictionPatches, weakPatches)
-                  : nullptr;
-
-    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
 
     IterationRegister iterationCount;
 
-    auto const report = [&](bool initial = false) {
-      if (writeData) {
-        dataWriter->reportSolution(programState, contactNetwork, globalFriction);
-        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")) {
-        std::vector<ScalarVector> stress(bodyCount);
-
-        for (size_t i=0; i<bodyCount; i++) {
-          const auto& body = contactNetwork.body(i);
-          body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
-                                           body->data()->getPoissonRatio(),
-                                           programState.u[i], stress[i]);
-
-        }
-
-        vtkWriter.write(programState.timeStep, programState.u, programState.v,
-                        programState.alpha, stress);
-      }
-    };
-    report(true);
+    ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, true);
 
     //DUNE_THROW(Dune::Exception, "Just need to stop here!");
 
@@ -465,7 +369,7 @@ int main(int argc, char *argv[]) {
 
       contactNetwork.setDeformation(programState.u);
 
-      report();
+      ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, false);
 
       if (programState.timeStep==timeSteps) {
         std::cout << "limit of timeSteps reached!" << std::endl;
diff --git a/src/multi-body-problem/multi-body-problem.cfg b/src/multi-body-problem/multi-body-problem.cfg
index 8a5362a1a389e38917a2a8ed8ad38963e67fe569..244aeb3bfbcd0d21824cefb96961ea6963fd128f 100644
--- a/src/multi-body-problem/multi-body-problem.cfg
+++ b/src/multi-body-problem/multi-body-problem.cfg
@@ -32,7 +32,7 @@ b               = 0.074    # [ ]
 
 
 [boundary.neumann]
-sigmaN          = 200.0      # [Pa]
+sigmaN          =  0.0     # 200.0 [Pa]
 
 
 [boundary.dirichlet]
@@ -47,20 +47,20 @@ relativeTau = 2e-4 # 1e-6
 
 [timeSteps]
 scheme = newmark
-timeSteps = 5000
+timeSteps = 5
 
 
 [problem]
 finalTime       = 100  # [s] #1000
-bodyCount       = 3
+bodyCount       = 4
 
 
 [io]
 data.write      = false
 printProgress   = true
 restarts.first  = 0
-restarts.spacing= 20
-restarts.write  = false #true
+restarts.spacing= 1 #20
+restarts.write  = true #true
 vtk.write       = true
 
 
diff --git a/src/strikeslip/CMakeLists.txt b/src/strikeslip/CMakeLists.txt
index ac71e9a9e908d56d8e6442e2f430e4287f709bc2..95dfa88f73b9d524482fba05ad01ebccf2bc527c 100644
--- a/src/strikeslip/CMakeLists.txt
+++ b/src/strikeslip/CMakeLists.txt
@@ -10,13 +10,13 @@ set(STRIKESLIP_SOURCE_FILES
   ../../dune/tectonic/data-structures/network/contactnetwork.cc
   ../../dune/tectonic/data-structures/enumparser.cc
   ../../dune/tectonic/factories/strikeslipfactory.cc
-  ../../dune/tectonic/io/vtk.cc
-  ../../dune/tectonic/io/hdf5/frictionalboundary-writer.cc
-  ../../dune/tectonic/io/hdf5/iteration-writer.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/io/hdf5/restart-io.cc
+  #../../dune/tectonic/io/hdf5/surface-writer.cc
+  #../../dune/tectonic/io/hdf5/time-writer.cc
   ../../dune/tectonic/problem-data/grid/simplexmanager.cc
   ../../dune/tectonic/problem-data/grid/trianglegeometry.cc
   ../../dune/tectonic/problem-data/grid/trianglegrids.cc
diff --git a/src/strikeslip/strikeslip.cc b/src/strikeslip/strikeslip.cc
index 412d10974efff94aa2e83010416f80245a5926a9..4cdf3f8195ae6b7f9fc2b6f5297002b4dbb7a1a6 100644
--- a/src/strikeslip/strikeslip.cc
+++ b/src/strikeslip/strikeslip.cc
@@ -90,9 +90,9 @@ size_t const dims = MY_DIM;
 
 Dune::ParameterTree getParameters(int argc, char *argv[]) {
   Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/strikeslip/strikeslip.cfg", parset);
+  Dune::ParameterTreeParser::readINITree("/home/joscha/software/dune/dune-tectonic/src/strikeslip/strikeslip.cfg", parset);
   Dune::ParameterTreeParser::readINITree(
-      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/strikeslip/strikeslip-%dD.cfg", dims), parset);
+      Dune::Fufem::formatString("/home/joscha/software/dune/dune-tectonic/src/strikeslip/strikeslip-%dD.cfg", dims), parset);
   Dune::ParameterTreeParser::readOptions(argc, argv, parset);
   return parset;
 }
@@ -130,8 +130,6 @@ int main(int argc, char *argv[]) {
 
     ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
 
-
-
     const size_t bodyCount = contactNetwork.nBodies();
 
     for (size_t i=0; i<contactNetwork.nLevels(); i++) {
@@ -175,10 +173,6 @@ int main(int argc, char *argv[]) {
     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",
@@ -246,18 +240,20 @@ int main(int argc, char *argv[]) {
     std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
     contactNetwork.frictionPatches(frictionPatches);
 
-    auto dataWriter =
-        writeData ? std::make_unique<
-                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
-                        *dataFile, vertexCoordinates, vertexBases,
-                        frictionPatches, weakPatches)
-                  : nullptr;
-
-    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
 
     IterationRegister iterationCount;
 
     auto const report = [&](bool initial = false) {
+        auto dataFile =
+            writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
+        auto dataWriter =
+            writeData ? std::make_unique<
+                            HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+                            *dataFile, vertexCoordinates, vertexBases,
+                            frictionPatches, weakPatches)
+                      : nullptr;
+
       if (writeData) {
         dataWriter->reportSolution(programState, contactNetwork, globalFriction);
         if (!initial)
@@ -288,6 +284,7 @@ int main(int argc, char *argv[]) {
 
         }
 
+        const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "../debug_print/vtk/");
         vtkWriter.write(programState.timeStep, programState.u, programState.v,
                         programState.alpha, stress);
       }