src/assemblers.hh

@@ -29,8 +29,9 @@ template <class GridView, int dimension> class MyAssembler {
   using CellBasis = P0Basis<GridView, double>;
   using VertexBasis = P1NodalBasis<GridView, double>;
 
-  CellBasis cellBasis;
-  VertexBasis vertexBasis;
+  CellBasis const cellBasis;
+  VertexBasis const vertexBasis;
+  GridView const &gridView;
 
 private:
   using Grid = typename GridView::Grid;
@@ -40,7 +41,6 @@ template <class GridView, int dimension> class MyAssembler {
   using LocalCellBasis = typename CellBasis::LocalFiniteElement;
   using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
-  GridView const &gridView;
   Assembler<CellBasis, CellBasis> cellAssembler;
   Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
@@ -49,41 +49,42 @@ template <class GridView, int dimension> class MyAssembler {
 
   void assembleFrictionalBoundaryMass(
       BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass);
+      ScalarMatrix &frictionalBoundaryMass) const;
 
   void assembleMass(Dune::VirtualFunction<
                         LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M);
+                    Matrix &M) const;
 
-  void assembleElasticity(double E, double nu, Matrix &A);
+  void assembleElasticity(double E, double nu, Matrix &A) const;
 
   void assembleViscosity(
       Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
           shearViscosity,
       Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
           bulkViscosity,
-      Matrix &C);
+      Matrix &C) const;
 
   void assembleBodyForce(
       Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f);
+      Vector &f) const;
 
   void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
                        Vector &f,
                        Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime);
+                       double relativeTime) const;
 
-  void assembleNormalStress(BoundaryPatch<GridView> const &frictionalBoundary,
-                            ScalarVector &normalStress, double youngModulus,
-                            double poissonRatio, Vector const &displacement);
+  void assembleWeightedNormalStress(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      ScalarVector &weightedNormalStress, double youngModulus,
+      double poissonRatio, Vector const &displacement) const;
 
   std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
       Config::FrictionModel frictionModel,
       BoundaryPatch<GridView> const &frictionalBoundary,
       GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &normalStress);
+      ScalarVector const &weightedNormalStress) const;
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress);
+                              Vector const &u, ScalarVector &stress) const;
 };
 #endif

src/boundary_writer.cc

-#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;
-  vertexCoordinateWriter.close();
-}
-
-template <class ScalarVector, class Vector>
-BoundaryWriter<ScalarVector, Vector>::~BoundaryWriter() {
-  displacementWriter.close();
-  velocityWriter.close();
-}
-
-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"

src/boundary_writer.hh

-#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);
-
-  virtual ~BoundaryWriter();
-
-  void writeKinetics(Vector const &u, Vector const &v);
-
-protected:
-  std::fstream displacementWriter;
-  std::fstream velocityWriter;
-  Dune::BitSetVector<1> const &boundaryNodes;
-
-  Projector projector_;
-};
-#endif

src/boundary_writer_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "explicitvectors.hh"
-
-template class BoundaryWriter<ScalarVector, Vector>;

src/coupledtimestepper.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "coupledtimestepper.hh"
-
-template <class Factory, class StateUpdater, class VelocityUpdater,
-          class ErrorNorm>
-CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
-    CoupledTimeStepper(double finalTime, Factory &factory,
-                       Dune::ParameterTree const &parset,
-                       std::shared_ptr<Nonlinearity> globalFriction,
-                       std::shared_ptr<StateUpdater> stateUpdater,
-                       std::shared_ptr<VelocityUpdater> velocityUpdater,
-                       ErrorNorm const &errorNorm,
-                       std::function<void(double, Vector &)> externalForces)
-    : finalTime_(finalTime),
-      factory_(factory),
-      parset_(parset),
-      globalFriction_(globalFriction),
-      stateUpdater_(stateUpdater),
-      velocityUpdater_(velocityUpdater),
-      externalForces_(externalForces),
-      errorNorm_(errorNorm) {}
-
-template <class Factory, class StateUpdater, class VelocityUpdater,
-          class ErrorNorm>
-FixedPointIterationCounter
-CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::step(
-    double relativeTime, double relativeTau) {
-  stateUpdater_->nextTimeStep();
-  velocityUpdater_->nextTimeStep();
-
-  auto const newRelativeTime = relativeTime + relativeTau;
-  Vector ell;
-  externalForces_(newRelativeTime, ell);
-
-  Matrix velocityMatrix;
-  Vector velocityRHS;
-  Vector velocityIterate;
-
-  auto const tau = relativeTau * finalTime_;
-  stateUpdater_->setup(tau);
-  velocityUpdater_->setup(ell, tau, newRelativeTime, velocityRHS,
-                          velocityIterate, velocityMatrix);
-  FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>
-      fixedPointIterator(factory_, parset_, globalFriction_, errorNorm_);
-  auto const iterations =
-      fixedPointIterator.run(stateUpdater_, velocityUpdater_, velocityMatrix,
-                             velocityRHS, velocityIterate);
-  return iterations;
-}
-
-#include "coupledtimestepper_tmpl.cc"

src/distances.hh → src/diameter.hh

-#ifndef SRC_DISTANCES_HH
-#define SRC_DISTANCES_HH
-
-#include <dune/tectonic/polyhedrondistance.hh>
-#include "sand-wedge-data/mygeometry.hh"
+#ifndef SRC_DIAMETER_HH
+#define SRC_DIAMETER_HH
 
 template <class Geometry> double diameter(Geometry const &geometry) {
   auto const numCorners = geometry.corners();
@@ -16,18 +13,4 @@ template <class Geometry> double diameter(Geometry const &geometry) {
   }
   return diameter;
 }
