src/time-stepping/rate/backward_euler.hh

 #ifndef SRC_TIME_STEPPING_RATE_BACKWARD_EULER_HH
 #define SRC_TIME_STEPPING_RATE_BACKWARD_EULER_HH
 
-template <class Vector, class Matrix, class Function, size_t dim>
-class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class BackwardEuler : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
 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,2> &_matrices,
+            const std::vector<Vector> &_u_initial,
+            const std::vector<Vector> &_v_initial,
+            const std::vector<Vector> &_a_initial,
+            const BoundaryNodes& _dirichletNodes,
+            const BoundaryFunctions& _dirichletFunctions);
 
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
+    virtual void setup(
+            const std::vector<Vector>&,
+            double,
+            double,
+            std::vector<Vector>&,
+            std::vector<Vector>&,
+            std::vector<Matrix>&) override;
 
-  std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
-      const override;
+    virtual void velocityObstacles(const Vector&, const Vector&, const Vector&, Vector&) override;
 
-private:
+    virtual void postProcess(const std::vector<Vector>&) override;
+
+    virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
 };
 #endif

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)
-    : RateUpdater<Vector, Matrix, Function, dim>(
-          _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
-          _dirichletFunction) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<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
-
-     Newmark means
-
-       a1 = 2/tau ( v1 - v0 ) - a0
-       u1 = tau/2 ( v1 + v0 ) + u0
-
-     in summary, we get at time t=1
-
-       M [2/tau ( u1 - u0 ) - a0] + C v1
-       + A [tau/2 ( v1 + v0 ) + u0] = ell
-
-     or
-
-       2/tau M v1 + C v1 + tau/2 A v1
-       = [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);
-  }
-  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 BoundaryFunctions, class BoundaryNodes>
+Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::Newmark(
+        const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+        const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+        const BoundaryNodes& _dirichletNodes,
+        const BoundaryFunctions& _dirichletFunctions) :
+    RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+        _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
+        _dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
+        const std::vector<Vector>& ell,
+        double _tau,
+        double relativeTime,
+        std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+        std::vector<Matrix>& AM) {
+
+    const size_t bodyCount = this->u.size();
+
+    rhs.resize(bodyCount);
+    iterate.resize(bodyCount);
+    AM.resize(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
+        this->tau = _tau;
+
+        /*  We start out with the formulation
+
+            M a + C v + A u = ell
+
+            Newmark means
+
+            a1 = 2/tau ( v1 - v0 ) - a0
+            u1 = tau/2 ( v1 + v0 ) + u0
+
+            in summary, we get at time t=1
+
+            M [2/tau ( u1 - u0 ) - a0] + C v1
+            + A [tau/2 ( v1 + v0 ) + u0] = ell
+
+            or
+
+            2/tau M v1 + C v1 + tau/2 A v1
+            = [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)
+        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;
+
+    for (size_t i=0; i<bodyCount; i++) {
+        auto& bodyIterate = iterate[i];
+
+        const auto& bodyDirichletNodes = this->dirichletNodes[i];
+        const auto& bodyDirichletFunctions = this->dirichletFunctions[i];
+        for (size_t bc=0; bc<bodyDirichletNodes.size(); ++bc) {
+            const auto& bcDirichletNodes = *bodyDirichletNodes[bc];
+            size_t dim = bcDirichletNodes[0].size();
+
+            double dirichletValue;
+            bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
+
+            std::cout << "dirichletValue: " << dirichletValue << std::endl;
+
+                  
+            for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
+                for (size_t j=0; j<dim; ++j) {
+                    if (bcDirichletNodes[k][j]) {
+                              std::cout << k << " " << j << std::endl;
+                        bodyIterate[k][j] = dirichletValue;
+                    }
+                }
+            }
+        }
+    }
+    print(iterate, "iterate:");
 }
 
-template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::postProcess(
-    Vector const &iterate) {
-  this->postProcessCalled = true;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) {
+    // v1 = 2/tau ( u1 - u0 ) - v0
+    v1Obstacles.resize(v0.size());
+    v1Obstacles = 0;
 
-  this->v = iterate;
+    Dune::MatrixVector::addProduct(v1Obstacles, 2.0 / this->tau, uObstacles);
+    Dune::MatrixVector::subtractProduct(v1Obstacles, 2.0 / this->tau, u0);
+    v1Obstacles -= v0;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::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);
+    // u1 = tau/2 ( v1 + v0 ) + u0
+    this->u = this->u_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>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>>
-Newmark<Vector, Matrix, Function, dim>::clone() const {
-  return std::make_shared<Newmark<Vector, Matrix, Function, dim>>(*this);
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::clone() const
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
+    return std::make_shared<Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(*this);
 }

src/time-stepping/rate/newmark.hh

 #ifndef SRC_TIME_STEPPING_RATE_NEWMARK_HH
 #define SRC_TIME_STEPPING_RATE_NEWMARK_HH
 
-template <class Vector, class Matrix, class Function, size_t dim>
-class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class Newmark : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
 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,2>& _matrices,
+            const std::vector<Vector>& _u_initial,
+            const std::vector<Vector>& _v_initial,
+            const std::vector<Vector>& _a_initial,
+            const BoundaryNodes& _dirichletNodes,
+            const BoundaryFunctions& _dirichletFunctions);
 
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
+    virtual void setup(
+            const std::vector<Vector>&,
+            double,
+            double,
+            std::vector<Vector>&,
+            std::vector<Vector>&,
+            std::vector<Matrix>&) override;
 
-  std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
-      const override;
+    virtual void velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) override;
+
+    virtual void postProcess(const std::vector<Vector>&) override;
+
+    virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
 };
 #endif

src/time-stepping/rate/rateupdater.cc

@@ -4,58 +4,59 @@
 
 #include "rateupdater.hh"
 
-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)
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::RateUpdater(
+    const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+    const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+    const BoundaryNodes& _dirichletNodes,
+    const BoundaryFunctions& _dirichletFunctions)
     : matrices(_matrices),
       u(_u_initial),
       v(_v_initial),
       a(_a_initial),
       dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
-  u_o = u;
-  v_o = v;
-  a_o = a;
-  postProcessCalled = false;
+      dirichletFunctions(_dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::nextTimeStep() {
+    u_o = u;
+    v_o = v;
+    a_o = a;
+    postProcessCalled = false;
 }
 
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(
-    Vector &displacement) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractDisplacement(std::vector<Vector>& displacements) const {
+    if (!postProcessCalled)
+        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+    displacements = u;
+}
 
-  displacement = u;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractOldDisplacement(std::vector<Vector>& displacements) const {
+    displacements = u_o;
 }
 
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
-    Vector &velocity) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractVelocity(std::vector<Vector>& velocity) const {
+    if (!postProcessCalled)
+        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 
-  velocity = v;
+    velocity = v;
 }
 
-template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(
-    Vector &oldVelocity) const {
-  oldVelocity = v_o;
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::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 {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractAcceleration(std::vector<Vector>& acceleration) const {
+    if (!postProcessCalled)
+        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 
-  acceleration = a;
+    acceleration = a;
 }
 
 #include "backward_euler.cc"

src/time-stepping/rate/rateupdater.hh

@@ -5,38 +5,46 @@
 
 #include <dune/common/bitsetvector.hh>
 
-#include "../../matrices.hh"
+#include "../../data-structures/matrices.hh"
 
-template <class Vector, class Matrix, class Function, size_t dim>
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
 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,2>& _matrices,
+            const std::vector<Vector>& _u_initial,
+            const std::vector<Vector>& _v_initial,
+            const std::vector<Vector>& _a_initial,
+            const BoundaryNodes& _dirichletNodes,
+            const BoundaryFunctions& _dirichletFunctions);
 
 public:
-  void nextTimeStep();
-  void virtual setup(Vector const &ell, double _tau, double relativeTime,
-                     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;
-
-  std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone()
-      const = 0;
+    void nextTimeStep();
+    void virtual setup(
+            const std::vector<Vector>& ell,
+            double _tau,
+            double relativeTime,
+            std::vector<Vector>& rhs,
+            std::vector<Vector>& iterate,
+            std::vector<Matrix>& AB) = 0;
+
+  void virtual postProcess(const std::vector<Vector>& iterate) = 0;
+  void virtual velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) = 0;
+  void extractDisplacement(std::vector<Vector>& displacements) const;
+  void extractOldDisplacement(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, BoundaryFunctions, BoundaryNodes>> virtual clone() const = 0;
 
 protected:
-  Matrices<Matrix> const &matrices;
-  Vector u, v, a;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
-  double dirichletValue;
+  const Matrices<Matrix,2>& matrices;
+  std::vector<Vector> u, v, a;
+  const BoundaryNodes& dirichletNodes;
+  const BoundaryFunctions& dirichletFunctions;
 
-  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/rateupdater_tmpl.cc

@@ -2,12 +2,14 @@
 #error MY_DIM unset
 #endif
 
-#include <dune/common/function.hh>
-
+#include "../../explicitgrid.hh"
 #include "../../explicitvectors.hh"
 
-using Function = Dune::VirtualFunction<double, double>;
+#include "../../data-structures/contactnetwork_tmpl.cc"
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
 
-template class RateUpdater<Vector, Matrix, Function, MY_DIM>;
-template class Newmark<Vector, Matrix, Function, MY_DIM>;
-template class BackwardEuler<Vector, Matrix, Function, MY_DIM>;
+template class RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;

src/time-stepping/rate_tmpl.cc

+#include "../data-structures/contactnetwork_tmpl.cc"
 #include "../explicitvectors.hh"
 
-#include <dune/common/function.hh>
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
 
-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);
+template
+auto initRateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+    Config::scheme scheme,
+    const BoundaryFunctions&  velocityDirichletFunctions,
+    const BoundaryNodes& velocityDirichletNodes,
+    const Matrices<Matrix, 2>& matrices,
+    const std::vector<Vector>& u_initial,
+    const std::vector<Vector>& v_initial,
+    const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>;

src/time-stepping/state.cc

@@ -6,20 +6,49 @@
 #include "state/ageinglawstateupdater.cc"
 #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) {
+template <class ScalarVector, class Vector, class ContactCoupling, class FrictionCouplingPair>
+auto initStateUpdater(
+        Config::stateModel model,
+        const std::vector<ScalarVector>& alpha_initial,
+        const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+        const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings) // contains frictionInfo
+-> std::shared_ptr<StateUpdater<ScalarVector, Vector>> {
+
+  assert(contactCouplings.size() == couplings.size());
+
+  auto stateUpdater = std::make_shared<StateUpdater<ScalarVector, Vector>>();
+
   switch (model) {
     case Config::AgeingLaw:
-      return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
-          alpha_initial, frictionalNodes, L, V0);
+      for (size_t i=0; i<couplings.size(); i++) {
+        const auto& coupling = couplings[i];
+        const auto nonmortarIdx = coupling->gridIdx_[0];
+
+        auto localUpdater =  std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
+                    alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+        localUpdater->setBodyIndex(nonmortarIdx);
+
+        stateUpdater->addLocalUpdater(localUpdater);
+      }
+      break;
     case Config::SlipLaw:
-      return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
-          alpha_initial, frictionalNodes, L, V0);
+      for (size_t i=0; i<couplings.size(); i++) {
+        const auto& coupling = couplings[i];
+        const auto nonmortarIdx = coupling->gridIdx_[0];
+
+        auto localUpdater =  std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
+                  alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+        localUpdater->setBodyIndex(nonmortarIdx);
+
+        stateUpdater->addLocalUpdater(localUpdater);
+      }
+      break;
     default:
       assert(false);
-  }
+      break;
+    }
+
+    return stateUpdater;
 }
 
 #include "state_tmpl.cc"

src/time-stepping/state.hh

@@ -3,16 +3,16 @@
 
 #include <memory>
 
-#include <dune/common/bitsetvector.hh>
-
-#include "../enums.hh"
+#include "../data-structures/enums.hh"
 #include "state/ageinglawstateupdater.hh"
 #include "state/sliplawstateupdater.hh"
 #include "state/stateupdater.hh"
 
-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);
-
+template <class ScalarVector, class Vector, class ContactCoupling, class FrictionCouplingPair>
+auto initStateUpdater(
+        Config::stateModel model,
+        const std::vector<ScalarVector>& alpha_initial,
+        const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+        const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings) // contains frictionInfo
+-> std::shared_ptr<StateUpdater<ScalarVector, Vector>>;
 #endif

src/time-stepping/state/ageinglawstateupdater.cc

 #include <cmath>
 
 #include "ageinglawstateupdater.hh"
-#include "../../tobool.hh"
+#include "../../utils/tobool.hh"
+#include "../../utils/debugutils.hh"
 
 template <class ScalarVector, class Vector>
 AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
-    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L, double _V0)
-    : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
+        const ScalarVector& alpha_initial,
+        const BitVector& nodes,
+        const double L,
+        const double V0) :
+    nodes_(nodes),
+    L_(L),
+    V0_(V0) {
+
+    localToGlobal_.resize(nodes_.count());
+    alpha_.resize(localToGlobal_.size());
+
+    size_t localIdx = 0;
+    for (size_t i=0; i<nodes_.size(); i++) {
+        if (not toBool(nodes_[i]))
+            continue;
+
+        localToGlobal_[localIdx] = i;
+        alpha_[localIdx] = alpha_initial[i];
+        localIdx++;
+    }
+}
 
 template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
-  alpha_o = alpha;
+  alpha_o_ = alpha_;
 }
 
 template <class ScalarVector, class Vector>
-void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
-  tau = _tau;
+void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double tau) {
+  tau_ = tau;
 }
 
 /*
   Compute [ 1-\exp(c*x) ] / x under the assumption that x is
   non-negative
 */
-double liftSingularity(double c, double x) {
+auto liftSingularity(
+        double c,
+        double x) {
+
   if (x <= 0)
     return -c;
   else
@@ -31,27 +53,32 @@ 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));
-  }
+void AgeingLawStateUpdater<ScalarVector, Vector>::solve(const Vector& velocity_field) {
+    for (size_t i=0; i<alpha_.size(); i++) {
+      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));
+    }
 }
 
 template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
