src/time-stepping/adaptivetimestepper.hh

+#ifndef SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
+#define SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
+
+#include <fstream>
+
+#include "coupledtimestepper.hh"
+
+struct IterationRegister {
+  void registerCount(FixedPointIterationCounter count);
+  void registerFinalCount(FixedPointIterationCounter count);
+
+  void reset();
+
+  FixedPointIterationCounter totalCount;
+  FixedPointIterationCounter finalCount;
+};
+
+template <class Factory, class Updaters, class ErrorNorm>
+class AdaptiveTimeStepper {
+  struct UpdatersWithCount {
+    Updaters updaters;
+    FixedPointIterationCounter count;
+  };
+
+  using Vector = typename Factory::Vector;
+  using ConvexProblem = typename Factory::ConvexProblem;
+  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+
+  using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
+
+public:
+  AdaptiveTimeStepper(Factory &factory, Dune::ParameterTree const &parset,
+                      std::shared_ptr<Nonlinearity> globalFriction,
+                      Updaters &current, double relativeTime,
+                      double relativeTau,
+                      std::function<void(double, Vector &)> externalForces,
+                      ErrorNorm const &errorNorm,
+                      std::function<bool(Updaters &, Updaters &)> mustRefine);
+
+  bool reachedEnd();
+  IterationRegister advance();
+
+  double relativeTime_;
+  double relativeTau_;
+
+private:
+  UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
+                         double rTau);
+
+  double finalTime_;
+  Factory &factory_;
+  Dune::ParameterTree const &parset_;
+  std::shared_ptr<Nonlinearity> globalFriction_;
+  Updaters &current_;
+  UpdatersWithCount R1_;
+  std::function<void(double, Vector &)> externalForces_;
+  std::function<bool(Updaters &, Updaters &)> mustRefine_;
+  ErrorNorm const &errorNorm_;
+
+  IterationRegister iterationRegister_;
+};
+#endif

src/fixedpointiterator_tmpl.cc → src/time-stepping/adaptivetimestepper_tmpl.cc

@@ -2,6 +2,9 @@
 #error MY_DIM unset
 #endif
 
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
 #include <dune/common/function.hh>
 
 #include <dune/solvers/norms/energynorm.hh>
@@ -10,12 +13,10 @@
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
-#include "explicitgrid.hh"
-#include "explicitvectors.hh"
-
-#include "solverfactory.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "rate/rateupdater.hh"
 #include "state/stateupdater.hh"
-#include "timestepping.hh"
+#include "updaters.hh"
 
 using Function = Dune::VirtualFunction<double, double>;
 using Factory = SolverFactory<
@@ -23,9 +24,9 @@ using Factory = SolverFactory<
     MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
     Grid>;
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyVelocityUpdater = TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
 
 using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
 
-template class FixedPointIterator<Factory, MyStateUpdater, MyVelocityUpdater,
-                                  ErrorNorm>;
+template class AdaptiveTimeStepper<Factory, MyUpdaters, ErrorNorm>;

src/time-stepping/coupledtimestepper.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "coupledtimestepper.hh"
+
+template <class Factory, class Updaters, class ErrorNorm>
+CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
+    double finalTime, Factory &factory, Dune::ParameterTree const &parset,
+    std::shared_ptr<Nonlinearity> globalFriction, Updaters updaters,
+    ErrorNorm const &errorNorm,
+    std::function<void(double, Vector &)> externalForces)
+    : finalTime_(finalTime),
+      factory_(factory),
+      parset_(parset),
+      globalFriction_(globalFriction),
+      updaters_(updaters),
+      externalForces_(externalForces),
+      errorNorm_(errorNorm) {}
+
+template <class Factory, class Updaters, class ErrorNorm>
+FixedPointIterationCounter
+CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
+                                                       double relativeTau) {
+  updaters_.state_->nextTimeStep();
+  updaters_.rate_->nextTimeStep();
+
+  auto const newRelativeTime = relativeTime + relativeTau;
+  Vector ell;
+  externalForces_(newRelativeTime, ell);
+
+  Matrix velocityMatrix;
+  Vector velocityRHS;
+  Vector velocityIterate;
+
+  auto const tau = relativeTau * finalTime_;
+  updaters_.state_->setup(tau);
+  updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS,
+                         velocityIterate, velocityMatrix);
+  FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
+      factory_, parset_, globalFriction_, errorNorm_);
+  auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
+                                                 velocityRHS, velocityIterate);
+  return iterations;
+}
+
+#include "coupledtimestepper_tmpl.cc"

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

