src/timestepping/backward_euler.hh → src/time-stepping/rate/backward_euler.hh

-#ifndef DUNE_TECTONIC_TIMESTEPPING_BACKWARD_EULER_HH
-#define DUNE_TECTONIC_TIMESTEPPING_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 TimeSteppingScheme<Vector, Matrix, Function, dim> {
+class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
 public:
-  BackwardEuler(Matrix const &_A, Matrix const &_M, Matrix const &_C,
-                Vector const &_u_initial, Vector const &_v_initial,
+  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);
 
-  void nextTimeStep() override;
   void setup(Vector const &, double, double, Vector &, Vector &,
              Matrix &) override;
   void postProcess(Vector const &) override;
-  void extractDisplacement(Vector &) const override;
-  void extractVelocity(Vector &) const override;
-  void extractOldVelocity(Vector &) const override;
 
-private:
-  Matrix const &A;
-  Matrix const &M;
-  Matrix const &C;
-  Vector u;
-  Vector v;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
-
-  Vector u_o;
-  Vector v_o;
+  std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
+      const override;
 
-  double tau;
-
-  bool postProcessCalled = false;
+private:
 };
 #endif

src/timestepping/newmark.cc → src/time-stepping/rate/newmark.cc

+#include <dune/solvers/common/arithmetic.hh>
+
+#include "newmark.hh"
+
 template <class Vector, class Matrix, class Function, size_t dim>
 Newmark<Vector, Matrix, Function, dim>::Newmark(
-    Matrix const &_A, Matrix const &_M, Matrix const &_C,
-    Vector const &_u_initial, Vector const &_v_initial,
-    Vector const &_a_initial, Dune::BitSetVector<dim> const &_dirichletNodes,
+    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
+    Vector const &_v_initial, Vector const &_a_initial,
+    Dune::BitSetVector<dim> const &_dirichletNodes,
     Function const &_dirichletFunction)