-    ScalarVector &_alpha) {
-  _alpha = alpha;
+        ScalarVector& alpha) {
+
+    std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
+    if (alpha.size() != nodes_.size()) {
+        alpha.resize(nodes_.size());
+    }
+
+    for (size_t i=0; i<localToGlobal_.size(); i++) {
+        alpha[localToGlobal_[i]] = alpha_[i];
+    }
 }
 
 template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-AgeingLawStateUpdater<ScalarVector, Vector>::clone() const {
+auto AgeingLawStateUpdater<ScalarVector, Vector>::clone() const
+ -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> {
+
   return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(*this);
 }

src/time-stepping/state/ageinglawstateupdater.hh

 #ifndef SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
 #define SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
 
+#include <dune/common/bitsetvector.hh>
+
 #include "stateupdater.hh"
 
 template <class ScalarVector, class Vector>
-class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
+class AgeingLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
 public:
-  AgeingLawStateUpdater(ScalarVector const &_alpha_initial,
-                        Dune::BitSetVector<1> const &_nodes, double _L,
-                        double _V0);
+    AgeingLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
 
   void nextTimeStep() override;
-  void setup(double _tau) override;
-  void solve(Vector const &velocity_field) override;
-  void extractAlpha(ScalarVector &) override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
 
-  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
 
 private:
-  ScalarVector alpha_o;
-  ScalarVector alpha;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double const V0;
-  double tau;
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
 };
 #endif

