From 1557d88dec65f13cde45c8e568f2ac002ad08617 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Mon, 26 Mar 2018 16:05:57 +0200
Subject: [PATCH] .

---
 src/hdf5/patchinfo-writer_tmpl.cc |  9 ++--
 src/hdf5/restart-io.cc            | 89 ++++++++++++++++++++++---------
 src/hdf5/restart-io.hh            | 29 ++++------
 src/hdf5/restart-io_tmpl.cc       |  1 +
 src/multi-body-problem.cc         | 17 +++---
 src/program_state.hh              | 14 ++---
 src/time-stepping/rate_tmpl.cc    |  2 +-
 7 files changed, 99 insertions(+), 62 deletions(-)

diff --git a/src/hdf5/patchinfo-writer_tmpl.cc b/src/hdf5/patchinfo-writer_tmpl.cc
index d2f3c014..16bbf512 100644
--- a/src/hdf5/patchinfo-writer_tmpl.cc
+++ b/src/hdf5/patchinfo-writer_tmpl.cc
@@ -4,11 +4,12 @@
 #include "../program_state.hh"
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<LeafGridView, double>;
+using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
 using MyFunction = BasisGridFunction<P1Basis, Vector>;
 
-template class GridEvaluator<LocalVector, LeafGridView>;
+template class GridEvaluator<LocalVector, DeformedGrid::LeafGridView>;
+
 template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, LeafGridView>::evaluate(MyFunction const &f) const;
+GridEvaluator<LocalVector, DeformedGrid::LeafGridView>::evaluate(MyFunction const &f) const;
 
-template class PatchInfoWriter<MyProgramState, P1Basis, LeafGridView>;
+template class PatchInfoWriter<MyProgramState, P1Basis, DeformedGrid::LeafGridView>;
diff --git a/src/hdf5/restart-io.cc b/src/hdf5/restart-io.cc
index 0d16b86c..d1940959 100644
--- a/src/hdf5/restart-io.cc
+++ b/src/hdf5/restart-io.cc
@@ -5,27 +5,66 @@
 #include "restart-io.hh"
 
 template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, const std::vector<size_t>& vertexCounts)
-    : displacementWriter_(group, "displacement", vertexCounts,
-                          Vector::block_type::dimension),
-      velocityWriter_(group, "velocity", vertexCounts,
-                      Vector::block_type::dimension),
-      accelerationWriter_(group, "acceleration", vertexCounts,
-                          Vector::block_type::dimension),
-      stateWriter_(group, "state", vertexCounts),
-      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCounts),
+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) {
-    // TODO
-  addEntry(displacementWriter_, programState.timeStep, programState.u[0]);
-  addEntry(velocityWriter_, programState.timeStep, programState.v[0]);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a[0]);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha[0]);
-  addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress[0]);
+  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,
@@ -34,15 +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;
-  // TODO
-  readEntry(displacementWriter_, timeStep, programState.u[0]);
-  readEntry(velocityWriter_, timeStep, programState.v[0]);
-  readEntry(accelerationWriter_, timeStep, programState.a[0]);
-  readEntry(stateWriter_, timeStep, programState.alpha[0]);
-  readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress[0]);
+
+  for (size_t i=0; i<bodyIO_.size(); i++) {
+    bodyIO_[i]->read(timeStep, programState);
+  }
+
   readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
   readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
 }
diff --git a/src/hdf5/restart-io.hh b/src/hdf5/restart-io.hh
index 1cf5167a..6b228632 100644
--- a/src/hdf5/restart-io.hh
+++ b/src/hdf5/restart-io.hh
@@ -5,43 +5,36 @@
 #include <dune/fufem/hdf5/sequenceio.hh>
 
 template <class ProgramState> class RestartBodyIO {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-
 public:
-  RestartBodyIO(HDF5::Grouplike &group, const std::vector<size_t>& vertexCounts);
+  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_;
-  HDF5::SequenceIO<0> relativeTimeWriter_;
-  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
 };
 
 template <class ProgramState> class RestartIO {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
 
 public:
-  RestartIO(HDF5::Grouplike &group, const std::vector<size_t>& vertexCounts);
+  RestartIO(HDF5::Grouplike& group, const std::vector<size_t>& vertexCounts);
+  ~RestartIO();
 
-  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:
-  HDF5::SequenceIO<2> displacementWriter_;
-  HDF5::SequenceIO<2> velocityWriter_;
-  HDF5::SequenceIO<2> accelerationWriter_;
-  HDF5::SequenceIO<1> stateWriter_;
-  HDF5::SequenceIO<1> weightedNormalStressWriter_;
+  std::vector<RestartBodyIO<ProgramState>* > bodyIO_;
+
   HDF5::SequenceIO<0> relativeTimeWriter_;
   HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
 };