-    : A(_A),
-      M(_M),
-      C(_C),
-      u(_u_initial),
-      v(_v_initial),
-      a(_a_initial),
-      dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::nextTimeStep() {
-  a_o = a;
-  v_o = v;
-  u_o = u;
-}
+    : 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,
@@ -26,9 +18,8 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
                                                    double relativeTime,
                                                    Vector &rhs, Vector &iterate,
                                                    Matrix &AM) {
-  postProcessCalled = false;
-
-  tau = _tau;
+  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
+  this->tau = _tau;
 
   /* We start out with the formulation
 
@@ -54,78 +45,58 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
 
   // set up LHS (for fixed tau, we'd only really have to do this once)
   {
-    Dune::MatrixIndexSet indices(A.N(), A.M());
-    indices.import(A);
-    indices.import(M);
-    indices.import(C);
+    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 / tau, M);
-  Arithmetic::addProduct(AM, 1.0, C);
-  Arithmetic::addProduct(AM, tau / 2.0, A);
+  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 / tau, M, v_o);
-    Arithmetic::addProduct(rhs, M, a_o);
-    Arithmetic::subtractProduct(rhs, tau / 2.0, A, v_o);
-    Arithmetic::subtractProduct(rhs, A, u_o);
+    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);
   }
 
-  // v_o makes a good initial iterate; we could use anything, though
-  iterate = 0.0;
+  iterate = this->v_o;
 
-  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
     for (size_t j = 0; j < dim; ++j)
-      if (dirichletNodes[i][j]) {
-        if (j == 0)
-          dirichletFunction.evaluate(relativeTime, iterate[i][j]);
-        else
-          iterate[i][j] = 0;
-      }
+      if (this->dirichletNodes[i][j])
+        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void Newmark<Vector, Matrix, Function, dim>::postProcess(
     Vector const &iterate) {
-  postProcessCalled = true;
+  this->postProcessCalled = true;
 
-  v = iterate;
+  this->v = iterate;
 
   // u1 = tau/2 ( v1 + v0 ) + u0
-  u = u_o;
-  Arithmetic::addProduct(u, tau / 2.0, v);
-  Arithmetic::addProduct(u, tau / 2.0, v_o);
+  this->u = this->u_o;
+  Arithmetic::addProduct(this->u, this->tau / 2.0, this->v);
+  Arithmetic::addProduct(this->u, this->tau / 2.0, this->v_o);
 
   // a1 = 2/tau ( v1 - v0 ) - a0
-  a = 0;
-  Arithmetic::addProduct(a, 2.0 / tau, v);
-  Arithmetic::subtractProduct(a, 2.0 / tau, v_o);
-  Arithmetic::subtractProduct(a, 1.0, a_o);
+  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);
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::extractDisplacement(
-    Vector &displacement) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  displacement = u;
-}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::extractVelocity(Vector &velocity)
-    const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  velocity = v;
-}
-
-template <class Vector, class MatrixType, class FunctionType, size_t dim>
-void Newmark<Vector, MatrixType, FunctionType, dim>::extractOldVelocity(
-    Vector &velocity) const {
-  velocity = v_o;
+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);
 }

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> {
+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);
+
+  void setup(Vector const &, double, double, Vector &, Vector &,
+             Matrix &) override;
+  void postProcess(Vector const &) override;
+
+  std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
+      const override;
+};
+#endif

src/time-stepping/rate/rateupdater.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#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)
+    : 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;
+}
+
+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!");
+
+  displacement = u;
+}
+
+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!");
+
+  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 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!");
+
+  acceleration = a;
+}
+
+#include "backward_euler.cc"
+#include "newmark.cc"
+#include "rateupdater_tmpl.cc"

src/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 "../../matrices.hh"
+
+template <class Vector, class Matrix, class Function, size_t dim>
+class RateUpdater {
+protected:
+  RateUpdater(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
+              Vector const &_v_initial, Vector const &_a_initial,
+              Dune::BitSetVector<dim> const &_dirichletNodes,
+              Function const &_dirichletFunction);
+
+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;
+
+protected:
+  Matrices<Matrix> const &matrices;
+  Vector u, v, a;
+  Dune::BitSetVector<dim> const &dirichletNodes;
+  Function const &dirichletFunction;
+  double dirichletValue;
+
+  Vector u_o, v_o, a_o;
+  double tau;
+
+  bool postProcessCalled = true;
+};
+#endif

src/time-stepping/rate/rateupdater_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include <dune/common/function.hh>
+
+#include "../../explicitvectors.hh"
+
+using Function = Dune::VirtualFunction<double, double>;
+
+template class RateUpdater<Vector, Matrix, Function, MY_DIM>;
+template class Newmark<Vector, Matrix, Function, MY_DIM>;
+template class BackwardEuler<Vector, Matrix, Function, MY_DIM>;

src/time-stepping/rate_tmpl.cc

+#include "../explicitvectors.hh"
+
+#include <dune/common/function.hh>
+
+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);

src/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>
+std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
+    Config::stateModel model, ScalarVector const &alpha_initial,
+    Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) {
+  switch (model) {
+    case Config::AgeingLaw:
+      return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
+          alpha_initial, frictionalNodes, L, V0);
+    case Config::SlipLaw:
+      return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
+          alpha_initial, frictionalNodes, L, V0);
+    default:
+      assert(false);
+  }
+}
+
+#include "state_tmpl.cc"

src/state.hh → src/time-stepping/state.hh

-#ifndef DUNE_TECTONIC_STATE_HH
-#define DUNE_TECTONIC_STATE_HH
+#ifndef SRC_TIME_STEPPING_STATE_HH
+#define SRC_TIME_STEPPING_STATE_HH
+
+#include <memory>
 
 #include <dune/common/bitsetvector.hh>
 
-#include "enums.hh"
+#include "../enums.hh"
+#include "state/ageinglawstateupdater.hh"
+#include "state/sliplawstateupdater.hh"
 #include "state/stateupdater.hh"
-#include "state/ruinastateupdater.hh"
-#include "state/dieterichstateupdater.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) {
-  switch (model) {
-    case Config::Dieterich:
-      return std::make_shared<DieterichStateUpdater<ScalarVector, Vector>>(
-          alpha_initial, frictionalNodes, L);
-    case Config::Ruina:
-      return std::make_shared<RuinaStateUpdater<ScalarVector, Vector>>(
-          alpha_initial, frictionalNodes, L);
-    default:
-      assert(false);
-  }
-}
+    Dune::BitSetVector<1> const &frictionalNodes, double L, double V0);
+
 #endif

src/time-stepping/state/ageinglawstateupdater.cc

+#include <cmath>
+
+#include "ageinglawstateupdater.hh"
+#include "../../tobool.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) {}
+
+template <class ScalarVector, class Vector>
+void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
+  alpha_o = alpha;
+}
+
+template <class ScalarVector, class Vector>
+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) {
+  if (x <= 0)
+    return -c;
+  else
+    return -std::expm1(c * x) / 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));
+  }
+}
+
+template <class ScalarVector, class Vector>
+void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
+    ScalarVector &_alpha) {
+  _alpha = alpha;
+}
+
+template <class ScalarVector, class Vector>
+std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+AgeingLawStateUpdater<ScalarVector, Vector>::clone() const {
+  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 "stateupdater.hh"
+
+template <class ScalarVector, class Vector>
+class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
+public:
+  AgeingLawStateUpdater(ScalarVector const &_alpha_initial,
+                        Dune::BitSetVector<1> const &_nodes, double _L,
+                        double _V0);
+
+  void nextTimeStep() override;
+  void setup(double _tau) override;
+  void solve(Vector const &velocity_field) override;
+  void extractAlpha(ScalarVector &) override;
+
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+
+private:
+  ScalarVector alpha_o;
+  ScalarVector alpha;
+  Dune::BitSetVector<1> const &nodes;
+  double const L;
+  double const V0;
+  double tau;
+};
+#endif

src/time-stepping/state/calculation.wxm

+/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
+/* [ Created with wxMaxima version 12.04.0 ] */
+
+/* [wxMaxima: input   start ] */
+d_eq    : 'diff(theta,t) = 1 - theta*(V/L);
+  ode2(d_eq, theta, t)$
+  subst(L/V[ref]*exp(alpha), theta, %)$
+  % / L * V[ref]$
+  ic1(%, t = 0, alpha = alpha[0])$
+d_alpha : log(expand(%));
+/* [wxMaxima: input   end   ] */
+
+/* [wxMaxima: input   start ] */
+r_eq : 'diff(theta,t) = -(V/L)*theta*log(theta*(V/L));
+  ode2(r_eq,theta,t)$
+  subst(L/V[ref] * exp(alpha),theta,%)$
+  ic1(%,t=0,alpha=alpha[0])$
+  map(exp,-(V/L)*%)$
+  block([logexpand : all], expand(% - log(V/V[ref])))$
+r_alpha : expand(logcontract(%));
+/* [wxMaxima: input   end   ] */
+
+/* Maxima can't load/batch files which end with a comment! */
+"Created with wxMaxima"$

src/time-stepping/state/explicit.tex

+\documentclass[preview]{standalone}
+\providecommand\Vref{V_\ast}
+\usepackage{amsmath}
+\usepackage{mathtools}
+\DeclarePairedDelimiter\paren{\lparen}{\rparen}
+\begin{document}
+\begin{equation*}
+  \alpha(t) =
+  \begin{cases*}
+    \log\paren*{e^{\alpha_0-tV/L} + \Vref \frac{1 - e^{-tV/L}}V}
+    &for the ageing law,\\
+    (e^{-\frac {tV}L} - 1) \log \frac V\Vref + \alpha_0 e^{-\frac {tV}L}
+    &for the slip law.
+  \end{cases*}
+\end{equation*}
+\end{document}

src/time-stepping/state/sliplawstateupdater.cc

+#include <cmath>
+
+#include "sliplawstateupdater.hh"
+#include "../../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) {}
+
+template <class ScalarVector, class Vector>
+void SlipLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
+  alpha_o = alpha;
+}
+
+template <class ScalarVector, class Vector>
+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);
+  }
+}
+
+template <class ScalarVector, class Vector>
+void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha(
+    ScalarVector &_alpha) {
+  _alpha = alpha;
+}
+
+template <class ScalarVector, class Vector>
+std::shared_ptr<StateUpdater<ScalarVector, Vector>>
+SlipLawStateUpdater<ScalarVector, Vector>::clone() const {
+  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 "stateupdater.hh"
+
+template <class ScalarVector, class Vector>
+class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
+public:
+  SlipLawStateUpdater(ScalarVector const &_alpha_initial,
+                      Dune::BitSetVector<1> const &_nodes, double _L,
+                      double _V0);
+
+  void nextTimeStep() override;
+  void setup(double _tau) override;
+  void solve(Vector const &velocity_field) override;
+  void extractAlpha(ScalarVector &) override;
+
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
+
+private:
+  ScalarVector alpha_o;
+  ScalarVector alpha;
+  Dune::BitSetVector<1> const &nodes;
+  double const L;
+  double const V0;
+  double tau;
+};
+
+#endif

