dune/tectonic/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 BoundaryFunctions, class BoundaryNodes>
+class BackwardEuler : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
+public:
+    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);
+
+    virtual void setup(
+            const std::vector<Vector>&,
+            double,
+            double,
+            std::vector<Vector>&,
+            std::vector<Vector>&,
+            std::vector<Matrix>&) override;
+
+    virtual void velocityObstacles(const Vector&, const Vector&, const Vector&, Vector&) override;
+
+    virtual void postProcess(const std::vector<Vector>&) override;
+
+    virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
+};
+#endif

dune/tectonic/time-stepping/rate/newmark.cc

+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include "newmark.hh"
+
+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 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;
+
+    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;
+
+    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 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);
+}

dune/tectonic/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 BoundaryFunctions, class BoundaryNodes>
+class Newmark : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
+public:
+    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);
+
+    virtual void setup(
+            const std::vector<Vector>&,
+            double,
+            double,
+            std::vector<Vector>&,
+            std::vector<Vector>&,
+            std::vector<Matrix>&) 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

dune/tectonic/time-stepping/rate/rateupdater.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "rateupdater.hh"
+
+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),
+      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 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;
+}
+
+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 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;
+}
+
+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 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;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+const Matrices<Matrix,2>& RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::getMatrices() const {
+    return matrices;
+}
+
+
+#include "backward_euler.cc"
+#include "newmark.cc"
+#include "rateupdater_tmpl.cc"

dune/tectonic/time-stepping/rate/rateupdater.hh

+#ifndef SRC_TIME_STEPPING_RATE_RATEUPDATER_HH
+#define SRC_TIME_STEPPING_RATE_RATEUPDATER_HH
+
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+
+#include "../../data-structures/matrices.hh"
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class RateUpdater {
+protected:
+    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(
+            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;
+
+  const Matrices<Matrix,2>& getMatrices() const;
+
+  std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> virtual clone() const = 0;
+
+protected:
+  const Matrices<Matrix,2>& matrices;
+  std::vector<Vector> u, v, a;
+  const BoundaryNodes& dirichletNodes;
+  const BoundaryFunctions& dirichletFunctions;
+
+  std::vector<Vector> u_o, v_o, a_o;
+  double tau;
+
+  bool postProcessCalled = true;
+};
+#endif

dune/tectonic/time-stepping/rate/rateupdater_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+
+#include "../../data-structures/network/contactnetwork.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
+template class RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;

dune/tectonic/time-stepping/rate_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+#include "../data-structures/network/contactnetwork.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+
+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>>;

dune/tectonic/time-stepping/state.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "state.hh"
+#include "state/ageinglawstateupdater.cc"
+#include "state/sliplawstateupdater.cc"
+
+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());
+
+  std::vector<size_t> leafVertexCounts(alpha_initial.size(), 0);
+  for (size_t i=0; i<leafVertexCounts.size(); i++) {
+    leafVertexCounts[i] = alpha_initial[i].size();
+  }
+
+  auto stateUpdater = std::make_shared<StateUpdater<ScalarVector, Vector>>(leafVertexCounts);
+
+  switch (model) {
+    case Config::AgeingLaw:
+      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:
+      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"

dune/tectonic/time-stepping/state.hh

+#ifndef SRC_TIME_STEPPING_STATE_HH
+#define SRC_TIME_STEPPING_STATE_HH
+
+#include <memory>
+
+#include "../data-structures/enums.hh"
+#include "state/ageinglawstateupdater.hh"
+#include "state/sliplawstateupdater.hh"
+#include "state/stateupdater.hh"
+
+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

dune/tectonic/time-stepping/state/CMakeLists.txt

+add_custom_target(tectonic_dune_time-stepping_state SOURCES
+  ageinglawstateupdater.hh
+  ageinglawstateupdater.cc
+  sliplawstateupdater.hh
+  sliplawstateupdater.cc
+  stateupdater.hh
+)
+
+#install headers
+install(FILES
+  ageinglawstateupdater.hh
+  sliplawstateupdater.hh
+  stateupdater.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

src/time-stepping/state/ageinglawstateupdater.cc → dune/tectonic/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,36 @@ 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) {
+        auto tangentVelocity = velocity_field[localToGlobal_[i]];
+        tangentVelocity[0] = 0.0;
+
+        double const V = tangentVelocity.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);
 }

dune/tectonic/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 LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
+public:
+    AgeingLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
+
+  void nextTimeStep() override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
+
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
+
+private:
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
+};
+#endif

dune/tectonic/time-stepping/state/explicit.aux

+\relax 