-
-template <class Geometry>
-double distanceToWeakeningRegion(
-    Geometry const &g,
-    ConvexPolyhedron<typename Geometry::GlobalCoordinate> const &weakPatch) {
-  using Coordinate = typename Geometry::GlobalCoordinate;
-
-  ConvexPolyhedron<Coordinate> bsg;
-  bsg.vertices.resize(g.corners());
-  for (size_t i = 0; i < bsg.vertices.size(); ++i)
-    bsg.vertices[i] = g.corner(i);
-
-  return distance(bsg, weakPatch, 1e-6 * MyGeometry::lengthScale);
-}
 #endif

src/enumparser.cc

@@ -9,8 +9,8 @@
 #include "enumparser.hh"
 
 template <class Enum>
-typename Dune::enable_if<
-    !Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value,
+typename std::enable_if<
+    !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
     std::istream &>::type
 operator>>(std::istream &lhs, Enum &e) {
   std::string s;
@@ -56,16 +56,13 @@ Config::stateModel StringToEnum<Config::stateModel>::convert(
   DUNE_THROW(Dune::Exception, "failed to parse enum");
 }
 
-Solver::VerbosityMode StringToEnum<Solver::VerbosityMode>::convert(
+Config::PatchType StringToEnum<Config::PatchType>::convert(
     std::string const &s) {
-  if (s == "full")
-    return Solver::FULL;
+  if (s == "Rectangular")
+    return Config::Rectangular;
 
-  if (s == "reduced")
-    return Solver::REDUCED;
-
-  if (s == "quiet")
-    return Solver::QUIET;
+  if (s == "Trapezoidal")
+    return Config::Trapezoidal;
 
   DUNE_THROW(Dune::Exception, "failed to parse enum");
 }
@@ -73,4 +70,4 @@ Solver::VerbosityMode StringToEnum<Solver::VerbosityMode>::convert(
 template std::istream &operator>>(std::istream &lhs, Config::FrictionModel &);
 template std::istream &operator>>(std::istream &lhs, Config::stateModel &);
 template std::istream &operator>>(std::istream &lhs, Config::scheme &);
-template std::istream &operator>>(std::istream &lhs, Solver::VerbosityMode &);
+template std::istream &operator>>(std::istream &lhs, Config::PatchType &);

src/enumparser.hh

@@ -3,7 +3,7 @@
 
 // Copyright Carsten Graeser 2012
 
-#include <dune/common/typetraits.hh>
+#include <type_traits>
 
 #include <dune/solvers/solvers/solver.hh>
 
@@ -23,13 +23,13 @@ template <> struct StringToEnum<Config::scheme> {
   static Config::scheme convert(std::string const &s);
 };
 
-template <> struct StringToEnum<Solver::VerbosityMode> {
-  static Solver::VerbosityMode convert(std::string const &s);
+template <> struct StringToEnum<Config::PatchType> {
+  static Config::PatchType convert(std::string const &s);
 };
 
 template <class Enum>
-typename Dune::enable_if<
-    !Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value,
+typename std::enable_if<
+    !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
     std::istream &>::type
 operator>>(std::istream &lhs, Enum &e);
 #endif

src/enums.hh

@@ -5,6 +5,7 @@ struct Config {
   enum FrictionModel { Truncated, Regularised };
   enum stateModel { AgeingLaw, SlipLaw };
   enum scheme { Newmark, BackwardEuler };
+  enum PatchType { Rectangular, Trapezoidal };
 };
 
 #endif

src/friction_writer.cc

-#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>
-FrictionWriter<ScalarVector, Vector>::~FrictionWriter() {
-  stateWriter.close();
-  coefficientWriter.close();
-}
-
-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"

src/friction_writer.hh

-#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);
-
-  ~FrictionWriter();
-
-  void writeOther(ScalarVector const &coefficient, ScalarVector const &alpha);
-
-private:
-  std::fstream coefficientWriter;
-  std::fstream stateWriter;
-  using BW::boundaryNodes;
-};
-
-#endif

src/friction_writer_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "explicitvectors.hh"
-
-template class FrictionWriter<ScalarVector, Vector>;

src/gridselector.hh

@@ -12,10 +12,7 @@
 #if !HAVE_ALUGRID
 #error ALUGRID was requested but not found
 #endif
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wignored-qualifiers"
 #include <dune/grid/alugrid.hh>
-#pragma clang diagnostic pop
 using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
 
 #elif WANT_GRID == WANT_UG

src/hdf5-writer.hh

+#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/file.hh>
+
+#include "hdf5/frictionalboundary-writer.hh"
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/patchinfo-writer.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(HDF5::Grouplike &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) {
+  }
+
+  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);
+  }
+
+private:
+  HDF5::Grouplike &file_;
+
+  IterationWriter iterationWriter_;
+  TimeWriter<ProgramState> timeWriter_;
+#if MY_DIM == 3
+  PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
+#endif
+  SurfaceWriter<ProgramState, GridView> surfaceWriter_;
+  FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
+};
+#endif

src/hdf5/frictionalboundary-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+#include "frictionalboundary-writer.hh"
+#include "restrict.hh"
+
+template <class ProgramState, class GridView>
+FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
+    HDF5::Grouplike &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);
+  HDF5::SingletonWriter<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"

src/hdf5/frictionalboundary-writer.hh

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

src/hdf5/frictionalboundary-writer_tmpl.cc

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

src/hdf5/iteration-writer.cc

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

src/hdf5/iteration-writer.hh

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

src/hdf5/patchinfo-writer.cc

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