diff --git a/src/hdf5/restart-io_tmpl.cc b/src/hdf5/restart-io_tmpl.cc
index da72bb42..717acace 100644
--- a/src/hdf5/restart-io_tmpl.cc
+++ b/src/hdf5/restart-io_tmpl.cc
@@ -4,4 +4,5 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
+template class RestartBodyIO<MyProgramState>;
 template class RestartIO<MyProgramState>;
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 38971c5b..e0b3a60f 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -432,13 +432,16 @@ int main(int argc, char *argv[]) {
     }
 
 
-    /*
-    Vector vertexCoordinates(leafVertexCount);
-    {
+
+    std::vector<Vector> vertexCoordinates(leafVertexCounts.size());
+    for (size_t i=0; i<vertexCoordinates.size(); i++) {
+      auto& vertexCoords = vertexCoordinates[i];
+      vertexCoords.resize(leafVertexCounts[i]);
+
       Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView, Dune::mcmgVertexLayout());
-      for (auto &&v : vertices(leafView))
-        vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
+          DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(*leafViews[i], Dune::mcmgVertexLayout());
+      for (auto &&v : vertices(*leafViews[i]))
+        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
     }
 
     using MyVertexBasis = typename Assembler::VertexBasis;
@@ -451,7 +454,7 @@ int main(int argc, char *argv[]) {
     MyVTKWriter<MyVertexBasis, typename Assembler::CellBasis> const vtkWriter(
         myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
-    */
+
     IterationRegister iterationCount;
     /*
     auto const report = [&](bool initial = false) {
diff --git a/src/program_state.hh b/src/program_state.hh
index b48103cf..1aea46f7 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -22,8 +22,6 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
   using Vector = VectorTEMPLATE;
   using ScalarVector = ScalarVectorTEMPLATE;
 
-  BodyState() {}
-
   BodyState(Vector * _u, Vector * _v, Vector * _a, ScalarVector * _alpha, ScalarVector * _weightedNormalStress)
     : u(_u),
       v(_v),
@@ -64,11 +62,11 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       alpha[i].resize(leafVertexCount),
       weightedNormalStress[i].resize(leafVertexCount),
 
-      bodies[i] = BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
+      bodies[i] = new BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
     }
   }
 
-  size_t size() {
+  size_t size() const {
       return bodyCount_;
   }
 
@@ -183,8 +181,11 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
                        parset.sub("a0.solver"));
   }
 
+private:
+  const size_t bodyCount_;
+
 public:
-  std::vector<const BodyState> bodies;
+  std::vector<BodyState* > bodies;
   std::vector<Vector> u;
   std::vector<Vector> v;
   std::vector<Vector> a;
@@ -194,7 +195,6 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
   double relativeTau;
   size_t timeStep;
 
-private:
-  const size_t bodyCount_;
 
+};
 #endif
diff --git a/src/time-stepping/rate_tmpl.cc b/src/time-stepping/rate_tmpl.cc
index 18593e15..c81bc044 100644
--- a/src/time-stepping/rate_tmpl.cc
+++ b/src/time-stepping/rate_tmpl.cc
@@ -7,7 +7,7 @@ using Function = Dune::VirtualFunction<double, double>;
 template std::shared_ptr<RateUpdater<Vector, Matrix, Function, MY_DIM>>
 initRateUpdater<Vector, Matrix, Function, MY_DIM>(
     Config::scheme scheme,
-    const std::vector<Function>& velocityDirichletFunctions,
+    const std::vector<const Function* >&  velocityDirichletFunctions,
     const std::vector<Dune::BitSetVector<MY_DIM>>& velocityDirichletNodes,
     const Matrices<Matrix>& matrices,
     const std::vector<Vector>& u_initial,
-- 
GitLab