src/time-stepping/rate_tmpl.cc

@@ -6,7 +6,10 @@ 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, Function const &velocityDirichletFunction,
-    Dune::BitSetVector<MY_DIM> const &velocityDirichletNodes,
-    Matrices<Matrix> const &matrices, Vector const &u_initial,
-    Vector const &v_initial, Vector const &a_initial);
+    Config::scheme scheme,
+    const std::vector<const Function* >&  velocityDirichletFunctions,
+    const std::vector<Dune::BitSetVector<MY_DIM>>& velocityDirichletNodes,
+    const Matrices<Matrix>& matrices,
+    const std::vector<Vector>& u_initial,
+    const std::vector<Vector>& v_initial,
+    const std::vector<Vector>& a_initial);

src/time-stepping/state.cc

@@ -7,9 +7,8 @@
 #include "state/sliplawstateupdater.cc"
 
 template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
-    Config::stateModel model, ScalarVector const &alpha_initial,
-    Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) {
+std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
+    const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0) {
   switch (model) {
     case Config::AgeingLaw:
       return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(

src/time-stepping/state.hh

@@ -12,7 +12,7 @@
 
 template <class ScalarVector, class Vector>
 std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
-    Config::stateModel model, ScalarVector const &alpha_initial,
-    Dune::BitSetVector<1> const &frictionalNodes, double L, double V0);
+    Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
+    const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);
 
 #endif

src/time-stepping/state/ageinglawstateupdater.cc

@@ -5,8 +5,10 @@
 
 template <class ScalarVector, class Vector>
 AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
-    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L, double _V0)
+    const std::vector<ScalarVector>& _alpha_initial,
+    const std::vector<Dune::BitSetVector<1>>& _nodes,
+    const std::vector<double>& _L,
+    const std::vector<double>& _V0)
     : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
 
 template <class ScalarVector, class Vector>
@@ -32,21 +34,28 @@ double liftSingularity(double c, double x) {
 
 template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::solve(
-    Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i) {
-    if (not toBool(nodes[i]))
-      continue;
-
-    double const V = velocity_field[i].two_norm();
-    double const mtoL = -tau / L;
-    alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) +
-                        V0 * liftSingularity(mtoL, V));
+  const std::vector<Vector>& velocity_field) {
+  for (size_t i = 0; i < alpha.size(); ++i) {
+    const auto& velocity_fieldi = velocity_field[i];
+    const auto& nodesi = nodes[i];
+    auto& alphai = alpha[i];
+    auto& alpha_oi = alpha_o[i];
+
+    for (size_t j = 0; j < nodesi.size(); ++j) {
+      if (not toBool(nodesi[j]))
+        continue;
+
+      double const V = velocity_fieldi[j].two_norm();
+      double const mtoL = -tau / L[i];
+      alphai[j] = std::log(std::exp(alpha_oi[j] + V * mtoL) +
+                        V0[i] * liftSingularity(mtoL, V));
+    }
   }
 }
 
 template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
-    ScalarVector &_alpha) {
+    std::vector<ScalarVector>& _alpha) {
   _alpha = alpha;
 }
 

src/time-stepping/state/ageinglawstateupdater.hh

@@ -6,23 +6,23 @@
 template <class ScalarVector, class Vector>
 class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
 public:
-  AgeingLawStateUpdater(ScalarVector const &_alpha_initial,
-                        Dune::BitSetVector<1> const &_nodes, double _L,
-                        double _V0);
+  AgeingLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
+                        const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
+                        const std::vector<double>& _V0);
 
   void nextTimeStep() override;
   void setup(double _tau) override;
-  void solve(Vector const &velocity_field) override;
-  void extractAlpha(ScalarVector &) override;
+  void solve(const std::vector<Vector>& velocity_field) override;
+  void extractAlpha(std::vector<ScalarVector> &) override;
 
   std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
 
 private:
-  ScalarVector alpha_o;
-  ScalarVector alpha;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double const V0;
+  std::vector<ScalarVector> alpha_o;
+  std::vector<ScalarVector> alpha;
+  const std::vector<Dune::BitSetVector<1>>& nodes;
+  const std::vector<double>& L;
+  const std::vector<double>& V0;
   double tau;
 };
 #endif

