src/time-stepping/coupledtimestepper.cc

@@ -26,17 +26,17 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
   updaters_.rate_->nextTimeStep();
 
   auto const newRelativeTime = relativeTime + relativeTau;
-  Vector ell;
+  std::vector<Vector> ell;
   externalForces_(newRelativeTime, ell);
 
-  Matrix velocityMatrix;
-  Vector velocityRHS;
-  Vector velocityIterate;
+  std::vector<Matrix> velocityMatrix;
+  std::vector<Vector> velocityRHS;
+  std::vector<Vector> velocityIterate;
 
   auto const tau = relativeTau * finalTime_;
-  updaters_.state_->setup(tau);
-  updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS,
-                         velocityIterate, velocityMatrix);
+  updaters_.state_->setup(tau); 
+  updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
+
   FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
       factory_, parset_, globalFriction_, errorNorm_);
   auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,

src/time-stepping/rate.cc

@@ -9,20 +9,22 @@
 template <class Vector, class Matrix, class Function, int dimension>
 std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
 initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial) {
+                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) {
   switch (scheme) {
     case Config::Newmark:
       return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     case Config::BackwardEuler:
       return std::make_shared<
           BackwardEuler<Vector, Matrix, Function, dimension>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     default:
       assert(false);
   }

src/time-stepping/rate.hh

@@ -9,8 +9,10 @@
 template <class Vector, class Matrix, class Function, int dimension>
 std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
 initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial);
+                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);
 #endif

src/time-stepping/rate/backward_euler.cc

-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
 
 #include "backward_euler.hh"
 
 template <class Vector, class Matrix, class Function, size_t dim>
 BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
-    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-    Vector const &_v_initial, Vector const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    Function const &_dirichletFunction)
+    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)
     : RateUpdater<Vector, Matrix, Function, dim>(
           _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
-          _dirichletFunction) {}
+          _dirichletFunctions) {}
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::setup(
-    Vector const &ell, double _tau, double relativeTime, Vector &rhs,
-    Vector &iterate, Matrix &AM) {
-  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
-  this->tau = _tau;
-
-  /* We start out with the formulation
-
-       M a + C v + A u = ell
-
-     Backward Euler means
-
-       a1 = 1.0/tau ( v1 - v0 )
-       u1 = tau v1 + u0
-
-     in summary, we get at time t=1
-
-       M [1.0/tau ( v1 - v0 )] + C v1
-       + A [tau v1 + u0] = ell
-
-     or
-
-       1.0/tau M v1 + C v1 + tau A v1
-       = [1.0/tau M + C + tau A] v1
-       = ell + 1.0/tau M v0 - A u0
-  */
-
-  // set up LHS (for fixed tau, we'd only really have to do this once)
-  {
-    Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
-                                 this->matrices.elasticity.M());
-    indices.import(this->matrices.elasticity);
-    indices.import(this->matrices.mass);
-    indices.import(this->matrices.damping);
-    indices.exportIdx(AM);
+void BackwardEuler<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
+                                                         double _tau,
+                                                         double relativeTime,
+                                                         std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+                                                         std::vector<Matrix>& AM) {
+  for (size_t i=0; i<this->u.size(); i++) {
+    this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
+    this->tau = _tau;
+
+    /* We start out with the formulation
+
+           M a + C v + A u = ell
+
+         Backward Euler means
+
+           a1 = 1.0/tau ( v1 - v0 )
+           u1 = tau v1 + u0
+
+         in summary, we get at time t=1
+
+           M [1.0/tau ( v1 - v0 )] + C v1
+           + A [tau v1 + u0] = ell
+
+         or
+
+           1.0/tau M v1 + C v1 + tau A v1
+           = [1.0/tau M + C + tau A] v1
+           = ell + 1.0/tau M v0 - A u0
+    */
+
+    // set up LHS (for fixed tau, we'd only really have to do this once)
+    Matrix& LHS = AM[i];
+    {
+      Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
+                                   this->matrices.elasticity[i]->M());
+      indices.import(*this->matrices.elasticity[i]);
+      indices.import(*this->matrices.mass[i]);
+      indices.import(*this->matrices.damping[i]);
+      indices.exportIdx(LHS);
+    }
+    LHS = 0.0;
+    Dune::MatrixVector::addProduct(LHS, 1.0 / this->tau, *this->matrices.mass[i]);
+    Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+    Dune::MatrixVector::addProduct(LHS, this->tau, *this->matrices.elasticity[i]);
+
+    // set up RHS
+    {
+      Vector& rhss = rhs[i];
+      rhss = ell[i];
+      Dune::MatrixVector::addProduct(rhss, 1.0 / this->tau, *this->matrices.mass[i],
+                           this->v_o[i]);
+      Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+    }
+
+    iterate = this->v_o;
+
+    const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
+    for (size_t k = 0; k < dirichletNodess.size(); ++k)
+      for (size_t j = 0; j < dim; ++j)
+        if (this->dirichletNodes[k][j])
+          iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
   }
