.

src/hdf5-writer.hh

@@ -20,19 +20,39 @@ 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>
@@ -40,11 +60,14 @@ class HDF5Writer {
                       // for the friction coefficient
                       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);
+    }
   }
 
   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,27 +36,24 @@ template <class ProgramState, class GridView>
 template <class Friction>
 void FrictionalBoundaryWriter<ProgramState, GridView>::write(
     ProgramState const &programState, Friction &friction) {
-    // TODO
+
   auto const frictionalBoundaryDisplacements =
-      restrictToSurface(programState.u[0], frictionalBoundary_);
+      restrictToSurface(programState.u[id_], frictionalBoundary_);
   addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
            frictionalBoundaryDisplacements);
 
   auto const frictionalBoundaryVelocities =
-          // TODO
-      restrictToSurface(programState.v[0], frictionalBoundary_);
+      restrictToSurface(programState.v[id_], frictionalBoundary_);
   addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
            frictionalBoundaryVelocities);
 
   auto const frictionalBoundaryStates =
-      // TODO
-      restrictToSurface(programState.alpha[0], frictionalBoundary_);
+      restrictToSurface(programState.alpha[id_], frictionalBoundary_);
   addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
            frictionalBoundaryStates);
 
-  // TODO
-  friction.updateAlpha(programState.alpha[0]);
-  auto const c = friction.coefficientOfFriction(programState.v[0]);
+  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

+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+using MyFriction = GlobalFriction<Matrix, Vector>;
+
+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_(
@@ -88,7 +90,7 @@ template <class ProgramState, class VertexBasis, class GridView>
 void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
     ProgramState const &programState) {
     // TODO
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[0]);
+  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/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,13 +25,12 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
 template <class ProgramState, class GridView>
 void SurfaceWriter<ProgramState, GridView>::write(
     ProgramState const &programState) {
-    // TODO
-  auto const surfaceDisplacements = restrictToSurface(programState.u[0], surface_);
+
+  auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_);
   addEntry(surfaceDisplacementWriter_, programState.timeStep,
            surfaceDisplacements);
 
-  // TODO
-  auto const surfaceVelocities = restrictToSurface(programState.v[0], 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, LeafGridView>;
+template class SurfaceWriter<MyProgramState, DeformedGrid::LeafGridView>;

src/multi-body-problem.cc

@@ -432,9 +432,16 @@ int main(int argc, char *argv[]) {
     }
 
 
-
+    using MyVertexBasis = typename Assembler::VertexBasis;
+    using MyCellBasis = typename Assembler::CellBasis;
     std::vector<Vector> vertexCoordinates(leafVertexCounts.size());
+    std::vector<const MyVertexBasis* > vertexBases(leafVertexCounts.size());
+    std::vector<const MyCellBasis* > cellBases(leafVertexCounts.size());
+
     for (size_t i=0; i<vertexCoordinates.size(); i++) {
+      vertexBases[i] = &(assemblers[i]->vertexBasis);
+      cellBases[i] = &(assemblers[i]->cellBasis);
+
       auto& vertexCoords = vertexCoordinates[i];
       vertexCoords.resize(leafVertexCounts[i]);
 
@@ -444,16 +451,15 @@ int main(int argc, char *argv[]) {
         vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
     }
 
-    using MyVertexBasis = typename Assembler::VertexBasis;
+
     auto dataWriter =
         writeData ? std::make_unique<
-                        HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
-                        *dataFile, vertexCoordinates, myAssembler.vertexBasis,
-                        surface, frictionalBoundary, weakPatch)
+                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+                        *dataFile, vertexCoordinates, vertexBases,
+                        surfaces, frictionalBoundaries, weakPatches)
                   : nullptr;