src/state/stateupdater.hh → src/time-stepping/state/stateupdater.hh

-#ifndef STATE_UPDATER_HH
-#define STATE_UPDATER_HH
+#ifndef SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
+#define SRC_TIME_STEPPING_STATE_STATEUPDATER_HH
 
-template <class ScalarVector, class Vector> class StateUpdater {
+template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
 public:
+  using ScalarVector = ScalarVectorTEMPLATE;
+
   void virtual nextTimeStep() = 0;
   void virtual setup(double _tau) = 0;
   void virtual solve(Vector const &velocity_field) = 0;
-  void virtual extractLogState(ScalarVector &logState) = 0;
+  void virtual extractAlpha(ScalarVector &alpha) = 0;
+
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> virtual clone() const = 0;
 };
 
 #endif

src/time-stepping/state_tmpl.cc

+#include "../explicitvectors.hh"
+
+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);

src/time-stepping/updaters.hh

+#ifndef SRC_TIME_STEPPING_UPDATERS_HH
+#define SRC_TIME_STEPPING_UPDATERS_HH
+
+template <class RateUpdaterTEMPLATE, class StateUpdaterTEMPLATE>
+struct Updaters {
+  using RateUpdater = RateUpdaterTEMPLATE;
+  using StateUpdater = StateUpdaterTEMPLATE;
+
+  Updaters() {}
+
+  Updaters(std::shared_ptr<RateUpdaterTEMPLATE> rateUpdater,
+           std::shared_ptr<StateUpdaterTEMPLATE> stateUpdater)
+      : rate_(rateUpdater), state_(stateUpdater) {}
+
+  bool operator==(Updaters const &other) const {
+    return other.rate_ == rate_ and other.state_ == state_;
+  }
+
+  Updaters<RateUpdater, StateUpdater> clone() const {
+    return Updaters<RateUpdater, StateUpdater>(rate_->clone(), state_->clone());
+  }
+
+  std::shared_ptr<RateUpdaterTEMPLATE> rate_;
+  std::shared_ptr<StateUpdaterTEMPLATE> state_;
+};
+
+#endif

