src/time-stepping/rate/rateupdater.hh

@@ -10,10 +10,10 @@
 template <class Vector, class Matrix, class Function, size_t dim>
 class RateUpdater {
 protected:
-  RateUpdater(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-              Vector const &_v_initial, Vector const &_a_initial,
-              Dune::BitSetVector<dim> const &_dirichletNodes,
-              Function const &_dirichletFunction);
+  RateUpdater(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+              const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+              const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+              const std::vector<Function>& _dirichletFunctions);
 
 public:
   void nextTimeStep();
@@ -21,22 +21,22 @@ class RateUpdater {
                      Vector &rhs, Vector &iterate, Matrix &AB) = 0;
 
   void virtual postProcess(Vector const &iterate) = 0;
-  void extractDisplacement(Vector &displacement) const;
-  void extractVelocity(Vector &velocity) const;
-  void extractOldVelocity(Vector &velocity) const;
-  void extractAcceleration(Vector &acceleration) const;
+  void extractDisplacement(std::vector<Vector>& displacements) const;
+  void extractVelocity(std::vector<Vector>& velocity) const;
+  void extractOldVelocity(std::vector<Vector>& velocity) const;
+  void extractAcceleration(std::vector<Vector>& acceleration) const;
 
   std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone()
       const = 0;
 
 protected:
-  Matrices<Matrix> const &matrices;
-  Vector u, v, a;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
+  const Matrices<Matrix>& matrices;
+  std::vector<Vector> u, v, a;
+  const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
+  const std::vector<Function>& dirichletFunctions;
   double dirichletValue;
 
-  Vector u_o, v_o, a_o;
+  std::vector<Vector> u_o, v_o, a_o;
   double tau;
 
   bool postProcessCalled = true;

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<Function>& velocityDirichletFunctions,
+    const std::vector<Dune::BitSetVector<dimension>>& 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<_Tp1>>& 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(mtVoL) * std::log(V / V0[i]) +
+                                               alpha_oi[j] * std::exp(mtVoL);
+      }
+   }
 }
 
 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<_Tp1>>& frictionalNodes, const std::vector<double>& L, const std::vector<double>& V0);