-    MyVTKWriter<MyVertexBasis, typename Assembler::CellBasis> const vtkWriter(
-        myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
+    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "body");
 
     IterationRegister iterationCount;
     /*

src/program_state.hh

@@ -66,6 +66,12 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
     }
   }
 
+  ~ProgramState() {
+    for (size_t i=0; i<bodyCount_; i++) {
+      delete bodies[i];
+    }
+  }
+
   size_t size() const {
       return bodyCount_;
   }
@@ -85,7 +91,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       const Body<Vector::block_type::dimension>& body) {
 
     using LocalVector = typename Vector::block_type;
-    using LocalMatrix = typename Matrix::block_type;
+    //using LocalMatrix = typename Matrix::block_type;
     auto constexpr dims = LocalVector::dimension;
 
 

src/vtk.cc

@@ -9,50 +9,81 @@
 #include "vtk.hh"
 
 template <class VertexBasis, class CellBasis>
-MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter(
-    CellBasis const &_cellBasis, VertexBasis const &_vertexBasis,
-    std::string _prefix)
-    : cellBasis(_cellBasis), vertexBasis(_vertexBasis), prefix(_prefix) {}
+BodyVTKWriter<VertexBasis, CellBasis>::BodyVTKWriter(
+    CellBasis const & cellBasis, VertexBasis const & vertexBasis,
+    std::string prefix)
+    : cellBasis_(cellBasis), vertexBasis_(vertexBasis), prefix_(prefix) {}
 
 template <class VertexBasis, class CellBasis>
 template <class Vector, class ScalarVector>
-void MyVTKWriter<VertexBasis, CellBasis>::write(
+void BodyVTKWriter<VertexBasis, CellBasis>::write(
     size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha,
     ScalarVector const &stress) const {
   Dune::VTKWriter<typename VertexBasis::GridView> writer(
-      vertexBasis.getGridView());
+      vertexBasis_.getGridView());
 
   auto const displacementPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
-          vertexBasis, u, "displacement");
+          vertexBasis_, u, "displacement");
   writer.addVertexData(displacementPointer);
 
   auto const velocityPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
-          vertexBasis, v, "velocity");
+          vertexBasis_, v, "velocity");
   writer.addVertexData(velocityPointer);
 
   auto const AlphaPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
-          vertexBasis, alpha, "Alpha");
+          vertexBasis_, alpha, "Alpha");
   writer.addVertexData(AlphaPointer);
 
   auto const stressPointer =
       std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(
-          cellBasis, stress, "stress");
+          cellBasis_, stress, "stress");
   writer.addCellData(stressPointer);
 
-  std::string const filename = prefix + std::to_string(record);
+  std::string const filename = prefix_ + "_" + std::to_string(record);
   writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
 template <class VertexBasis, class CellBasis>
-void MyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
+void BodyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
   Dune::VTKWriter<typename VertexBasis::GridView> writer(
-      vertexBasis.getGridView());
+      vertexBasis_.getGridView());
 
-  std::string const filename = prefix + "_grid";
+  std::string const filename = prefix_ + "_grid";
   writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
+template <class VertexBasis, class CellBasis>
+MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
+  std::string prefix): bodyVTKWriters_(vertexBases.size()), prefix_(prefix) {
+
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    bodyVTKWriters_[i] = new BodyVTKWriter<VertexBasis, CellBasis>(*cellBases[i], *vertexBases[i], prefix_+std::to_string(i));
+  }
+}
+
+template <class VertexBasis, class CellBasis>
+MyVTKWriter<VertexBasis, CellBasis>::~MyVTKWriter() {
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    delete bodyVTKWriters_[i];
+  }
+}
+
+template <class VertexBasis, class CellBasis>
+template <class Vector, class ScalarVector>
+void MyVTKWriter<VertexBasis, CellBasis>::write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+                        const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const {
+
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    bodyVTKWriters_[i]->write(record, u[i], v[i], alpha[i], stress[i]);
+  }
+}
+template <class VertexBasis, class CellBasis>
+void MyVTKWriter<VertexBasis, CellBasis>::writeGrids() const {
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    bodyVTKWriters_[i]->writeGrid();
+  }
+}
 #include "vtk_tmpl.cc"

src/vtk.hh

@@ -3,13 +3,14 @@
 
 #include <string>
 
-template <class VertexBasis, class CellBasis> class MyVTKWriter {
-  CellBasis const &cellBasis;
-  VertexBasis const &vertexBasis;
-  std::string const prefix;
+template <class VertexBasis, class CellBasis> class BodyVTKWriter {
+private:
+  CellBasis const &cellBasis_;
+  VertexBasis const &vertexBasis_;
+  std::string const prefix_;
 
 public:
-  MyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis,
+  BodyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis,
               std::string prefix);
 
   template <class Vector, class ScalarVector>
@@ -18,4 +19,22 @@ template <class VertexBasis, class CellBasis> class MyVTKWriter {
 
   void writeGrid() const;
 };
+
+template <class VertexBasis, class CellBasis> class MyVTKWriter {
+private:
+  std::vector<BodyVTKWriter<VertexBasis, CellBasis>* > bodyVTKWriters_;
+  std::string const prefix_;
+
+public:
+  MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
+              std::string prefix);
+
+  ~MyVTKWriter();
+
+  template <class Vector, class ScalarVector>
+  void write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+             const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;
+
+  void writeGrids() const;
+};
 #endif

src/vtk_tmpl.cc

@@ -8,11 +8,11 @@
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
-using MyP0Basis = P0Basis<LeafGridView, double>;
-using P1Basis = P1NodalBasis<LeafGridView, double>;
+using MyP0Basis = P0Basis<DeformedGrid::LeafGridView, double>;
+using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
 
 template class MyVTKWriter<P1Basis, MyP0Basis>;
 
 template void MyVTKWriter<P1Basis, MyP0Basis>::write<Vector, ScalarVector>(
-    size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha,
-    ScalarVector const &stress) const;
+    size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+        const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;