-  AM = 0.0;
-  Arithmetic::addProduct(AM, 1.0 / this->tau, this->matrices.mass);
-  Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
-  Arithmetic::addProduct(AM, this->tau, this->matrices.elasticity);
-
-  // set up RHS
-  {
-    rhs = ell;
-    Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
-                           this->v_o);
-    Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
-  }
-
-  iterate = this->v_o;
-
-  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
-    for (size_t j = 0; j < dim; ++j)
-      if (this->dirichletNodes[i][j])
-        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
-    Vector const &iterate) {
+    const std::vector<Vector>& iterate) {
   this->postProcessCalled = true;
 
   this->v = iterate;
 
   this->u = this->u_o;
-  Arithmetic::addProduct(this->u, this->tau, this->v);
 
-  this->a = this->v;
-  this->a -= this->v_o;
-  this->a /= this->tau;
+  for (size_t i=0; i<this->u.size(); i++) {
+    Dune::MatrixVector::addProduct(this->u[i], this->tau, this->v[i]);
+
+    Vector& ai = this->a[i];
+    ai = this->v[i];
+    ai -= this->v_o[i];
+    ai /= this->tau;
+  }
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>

src/time-stepping/rate/backward_euler.hh

@@ -4,14 +4,15 @@
 template <class Vector, class Matrix, class Function, size_t dim>
 class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
 public:
-  BackwardEuler(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-                Vector const &_v_initial, Vector const &_a_initial,
-                Dune::BitSetVector<dim> const &_dirichletNodes,
-                Function const &_dirichletFunction);
+  BackwardEuler(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);
 
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
+  void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
+             std::vector<Matrix>&) override;
+
+  void postProcess(const std::vector<Vector>&) override;
 
   std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
       const override;

src/time-stepping/rate/newmark.cc

-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
 
 #include "newmark.hh"
 
 template <class Vector, class Matrix, class Function, size_t dim>
 Newmark<Vector, Matrix, Function, dim>::Newmark(
-    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-    Vector const &_v_initial, Vector const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    Function const &_dirichletFunction)
+    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)
     : RateUpdater<Vector, Matrix, Function, dim>(
           _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
-          _dirichletFunction) {}
+          _dirichletFunctions) {}
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
+void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
                                                    double _tau,
                                                    double relativeTime,
-                                                   Vector &rhs, Vector &iterate,
-                                                   Matrix &AM) {
-  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
-  this->tau = _tau;
+                                                   std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+                                                   std::vector<Matrix>& AM) {
+  for (size_t i=0; i<this->u.size(); i++) {
+    this->dirichletFunctions[i].evaluate(relativeTime, this->dirichletValues[i]);
+    this->tau = _tau;
 
-  /* We start out with the formulation
+    /* We start out with the formulation
 
        M a + C v + A u = ell
 
@@ -41,58 +43,64 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
        = [2/tau M + C + tau/2 A] v1
        = ell + 2/tau M v0 + M a0
        - tau/2 A v0 - A u0
-  */
-
-  // set up LHS (for fixed tau, we'd only really have to do this once)
-  {
-    Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
-                                 this->matrices.elasticity.M());
-    indices.import(this->matrices.elasticity);
-    indices.import(this->matrices.mass);
-    indices.import(this->matrices.damping);
-    indices.exportIdx(AM);
+    */
+
+    // set up LHS (for fixed tau, we'd only really have to do this once)
+    Matrix& LHS = AM[i];
+    {
+      Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
+                                 this->matrices.elasticity[i]->M());
+      indices.import(*this->matrices.elasticity[i]);
+      indices.import(*this->matrices.mass[i]);
+      indices.import(*this->matrices.damping[i]);
+      indices.exportIdx(LHS);
+    }
+    LHS = 0.0;
+    Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
+    Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+    Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
+
+    // set up RHS
+    {
+      Vector& rhss = rhs[i];
+      rhss = ell[i];
+      Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
+                           this->v_o[i]);
+      Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
+      Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
+                                this->v_o[i]);
+      Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+    }
+
+    iterate = this->v_o;
+
+    const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
+    for (size_t k = 0; k < dirichletNodess.size(); ++k)
+      for (size_t j = 0; j < dim; ++j)
+        if (this->dirichletNodess[k][j])
+          iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
   }