src/timestepping.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/istl/matrixindexset.hh>
-#include <dune/fufem/arithmetic.hh>
-
-#include "timestepping.hh"
-
-#include "timestepping/backward_euler.cc"
-#include "timestepping/newmark.cc"
-
-#include "timestepping_tmpl.cc"

src/timestepping/newmark.hh

-#ifndef DUNE_TECTONIC_TIMESTEPPING_NEWMARK_HH
-#define DUNE_TECTONIC_TIMESTEPPING_NEWMARK_HH
-
-template <class Vector, class Matrix, class Function, size_t dim>
-class Newmark : public TimeSteppingScheme<Vector, Matrix, Function, dim> {
-public:
-  Newmark(Matrix const &_A, Matrix const &_M, Matrix const &_C,
-          Vector const &_u_initial, Vector const &_v_initial,
-          Vector const &_a_initial,
-          Dune::BitSetVector<dim> const &_dirichletNodes,
-          Function const &_dirichletFunction);
-
-  void nextTimeStep() override;
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
-  void extractDisplacement(Vector &) const override;
-  void extractVelocity(Vector &) const override;
-  void extractOldVelocity(Vector &) const override;
-
-private:
-  Matrix const &A;
-  Matrix const &M;
-  Matrix const &C;
-  Vector u;
-  Vector v;
-  Vector a;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
-
-  Vector u_o;
-  Vector v_o;
-  Vector a_o;
-
-  double tau;
-
-  bool postProcessCalled = false;
-};
-#endif