dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -5,4 +5,4 @@
 Module: dune-tectonic
 Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Depends: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/globalratestatefriction.hh

@@ -32,7 +32,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       : restrictions(), localToGlobal(), zeroFriction() {
     auto const gridView = frictionalBoundary.gridView();
     Dune::MultipleCodimMultipleGeomTypeMapper<
-        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
     for (auto it = gridView.template begin<block_size>();
          it != gridView.template end<block_size>(); ++it) {
       auto const i = vertexMapper.index(*it);
@@ -44,8 +44,8 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       restrictions.emplace_back(weights[i], weightedNormalStress[i],
                                 frictionInfo(geoToPoint(it->geometry())));
     }
-    assert(restrictions.size() == frictionalBoundary.numVertices());
-    assert(localToGlobal.size() == frictionalBoundary.numVertices());
+    assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
   }
 
   void updateAlpha(ScalarVector const &alpha) override {

dune/tectonic/myblockproblem.hh

@@ -57,7 +57,7 @@ class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
       : Base(parset, problem),
         maxEigenvalues_(problem.f.size()),
-        localBisection(0.0, 1.0, 0.0, true, 0.0) {
+        localBisection(0.0, 1.0, 0.0, 0.0) {
     for (size_t i = 0; i < problem.f.size(); ++i) {
       LocalVectorType eigenvalues;
       Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);

src/CMakeLists.txt

@@ -19,6 +19,29 @@ set(SW_SOURCE_FILES
   time-stepping/state.cc
   vtk.cc
 )
+
+set(MSW_SOURCE_FILES
+  assemblers.cc
+  enumparser.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
+  multi-body-problem-data/cuboidgeometry.cc
+  multi-body-problem-data/mygrids.cc
+  multi-body-problem.cc
+  spatial-solving/fixedpointiterator.cc
+  spatial-solving/solverfactory.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+  vtk.cc
+)
+
 set(UGW_SOURCE_FILES
   assemblers.cc # FIXME
   one-body-problem-data/mygrid.cc
@@ -28,15 +51,20 @@ set(UGW_SOURCE_FILES
 
 foreach(_dim 2 3)
   set(_sw_target one-body-problem-${_dim}D)
+  set(_msw_target multi-body-problem-${_dim}D)
   set(_ugw_target uniform-grid-writer-${_dim}D)
 
   add_executable(${_sw_target} ${SW_SOURCE_FILES})
+  add_executable(${_msw_target} ${MSW_SOURCE_FILES})
   add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
   add_dune_ug_flags(${_sw_target})
+  add_dune_ug_flags(${_msw_target})
   add_dune_ug_flags(${_ugw_target})
 
   add_dune_hdf5_flags(${_sw_target})
+  add_dune_hdf5_flags(${_msw_target})
 
   set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
   set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/assemblers_tmpl.cc

@@ -4,4 +4,4 @@
 
 #include "explicitgrid.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+template class MyAssembler<DeformedGrid::LeafGridView, MY_DIM>;

src/explicitgrid.hh

@@ -2,6 +2,13 @@
 #define SRC_EXPLICITGRID_HH
 
 #include "gridselector.hh"
+#include "explicitvectors.hh"
 
-using GridView = Grid::LeafGridView;
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
+using LeafGridView = Grid::LeafGridView;
+using LevelGridView = Grid::LevelGridView;
+
+using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
+using DeformedGrid = DeformedGridComplex::DeformedGridType;
 #endif

src/hdf5-writer.hh

@@ -20,31 +20,54 @@ class HDF5Writer {
   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),
+  HDF5Writer(HDF5::Grouplike &file,
+             const std::vector<Vector>& vertexCoordinates,
+             const std::vector<const VertexBasis* >& vertexBases,
+             const std::vector<Patch>& surfaces,
+             const std::vector<Patch>& frictionalBoundaries,
+             const std::vector<ConvexPolyhedron<LocalVector>>& weakPatches)
+      : bodyCount_(vertexCoordinates.size()),
+        file_(file),
         iterationWriter_(file_),
         timeWriter_(file_),
 #if MY_DIM == 3
-        patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
+        patchInfoWriters_(bodyCount_),
 #endif
-        surfaceWriter_(file_, vertexCoordinates, surface),
-        frictionalBoundaryWriter_(file_, vertexCoordinates,
-                                  frictionalBoundary) {
+        surfaceWriters_(bodyCount_),
+        frictionalBoundaryWriters_(bodyCount_) {
+
+    for (size_t i=0; i<bodyCount_; i++) {
+#if MY_DIM == 3
+      patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, *vertexBases[i], frictionalBoundaries[i], weakPatches[i], i);
+#endif
+      surfaceWriters_[i] = new SurfaceWriter<ProgramState, GridView>(file_, vertexCoordinates[i], surfaces[i], i);
+      frictionalBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates[i], frictionalBoundaries[i], i);
+    }
+  }
+
+  ~HDF5Writer() {
+    for (size_t i=0; i<bodyCount_; i++) {
+#if MY_DIM == 3
+      delete patchInfoWriters_[i];
+#endif
+      delete surfaceWriters_[i];
+      delete frictionalBoundaryWriters_[i];
+    }
   }
 
   template <class Friction>
   void reportSolution(ProgramState const &programState,
                       // for the friction coefficient
-                      Friction &friction) {
+                      std::vector<std::shared_ptr<Friction>>& friction) {
     timeWriter_.write(programState);
+
+    for (size_t i=0; i<bodyCount_; i++) {
 #if MY_DIM == 3
-    patchInfoWriter_.write(programState);
+      patchInfoWriters_[i]->write(programState);
 #endif
-    surfaceWriter_.write(programState);
-    frictionalBoundaryWriter_.write(programState, friction);
+      surfaceWriters_[i]->write(programState);
+      frictionalBoundaryWriters_[i]->write(programState, *friction[i]);
+    }
   }
 
   void reportIterations(ProgramState const &programState,
@@ -53,14 +76,16 @@ class HDF5Writer {
   }
 
 private:
+  const size_t bodyCount_;
+
   HDF5::Grouplike &file_;
 
   IterationWriter iterationWriter_;
   TimeWriter<ProgramState> timeWriter_;
 #if MY_DIM == 3
-  PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
+  std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_;
 #endif
-  SurfaceWriter<ProgramState, GridView> surfaceWriter_;
-  FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
+  std::vector<SurfaceWriter<ProgramState, GridView>* > surfaceWriters_;
+  std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionalBoundaryWriters_;
 };
 #endif

src/hdf5/frictionalboundary-writer.cc

@@ -10,8 +10,9 @@
 template <class ProgramState, class GridView>
 FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
     HDF5::Grouplike &file, Vector const &vertexCoordinates,
-    Patch const &frictionalBoundary)
-    : group_(file, "frictionalBoundary"),
+    Patch const &frictionalBoundary, size_t const id)
+    : id_(id),
+      group_(file, "frictionalBoundary"+std::to_string(id_)),
       frictionalBoundary_(frictionalBoundary),
       frictionalBoundaryDisplacementWriter_(group_, "displacement",
                                             frictionalBoundary.numVertices(),
@@ -35,23 +36,24 @@ 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_);
+      restrictToSurface(programState.u[id_], frictionalBoundary_);
   addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
            frictionalBoundaryDisplacements);
 
   auto const frictionalBoundaryVelocities =
-      restrictToSurface(programState.v, frictionalBoundary_);
+      restrictToSurface(programState.v[id_], frictionalBoundary_);
   addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
            frictionalBoundaryVelocities);
 
   auto const frictionalBoundaryStates =