src/time-stepping/state/sliplawstateupdater.cc

 #include <cmath>
 
 #include "sliplawstateupdater.hh"
-#include "../../tobool.hh"
+#include "../../utils/tobool.hh"
 
 template <class ScalarVector, class Vector>
 SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater(
-    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L, double _V0)
-    : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
+        const ScalarVector& alpha_initial,
+        const BitVector& nodes,
+        const double L,
+        const double V0) :
+    nodes_(nodes),
+    L_(L),
+    V0_(V0) {
+
+    localToGlobal_.resize(nodes_.count());
+    alpha_.resize(localToGlobal_.size());
+
+    size_t localIdx = 0;
+    for (size_t i=0; i<nodes_.size(); i++) {
+        if (not toBool(nodes_[i]))
+            continue;
+
+        localToGlobal_[localIdx] = i;
+        alpha_[localIdx] = alpha_initial[i];
+        localIdx++;
+    }
+}
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
-  alpha_o = alpha;
+    alpha_o_ = alpha_;
 }
 
 template <class ScalarVector, class Vector>
-void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
-  tau = _tau;
+void SlipLawStateUpdater<ScalarVector, Vector>::setup(double tau) {
+    tau_ = 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 Vector& velocity_field) {
+
+    for (size_t i=0; i<alpha_.size(); i++) {
+        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);
+    }
 }
 
 template <class ScalarVector, class Vector>
 void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
