diff --git a/CMakeLists.txt b/CMakeLists.txt
index 653046b1a43ea8cd9045824d9c88b3996622e602..f8e6fc720cff31d243c8f2b14d023998d5403de9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,6 +18,9 @@ include(DuneMacros)
 # start a dune project with information from dune.module
 dune_project()
 
+find_package(Boost REQUIRED system filesystem)
+find_package(HDF5 REQUIRED)
+
 add_subdirectory("src")
 add_subdirectory("dune")
 add_subdirectory("doc")
diff --git a/dune/tectonic/globalfriction.hh b/dune/tectonic/globalfriction.hh
index 3fb7c39999db0cc34b21e08e1034c82918156faf..78ba058d9f1f65addd632018bee30497d3d43fe1 100644
--- a/dune/tectonic/globalfriction.hh
+++ b/dune/tectonic/globalfriction.hh
@@ -82,10 +82,11 @@ template <class Matrix, class Vector> class GlobalFriction {
     return res->regularity(x);
   }
 
-  void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) const {
-    coefficient.resize(x.size());
+  ScalarVector coefficientOfFriction(Vector const &x) const {
+    ScalarVector ret(x.size());
     for (size_t i = 0; i < x.size(); ++i)
-      coefficient[i] = restriction(i)->coefficientOfFriction(x[i]);
+      ret[i] = restriction(i)->coefficientOfFriction(x[i]);
+    return ret;
   }
 
   void virtual updateAlpha(ScalarVector const &alpha) = 0;
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7da8041f580388de467d5bf4b2657e87fa9b25a2..ad428776a1170c39d2ac9c3b2629585f770293ca 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,11 +1,15 @@
 set(SOURCE_FILES
   adaptivetimestepper.cc
   assemblers.cc
-  boundary_writer.cc
   coupledtimestepper.cc
   enumparser.cc
   fixedpointiterator.cc
-  friction_writer.cc
+  hdf5/frictionalboundary-writer.cc
+  hdf5/iteration-writer.cc
+  hdf5/patchinfo-writer.cc
+  hdf5/restart-io.cc
+  hdf5/surface-writer.cc
+  hdf5/time-writer.cc
   rate.cc
   rate/rateupdater.cc
   sand-wedge.cc
@@ -22,11 +26,6 @@ dune_symlink_to_source_files("sand-wedge-data/parset.cfg")
 dune_symlink_to_source_files("sand-wedge-data/parset-2D.cfg")
 dune_symlink_to_source_files("sand-wedge-data/parset-3D.cfg")
 
-find_package(Boost REQUIRED system filesystem serialization)
-
-# dataio.hh expects this to be set
-set_property(DIRECTORY APPEND PROPERTY COMPILE_DEFINITIONS "HAVE_BOOST_SERIALIZATION")
-
 foreach(_dim 2 3)
   set(_target sand-wedge-${_dim}D)
   add_executable(${_target} ${SOURCE_FILES})
@@ -36,8 +35,14 @@ foreach(_dim 2 3)
   set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS})
 
   target_link_libraries(${_target} ${Boost_FILESYSTEM_LIBRARY})
-  target_link_libraries(${_target} ${Boost_SERIALIZATION_LIBRARY})
   target_link_libraries(${_target} ${Boost_SYSTEM_LIBRARY})
 
+  set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIR})
+  target_link_libraries(${_target} ${HDF5_LIBRARIES})
+  target_link_libraries(${_target} "dl") # FIXME: missing from HDF5_LIBRARIES
+
   set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()
+
+# Mark UG as required
+find_package(UG REQUIRED)
diff --git a/src/adaptivetimestepper.cc b/src/adaptivetimestepper.cc
index 1be9b39924b4c5f351345b3b5dd831d468a954fa..8059e28e3b7f0b28e89fd1e80700e461002afd03 100644
--- a/src/adaptivetimestepper.cc
+++ b/src/adaptivetimestepper.cc
@@ -4,6 +4,19 @@
 
 #include "adaptivetimestepper.hh"
 
+void IterationRegister::registerCount(FixedPointIterationCounter count) {
+  totalCount += count;
+}
+
+void IterationRegister::registerFinalCount(FixedPointIterationCounter count) {
+  finalCount = count;
+}
+
+void IterationRegister::reset() {
+  totalCount = FixedPointIterationCounter();
+  finalCount = FixedPointIterationCounter();
+}
+
 template <class Factory, class Updaters, class ErrorNorm>
 AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
     Factory &factory, Dune::ParameterTree const &parset,
@@ -19,17 +32,9 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
       parset_(parset),
       globalFriction_(globalFriction),
       current_(current),
-      R1_(current_.clone()),
       externalForces_(externalForces),
       mustRefine_(mustRefine),
-      iterationWriter_("iterations", std::fstream::out),
-      errorNorm_(errorNorm) {
-  MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                          globalFriction_, R1_, errorNorm,
-                                          externalForces_);
-  stepAndReport("R1", coupledTimeStepper, relativeTime_, relativeTau_);
-  iterationWriter_ << std::endl;
-}
+      errorNorm_(errorNorm) {}
 
 template <class Factory, class Updaters, class ErrorNorm>
 bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
@@ -39,28 +44,15 @@ bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
 template <class Factory, class Updaters, class ErrorNorm>
 bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::coarsen() {
   bool didCoarsen = false;
-
   while (relativeTime_ + relativeTau_ < 1.0) {
-    R2_ = R1_.clone();
-    {
-      MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                              globalFriction_, R2_, errorNorm_,
-                                              externalForces_);
-      stepAndReport("R2", coupledTimeStepper, relativeTime_ + relativeTau_,
-                    relativeTau_);
-    }
-
-    Updaters C = current_.clone();
-    {
-      MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                              globalFriction_, C, errorNorm_,
-                                              externalForces_);
+    R2_ = { R1_.updaters.clone() };
+    step(R2_, relativeTime_ + relativeTau_, relativeTau_);
 
-      stepAndReport("C", coupledTimeStepper, relativeTime_, 2.0 * relativeTau_);
-    }
+    UpdatersWithCount C{ current_.clone() };
+    step(C, relativeTime_, 2.0 * relativeTau_);
 