-      restrictToSurface(programState.alpha, frictionalBoundary_);
+      restrictToSurface(programState.alpha[id_], frictionalBoundary_);
   addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
            frictionalBoundaryStates);
 
-  friction.updateAlpha(programState.alpha);
-  auto const c = friction.coefficientOfFriction(programState.v);
+  friction.updateAlpha(programState.alpha[id_]);
+  auto const c = friction.coefficientOfFriction(programState.v[id_]);
   auto const frictionalBoundaryCoefficient =
       restrictToSurface(c, frictionalBoundary_);
   addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,

src/hdf5/frictionalboundary-writer.hh

@@ -12,12 +12,14 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
 
 public:
   FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                           Patch const &frictionalBoundary);
+                           Patch const &frictionalBoundary, size_t const id);
 
   template <class Friction>
   void write(ProgramState const &programState, Friction &friction);
 
 private:
+  size_t const id_;
+
   HDF5::Group group_;
 
   Patch const &frictionalBoundary_;

src/hdf5/frictionalboundary-writer_tmpl.cc

@@ -6,6 +6,8 @@
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 using MyFriction = GlobalFriction<Matrix, Vector>;
 
-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
+template class FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>;
+template void FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>::write(
     MyProgramState const &programState, MyFriction &friction);
+
+

src/hdf5/patchinfo-writer.cc