src/time-stepping/state/sliplawstateupdater.cc

@@ -5,8 +5,9 @@
 
 template <class ScalarVector, class Vector>
 SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
-    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L, double _V0)
+    const std::vector<ScalarVector>& _alpha_initial,
+    const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
+    const std::vector<double>& _V0)
     : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
 
 template <class ScalarVector, class Vector>
@@ -21,21 +22,29 @@ void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::solve(
-    Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i) {
-    if (not toBool(nodes[i]))
-      continue;
-
-    double const V = velocity_field[i].two_norm();
-    double const mtVoL = -tau * V / L;
-    alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) +
-                                           alpha_o[i] * std::exp(mtVoL);
-  }
+  const std::vector<Vector>& velocity_field) {
+
+  for (size_t i = 0; i < alpha.size(); ++i) {
+      const auto& velocity_fieldi = velocity_field[i];
+      const auto& nodesi = nodes[i];
+      auto& alphai = alpha[i];
+      auto& alpha_oi = alpha_o[i];
+
+      for (size_t j = 0; j < nodesi.size(); ++j) {
+        if (not toBool(nodesi[j]))
+          continue;
+
+        double const V = velocity_fieldi[j].two_norm();
+        double const mtoL = -tau * V / L[i];
+        alphai[j] = (V <= 0) ? alpha_oi[j] : std::expm1(mtoL) * std::log(V / V0[i]) +
+                                               alpha_oi[j] * std::exp(mtoL);
+      }
+   }
 }
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
-    ScalarVector &_alpha) {
+    std::vector<ScalarVector>& _alpha) {
   _alpha = alpha;
 }
 

src/time-stepping/state/sliplawstateupdater.hh

@@ -6,24 +6,24 @@
 template <class ScalarVector, class Vector>
 class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
 public:
-  SlipLawStateUpdater(ScalarVector const &_alpha_initial,
-                      Dune::BitSetVector<1> const &_nodes, double _L,
-                      double _V0);
+    SlipLawStateUpdater(const std::vector<ScalarVector>& _alpha_initial,
+                          const std::vector<Dune::BitSetVector<1>>& _nodes, const std::vector<double>& _L,
+                          const std::vector<double>& _V0);
 
-  void nextTimeStep() override;
-  void setup(double _tau) override;
-  void solve(Vector const &velocity_field) override;
-  void extractAlpha(ScalarVector &) override;
+    void nextTimeStep() override;
+    void setup(double _tau) override;
+    void solve(const std::vector<Vector>& velocity_field) override;
+    void extractAlpha(std::vector<ScalarVector> &) override;
 
-  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+    std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
 
 private:
-  ScalarVector alpha_o;
-  ScalarVector alpha;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double const V0;
-  double tau;
+    std::vector<ScalarVector> alpha_o;
+    std::vector<ScalarVector> alpha;
+    const std::vector<Dune::BitSetVector<1>>& nodes;
+    const std::vector<double>& L;
+    const std::vector<double>& V0;
+    double tau;
 };
 
 #endif

src/time-stepping/state/stateupdater.hh

@@ -7,8 +7,8 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
 
   void virtual nextTimeStep() = 0;
   void virtual setup(double _tau) = 0;
-  void virtual solve(Vector const &velocity_field) = 0;
-  void virtual extractAlpha(ScalarVector &alpha) = 0;
+  void virtual solve(const std::vector<Vector>& velocity_field) = 0;
+  void virtual extractAlpha(std::vector<ScalarVector>& alpha) = 0;
 
   std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
 };

src/time-stepping/state_tmpl.cc

@@ -2,5 +2,5 @@
 
 template std::shared_ptr<StateUpdater<ScalarVector, Vector>>
 initStateUpdater<ScalarVector, Vector>(
-    Config::stateModel model, ScalarVector const &alpha_initial,
-    Dune::BitSetVector<1> const &frictionalNodes, double L, double V0);
+    Config::stateModel model, const std::vector<ScalarVector>& alpha_initial,
+        const std::vector<Dune::BitSetVector<1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);

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<GridView, double>;
-using P1Basis = P1NodalBasis<GridView, 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;