-#ifndef SRC_COUPLEDTIMESTEPPER_HH
-#define SRC_COUPLEDTIMESTEPPER_HH
+#ifndef SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
+#define SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
 
 #include <functional>
 #include <memory>
 
 #include <dune/common/parametertree.hh>
-#include "fixedpointiterator.hh"
 
-template <class Factory, class StateUpdater, class VelocityUpdater,
-          class ErrorNorm>
+#include "../spatial-solving/fixedpointiterator.hh"
+
+template <class Factory, class Updaters, class ErrorNorm>
 class CoupledTimeStepper {
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
@@ -19,9 +19,7 @@ class CoupledTimeStepper {
   CoupledTimeStepper(double finalTime, Factory &factory,
                      Dune::ParameterTree const &parset,
                      std::shared_ptr<Nonlinearity> globalFriction,
-                     std::shared_ptr<StateUpdater> stateUpdater,
-                     std::shared_ptr<VelocityUpdater> velocityUpdater,
-                     ErrorNorm const &errorNorm,
+                     Updaters updaters, ErrorNorm const &errorNorm,
                      std::function<void(double, Vector &)> externalForces);
 
   FixedPointIterationCounter step(double relativeTime, double relativeTau);
@@ -31,8 +29,7 @@ class CoupledTimeStepper {
   Factory &factory_;
   Dune::ParameterTree const &parset_;
   std::shared_ptr<Nonlinearity> globalFriction_;
-  std::shared_ptr<StateUpdater> stateUpdater_;
-  std::shared_ptr<VelocityUpdater> velocityUpdater_;
+  Updaters updaters_;
   std::function<void(double, Vector &)> externalForces_;
   ErrorNorm const &errorNorm_;
 };

src/coupledtimestepper_tmpl.cc → src/time-stepping/coupledtimestepper_tmpl.cc

@@ -2,6 +2,9 @@
 #error MY_DIM unset
 #endif
 
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
 #include <dune/common/function.hh>
 
 #include <dune/solvers/norms/energynorm.hh>
@@ -10,12 +13,10 @@
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
-#include "explicitgrid.hh"
-#include "explicitvectors.hh"
-
-#include "solverfactory.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "rate/rateupdater.hh"
 #include "state/stateupdater.hh"
-#include "timestepping.hh"
+#include "updaters.hh"
 
 using Function = Dune::VirtualFunction<double, double>;
 using Factory = SolverFactory<
@@ -23,9 +24,9 @@ using Factory = SolverFactory<
     MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
     Grid>;
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyVelocityUpdater = TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
 
 using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
 
-template class CoupledTimeStepper<Factory, MyStateUpdater, MyVelocityUpdater,
-                                  ErrorNorm>;
+template class CoupledTimeStepper<Factory, MyUpdaters, ErrorNorm>;

src/time-stepping/rate.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "rate.hh"
+#include "rate/backward_euler.hh"
+#include "rate/newmark.hh"
+
+template <class Vector, class Matrix, class Function, int dimension>
+std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
+initRateUpdater(Config::scheme scheme,
+                Function const &velocityDirichletFunction,
+                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
+                Matrices<Matrix> const &matrices, Vector const &u_initial,
+                Vector const &v_initial, Vector const &a_initial) {
+  switch (scheme) {
+    case Config::Newmark:
+      return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
+          matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
+          velocityDirichletFunction);
+    case Config::BackwardEuler:
+      return std::make_shared<
+          BackwardEuler<Vector, Matrix, Function, dimension>>(
+          matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
+          velocityDirichletFunction);
+    default:
+      assert(false);
+  }
+}
+
+#include "rate_tmpl.cc"

src/time-stepping/rate.hh

+#ifndef SRC_TIME_STEPPING_RATE_HH
+#define SRC_TIME_STEPPING_RATE_HH
+
+#include <memory>
+
+#include "../enums.hh"
+#include "rate/rateupdater.hh"
+
+template <class Vector, class Matrix, class Function, int dimension>
+std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
+initRateUpdater(Config::scheme scheme,
+                Function const &velocityDirichletFunction,
+                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
+                Matrices<Matrix> const &matrices, Vector const &u_initial,
+                Vector const &v_initial, Vector const &a_initial);
+#endif

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

 #include <dune/solvers/common/arithmetic.hh>
 
+#include "backward_euler.hh"
+
 template <class Vector, class Matrix, class Function, size_t dim>
 BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
-    Matrix const &_A, Matrix const &_M, Matrix const &_C,
-    Vector const &_u_initial, Vector const &_v_initial,
+    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),
-      dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::nextTimeStep() {
-  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 BackwardEuler<Vector, Matrix, Function, dim>::setup(
     Vector const &ell, double _tau, double relativeTime, Vector &rhs,
     Vector &iterate, Matrix &AM) {
-  postProcessCalled = false;
-  dirichletFunction.evaluate(relativeTime, dirichletValue);
-  tau = _tau;
+  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
+  this->tau = _tau;
 
   /* We start out with the formulation
 
@@ -51,69 +42,51 @@ void BackwardEuler<Vector, Matrix, Function, dim>::setup(
 
   // 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, 1.0 / tau, M);
-  Arithmetic::addProduct(AM, 1.0, C);
-  Arithmetic::addProduct(AM, tau, A);
+  Arithmetic::addProduct(AM, 1.0 / this->tau, this->matrices.mass);
+  Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
+  Arithmetic::addProduct(AM, this->tau, this->matrices.elasticity);
 
   // set up RHS
   {
     rhs = ell;
-    Arithmetic::addProduct(rhs, 1.0 / tau, M, v_o);
-    Arithmetic::subtractProduct(rhs, A, u_o);
+    Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
+                           this->v_o);
+    Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
   }
 
-  iterate = v_o;
+  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])
-        iterate[i][j] = (j == 0) ? dirichletValue : 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 BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
     Vector const &iterate) {
-  postProcessCalled = true;
-
-  v = iterate;
-
-  u = u_o;
-  Arithmetic::addProduct(u, tau, v);
-}
+  this->postProcessCalled = true;
 
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::extractDisplacement(
-    Vector &displacement) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  displacement = u;
-}
+  this->v = iterate;
 
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::extractVelocity(
-    Vector &velocity) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+  this->u = this->u_o;
+  Arithmetic::addProduct(this->u, this->tau, this->v);
 
-  velocity = v;
-}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::extractOldVelocity(
-    Vector &velocity) const {
-  velocity = v_o;
+  this->a = this->v;
+  this->a -= this->v_o;
+  this->a /= this->tau;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dim>>
+std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>>
 BackwardEuler<Vector, Matrix, Function, dim>::clone() const {
   return std::make_shared<BackwardEuler<Vector, Matrix, Function, dim>>(*this);
 }

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> {
+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);
+
+  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;
+
+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,
@@ -28,9 +18,8 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
                                                    double relativeTime,
                                                    Vector &rhs, Vector &iterate,
                                                    Matrix &AM) {
-  postProcessCalled = false;
-  dirichletFunction.evaluate(relativeTime, dirichletValue);
-  tau = _tau;
+  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
+  this->tau = _tau;
 
   /* We start out with the formulation
 
@@ -56,79 +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);
   }
 
-  iterate = v_o;
+  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])
-        iterate[i][j] = (j == 0) ? dirichletValue : 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);
-}
-
-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 Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::extractOldVelocity(
-    Vector &velocity) const {
-  velocity = v_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>
-std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, 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);
 }

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/timestepping_tmpl.cc → src/time-stepping/rate/rateupdater_tmpl.cc

@@ -4,9 +4,10 @@
 
 #include <dune/common/function.hh>
 
-#include "explicitvectors.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/state.hh → src/time-stepping/state.hh

-#ifndef SRC_STATE_HH
-#define SRC_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 "state/stateupdater.hh"
+#include "../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(

src/state/ageinglawstateupdater.cc → src/time-stepping/state/ageinglawstateupdater.cc

+#include <cmath>
+
 #include "ageinglawstateupdater.hh"
-#include "../tobool.hh"
+#include "../../tobool.hh"
 
 template <class ScalarVector, class Vector>
 AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(

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

-#ifndef SRC_STATE_AGEINGLAWSTATEUPDATER_HH
-#define SRC_STATE_AGEINGLAWSTATEUPDATER_HH
+#ifndef SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
+#define SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
 
 #include "stateupdater.hh"
 
@@ -15,7 +15,7 @@ class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> {
   void solve(Vector const &velocity_field) override;
   void extractAlpha(ScalarVector &) override;
 
-  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const;
+  std::shared_ptr<StateUpdater<ScalarVector, Vector>> clone() const override;
 
 private:
   ScalarVector alpha_o;