-    ScalarVector &_alpha) {
-  _alpha = alpha;
+        ScalarVector& alpha) {
+
+    assert(alpha.size() == nodes_.size());
+
+    for (size_t i=0; i<localToGlobal_.size(); i++) {
+        alpha[localToGlobal_[i]] = alpha_[i];
+    }
 }
 
 template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-SlipLawStateUpdater<ScalarVector, Vector>::clone() const {
+auto SlipLawStateUpdater<ScalarVector, Vector>::clone() const
+-> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> {
   return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(*this);
-}
+}                                                      

src/time-stepping/state/sliplawstateupdater.hh

 #ifndef SRC_TIME_STEPPING_STATE_SLIPLAWSTATEUPDATER_HH
 #define SRC_TIME_STEPPING_STATE_SLIPLAWSTATEUPDATER_HH
 
+#include <dune/common/bitsetvector.hh>
+
 #include "stateupdater.hh"
 
 template <class ScalarVector, class Vector>
-class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
+class SlipLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
 public:
-  SlipLawStateUpdater(ScalarVector const &_alpha_initial,
-                      Dune::BitSetVector<1> const &_nodes, double _L,
-                      double _V0);
+    SlipLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
 
   void nextTimeStep() override;
-  void setup(double _tau) override;
-  void solve(Vector const &velocity_field) override;
-  void extractAlpha(ScalarVector &) override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
 