-    if (!mustRefine_(C, R2_)) {
-      R2_ = { nullptr, nullptr };
+    if (!mustRefine_(C.updaters, R2_.updaters)) {
+      R2_ = {};
       R1_ = C;
       relativeTau_ *= 2.0;
       didCoarsen = true;
@@ -74,21 +66,12 @@ bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::coarsen() {
 template <class Factory, class Updaters, class ErrorNorm>
 void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::refine() {
   while (true) {
-    Updaters F1 = current_.clone();
-    Updaters F2;
-    {
-      MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                              globalFriction_, F1, errorNorm_,
-                                              externalForces_);
-      stepAndReport("F1", coupledTimeStepper, relativeTime_,
-                    relativeTau_ / 2.0);
-
-      F2 = F1.clone();
-      stepAndReport("F2", coupledTimeStepper,
-                    relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
-    }
+    UpdatersWithCount F1{ current_.clone() };
+    step(F1, relativeTime_, relativeTau_ / 2.0);
+    UpdatersWithCount F2{ F1.updaters.clone() };
+    step(F2, relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
 
-    if (!mustRefine_(R1_, F2)) {
+    if (!mustRefine_(R1_.updaters, F2.updaters)) {
       break;
     } else {
       R1_ = F1;
@@ -99,23 +82,33 @@ void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::refine() {
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
+IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
+  if (R2_.updaters == Updaters()) {
+    R1_ = { current_.clone() };
+    step(R1_, relativeTime_, relativeTau_); // counting again upon restart
+  } else {
+    R1_ = R2_;
+  }
+
+  iterationRegister_.reset();
   if (!coarsen())
     refine();
 
-  iterationWriter_ << std::endl;
-
-  current_ = R1_;
-  R1_ = R2_;
+  current_ = R1_.updaters;
+  iterationRegister_.registerFinalCount(R1_.count);
   relativeTime_ += relativeTau_;
+
+  return iterationRegister_;
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::stepAndReport(
-    std::string type, MyCoupledTimeStepper &stepper, double rTime,
-    double rTau) {
-  iterationWriter_ << type << " " << stepper.step(rTime, rTau) << " "
-                   << std::flush;
+void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
+    UpdatersWithCount &updatersAndCount, double rTime, double rTau) {
+  updatersAndCount.count =
+      MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
+                           updatersAndCount.updaters, errorNorm_,
+                           externalForces_).step(rTime, rTau);
+  iterationRegister_.registerCount(updatersAndCount.count);
 }
 
 #include "adaptivetimestepper_tmpl.cc"
diff --git a/src/adaptivetimestepper.hh b/src/adaptivetimestepper.hh
index 7fa6298c33cb6848948ab131432b3d347c7a83e0..dcda38d0e7ab13abbe7c4693a05542fd1bf47031 100644
--- a/src/adaptivetimestepper.hh
+++ b/src/adaptivetimestepper.hh
@@ -5,8 +5,23 @@
 
 #include "coupledtimestepper.hh"
 
+struct IterationRegister {
+  void registerCount(FixedPointIterationCounter count);
+  void registerFinalCount(FixedPointIterationCounter count);
+
+  void reset();
+
+  FixedPointIterationCounter totalCount;
+  FixedPointIterationCounter finalCount;
+};
+
 template <class Factory, class Updaters, class ErrorNorm>
 class AdaptiveTimeStepper {
+  struct UpdatersWithCount {
+    Updaters updaters;
+    FixedPointIterationCounter count;
+  };
+
   using Vector = typename Factory::Vector;
   using ConvexProblem = typename Factory::ConvexProblem;
   using Nonlinearity = typename ConvexProblem::NonlinearityType;
@@ -25,25 +40,25 @@ class AdaptiveTimeStepper {
   bool reachedEnd();
   bool coarsen();
   void refine();
-  void advance();
+  IterationRegister advance();
 
   double relativeTime_;
   double relativeTau_;
 
 private:
-  void stepAndReport(std::string type, MyCoupledTimeStepper &stepper,
-                     double rTime, double rTau);
+  void step(UpdatersWithCount &updatersWithCount, double rTime, double rTau);
 
   double finalTime_;
   Factory &factory_;
   Dune::ParameterTree const &parset_;
   std::shared_ptr<Nonlinearity> globalFriction_;
   Updaters &current_;
-  Updaters R1_;
-  Updaters R2_;
+  UpdatersWithCount R1_;
+  UpdatersWithCount R2_;
   std::function<void(double, Vector &)> externalForces_;
   std::function<bool(Updaters &, Updaters &)> mustRefine_;
-  std::fstream iterationWriter_;
   ErrorNorm const &errorNorm_;
+
+  IterationRegister iterationRegister_;
 };
 #endif
diff --git a/src/boundary_writer.cc b/src/boundary_writer.cc
deleted file mode 100644
index 82d5e2fbbdd33f0fa53f809f9d00f98dd56f58d0..0000000000000000000000000000000000000000
--- a/src/boundary_writer.cc
+++ /dev/null
@@ -1,39 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "boundary_writer.hh"
-#include "tobool.hh"
-
-template <class ScalarVector, class Vector>
-BoundaryWriter<ScalarVector, Vector>::BoundaryWriter(
-    Vector const &vertexCoordinates,
-    Dune::BitSetVector<1> const &_boundaryNodes, std::string prefix,
-    Projector projector)
-    : displacementWriter(prefix + "Displacements", std::fstream::out),
-      velocityWriter(prefix + "Velocities", std::fstream::out),
-      boundaryNodes(_boundaryNodes),
-      projector_(projector) {
-  std::fstream vertexCoordinateWriter(prefix + "Coordinates",
-                                      std::fstream::out);
-  for (size_t i = 0; i < boundaryNodes.size(); ++i)
-    if (toBool(boundaryNodes[i]))
-      vertexCoordinateWriter << vertexCoordinates[i] << std::endl;
-}
-
-template <class ScalarVector, class Vector>
-void BoundaryWriter<ScalarVector, Vector>::writeKinetics(Vector const &u,
-                                                         Vector const &v) {
-  for (size_t i = 0; i < boundaryNodes.size(); ++i) {
-    if (!toBool(boundaryNodes[i]))
-      continue;
-
-    displacementWriter << projector_(u[i]) << " ";
-    velocityWriter << projector_(v[i]) << " ";
-  }
-
-  displacementWriter << std::endl;
-  velocityWriter << std::endl;
-}
-
-#include "boundary_writer_tmpl.cc"
diff --git a/src/boundary_writer.hh b/src/boundary_writer.hh
deleted file mode 100644
index 9f0680f3905d343d612724d44355c4f5a146fe34..0000000000000000000000000000000000000000
--- a/src/boundary_writer.hh
+++ /dev/null
@@ -1,27 +0,0 @@
-#ifndef SRC_BOUNDARY_WRITER_HH
-#define SRC_BOUNDARY_WRITER_HH
-
-#include <fstream>
-#include <string>
-
-#include <dune/common/bitsetvector.hh>
-
-template <class ScalarVector, class Vector> class BoundaryWriter {
-protected:
-  using Projector = std::function<double(typename Vector::block_type const &)>;
-
-public:
-  BoundaryWriter(Vector const &vertexCoordinates,
-                 Dune::BitSetVector<1> const &_boundaryNodes,
-                 std::string prefix, Projector projector);
-
-  void writeKinetics(Vector const &u, Vector const &v);
-
-protected:
-  std::fstream displacementWriter;
-  std::fstream velocityWriter;
-  Dune::BitSetVector<1> const &boundaryNodes;
-
-  Projector projector_;
-};
-#endif
diff --git a/src/boundary_writer_tmpl.cc b/src/boundary_writer_tmpl.cc
deleted file mode 100644
index 378584de089e6f52cb492bcf131dcae58837a2f2..0000000000000000000000000000000000000000
--- a/src/boundary_writer_tmpl.cc
+++ /dev/null
@@ -1,7 +0,0 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "explicitvectors.hh"
-
-template class BoundaryWriter<ScalarVector, Vector>;
diff --git a/src/fixedpointiterator.cc b/src/fixedpointiterator.cc
index cff97d961a7c7965d00ba6ea77073c42648dc97e..68dfd4febbc91337e29e32afffeb763583becbb7 100644
--- a/src/fixedpointiterator.cc
+++ b/src/fixedpointiterator.cc
@@ -13,6 +13,12 @@
 
 #include "fixedpointiterator.hh"
 
+void FixedPointIterationCounter::operator+=(
+    FixedPointIterationCounter const &other) {
+  iterations += other.iterations;
+  multigridIterations += other.multigridIterations;
+}
+
 template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
     Factory &factory, Dune::ParameterTree const &parset,
diff --git a/src/fixedpointiterator.hh b/src/fixedpointiterator.hh
index 9a4a336fc5087fd026db4d7a20d35d90e801c483..4833b496bb9d91e346ec972f80b074f895c1926e 100644
--- a/src/fixedpointiterator.hh
+++ b/src/fixedpointiterator.hh
@@ -11,6 +11,8 @@
 struct FixedPointIterationCounter {
   size_t iterations;
   size_t multigridIterations;
+
+  void operator+=(FixedPointIterationCounter const &other);
 };
 
 std::ostream &operator<<(std::ostream &stream,
diff --git a/src/friction_writer.cc b/src/friction_writer.cc
deleted file mode 100644
index ce71828dd4a2fd5117e21476baa0a927ddb94596..0000000000000000000000000000000000000000
--- a/src/friction_writer.cc
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "friction_writer.hh"
-#include "tobool.hh"
-
-template <class ScalarVector, class Vector>
-FrictionWriter<ScalarVector, Vector>::FrictionWriter(
-    Vector const &vertexCoordinates,
-    Dune::BitSetVector<1> const &_boundaryNodes, std::string prefix,
-    typename BW::Projector projector)
-    : BW(vertexCoordinates, _boundaryNodes, prefix, projector),
-      coefficientWriter(prefix + "Coefficients", std::fstream::out),
-      stateWriter(prefix + "Alpha", std::fstream::out) {}
-
-template <class ScalarVector, class Vector>
-void FrictionWriter<ScalarVector, Vector>::writeOther(
-    ScalarVector const &coefficient, ScalarVector const &alpha) {
-  for (size_t i = 0; i < boundaryNodes.size(); ++i) {
-    if (!toBool(boundaryNodes[i]))
-      continue;
-
-    coefficientWriter << coefficient[i] << " ";
-    stateWriter << alpha[i] << " ";
-  }
-
-  stateWriter << std::endl;
-  coefficientWriter << std::endl;
-}
-
-#include "friction_writer_tmpl.cc"
diff --git a/src/friction_writer.hh b/src/friction_writer.hh
deleted file mode 100644
index 957f5f9f3757d7c60792576cfd1456508088ff43..0000000000000000000000000000000000000000
--- a/src/friction_writer.hh
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef SRC_FRICTION_WRITER_HH
-#define SRC_FRICTION_WRITER_HH
-
-#include <fstream>
-#include <string>
-
-#include <dune/common/bitsetvector.hh>
-
-#include "boundary_writer.hh"
-
-template <class ScalarVector, class Vector>
-class FrictionWriter : public BoundaryWriter<ScalarVector, Vector> {
-  using BW = BoundaryWriter<ScalarVector, Vector>;
-
-public:
-  FrictionWriter(Vector const &vertexCoordinates,
-                 Dune::BitSetVector<1> const &_boundaryNodes,
-                 std::string prefix, typename BW::Projector projector);
-
-  void writeOther(ScalarVector const &coefficient, ScalarVector const &alpha);
-
-private:
-  std::fstream coefficientWriter;
-  std::fstream stateWriter;
-  using BW::boundaryNodes;
-};
-
-#endif
diff --git a/src/friction_writer_tmpl.cc b/src/friction_writer_tmpl.cc
deleted file mode 100644
index 5ac5e7d55873569ad0a84eb55fdd3c2277d44183..0000000000000000000000000000000000000000
--- a/src/friction_writer_tmpl.cc
+++ /dev/null
@@ -1,7 +0,0 @@
-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "explicitvectors.hh"
-
-template class FrictionWriter<ScalarVector, Vector>;
diff --git a/src/hdf5-writer.hh b/src/hdf5-writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e4e5358e5a78ef582cf9f9ca47c29a425d3f7296
--- /dev/null
+++ b/src/hdf5-writer.hh
@@ -0,0 +1,73 @@
+#ifndef SRC_HDF5_WRITER_HH
+#define SRC_HDF5_WRITER_HH
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/hdf5/hdf5file.hh>
+
+#include "hdf5/frictionalboundary-writer.hh"
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/patchinfo-writer.hh"
+#include "hdf5/restart-io.hh"
+#include "hdf5/surface-writer.hh"
+#include "hdf5/time-writer.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5Writer {
+private:
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+  using LocalVector = typename Vector::block_type;
+
+public:
+  HDF5Writer(HDF5Grouplike &file, Vector const &vertexCoordinates,
+             VertexBasis const &vertexBasis, Patch const &surface,
+             Patch const &frictionalBoundary,
+             ConvexPolyhedron<LocalVector> const &weakPatch)
+      : file_(file),
+        iterationWriter_(file_),
+        timeWriter_(file_),
+#if MY_DIM == 3
+        patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
+#endif
+        surfaceWriter_(file_, vertexCoordinates, surface),
+        frictionalBoundaryWriter_(file_, vertexCoordinates, frictionalBoundary),
+        restartWriter_(file_, vertexCoordinates.size()) {
+  }
+
+  template <class Friction>
+  void reportSolution(ProgramState const &programState,
+                      // for the friction coefficient
+                      Friction &friction) {
+    timeWriter_.write(programState);
+#if MY_DIM == 3
+    patchInfoWriter_.write(programState);
+#endif
+    surfaceWriter_.write(programState);
+    frictionalBoundaryWriter_.write(programState, friction);
+  }
+
+  void reportIterations(ProgramState const &programState,
+                        IterationRegister const &iterationCount) {
+    iterationWriter_.write(programState.timeStep, iterationCount);
+  }
+
+  void writeRestart(ProgramState const &programState) {
+    restartWriter_.write(programState);
+  }
+
+private:
+  HDF5Grouplike &file_;
+
+  IterationWriter iterationWriter_;
+  TimeWriter<ProgramState> timeWriter_;
+#if MY_DIM == 3
+  PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
+#endif
+  SurfaceWriter<ProgramState, GridView> surfaceWriter_;
+  FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
+
+  RestartIO<ProgramState> restartWriter_;
+};
+#endif
diff --git a/src/hdf5/frictionalboundary-writer.cc b/src/hdf5/frictionalboundary-writer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b98290dfc103da2ded9f903c1d04a5ea1cf128d2
--- /dev/null
+++ b/src/hdf5/frictionalboundary-writer.cc
@@ -0,0 +1,59 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "frictionalboundary-writer.hh"
+#include "restrict.hh"
+
+template <class ProgramState, class GridView>
+FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
+    HDF5Grouplike &file, Vector const &vertexCoordinates,
+    Patch const &frictionalBoundary)
+    : group_(file, "frictionalBoundary"),
+      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);
+  HDF5SingletonWriter<2> frictionalBoundaryCoordinateWriter(
+      group_, "coordinates", frictionalBoundaryCoordinates.size(),
+      Vector::block_type::dimension);
+  setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
+}
+
+template <class ProgramState, class GridView>
+template <class Friction>
+void FrictionalBoundaryWriter<ProgramState, GridView>::write(
+    ProgramState const &programState, Friction &friction) {
+  auto const frictionalBoundaryDisplacements =
+      restrictToSurface(programState.u, frictionalBoundary_);
+  addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
+           frictionalBoundaryDisplacements);
+
+  auto const frictionalBoundaryVelocities =
+      restrictToSurface(programState.v, frictionalBoundary_);
+  addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
+           frictionalBoundaryVelocities);
+
+  auto const frictionalBoundaryStates =
+      restrictToSurface(programState.alpha, frictionalBoundary_);
+  addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
+           frictionalBoundaryStates);
+
+  friction.updateAlpha(programState.alpha);
+  auto const c = friction.coefficientOfFriction(programState.v);
+  auto const frictionalBoundaryCoefficient =
+      restrictToSurface(c, frictionalBoundary_);
+  addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
+           frictionalBoundaryCoefficient);
+}
+
+#include "frictionalboundary-writer_tmpl.cc"
diff --git a/src/hdf5/frictionalboundary-writer.hh b/src/hdf5/frictionalboundary-writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..772b82df1e0ea82cca3ae8a30d820480bd7c0aa6
--- /dev/null
+++ b/src/hdf5/frictionalboundary-writer.hh
@@ -0,0 +1,31 @@
+#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
+#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/fufem/hdf5/hdf5-sequence-io.hh>
+#include <dune/fufem/hdf5/hdf5-singleton-writer.hh>
+
+template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
+  using ScalarVector = typename ProgramState::ScalarVector;
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+public:
+  FrictionalBoundaryWriter(HDF5Grouplike &file, Vector const &vertexCoordinates,
+                           Patch const &frictionalBoundary);
+
+  template <class Friction>
+  void write(ProgramState const &programState, Friction &friction);
+
+private:
+  HDF5Group group_;
+
+  Patch const &frictionalBoundary_;
+
+  HDF5SequenceIO<2> frictionalBoundaryDisplacementWriter_;
+  HDF5SequenceIO<2> frictionalBoundaryVelocityWriter_;
+  HDF5SequenceIO<1> frictionalBoundaryStateWriter_;
+  HDF5SequenceIO<1> frictionalBoundaryCoefficientWriter_;
+};
+#endif
diff --git a/src/hdf5/frictionalboundary-writer_tmpl.cc b/src/hdf5/frictionalboundary-writer_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e39ffe0fbe2388ea1df00157a7f18458daa4a200
--- /dev/null
+++ b/src/hdf5/frictionalboundary-writer_tmpl.cc
@@ -0,0 +1,11 @@
+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+using MyFriction = GlobalFriction<Matrix, Vector>;
+
+template class FrictionalBoundaryWriter<MyProgramState, GridView>;
+template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
+    MyProgramState const &programState, MyFriction &friction);
diff --git a/src/hdf5/iteration-writer.cc b/src/hdf5/iteration-writer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e48f091b970f2e4f657241d95717fca99e8e34b4
--- /dev/null
+++ b/src/hdf5/iteration-writer.cc
@@ -0,0 +1,26 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "iteration-writer.hh"
+
+IterationWriter::IterationWriter(HDF5Grouplike &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/src/hdf5/iteration-writer.hh b/src/hdf5/iteration-writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cdcd80d3f3097ade55064694541188007b432977
--- /dev/null
+++ b/src/hdf5/iteration-writer.hh
@@ -0,0 +1,24 @@
+#ifndef SRC_HDF_ITERATION_WRITER_HH
+#define SRC_HDF_ITERATION_WRITER_HH
+
+#include <dune/fufem/hdf5/hdf5-sequence-io.hh>
+
+#include "../adaptivetimestepper.hh"
+
+class IterationWriter {
+public:
+  IterationWriter(HDF5Grouplike &file);
+
+  void write(size_t timeStep, IterationRegister const &iterationCount);
+
+private:
+  HDF5Group group_;
+  HDF5Group fpiSubGroup_;
+  HDF5Group mgSubGroup_;
+
+  HDF5SequenceIO<0, size_t> finalMGIterationWriter_;
+  HDF5SequenceIO<0, size_t> finalFPIIterationWriter_;
+  HDF5SequenceIO<0, size_t> totalMGIterationWriter_;
+  HDF5SequenceIO<0, size_t> totalFPIIterationWriter_;
+};
+#endif
diff --git a/src/hdf5/patchinfo-writer.cc b/src/hdf5/patchinfo-writer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f0195b1a927a24adcddeaf543bc82ed6a696af9a
--- /dev/null
+++ b/src/hdf5/patchinfo-writer.cc
@@ -0,0 +1,94 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/grid/hierarchic-approximation.hh>
+
+#include "patchinfo-writer.hh"
+
+template <class LocalVector, class GridView>
+GridEvaluator<LocalVector, GridView>::GridEvaluator(
+    ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
+  double const bufferWidth = 0.05 * MyGeometry::lengthScale;
+  auto const xminmax = std::minmax_element(
+      weakPatch.vertices.begin(), weakPatch.vertices.end(),
+      [](LocalVector const &a, LocalVector const &b) { return a[0] < b[0]; });
+  double const xmin = (*xminmax.first)[0] - bufferWidth;
+  double const xspan = (*xminmax.second)[0] + bufferWidth - xmin;
+
+  double const zmin = -MyGeometry::depth / 2.0;
+  double const zspan = MyGeometry::depth;
+
+  double spatialResolution = 0.01 * MyGeometry::lengthScale;
+  size_t const xsteps = std::round(xspan / spatialResolution);
+  size_t const zsteps = std::round(zspan / spatialResolution);
+
+  xCoordinates.resize(xsteps + 1);
+  for (size_t xi = 0; xi <= xsteps; xi++)
+    xCoordinates[xi] = xmin + xi * xspan / xsteps;
+  zCoordinates.resize(zsteps + 1);
+  for (size_t zi = 0; zi <= zsteps; zi++)
+    zCoordinates[zi] = zmin + zi * zspan / zsteps;
+
+  HierarchicApproximation<typename GridView::Grid, GridView> const
+      hApproximation(gridView.grid(), gridView, 1e-6 * MyGeometry::lengthScale);
+
+  LocalVector global(0);
+  localInfo.resize(xsteps + 1);
+  for (size_t xi = 0; xi < xCoordinates.size(); ++xi) {
+    localInfo[xi].resize(zsteps + 1);
+    for (size_t zi = 0; zi < zCoordinates.size(); ++zi) {
+      global[0] = xCoordinates[xi];
+      global[2] = zCoordinates[zi];
+      auto const element = hApproximation.findEntity(global);
+      localInfo[xi][zi] =
+          std::make_pair(element, element->geometry().local(global));
+    }
+  }
+}
+
+template <class LocalVector, class GridView>
+template <class Function>
+Dune::Matrix<typename Function::RangeType>
+GridEvaluator<LocalVector, GridView>::evaluate(Function const &f) const {
+  Dune::Matrix<typename Function::RangeType> ret(xCoordinates.size(),
+                                                 zCoordinates.size());
+  for (size_t xi = 0; xi < localInfo.size(); ++xi) {
+    auto const &localInfoX = localInfo[xi];
+    for (size_t zi = 0; zi < localInfoX.size(); ++zi) {
+      auto const &localInfoXZ = localInfoX[zi];
+      f.evaluateLocal(localInfoXZ.first, localInfoXZ.second, ret[xi][zi]);
+    }
+  }
+  return ret;
+}
+
+template <class ProgramState, class VertexBasis, class GridView>
+PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
+    HDF5Grouplike &file, VertexBasis const &vertexBasis,
+    Patch const &frictionalBoundary,
+    ConvexPolyhedron<LocalVector> const &weakPatch)
+    : group_(file, "weakPatchGrid"),
+      vertexBasis_(vertexBasis),
+      gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
+      weakPatchGridVelocityWriter_(
+          group_, "velocity", gridEvaluator_.xCoordinates.size(),
+          gridEvaluator_.zCoordinates.size(), Vector::block_type::dimension) {
+  HDF5SingletonWriter<1> weakPatchGridXCoordinateWriter(
+      group_, "xCoordinates", gridEvaluator_.xCoordinates.size());
+  setEntry(weakPatchGridXCoordinateWriter, gridEvaluator_.xCoordinates);
+
+  HDF5SingletonWriter<1> weakPatchGridZCoordinateWriter(
+      group_, "zCoordinates", gridEvaluator_.zCoordinates.size());
+  setEntry(weakPatchGridZCoordinateWriter, gridEvaluator_.zCoordinates);
+}
+
+template <class ProgramState, class VertexBasis, class GridView>
+void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
+    ProgramState const &programState) {
+  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
+  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
+}
+
+#include "patchinfo-writer_tmpl.cc"
diff --git a/src/hdf5/patchinfo-writer.hh b/src/hdf5/patchinfo-writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c091f4491f92864f03f64b1d6487e3a3a5117eac
--- /dev/null
+++ b/src/hdf5/patchinfo-writer.hh
@@ -0,0 +1,54 @@
+#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/hdf5-sequence-io.hh>
+#include <dune/fufem/hdf5/hdf5-singleton-writer.hh>
+
+#include "../sand-wedge-data/mygeometry.hh"
+
+template <class LocalVector, class GridView> class GridEvaluator {
+  using Element = typename GridView::Grid::template Codim<0>::Entity;
+  using ElementPointer =
+      typename GridView::Grid::template Codim<0>::EntityPointer;
+
+public:
+  GridEvaluator(ConvexPolyhedron<LocalVector> const &weakPatch,
+                GridView const &gridView);
+
+  template <class Function>
+  Dune::Matrix<typename Function::RangeType> evaluate(Function const &f) const;
+
+  Dune::BlockVector<Dune::FieldVector<double, 1>> xCoordinates;
+  Dune::BlockVector<Dune::FieldVector<double, 1>> zCoordinates;
+
+private:
+  std::vector<std::vector<std::pair<ElementPointer, 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(HDF5Grouplike &file, VertexBasis const &vertexBasis,
+                  Patch const &frictionalBoundary,
+                  ConvexPolyhedron<LocalVector> const &weakPatch);
+
+  void write(ProgramState const &programState);
+
+private:
+  HDF5Group group_;
+
+  VertexBasis const &vertexBasis_;
+
+  GridEvaluator<LocalVector, GridView> const gridEvaluator_;
+  HDF5SequenceIO<3> weakPatchGridVelocityWriter_;
+};
+#endif
diff --git a/src/hdf5/patchinfo-writer_tmpl.cc b/src/hdf5/patchinfo-writer_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ef346e105720a488029f2139cf5c4466af2f29fa
--- /dev/null
+++ b/src/hdf5/patchinfo-writer_tmpl.cc
@@ -0,0 +1,14 @@
+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+using P1Basis = P1NodalBasis<GridView, double>;
+using MyFunction = BasisGridFunction<P1Basis, Vector>;
+
+template class GridEvaluator<LocalVector, GridView>;
+template Dune::Matrix<typename MyFunction::RangeType>
+GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
+
+template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;
diff --git a/src/hdf5/restart-io.cc b/src/hdf5/restart-io.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c7fc74222443b8516a5c610297d48e41e6e61110
--- /dev/null
+++ b/src/hdf5/restart-io.cc
@@ -0,0 +1,49 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "restart-io.hh"
+
+template <class ProgramState>
+RestartIO<ProgramState>::RestartIO(HDF5Grouplike &file, size_t vertexCount)
+    : group_(file, "restarts"),
+      displacementWriter_(group_, "displacement", vertexCount,
+                          Vector::block_type::dimension),
+      velocityWriter_(group_, "velocity", vertexCount,
+                      Vector::block_type::dimension),
+      accelerationWriter_(group_, "acceleration", vertexCount,
+                          Vector::block_type::dimension),
+      stateWriter_(group_, "state", vertexCount),
+      weightedNormalStressWriter_(group_, "weightedNormalStress", vertexCount),
+      relativeTimeWriter_(group_, "relativeTime"),
+      relativeTimeIncrementWriter_(group_, "relativeTimeIncrement") {}
+
+template <class ProgramState>
+void RestartIO<ProgramState>::write(ProgramState const &programState) {
+  addEntry(displacementWriter_, programState.timeStep, programState.u);
+  addEntry(velocityWriter_, programState.timeStep, programState.v);
+  addEntry(accelerationWriter_, programState.timeStep, programState.u);
+  addEntry(stateWriter_, programState.timeStep, programState.alpha);
+  addEntry(weightedNormalStressWriter_, programState.timeStep,
+           programState.weightedNormalStress);
+  addEntry(relativeTimeWriter_, programState.timeStep,
+           programState.relativeTime);
+  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
+           programState.relativeTau);
+}
+
+template <class ProgramState>
+void RestartIO<ProgramState>::read(size_t timeStep,
+                                   ProgramState &programState) {
+  programState.timeStep = timeStep;
+  readEntry(displacementWriter_, timeStep, programState.u);
+  readEntry(velocityWriter_, timeStep, programState.v);
+  readEntry(accelerationWriter_, timeStep, programState.u);
+  readEntry(stateWriter_, timeStep, programState.alpha);
+  readEntry(weightedNormalStressWriter_, timeStep,
+            programState.weightedNormalStress);
+  readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
+  readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
+}
+
+#include "restart-io_tmpl.cc"
diff --git a/src/hdf5/restart-io.hh b/src/hdf5/restart-io.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c42b102279a443aefe4ad1dfc700385c2b82e29b
--- /dev/null
+++ b/src/hdf5/restart-io.hh
@@ -0,0 +1,29 @@
+#ifndef SRC_HDF_RESTART_HDF_HH
+#define SRC_HDF_RESTART_HDF_HH
+
+#include <dune/fufem/hdf5/hdf5file.hh>
+#include <dune/fufem/hdf5/hdf5-sequence-io.hh>
+
+template <class ProgramState> class RestartIO {
+  using ScalarVector = typename ProgramState::ScalarVector;
+  using Vector = typename ProgramState::Vector;
+
+public:
+  RestartIO(HDF5Grouplike &file, size_t vertexCount);
+
+  void write(ProgramState const &programState);
+
+  void read(size_t timeStep, ProgramState &programState);
+
+private:
+  HDF5Group group_;
+
+  HDF5SequenceIO<2> displacementWriter_;
+  HDF5SequenceIO<2> velocityWriter_;
+  HDF5SequenceIO<2> accelerationWriter_;
+  HDF5SequenceIO<1> stateWriter_;
+  HDF5SequenceIO<1> weightedNormalStressWriter_;
+  HDF5SequenceIO<0> relativeTimeWriter_;
+  HDF5SequenceIO<0> relativeTimeIncrementWriter_;
+};
+#endif
diff --git a/src/hdf5/restart-io_tmpl.cc b/src/hdf5/restart-io_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..da72bb4262290b5f318d897739cd9961e7194227
--- /dev/null
+++ b/src/hdf5/restart-io_tmpl.cc
@@ -0,0 +1,7 @@
+#include "../explicitvectors.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class RestartIO<MyProgramState>;
diff --git a/src/hdf5/restrict.hh b/src/hdf5/restrict.hh
new file mode 100644
index 0000000000000000000000000000000000000000..497c5890f834a1a3e747fd265d00291ed958192a
--- /dev/null
+++ b/src/hdf5/restrict.hh
@@ -0,0 +1,23 @@
+#ifndef SRC_HDF5_RESTRICT_HH
+#define SRC_HDF5_RESTRICT_HH
+
+#include <cassert>
+
+#include <dune/common/bitsetvector.hh>
+
+#include "../tobool.hh"
+
+template <class Vector, class Patch>
+Vector restrictToSurface(Vector const &v1, Patch const &patch) {
+  auto const &vertices = *patch.getVertices();
+  assert(vertices.size() == v1.size());
+
+  Vector ret(vertices.count());
+  auto target = ret.begin();
+  for (size_t i = 0; i < v1.size(); ++i)
+    if (toBool(vertices[i]))
+      *(target++) = v1[i];
+  assert(target == ret.end());
+  return ret;
+}
+#endif
diff --git a/src/hdf5/surface-writer.cc b/src/hdf5/surface-writer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..51b9508a98970715173b3941aa87f7630d494d19
--- /dev/null
+++ b/src/hdf5/surface-writer.cc
@@ -0,0 +1,35 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "restrict.hh"
+#include "surface-writer.hh"
+
+template <class ProgramState, class GridView>
+SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
+    HDF5Grouplike &file, Vector const &vertexCoordinates, Patch const &surface)
+    : group_(file, "surface"),
+      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);
+  HDF5SingletonWriter<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, surface_);
+  addEntry(surfaceDisplacementWriter_, programState.timeStep,
+           surfaceDisplacements);
+
+  auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
+  addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
+}
+
+#include "surface-writer_tmpl.cc"
diff --git a/src/hdf5/surface-writer.hh b/src/hdf5/surface-writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7ac92a4c96cee5c7ed194a9414660c812ea4be92
--- /dev/null
+++ b/src/hdf5/surface-writer.hh
@@ -0,0 +1,27 @@
+#ifndef SRC_HDF_SURFACE_WRITER_HH
+#define SRC_HDF_SURFACE_WRITER_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/fufem/hdf5/hdf5-sequence-io.hh>
+#include <dune/fufem/hdf5/hdf5-singleton-writer.hh>
+
+template <class ProgramState, class GridView> class SurfaceWriter {
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+public:
+  SurfaceWriter(HDF5Grouplike &file, Vector const &vertexCoordinates,
+                Patch const &surface);
+
+  void write(ProgramState const &programState);
+
+private:
+  HDF5Group group_;
+
+  Patch const &surface_;
+
+  HDF5SequenceIO<2> surfaceDisplacementWriter_;
+  HDF5SequenceIO<2> surfaceVelocityWriter_;
+};
+#endif
diff --git a/src/hdf5/surface-writer_tmpl.cc b/src/hdf5/surface-writer_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e4570aa79c09b1113bfd0ea29a670c571ee4c55c
--- /dev/null
+++ b/src/hdf5/surface-writer_tmpl.cc
@@ -0,0 +1,8 @@
+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class SurfaceWriter<MyProgramState, GridView>;
diff --git a/src/hdf5/time-writer.cc b/src/hdf5/time-writer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1c1419461347592b48dcbdf57848abf4afcf5bcf
--- /dev/null
+++ b/src/hdf5/time-writer.cc
@@ -0,0 +1,21 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "time-writer.hh"
+
+template <class ProgramState>
+TimeWriter<ProgramState>::TimeWriter(HDF5Grouplike &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/src/hdf5/time-writer.hh b/src/hdf5/time-writer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e8f7241790f0466487630cac1e209e2950d70de1
--- /dev/null
+++ b/src/hdf5/time-writer.hh
@@ -0,0 +1,17 @@
+#ifndef SRC_HDF_TIME_WRITER_HH
+#define SRC_HDF_TIME_WRITER_HH
+
+#include <dune/fufem/hdf5/hdf5file.hh>
+#include <dune/fufem/hdf5/hdf5-sequence-io.hh>
+
+template <class ProgramState> class TimeWriter {
+public:
+  TimeWriter(HDF5Grouplike &file);
+  void write(ProgramState const &programState);
+
+private:
+  HDF5Grouplike &file_;
+  HDF5SequenceIO<0> relativeTimeWriter_;
+  HDF5SequenceIO<0> relativeTimeIncrementWriter_;
+};
+#endif
diff --git a/src/hdf5/time-writer_tmpl.cc b/src/hdf5/time-writer_tmpl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7b14edeb8eb05e7ebff3d22b61f8ef249c262ada
--- /dev/null
+++ b/src/hdf5/time-writer_tmpl.cc
@@ -0,0 +1,7 @@
+#include "../explicitvectors.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class TimeWriter<MyProgramState>;
diff --git a/src/program_state.hh b/src/program_state.hh
index c8b02248babfe4c0173bd42c600a2f921a96ca44..a1734168d53414f303eb4ddcb0ed5290892c1c6b 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -4,7 +4,6 @@
 #include <dune/common/parametertree.hh>
 
 #include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/dunedataio.hh>
 
 #include <dune/tectonic/body.hh>
 
@@ -12,8 +11,11 @@
 #include "matrices.hh"
 #include "solverfactory.hh"
 
-template <class Vector, class ScalarVector> class ProgramState {
+template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 public:
+  using Vector = VectorTEMPLATE;
+  using ScalarVector = ScalarVectorTEMPLATE;
+
   ProgramState(int leafVertexCount)
       : u(leafVertexCount),
         v(leafVertexCount),
@@ -112,21 +114,4 @@ template <class Vector, class ScalarVector> class ProgramState {
   size_t timeStep;
 };
 
-namespace boost {
-namespace serialization {
-  template <class Archive, class Vector, class ScalarVector>
-  void serialize(Archive &ar, ProgramState<Vector, ScalarVector> &ps,
-                 const unsigned int version) {
-    ar &ps.u;
-    ar &ps.v;
-    ar &ps.a;
-    ar &ps.alpha;
-    ar &ps.weightedNormalStress;
-    ar &ps.relativeTime;
-    ar &ps.relativeTau;
-    ar &ps.timeStep;
-  }
-}
-}
-
 #endif
diff --git a/src/sand-wedge-data/mygeometry.cc b/src/sand-wedge-data/mygeometry.cc
index ede4a864b4d1dce08158043222006afadcbfcf92..fef8452a2c698e745ba5e8545df0a07188b02270 100644
--- a/src/sand-wedge-data/mygeometry.cc
+++ b/src/sand-wedge-data/mygeometry.cc
@@ -12,13 +12,6 @@
 
 #include "mygeometry.hh"
 
-double MyGeometry::verticalProjection(LocalVector const &x) {
-  return zenith * x;
-}
-double MyGeometry::horizontalProjection(LocalVector const &x) {
-  return orthogonalZenith * x;
-}
-
 void MyGeometry::write() {
   std::fstream writer("geometry", std::fstream::out);
   writer << "A = " << A << std::endl;
diff --git a/src/sand-wedge-data/mygeometry.hh b/src/sand-wedge-data/mygeometry.hh
index b6eac051dc81f681be515a9f8ac519087e2bc93d..720a5bf3b30b959ea1b50f0585b088a7a8fd2ef9 100644
--- a/src/sand-wedge-data/mygeometry.hh
+++ b/src/sand-wedge-data/mygeometry.hh
@@ -65,7 +65,6 @@ namespace reference {
   LocalVector2D const I = createVector(Y[0] + G[0], Y[1] + G[1]);
 
   LocalVector2D const zenith = createVector(0, 1);
-  LocalVector2D const orthogonalZenith = createVector(-1, 0);
 
   LocalMatrix2D const rotation =
       createMatrix(std::cos(leftAngle), -std::sin(leftAngle),
@@ -102,10 +101,6 @@ LocalVector const Y = rotate(reference::Y);
 LocalVector const Z = rotate(reference::Z);
 
 LocalVector const zenith = rotate(reference::zenith);
-LocalVector const orthogonalZenith = rotate(reference::orthogonalZenith);
-
-double verticalProjection(LocalVector const &);
-double horizontalProjection(LocalVector const &);
 
 void write();
 
diff --git a/src/sand-wedge-data/special_writer.hh b/src/sand-wedge-data/special_writer.hh
deleted file mode 100644
index efd476f476f63c7de49fd5ceeed00c54bf41925d..0000000000000000000000000000000000000000
--- a/src/sand-wedge-data/special_writer.hh
+++ /dev/null
@@ -1,83 +0,0 @@
-#ifndef SRC_SPECIAL_WRITER_HH
-#define SRC_SPECIAL_WRITER_HH
-
-#include <fstream>
-#include <utility>
-
-#include <dune/common/fvector.hh>
-#include <dune/grid/utility/hierarchicsearch.hh>
-#include <dune/fufem/functions/virtualgridfunction.hh>
-
-#include "mygeometry.hh"
-
-template <class GridView, int dimension> class SpecialWriter {
-  using LocalVector = Dune::FieldVector<double, dimension>;
-  using Element = typename GridView::Grid::template Codim<0>::Entity;
-  using ElementPointer =
-      typename GridView::Grid::template Codim<0>::EntityPointer;
-
-  void writeHorizontal(LocalVector const &v) {
-    writer_ << MyGeometry::horizontalProjection(v) << " ";
-  }
-  void writeVertical(LocalVector const &v) {
-    writer_ << MyGeometry::verticalProjection(v) << " ";
-  }
-
-  std::pair<ElementPointer, LocalVector> globalToLocal(
-      LocalVector const &x) const {
-    auto const element = hsearch_.findEntity(x);
-    return std::make_pair(element, element->geometry().local(x));
-  }
-
-  std::fstream writer_;
-
-  Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
-
-  std::pair<ElementPointer, LocalVector> const G;
-  std::pair<ElementPointer, LocalVector> const H;
-  std::pair<ElementPointer, LocalVector> const J;
-  std::pair<ElementPointer, LocalVector> const I;
-  std::pair<ElementPointer, LocalVector> const U;
-  std::pair<ElementPointer, LocalVector> const Z;
-
-public:
-  SpecialWriter(std::string filename, GridView const &gridView)
-      : writer_(filename, std::fstream::out),
-        hsearch_(gridView.grid(), gridView),
-        G(globalToLocal(MyGeometry::G)),
-        H(globalToLocal(MyGeometry::H)),
-        J(globalToLocal(MyGeometry::J)),
-        I(globalToLocal(MyGeometry::I)),
-        U(globalToLocal(MyGeometry::U)),
-        Z(globalToLocal(MyGeometry::Z)) {
-    writer_ << "Gh Hh Jh Ih Uv Uh Zv Zh" << std::endl;
-  }
-
-  void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
-                 specialField) {
-    LocalVector value;
-
-    specialField.evaluateLocal(*G.first, G.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*H.first, H.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*J.first, J.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*I.first, I.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*U.first, U.second, value);
-    writeVertical(value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*Z.first, Z.second, value);
-    writeVertical(value);
-    writeHorizontal(value);
-
-    writer_ << std::endl;
-  }
-};
-#endif
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 3b7aa548b7a879dbeae11c300b41b67f4646e4dd..adc33009ad7323197a2ef61fdb78b20860b08855 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -39,9 +39,7 @@
 #include <dune/fufem/boundarypatch.hh>
 #pragma clang diagnostic pop
 
-#include <dune/fufem/dataio.hh>
 #include <dune/fufem/dunepython.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/sharedpointermap.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
@@ -52,14 +50,16 @@
 
 #include <dune/tectonic/myblockproblem.hh>
 #include <dune/tectonic/globalfriction.hh>
+#include <dune/fufem/hdf5/hdf5file.hh>
 
 #include "adaptivetimestepper.hh"
 #include "assemblers.hh"
 #include "diameter.hh"
 #include "enumparser.hh"
 #include "enums.hh"
-#include "friction_writer.hh"
 #include "gridselector.hh"
+#include "hdf5-writer.hh"
+#include "hdf5/restart-io.hh"
 #include "matrices.hh"
 #include "program_state.hh"
 #include "rate.hh"
@@ -67,7 +67,6 @@
 #include "sand-wedge-data/mygeometry.hh"
 #include "sand-wedge-data/myglobalfrictiondata.hh"
 #include "sand-wedge-data/mygrid.hh"
-#include "sand-wedge-data/special_writer.hh"
 #include "sand-wedge-data/weakpatch.hh"
 #include "solverfactory.hh"
 #include "state.hh"
@@ -148,18 +147,9 @@ int main(int argc, char *argv[]) {
 
     auto myFaces = gridConstructor.constructFaces(leafView);
 
-    // Neumann boundary
     BoundaryPatch<GridView> const neumannBoundary(leafView);
-
-    // Frictional Boundary
     BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    Dune::BitSetVector<1> frictionalNodes(leafVertexCount);
-    frictionalBoundary.getVertices(frictionalNodes);
-
-    // Surface
     BoundaryPatch<GridView> const &surface = myFaces.upper;
-    Dune::BitSetVector<1> surfaceNodes(leafVertexCount);
-    surface.getVertices(surfaceNodes);
 
     // Dirichlet Boundary
     Dune::BitSetVector<dims> noNodes(leafVertexCount);
@@ -228,14 +218,14 @@ int main(int argc, char *argv[]) {
     auto const restartTemplate = parset.get<std::string>("restarts.template");
     auto const restartDirectory = fs::path(restartTemplate).parent_path();
 
+    HDF5File ioFile("output.h5");
     if (firstRestart != 0)
-      DataIO::loadData(programState,
-                       str(boost::format(restartTemplate) % firstRestart));
+      RestartIO<ProgramState<Vector, ScalarVector>>(ioFile, leafVertexCount)
+          .read(firstRestart, programState);
     else
       programState.setupInitialConditions(parset, computeExternalForces,
                                           matrices, myAssembler, dirichletNodes,
                                           noNodes, frictionalBoundary, body);
-    assert(programState.timeStep == firstRestart);
 
     MyGlobalFrictionData<LocalVector> frictionInfo(
         parset.sub("boundary.friction"), weakPatch);
@@ -255,54 +245,24 @@ int main(int argc, char *argv[]) {
       }
     }
 
-    // Set up writers
-    FrictionWriter<ScalarVector, Vector> frictionWriter(
-        vertexCoordinates, frictionalNodes, "friction",
-        MyGeometry::horizontalProjection);
-    BoundaryWriter<ScalarVector, Vector> verticalSurfaceWriter(
-        vertexCoordinates, surfaceNodes, "verticalSurface",
-        MyGeometry::verticalProjection);
-    BoundaryWriter<ScalarVector, Vector> horizontalSurfaceWriter(
-        vertexCoordinates, surfaceNodes, "horizontalSurface",
-        MyGeometry::horizontalProjection);
-    SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities",
-                                                        leafView);
-    SpecialWriter<GridView, dims> specialDisplacementWriter(
-        "specialDisplacements", leafView);
-
-    std::fstream timeStepWriter("timeSteps", std::fstream::out);
-    timeStepWriter << std::fixed << std::setprecision(10);
-
+    HDF5Writer<ProgramState<Vector, ScalarVector>,
+               typename MyAssembler::VertexBasis, GridView>
+        hdf5Writer(ioFile, vertexCoordinates, myAssembler.vertexBasis, surface,
+                   frictionalBoundary, weakPatch);
     MyVTKWriter<typename MyAssembler::VertexBasis,
                 typename MyAssembler::CellBasis> const
         vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
-    auto const report = [&]() {
-      if (programState.timeStep != 0)
-        timeStepWriter << programState.relativeTime << " "
-                       << programState.relativeTau << std::endl;
-
-      horizontalSurfaceWriter.writeKinetics(programState.u, programState.v);
-      verticalSurfaceWriter.writeKinetics(programState.u, programState.v);
-
-      ScalarVector c;
-      myGlobalFriction->coefficientOfFriction(programState.v, c);
-      frictionWriter.writeKinetics(programState.u, programState.v);
-      frictionWriter.writeOther(c, programState.alpha);
-
-      BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
-          myAssembler.vertexBasis, programState.v);
-      BasisGridFunction<typename MyAssembler::VertexBasis, Vector> displacement(
-          myAssembler.vertexBasis, programState.u);
-      specialVelocityWriter.write(velocity);
-      specialDisplacementWriter.write(displacement);
-
-      if (programState.timeStep % restartSpacing == 0) {
-        if (!fs::is_directory(restartDirectory))
-          fs::create_directories(restartDirectory);
-        DataIO::writeData(programState, str(boost::format(restartTemplate) %
-                                            programState.timeStep));
-      }
+    IterationRegister iterationCount;
+    auto const report = [&](bool initial = false) {
+      hdf5Writer.reportSolution(programState, *myGlobalFriction);
+      if (!initial)
+        hdf5Writer.reportIterations(programState, iterationCount);
+
+      if (!initial and programState.timeStep % restartSpacing == 0)
+        hdf5Writer.writeRestart(programState);
+
+      ioFile.flush();
 
       if (parset.get<bool>("io.writeVTK")) {
         ScalarVector stress;
@@ -313,7 +273,7 @@ int main(int argc, char *argv[]) {
                         programState.alpha, stress);
       }
     };
-    report();
+    report(true);
 
     // Set up TNNMG solver
     using NonlinearFactory = SolverFactory<
@@ -330,7 +290,7 @@ int main(int argc, char *argv[]) {
                         programState.u, programState.v, programState.a),
         initStateUpdater<ScalarVector, Vector>(
             parset.get<Config::stateModel>("boundary.friction.stateModel"),
-            programState.alpha, frictionalNodes,
+            programState.alpha, *frictionalBoundary.getVertices(),
             parset.get<double>("boundary.friction.L"),
             parset.get<double>("boundary.friction.V0")));
 
@@ -355,7 +315,7 @@ int main(int argc, char *argv[]) {
     while (!adaptiveTimeStepper.reachedEnd()) {
       programState.timeStep++;
 
-      adaptiveTimeStepper.advance();
+      iterationCount = adaptiveTimeStepper.advance();
 
       programState.relativeTime = adaptiveTimeStepper.relativeTime_;
       programState.relativeTau = adaptiveTimeStepper.relativeTau_;
diff --git a/src/updaters.hh b/src/updaters.hh
index 6d27f445b301eacb9d8b73eb48fa9f6cb4764da7..f72c4a718c5980684f91294e6666827625311070 100644
--- a/src/updaters.hh
+++ b/src/updaters.hh
@@ -12,6 +12,10 @@ struct Updaters {
            std::shared_ptr<StateUpdaterTEMPLATE> stateUpdater)
       : rate_(rateUpdater), state_(stateUpdater) {}
 
+  bool operator==(Updaters const &other) const {
+    return other.rate_ == rate_ and other.state_ == state_;
+  }
+
   Updaters<RateUpdater, StateUpdater> clone() const {
     return Updaters<RateUpdater, StateUpdater>(rate_->clone(), state_->clone());
   }