-  AM = 0.0;
-  Arithmetic::addProduct(AM, 2.0 / this->tau, this->matrices.mass);
-  Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
-  Arithmetic::addProduct(AM, this->tau / 2.0, this->matrices.elasticity);
-
-  // set up RHS
-  {
-    rhs = ell;
-    Arithmetic::addProduct(rhs, 2.0 / this->tau, this->matrices.mass,
-                           this->v_o);
-    Arithmetic::addProduct(rhs, this->matrices.mass, this->a_o);
-    Arithmetic::subtractProduct(rhs, this->tau / 2.0, this->matrices.elasticity,
-                                this->v_o);
-    Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
-  }
-
-  iterate = this->v_o;
-
-  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
-    for (size_t j = 0; j < dim; ++j)
-      if (this->dirichletNodes[i][j])
-        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::postProcess(
-    Vector const &iterate) {
+void Newmark<Vector, Matrix, Function, dim>::postProcess(const std::vector<Vector>& iterate) {
   this->postProcessCalled = true;
 
   this->v = iterate;
 
   // u1 = tau/2 ( v1 + v0 ) + u0
   this->u = this->u_o;
-  Arithmetic::addProduct(this->u, this->tau / 2.0, this->v);
-  Arithmetic::addProduct(this->u, this->tau / 2.0, this->v_o);
-
-  // a1 = 2/tau ( v1 - v0 ) - a0
-  this->a = 0.0;
-  Arithmetic::addProduct(this->a, 2.0 / this->tau, this->v);
-  Arithmetic::subtractProduct(this->a, 2.0 / this->tau, this->v_o);
-  Arithmetic::subtractProduct(this->a, 1.0, this->a_o);
+
+  for (size_t i=0; i<this->u.size(); i++) {
+    Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
+    Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
+
+    // a1 = 2/tau ( v1 - v0 ) - a0
+    this->a[i] = 0.0;
+    Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
+    Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
+    Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
+  }
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>

src/time-stepping/rate/newmark.hh

@@ -4,14 +4,15 @@
 template <class Vector, class Matrix, class Function, size_t dim>
 class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
 public:
-  Newmark(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-          Vector const &_v_initial, Vector const &_a_initial,
-          Dune::BitSetVector<dim> const &_dirichletNodes,
-          Function const &_dirichletFunction);
+  Newmark(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);
 
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
+  void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
+             std::vector<Matrix>&) override;
+
+  void postProcess(const std::vector<Vector>&) override;
 
   std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
       const override;

src/time-stepping/rate/rateupdater.cc

@@ -6,16 +6,16 @@
 
 template <class Vector, class Matrix, class Function, size_t dim>
 RateUpdater<Vector, Matrix, Function, dim>::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)
+    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)
     : matrices(_matrices),
       u(_u_initial),
       v(_v_initial),
       a(_a_initial),
       dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
+      dirichletFunctions(_dirichletFunctions) {}
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
@@ -26,17 +26,15 @@ void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(
-    Vector &displacement) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(std::vector<Vector>& displacements) const {
   if (!postProcessCalled)
     DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 
-  displacement = u;
+  displacements = u;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
-    Vector &velocity) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(std::vector<Vector>& velocity) const {
   if (!postProcessCalled)
     DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 
@@ -44,14 +42,12 @@ void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(
-    Vector &oldVelocity) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
   oldVelocity = v_o;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(
-    Vector &acceleration) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(std::vector<Vector>& acceleration) const {
   if (!postProcessCalled)
     DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 

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);