-  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
 
 private:
-  ScalarVector alpha_o;
-  ScalarVector alpha;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double const V0;
-  double tau;
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
 };
 
 #endif

src/time-stepping/state/stateupdater.hh

 #ifndef SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
 #define SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
 
-template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
+#include <memory>
+#include <vector>
+
+// state updater for each coupling
+template <class ScalarVectorTEMPLATE, class Vector> class LocalStateUpdater {
 public:
   using ScalarVector = ScalarVectorTEMPLATE;
 
+  void setBodyIndex(const size_t bodyIdx) {
+      bodyIdx_ = bodyIdx;
+  }
+
+  auto bodyIndex() const {
+      return bodyIdx_;
+  }
+
   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 Vector& velocity_field) = 0;
+  void virtual extractAlpha(ScalarVector& alpha) = 0;
+
+  std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
 
-  std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
+private:
+  size_t bodyIdx_;
 };
 
+
+template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
+public:
+  using ScalarVector = ScalarVectorTEMPLATE;
+  using LocalStateUpdater = LocalStateUpdater<ScalarVector, Vector>;
+
+  void addLocalUpdater(std::shared_ptr<LocalStateUpdater> localStateUpdater) {
+    localStateUpdaters_.emplace_back(localStateUpdater);
+  }
+
+  void nextTimeStep() {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        localStateUpdaters_[i]->nextTimeStep();
+    }
+  }
+
+  void setup(double tau) {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        localStateUpdaters_[i]->setup(tau);
+    }
+  }
+
+  void solve(const std::vector<Vector>& velocity_field) {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        auto& localStateUpdater = localStateUpdaters_[i];
+        localStateUpdater->solve(velocity_field[localStateUpdater->bodyIndex()]);
+    }
+  }
+
+  void extractAlpha(std::vector<ScalarVector>& alpha) {
+    for (size_t i=0; i<localStateUpdaters_.size(); i++) {
+        auto& localStateUpdater = localStateUpdaters_[i];
+        localStateUpdater->extractAlpha(alpha[localStateUpdater->bodyIndex()]);
+    }
+  }
+
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const {
+      return std::make_shared<StateUpdater<ScalarVector, Vector>>(*this);
+  }
+
+private:
+  std::vector<std::shared_ptr<LocalStateUpdater>> localStateUpdaters_;
+};
 #endif

src/time-stepping/state_tmpl.cc

 #include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include <dune/common/promotiontraits.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
+#include "../frictioncouplingpair.hh"
+
+using field_type = typename Dune::PromotionTraits<typename Vector::field_type,
+                                            typename DeformedGrid::ctype>::PromotedType;
+
+using MyContactCoupling = Dune::Contact::DualMortarCoupling<field_type, DeformedGrid>;
+using MyFrictionCouplingPair = FrictionCouplingPair<DeformedGrid, LocalVector, field_type>;
 
 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);
+initStateUpdater<ScalarVector, Vector, MyContactCoupling, MyFrictionCouplingPair>(
+    Config::stateModel model,
+    const std::vector<ScalarVector>& alpha_initial,
+    const std::vector<std::shared_ptr<MyContactCoupling>>& contactCouplings,
+    const std::vector<std::shared_ptr<MyFrictionCouplingPair>>& couplings);

src/utils/almostequal.hh

+#ifndef SRC_ALMOSTEQUAL_HH
+#define SRC_ALMOSTEQUAL_HH
+
+#include <type_traits>
+#include <limits>
+#include <math.h>
+
+template <typename ctype>
+typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type almost_equal(ctype x, ctype y, int ulp) {
+    return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp || std::abs(x-y) < std::numeric_limits<ctype>::min();
+}
+
+#endif