@@ -68,8 +68,10 @@ 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"),
+    ConvexPolyhedron<LocalVector> const &weakPatch,
+    size_t const id)
+    : id_(id),
+      group_(file, "weakPatchGrid"+std::to_string(id_)),
       vertexBasis_(vertexBasis),
       gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
       weakPatchGridVelocityWriter_(
@@ -87,7 +89,8 @@ PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
 template <class ProgramState, class VertexBasis, class GridView>
 void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
     ProgramState const &programState) {
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+    // TODO
+  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]);
   auto const gridVelocity = gridEvaluator_.evaluate(velocity);
   addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
 }

src/hdf5/patchinfo-writer.hh

@@ -36,11 +36,14 @@ class PatchInfoWriter {
 public:
   PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
                   Patch const &frictionalBoundary,
-                  ConvexPolyhedron<LocalVector> const &weakPatch);
+                  ConvexPolyhedron<LocalVector> const &weakPatch,
+                  size_t const id);
 
   void write(ProgramState const &programState);
 
 private:
+  const size_t id_;
+
   HDF5::Group group_;
 
   VertexBasis const &vertexBasis_;

src/hdf5/patchinfo-writer_tmpl.cc

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

src/hdf5/restart-io.cc

@@ -5,26 +5,66 @@
 #include "restart-io.hh"
 
 template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
-    : displacementWriter_(group, "displacement", vertexCount,
-                          Vector::block_type::dimension),
-      velocityWriter_(group, "velocity", vertexCount,
-                      Vector::block_type::dimension),
-      accelerationWriter_(group, "acceleration", vertexCount,
-                          Vector::block_type::dimension),
-      stateWriter_(group, "state", vertexCount),
-      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
+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") {}
+      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) {
-  addEntry(displacementWriter_, programState.timeStep, programState.u);
-  addEntry(velocityWriter_, programState.timeStep, programState.v);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha);
-  addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress);
+  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,
@@ -33,14 +73,15 @@ void RestartIO<ProgramState>::write(ProgramState const &programState) {
 
 template <class ProgramState>
 void RestartIO<ProgramState>::read(size_t timeStep,
-                                   ProgramState &programState) {
+                                   ProgramState& programState) {
+  assert(programState.size() == bodyIO_.size());
+
   programState.timeStep = timeStep;
-  readEntry(displacementWriter_, timeStep, programState.u);
-  readEntry(velocityWriter_, timeStep, programState.v);
-  readEntry(accelerationWriter_, timeStep, programState.a);
-  readEntry(stateWriter_, timeStep, programState.alpha);
-  readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress);
+
+  for (size_t i=0; i<bodyIO_.size(); i++) {
+    bodyIO_[i]->read(timeStep, programState);
+  }
+
   readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
   readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
 }

src/hdf5/restart-io.hh

@@ -4,23 +4,37 @@
 #include <dune/fufem/hdf5/file.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>
 
-template <class ProgramState> class RestartIO {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-
+template <class ProgramState> class RestartBodyIO {
 public:
-  RestartIO(HDF5::Grouplike &group, size_t vertexCount);
+  RestartBodyIO(HDF5::Grouplike& group, const size_t vertexCount, const size_t id);
 
-  void write(ProgramState const &programState);
+  void write(const ProgramState& programState);
 
-  void read(size_t timeStep, 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_;
+};
+
+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_;
 };

src/hdf5/restart-io_tmpl.cc

@@ -4,4 +4,5 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
+template class RestartBodyIO<MyProgramState>;
 template class RestartIO<MyProgramState>;

src/hdf5/surface-writer.cc

@@ -7,8 +7,9 @@
 
 template <class ProgramState, class GridView>
 SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
-    HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface)
-    : group_(file, "surface"),
+    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),
@@ -24,11 +25,12 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
 template <class ProgramState, class GridView>
 void SurfaceWriter<ProgramState, GridView>::write(
     ProgramState const &programState) {
-  auto const surfaceDisplacements = restrictToSurface(programState.u, surface_);
+
+  auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_);
   addEntry(surfaceDisplacementWriter_, programState.timeStep,
            surfaceDisplacements);
 
-  auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
+  auto const surfaceVelocities = restrictToSurface(programState.v[id_], surface_);
   addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
 }
 

src/hdf5/surface-writer.hh

@@ -12,11 +12,13 @@ template <class ProgramState, class GridView> class SurfaceWriter {
 
 public:
   SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                Patch const &surface);
+                Patch const &surface, size_t const id);
 
   void write(ProgramState const &programState);
 
 private:
+  size_t const id_;
+
   HDF5::Group group_;
 
   Patch const &surface_;

src/hdf5/surface-writer_tmpl.cc

@@ -5,4 +5,4 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
-template class SurfaceWriter<MyProgramState, GridView>;
+template class SurfaceWriter<MyProgramState, DeformedGrid::LeafGridView>;