dune/tectonic/time-stepping/state/explicit.log

+This is pdfTeX, Version 3.14159265-2.6-1.40.17 (TeX Live 2016/Debian) (preloaded format=pdflatex 2018.9.22)  2 SEP 2019 17:19
+entering extended mode
+ restricted \write18 enabled.
+ %&-line parsing enabled.
+**explicit.tex
+(./explicit.tex
+LaTeX2e <2017/01/01> patch level 3
+Babel <3.9r> and hyphenation patterns for 83 language(s) loaded.
+(/usr/share/texlive/texmf-dist/tex/latex/standalone/standalone.cls
+Document Class: standalone 2015/07/15 v1.2 Class to compile TeX sub-files stand
+alone
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifluatex.sty
+Package: ifluatex 2016/05/16 v1.4 Provides the ifluatex switch (HO)
+Package ifluatex Info: LuaTeX not detected.
+)
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifpdf.sty
+Package: ifpdf 2016/05/14 v3.1 Provides the ifpdf switch
+)
+(/usr/share/texlive/texmf-dist/tex/generic/ifxetex/ifxetex.sty
+Package: ifxetex 2010/09/12 v0.6 Provides ifxetex conditional
+)
+(/usr/share/texlive/texmf-dist/tex/latex/xkeyval/xkeyval.sty
+Package: xkeyval 2014/12/03 v2.7a package option processing (HA)
+
+(/usr/share/texlive/texmf-dist/tex/generic/xkeyval/xkeyval.tex
+(/usr/share/texlive/texmf-dist/tex/generic/xkeyval/xkvutils.tex
+\XKV@toks=\toks14
+\XKV@tempa@toks=\toks15
+
+(/usr/share/texlive/texmf-dist/tex/generic/xkeyval/keyval.tex))
+\XKV@depth=\count79
+File: xkeyval.tex 2014/12/03 v2.7a key=value parser (HA)
+))
+\sa@internal=\count80
+\c@sapage=\count81
+
+(/usr/share/texlive/texmf-dist/tex/latex/standalone/standalone.cfg
+File: standalone.cfg 2015/07/15 v1.2 Default configuration file for 'standalone
+' class
+)
+(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls
+Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
+(/usr/share/texlive/texmf-dist/tex/latex/base/size10.clo
+File: size10.clo 2014/09/29 v1.4h Standard LaTeX file (size option)
+)
+\c@part=\count82
+\c@section=\count83
+\c@subsection=\count84
+\c@subsubsection=\count85
+\c@paragraph=\count86
+\c@subparagraph=\count87
+\c@figure=\count88
+\c@table=\count89
+\abovecaptionskip=\skip41
+\belowcaptionskip=\skip42
+\bibindent=\dimen102
+)
+(/usr/share/texmf/tex/latex/preview/preview.sty
+Package: preview 2010/02/14 11.90 (AUCTeX/preview-latex)
+
+(/usr/share/texmf/tex/latex/preview/prtightpage.def
+\PreviewBorder=\dimen103
+)
+\pr@snippet=\count90
+\pr@box=\box26
+\pr@output=\toks16
+))
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty
+Package: amsmath 2016/11/05 v2.16a AMS math features
+\@mathmargin=\skip43
+
+For additional information on amsmath, use the `?' option.
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty
+Package: amstext 2000/06/29 v2.01 AMS text
+
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty
+File: amsgen.sty 1999/11/30 v2.0 generic functions
+\@emptytoks=\toks17
+\ex@=\dimen104
+))
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty
+Package: amsbsy 1999/11/29 v1.2d Bold Symbols
+\pmbraise@=\dimen105
+)
+(/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty
+Package: amsopn 2016/03/08 v2.02 operator names
+)
+\inf@bad=\count91
+LaTeX Info: Redefining \frac on input line 213.
+\uproot@=\count92
+\leftroot@=\count93
+LaTeX Info: Redefining \overline on input line 375.
+\classnum@=\count94
+\DOTSCASE@=\count95
+LaTeX Info: Redefining \ldots on input line 472.
+LaTeX Info: Redefining \dots on input line 475.
+LaTeX Info: Redefining \cdots on input line 596.
+\Mathstrutbox@=\box27
+\strutbox@=\box28
+\big@size=\dimen106
+LaTeX Font Info:    Redeclaring font encoding OML on input line 712.
+LaTeX Font Info:    Redeclaring font encoding OMS on input line 713.
+\macc@depth=\count96
+\c@MaxMatrixCols=\count97
+\dotsspace@=\muskip10
+\c@parentequation=\count98
+\dspbrk@lvl=\count99
+\tag@help=\toks18
+\row@=\count100
+\column@=\count101
+\maxfields@=\count102
+\andhelp@=\toks19
+\eqnshift@=\dimen107
+\alignsep@=\dimen108
+\tagshift@=\dimen109
+\tagwidth@=\dimen110
+\totwidth@=\dimen111
+\lineht@=\dimen112
+\@envbody=\toks20
+\multlinegap=\skip44
+\multlinetaggap=\skip45
+\mathdisplay@stack=\toks21
+LaTeX Info: Redefining \[ on input line 2817.
+LaTeX Info: Redefining \] on input line 2818.
+)
+(/usr/share/texlive/texmf-dist/tex/latex/mathtools/mathtools.sty
+Package: mathtools 2015/11/12 v1.18 mathematical typesetting tools
+
+(/usr/share/texlive/texmf-dist/tex/latex/tools/calc.sty
+Package: calc 2014/10/28 v4.3 Infix arithmetic (KKT,FJ)
+\calc@Acount=\count103
+\calc@Bcount=\count104
+\calc@Adimen=\dimen113
+\calc@Bdimen=\dimen114
+\calc@Askip=\skip46
+\calc@Bskip=\skip47
+LaTeX Info: Redefining \setlength on input line 80.
+LaTeX Info: Redefining \addtolength on input line 81.
+\calc@Ccount=\count105
+\calc@Cskip=\skip48
+)
+(/usr/share/texlive/texmf-dist/tex/latex/mathtools/mhsetup.sty
+Package: mhsetup 2010/01/21 v1.2a programming setup (MH)
+)
+LaTeX Info: Thecontrolsequence`\('isalreadyrobust on input line 129.
+LaTeX Info: Thecontrolsequence`\)'isalreadyrobust on input line 129.
+LaTeX Info: Thecontrolsequence`\['isalreadyrobust on input line 129.
+LaTeX Info: Thecontrolsequence`\]'isalreadyrobust on input line 129.
+\g_MT_multlinerow_int=\count106
+\l_MT_multwidth_dim=\dimen115
+\origjot=\skip49
+\l_MT_shortvdotswithinadjustabove_dim=\dimen116
+\l_MT_shortvdotswithinadjustbelow_dim=\dimen117
+\l_MT_above_intertext_sep=\dimen118
+\l_MT_below_intertext_sep=\dimen119
+\l_MT_above_shortintertext_sep=\dimen120
+\l_MT_below_shortintertext_sep=\dimen121
+)
+No file explicit.aux.
+\openout1 = `explicit.aux'.
+
+LaTeX Font Info:    Checking defaults for OML/cmm/m/it on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for T1/cmr/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for OT1/cmr/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for OMS/cmsy/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for OMX/cmex/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+LaTeX Font Info:    Checking defaults for U/cmr/m/n on input line 6.
+LaTeX Font Info:    ... okay on input line 6.
+Preview: Fontsize 10pt
+Preview: PDFoutput 1
+(/usr/share/texlive/texmf-dist/tex/latex/graphics/graphicx.sty
+Package: graphicx 2014/10/28 v1.0g Enhanced LaTeX Graphics (DPC,SPQR)
+
+(/usr/share/texlive/texmf-dist/tex/latex/graphics/graphics.sty
+Package: graphics 2016/10/09 v1.0u Standard LaTeX Graphics (DPC,SPQR)
+
+(/usr/share/texlive/texmf-dist/tex/latex/graphics/trig.sty
+Package: trig 2016/01/03 v1.10 sin cos tan (DPC)
+)
+(/usr/share/texlive/texmf-dist/tex/latex/graphics-cfg/graphics.cfg
+File: graphics.cfg 2016/06/04 v1.11 sample graphics configuration
+)
+Package graphics Info: Driver file: pdftex.def on input line 99.
+
+(/usr/share/texlive/texmf-dist/tex/latex/graphics-def/pdftex.def
+File: pdftex.def 2017/01/12 v0.06k Graphics/color for pdfTeX
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/infwarerr.sty
+Package: infwarerr 2016/05/16 v1.4 Providing info/warning/error messages (HO)
+)
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ltxcmds.sty
+Package: ltxcmds 2016/05/16 v1.23 LaTeX kernel commands for general use (HO)
+)
+\Gread@gobject=\count107
+
+(/usr/share/texlive/texmf-dist/tex/context/base/mkii/supp-pdf.mkii
+[Loading MPS to PDF converter (version 2006.09.02).]
+\scratchcounter=\count108
+\scratchdimen=\dimen122
+\scratchbox=\box29
+\nofMPsegments=\count109
+\nofMParguments=\count110
+\everyMPshowfont=\toks22
+\MPscratchCnt=\count111
+\MPscratchDim=\dimen123
+\MPnumerator=\count112
+\makeMPintoPDFobject=\count113
+\everyMPtoPDFconversion=\toks23
+))) (/usr/share/texlive/texmf-dist/tex/generic/oberdiek/pdftexcmds.sty
+Package: pdftexcmds 2016/05/21 v0.22 Utility functions of pdfTeX for LuaTeX (HO
+)
+Package pdftexcmds Info: LuaTeX not detected.
+Package pdftexcmds Info: \pdf@primitive is available.
+Package pdftexcmds Info: \pdf@ifprimitive is available.
+Package pdftexcmds Info: \pdfdraftmode found.
+)
+(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/epstopdf-base.sty
+Package: epstopdf-base 2016/05/15 v2.6 Base part for package epstopdf
+
+(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/grfext.sty
+Package: grfext 2016/05/16 v1.2 Manage graphics extensions (HO)
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/kvdefinekeys.sty
+Package: kvdefinekeys 2016/05/16 v1.4 Define keys (HO)
+))
+(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/kvoptions.sty
+Package: kvoptions 2016/05/16 v3.12 Key value format for package options (HO)
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/kvsetkeys.sty
+Package: kvsetkeys 2016/05/16 v1.17 Key value parser (HO)
+
+(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/etexcmds.sty
+Package: etexcmds 2016/05/16 v1.6 Avoid name clashes with e-TeX commands (HO)
+Package etexcmds Info: Could not find \expanded.
+(etexcmds)             That can mean that you are not using pdfTeX 1.50 or
+(etexcmds)             that some package has redefined \expanded.
+(etexcmds)             In the latter case, load this package earlier.
+)))
+Package epstopdf-base Info: Redefining graphics rule for `.eps' on input line 4
+38.
+Package grfext Info: Graphics extension search list:
+(grfext)             [.png,.pdf,.jpg,.mps,.jpeg,.jbig2,.jb2,.PNG,.PDF,.JPG,.JPE
+G,.JBIG2,.JB2,.eps]
+(grfext)             \AppendGraphicsExtensions on input line 456.
+
+(/usr/share/texlive/texmf-dist/tex/latex/latexconfig/epstopdf-sys.cfg
+File: epstopdf-sys.cfg 2010/07/13 v1.3 Configuration of (r)epstopdf for TeX Liv
+e
+))
+\Gin@req@height=\dimen124
+\Gin@req@width=\dimen125
+)
+Preview: Tightpage 0 0 0 0
+[1{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map}] (./explicit.aux) ) 
+Here is how much of TeX's memory you used:
+ 3166 strings out of 493013
+ 46673 string characters out of 6135682
+ 112782 words of memory out of 5000000
+ 6679 multiletter control sequences out of 15000+600000
+ 4094 words of font info for 16 fonts, out of 8000000 for 9000
+ 1141 hyphenation exceptions out of 8191
+ 53i,12n,56p,255b,147s stack positions out of 5000i,500n,10000p,200000b,80000s
+</usr/
+share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmex10.pfb></usr/share/
+texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi10.pfb></usr/share/texliv
+e/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi5.pfb></usr/share/texlive/texmf
+-dist/fonts/type1/public/amsfonts/cm/cmmi7.pfb></usr/share/texlive/texmf-dist/f
+onts/type1/public/amsfonts/cm/cmr10.pfb></usr/share/texlive/texmf-dist/fonts/ty
+pe1/public/amsfonts/cm/cmr5.pfb></usr/share/texlive/texmf-dist/fonts/type1/publ
+ic/amsfonts/cm/cmr7.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfo
+nts/cm/cmsy10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm
+/cmsy5.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmsy7.
+pfb>
+Output written on explicit.pdf (1 page, 82543 bytes).
+PDF statistics:
+ 48 PDF objects out of 1000 (max. 8388607)
+ 34 compressed objects within 1 object stream
+ 0 named destinations out of 1000 (max. 500000)
+ 1 words of extra memory for PDF output out of 10000 (max. 10000000)
+

src/time-stepping/state/sliplawstateupdater.cc → dune/tectonic/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) {
+        auto tangentVelocity = velocity_field[localToGlobal_[i]];
+        tangentVelocity[0] = 0.0;
+
+        double const V = tangentVelocity.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);
-}
+}                                                      

dune/tectonic/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 LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
+public:
+    SlipLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
+
+  void nextTimeStep() override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
+
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
+
+private:
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
+};
+
+#endif