src/utils/debugutils.hh

+#ifndef DEBUG_UTILS_
+#define DEBUG_UTILS_
+
+#include <string>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+//using namespace std;
+
+
+    template <int s>
+    void print(const Dune::BitSetVector<s>& x, std::string message){
+        std::cout << message << std::endl;
+        for (size_t i=0; i<x.size(); i++) {
+            std::cout << "block " << i << ": ";
+            for (size_t j=0; j<s; j++) {
+                 std::cout << x[i][j] << " ";
+            }
+            std::cout << std::endl;
+        }
+        std::cout << std::endl << std::endl;
+    }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::FieldVector<ctype, dim>& x, std::string message){
+       if (message != "")
+        std::cout << message << std::endl;
+
+       for (size_t i=0; i<x.size(); i++) {
+           std::cout << x[i] << ";" << std::endl;
+       }
+
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::BlockVector<Dune::FieldVector<ctype, dim>>& x, std::string message){
+       std::cout << message  << " " << "size " << x.size() << std::endl;
+       for (size_t i=0; i<x.size(); i++) {
+           print(x[i], "");
+       }
+       std::cout << std::endl << std::endl;
+   }
+
+   template <int dim, typename ctype=double>
+   void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, bool endOfLine = true){
+        if (message != "")
+            std::cout << message << std::endl;
+
+        for (size_t i=0; i<mat.N(); i++) {
+            for (size_t j=0; j<mat.M(); j++) {
+                if (mat.exists(i,j))
+                    std::cout << mat[i][j] << " ";
+                else
+                    std::cout << "0 ";
+            }
+
+            if (endOfLine)
+                std::cout << ";"<< std::endl;
+        }
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, const size_t line, bool endOfLine = true){
+        if (message != "")
+            std::cout << message << std::endl;
+
+        assert(line<mat.N());
+        for (size_t j=0; j<mat.M(); j++) {
+            if (mat.exists(line,j))
+                    std::cout << mat[line][j] << " ";
+                else
+                    std::cout << "0 ";
+            }
+
+            if (endOfLine)
+                std::cout << ";"<< std::endl;
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::BCRSMatrix<Dune::FieldMatrix<ctype, dim, dim>>& mat, std::string message){
+        std::cout << message << " " << "size (" << mat.N() << ", " << mat.M() << ")" <<std::endl;
+
+        Dune::FieldMatrix<ctype, dim, dim> zeroEntry = 0;
+
+        for (size_t i=0; i<mat.N(); i++) {
+            for (size_t d=0; d<dim; d++) {
+                for (size_t j=0; j<mat.M(); j++) {
+                    if (mat.exists(i,j))
+                        print(mat[i][j], "", d, j==mat.M()-1);
+                    else
+                        print(zeroEntry, "", d, j==mat.M()-1);
+                }
+            }
+       }
+       std::cout << std::endl;
+   }
+
+   template <class T=Dune::FieldMatrix<double,1,1>>
+   void print(const Dune::Matrix<T>& mat, std::string message){
+       std::cout << message << std::endl;
+       for (size_t i=0; i<mat.N(); i++) {
+           for (size_t j=0; j<mat.M(); j++) {
+                   std::cout << mat[i][j] << " ";
+           }
+           std::cout << std::endl;
+       }
+       std::cout << std::endl;
+   }
+
+   template <class T>
+   void print(const std::vector<T>& x, std::string message){
+      std::cout << message << std::endl;
+      for (size_t i=0; i<x.size(); i++) {
+          std::cout << x[i] << std::endl;
+      }
+      std::cout << std::endl << std::endl;
+   }
+
+ /*  void print(const std::vector<Dune::FieldVector<double,1>>& x, std::string message){
+      std::cout << message << std::endl;
+      for (size_t i=0; i<x.size(); i++) {
+          std::cout << x[i] << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }*/
+
+   template <class T>
+   void print(const std::set<T>& x, std::string message){
+      std::cout << message << std::endl;
+
+      for (const auto& c : x) {
+          std::cout << c << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }
+/*
+   void print(const std::set<int>& x, std::string message){
+      std::cout << message << std::endl;
+      std::set<int>::iterator dofIt = x.begin();
+      std::set<int>::iterator dofEndIt = x.end();
+      for (; dofIt != dofEndIt; ++dofIt) {
+          std::cout << *dofIt << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }
+
+   void step(const double stepNumber, std::string message=""){
+      std::cout << message << " Step " << stepNumber << "!" << std::endl;
+   }*/
+
+   template <class GridView, class VectorType>
+   void writeToVTK(const GridView& gridView, const VectorType& x, const std::string path, const std::string name) {
+       Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf(lStream.rdbuf());
+
+       vtkwriter.addVertexData(x,"data");
+       vtkwriter.pwrite(name, path, path);
+
+       std::cout.rdbuf( lBufferOld );
+   }
+
+   template <class GridView>
+   void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
+       Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf(lStream.rdbuf());
+
+       vtkwriter.pwrite(name, path, path);
+
+       std::cout.rdbuf( lBufferOld );
+   }
+
+   /*
+   template <class VectorType, class DGBasisType>
+   void writeToVTK(const DGBasisType& basis, const VectorType& x, const std::string path, const std::string name) {
+       Dune::VTKWriter<typename DGBasisType::GridView> vtkwriter(basis.getGridView());
+       VectorType toBePlotted(x);
+
+       toBePlotted.resize(basis.getGridView().indexSet().size(DGBasisType::GridView::Grid::dimension));
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf( lStream.rdbuf() );
+
+       vtkwriter.addVertexData(toBePlotted,"data");
+       vtkwriter.pwrite(name, path, path);
+
+       std::cout.rdbuf( lBufferOld );
+   }*/
+
+   template <class Intersection, class Basis>
+   void printIntersection(const Intersection& it, const Basis& basis, std::string message = "") {
+       if (message != "") {
+           std::cout << message << std::endl;
+       }
+
+       const auto& inside = it.inside();
+       const auto& localCoefficients = basis.getLocalFiniteElement(inside).localCoefficients();
+
+       std::cout << "dofs: ";
+       for (size_t i=0; i<localCoefficients.size(); i++) {
+           unsigned int entity = localCoefficients.localKey(i).subEntity();
+           unsigned int codim  = localCoefficients.localKey(i).codim();
+
+           if (it.containsInsideSubentity(entity, codim)) {
+               int dofIndex = basis.index(it.inside(), i);
+               std::cout << dofIndex << " ";
+               break;
+           }
+       }
+       std::cout << std::endl;
+   }
+
+   template <class BoundaryPatch, class Basis>
+   void printBoundary(const BoundaryPatch& patch, const Basis& basis, std::string message = "") {
+       if (message != "") {
+           std::cout << message << std::endl;
+       }
+
+       for (auto bIt = patch.begin(); bIt != patch.end(); ++bIt) {
+            printIntersection(*bIt, basis, "");
+       }
+
+       std::cout << std::endl;
+   }
+
+   template <class BasisType, typename ctype=double>
+   void printBasisDofLocation(const BasisType& basis) {
+     /*  typedef typename BasisType::GridView GridView;
+
+       const int dim = GridView::dimension;
+
+       std::map<int, int> indexTransformation;
+       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+       const GridView& gridView = basis.getGridView();
+       const int gridN = std::pow(gridView.indexSet().size(dim), 1.0/dim)-1;
+
+       for (auto& it : elements(gridView)) {
+           const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+           const auto& geometry = it.geometry();
+
+           for(int i=0; i<geometry.corners(); ++i) {
+               const auto& vertex = geometry.corner(i);
+               const auto& local = geometry.local(vertex);
+
+               int vertexIndex = vertex[0]*gridN;
+               for (int j=1; j<dim; ++j){
+                   vertexIndex += vertex[j]*gridN*std::pow(gridN+1, j);
+               }
+
+               const int localBasisSize = tFE.localBasis().size();
+               std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+               tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+               for(int k=0; k<localBasisSize; ++k) {
+                   if (localBasisRep[k]==1) {
+                       int dofIndex = basis.index(it, k);
+
+                       if (indexTransformation.count(dofIndex)==0) {
+                           indexTransformation[dofIndex] = vertexIndex;
+                           indexCoords[dofIndex] = vertex;
+                       }
+
+                       break;
+                   }
+               }
+           }
+       }
+
+       std::cout << std::endl;
+       std::cout << "Coarse level: Geometric vertex indices vs. associated basis dof indices " << indexTransformation.size() << std::endl;
+       std::map<int,int>::iterator mapIt = indexTransformation.begin();
+       std::map<int,int>::iterator mapEndIt = indexTransformation.end();
+       for (; mapIt!=mapEndIt; ++mapIt) {
+           std::cout << mapIt->second << " => " << mapIt->first << " Coords: ";
+
+           const Dune::FieldVector<ctype, dim>& coords = indexCoords[mapIt->first];
+           for (size_t i=0; i<coords.size(); i++) {
+               std::cout << coords[i] << " ";
+           }
+           std::cout << std::endl;
+       }
+       std::cout << std::endl;*/
+
+       const int dim = BasisType::GridView::dimension;
+
+       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+       const auto& gridView = basis.getGridView();
+
+       for (auto& it : elements(gridView)) {
+            const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+            const auto& geometry = it.geometry();
+
+            for(int i=0; i<geometry.corners(); ++i) {
+                const auto& vertex = geometry.corner(i);
+                const auto& local = geometry.local(vertex);
+
+                const int localBasisSize = tFE.localBasis().size();
+                std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+                tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+                for(int k=0; k<localBasisSize; ++k) {
+                    if (1.0 - localBasisRep[k]< 1e-14) {
+                        int dofIndex = basis.index(it, k);
+                        indexCoords[dofIndex] = vertex;
+                        break;
+                    }
+                }
+              }
+       }
+
+       std::cout << "Coords of basis dofs: " << indexCoords.size() << std::endl;
+       auto&& mapIt = indexCoords.begin();
+       const auto& mapEndIt = indexCoords.end();
+       for (; mapIt!=mapEndIt; ++mapIt) {
+            std::cout << mapIt->first << " Coords: ";
+
+            const auto& coords = mapIt->second;
+            for (size_t i=0; i<coords.size(); i++) {
+                std::cout << coords[i] << " ";
+            }
+            std::cout << std::endl;
+       }
+
+       std::cout << std::endl;
+   }
+
+   template <class GridView>
+   void printDofLocation(const GridView& gridView) {
+
+       std::cout << "-- GridView vertices --" << std::endl;
+       std::cout << "Index   Coords         " << std::endl;
+       std::cout << "-----------------------" << std::endl;
+
+       Dune::MultipleCodimMultipleGeomTypeMapper<
+           GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+       for (auto&& it : vertices(gridView)) {
+            const auto i = vertexMapper.index(it);
+            const auto& geometry = it.geometry();
+            const auto& corner = geometry.corner(0);
+
+            std::cout << std::setw(5) << i << "   ";
+            for(size_t i=0; i<corner.size(); ++i) {
+                std::cout << std::setw(5) << std::setprecision(3) << corner[i] << " ";
+            }
+            std::cout << std::endl;
+       }
+       std::cout << std::endl;
+   }
+
+   template <class Vector, class NBodyAssembler>
+   void printMortarBasis(const NBodyAssembler& nBodyAssembler) {
+       std::cout << "-- Mortar basis in canonical coords --" << std::endl;
+       std::cout << "--------------------------------------" << std::endl;
+
+       const auto& dim = nBodyAssembler.getTransformationMatrix().N();
+
+       Vector mortarBasisVector(dim);
+       std::vector<Vector> canonicalBasisVector(nBodyAssembler.nGrids());
+       for (size_t i=0; i<dim; i++) {
+            mortarBasisVector = 0;
+            mortarBasisVector[i] = 1;
+
+            nBodyAssembler.postprocess(mortarBasisVector, canonicalBasisVector);
+            print(canonicalBasisVector, "canonicalBasisVector " + std::to_string(i));
+       }
+       std::cout << std::endl;
+   }
+#endif

src/tobool.hh → src/utils/tobool.hh

@@ -4,7 +4,7 @@
 #include <dune/common/bitsetvector.hh>
 
 template <class Alloc>
-bool toBool(Dune::BitSetVectorConstReference<1, Alloc> x) {
+bool toBool(const Dune::BitSetVectorConstReference<1, Alloc> x) {
   return x[0];
 }
 #endif

todo.txt

+------------
+--  ToDo  --
+------------
+
+
+1. LevelContactNetwork
+
+    
+1. build n-body system 
+    Data structure: LevelContactNetwork
+    Factories:      StackedBlocksFactory
+    
+    - extend to multilevel LevelContactNetwork
+    - write new multilevel Cantor network factory
+    
+2. initialize/set up program state
+    Data structure: ProgramState
+    
+    - test setupInitialConditions()
+
+3. assemble RSD friction
+
+4. set up TNNMG solver
+    - rate updater
+    - state updater
+    
+5. adaptive time stepper
+
+