src/spatial-solving/solverfactory_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/myblockproblem.hh>
+
+template class SolverFactory<
+    MY_DIM,
+    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+    Grid>;
+template class SolverFactory<
+    MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
+                ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+    Grid>;

src/state/compute_state_dieterich_euler.cc

-#include <cassert>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tnnmg/problem-classes/directionalconvexfunction.hh>
-
-#include "compute_state_dieterich_euler.hh"
-
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/fufem/interval.hh>
-
-/* Dieterich's law:
-
-   \f$ \dot \theta = 1 - V \theta/L \f$
-
-   Transformed using \alpha = \log \theta:
-
-   \f$ \dot \alpha = \exp(-\alpha) - V/L \f$
-
-   discretised in time (using backward euler):
-
-   \f$\delta \alpha_1 = \exp(-\alpha_1) - V/L\f$
-
-   or
-
-   \f$ 0 = ( \alpha_1 - \alpha_0 )/\tau - ( \exp(-\alpha_1) - V/L ) \f$
-
-   this is obviously monotone in \f$\alpha_1\f$; we can thus see if as
-   a convex minimisation problem.
-
-   find the antiderivative:
-
-   \f$ J(\alpha_1)
-   = 0.5/\tau * \alpha_1^2 - ( \alpha_0/\tau - V/L ) \alpha_1 +
-   \exp(-\alpha_1)\f$
-
-   Consequence: linear part is \f$ \alpha_0/\tau - V/L \f$; nonlinear
-   part is \f$\exp(-\alpha_1)\f$ which has derivative \f$
-   -\exp(-\alpha_1) \f$.
-*/
-
-//! The nonlinearity \f$f(x) = exp(-x)\f$
-class DieterichNonlinearity {
-public:
-  using VectorType = Dune::FieldVector<double, 1>;
-  using MatrixType = Dune::FieldMatrix<double, 1, 1>;
-
-  DieterichNonlinearity() {}
-
-  void directionalSubDiff(VectorType const &u, VectorType const &v,
-                          Interval<double> &D) const {
-    D[0] = D[1] = v[0] * (-std::exp(-u[0]));
-  }
-
-  void directionalDomain(VectorType const &, VectorType const &,
-                         Interval<double> &dom) const {
-    dom[0] = -std::numeric_limits<double>::max();
-    dom[1] = std::numeric_limits<double>::max();
-  }
-};
-
-double state_update_dieterich_euler(double tau, double VoL, double logState_o) {
-  DieterichNonlinearity::VectorType const start(0);
-  DieterichNonlinearity::VectorType const direction(1);
-  DieterichNonlinearity phi;
-  DirectionalConvexFunction<DieterichNonlinearity> const J(
-      1.0 / tau, logState_o / tau - VoL, phi, start, direction);
-
-  int bisectionsteps = 0;
-  Bisection const bisection;
-  return bisection.minimize(J, 0.0, 0.0, bisectionsteps);
-}

src/state/compute_state_dieterich_euler.hh

-#ifndef COMPUTE_STATE_DIETERICH_EULER_HH
-#define COMPUTE_STATE_DIETERICH_EULER_HH
-
-double state_update_dieterich_euler(double h, double VoL, double logState_o);
-#endif

src/state/compute_state_ruina.cc

-#include <cassert>
-
-#include <cmath>
-
-#include "compute_state_ruina.hh"
-
-double state_update_ruina(double tau, double uol, double logState_o) {
-  if (uol == 0)
-    return logState_o;
-
-  double const ret = logState_o - uol * std::log(uol / tau);
-  return ret / (1 + uol);
-}

src/state/compute_state_ruina.hh

-#ifndef COMPUTE_STATE_RUINA_HH
-#define COMPUTE_STATE_RUINA_HH
-
-double state_update_ruina(double h, double uol, double logState_o);
-#endif

src/state/dieterichstateupdater.hh

-#ifndef DIETERICH_STATE_UPDATER_HH
-#define DIETERICH_STATE_UPDATER_HH
-
-#include "compute_state_dieterich_euler.hh"
-#include "stateupdater.hh"
-
-template <class SingletonVector, class Vector>
-class DieterichStateUpdater : public StateUpdater<SingletonVector, Vector> {
-public:
-  DieterichStateUpdater(SingletonVector _logState_initial,
-                        Dune::BitSetVector<1> const &_nodes, double _L);
-
-  virtual void nextTimeStep();
-  virtual void setup(double _tau);
-  virtual void solve(Vector const &velocity_field);
-  virtual void extractLogState(SingletonVector &);
-
-private:
-  SingletonVector logState_o;
-  SingletonVector logState;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double tau;
-};
-
-template <class SingletonVector, class Vector>
-DieterichStateUpdater<SingletonVector, Vector>::DieterichStateUpdater(
-    SingletonVector _logState_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L)
-    : logState(_logState_initial), nodes(_nodes), L(_L) {}
-
-template <class SingletonVector, class Vector>
-void DieterichStateUpdater<SingletonVector, Vector>::nextTimeStep() {
-  logState_o = logState;
-}
-
-template <class SingletonVector, class Vector>
-void DieterichStateUpdater<SingletonVector, Vector>::setup(double _tau) {
-  tau = _tau;
-}
-
-template <class SingletonVector, class Vector>
-void DieterichStateUpdater<SingletonVector, Vector>::solve(
-    Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i)
-    if (nodes[i][0]) {
-      double const V = velocity_field[i].two_norm();
-      logState[i] = state_update_dieterich_euler(tau, V / L, logState_o[i]);
-    }
-}
-
-template <class SingletonVector, class Vector>
-void DieterichStateUpdater<SingletonVector, Vector>::extractLogState(
-    SingletonVector &_logState) {
-  _logState = logState;
-}
-
-#endif

src/state/ruinastateupdater.hh

-#ifndef RUINA_STATE_UPDATER_HH
-#define RUINA_STATE_UPDATER_HH
-
-#include "compute_state_ruina.hh"
-#include "stateupdater.hh"
-
-template <class SingletonVector, class Vector>
-class RuinaStateUpdater : public StateUpdater<SingletonVector, Vector> {
-public:
-  RuinaStateUpdater(SingletonVector _logState_initial,
-                    Dune::BitSetVector<1> const &_nodes, double _L);
-
-  virtual void nextTimeStep();
-  virtual void setup(double _tau);
-  virtual void solve(Vector const &velocity_field);
-  virtual void extractLogState(SingletonVector &);
-
-private:
-  SingletonVector logState_o;
-  SingletonVector logState;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double tau;
-};
-
-template <class SingletonVector, class Vector>
-RuinaStateUpdater<SingletonVector, Vector>::RuinaStateUpdater(
-    SingletonVector _logState_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L)
-    : logState(_logState_initial), nodes(_nodes), L(_L) {}
-
-template <class SingletonVector, class Vector>
-void RuinaStateUpdater<SingletonVector, Vector>::nextTimeStep() {
-  logState_o = logState;
-}
-
-template <class SingletonVector, class Vector>
-void RuinaStateUpdater<SingletonVector, Vector>::setup(double _tau) {
-  tau = _tau;
-}
-
-template <class SingletonVector, class Vector>
-void RuinaStateUpdater<SingletonVector, Vector>::solve(
-    Vector const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i)
-    if (nodes[i][0]) {
-      double const V = velocity_field[i].two_norm();
-      logState[i] = state_update_ruina(tau, V * tau / L, logState_o[i]);
-    }
-}
-
-template <class SingletonVector, class Vector>
-void RuinaStateUpdater<SingletonVector, Vector>::extractLogState(
-    SingletonVector &_logState) {
-  _logState = logState;
-}
-
-#endif

src/state/stateupdater.hh

-#ifndef STATE_UPDATER_HH
-#define STATE_UPDATER_HH
-
-template <class SingletonVector, class Vector> class StateUpdater {
-public:
-  virtual void nextTimeStep() = 0;
-  virtual void setup(double _tau) = 0;
-  virtual void solve(Vector const &velocity_field) = 0;
-  virtual void extractLogState(SingletonVector &logState) = 0;
-};
-
-#endif

src/time-stepping/adaptivetimestepper.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "adaptivetimestepper.hh"
+
+void IterationRegister::registerCount(FixedPointIterationCounter count) {
+  totalCount += count;
+}
+
+void IterationRegister::registerFinalCount(FixedPointIterationCounter count) {
+  finalCount = count;
+}
+
+void IterationRegister::reset() {
+  totalCount = FixedPointIterationCounter();
+  finalCount = FixedPointIterationCounter();
+}
+
+template <class Factory, class Updaters, class ErrorNorm>
+AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::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)
+    : relativeTime_(relativeTime),
+      relativeTau_(relativeTau),
+      finalTime_(parset.get<double>("problem.finalTime")),
+      factory_(factory),
+      parset_(parset),
+      globalFriction_(globalFriction),
+      current_(current),
+      R1_(),
+      externalForces_(externalForces),
+      mustRefine_(mustRefine),
+      errorNorm_(errorNorm) {}
+
+template <class Factory, class Updaters, class ErrorNorm>
+bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
+  return relativeTime_ >= 1.0;
+}
+
+template <class Factory, class Updaters, class ErrorNorm>
+IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
+  /*
+    |     C     | We check here if making the step R1 of size tau is a
+    |  R1 | R2  | good idea. To check if we can coarsen, we compare
+    |F1|F2|     | the result of (R1+R2) with C, i.e. two steps of size
+                  tau with one of size 2*tau. To check if we need to
+    refine, we compare the result of (F1+F2) with R1, i.e. two steps
+    of size tau/2 with one of size tau. The method makes multiple
+    coarsening/refining attempts, with coarsening coming first. */
+  if (R1_.updaters == Updaters())
+    R1_ = step(current_, relativeTime_, relativeTau_);
+
+  bool didCoarsen = false;
+  iterationRegister_.reset();
+  UpdatersWithCount R2;
+  UpdatersWithCount C;
+  while (relativeTime_ + relativeTau_ <= 1.0) {
+    R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
+    C = step(current_, relativeTime_, 2 * relativeTau_);
+    if (mustRefine_(R2.updaters, C.updaters))
+      break;
+
+    didCoarsen = true;
+    relativeTau_ *= 2;
+    R1_ = C;
+  }
+  UpdatersWithCount F1;
+  UpdatersWithCount F2;
+  if (!didCoarsen) {
+    while (true) {
+      F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
+      F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
+                relativeTau_ / 2.0);
+      if (!mustRefine_(F2.updaters, R1_.updaters))
+        break;
+
+      relativeTau_ /= 2.0;
+      R1_ = F1;
+      R2 = F2;
+    }
+  }
+
+  iterationRegister_.registerFinalCount(R1_.count);
+  relativeTime_ += relativeTau_;
+  current_ = R1_.updaters;
+  R1_ = R2;
+
+  return iterationRegister_;
+}
+
+template <class Factory, class Updaters, class ErrorNorm>
+typename AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
+    Updaters const &oldUpdaters, double rTime, double rTau) {
+  UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
+  newUpdatersAndCount.count =
+      MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
+                           newUpdatersAndCount.updaters, errorNorm_,
+                           externalForces_)
+          .step(rTime, rTau);
+  iterationRegister_.registerCount(newUpdatersAndCount.count);
+  return newUpdatersAndCount;
+}
+
+#include "adaptivetimestepper_tmpl.cc"

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/time-stepping/adaptivetimestepper_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/common/function.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/myblockproblem.hh>
+
+#include "../spatial-solving/solverfactory.hh"
+#include "rate/rateupdater.hh"
+#include "state/stateupdater.hh"
+#include "updaters.hh"
+
+using Function = Dune::VirtualFunction<double, double>;
+using Factory = SolverFactory<
+    MY_DIM,
+    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+    Grid>;
+using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+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/time-stepping/coupledtimestepper.hh

+#ifndef SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
+#define SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
+
+#include <functional>
+#include <memory>
+
+#include <dune/common/parametertree.hh>
+
+#include "../spatial-solving/fixedpointiterator.hh"
+
+template <class Factory, class Updaters, class ErrorNorm>
+class CoupledTimeStepper {
+  using Vector = typename Factory::Vector;
+  using Matrix = typename Factory::Matrix;
+  using ConvexProblem = typename Factory::ConvexProblem;
+  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+
+public:
+  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);
+
+  FixedPointIterationCounter step(double relativeTime, double relativeTau);
+
+private:
+  double finalTime_;
+  Factory &factory_;
+  Dune::ParameterTree const &parset_;
+  std::shared_ptr<Nonlinearity> globalFriction_;
+  Updaters updaters_;
+  std::function<void(double, Vector &)> externalForces_;
+  ErrorNorm const &errorNorm_;
+};
+#endif

src/time-stepping/coupledtimestepper_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/common/function.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/myblockproblem.hh>
+
+#include "../spatial-solving/solverfactory.hh"
+#include "rate/rateupdater.hh"
+#include "state/stateupdater.hh"
+#include "updaters.hh"
+
+using Function = Dune::VirtualFunction<double, double>;
+using Factory = SolverFactory<
+    MY_DIM,
+    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+    Grid>;
+using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+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;
-
-  tau = _tau;
+  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
+  this->tau = _tau;
 
   /* We start out with the formulation
 
@@ -49,75 +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);
   }
 
-  // v_o makes a good initial iterate; we could use anything, though
-  iterate = 0.0;
-
-  for (size_t i = 0; i < dirichletNodes.size(); ++i)
-    switch (dirichletNodes[i].count()) {
-      case 0:
-        continue;
-      case dim:
-        iterate[i] = 0;
-        dirichletFunction.evaluate(relativeTime, iterate[i][0]);
-        break;
-      case 1:
-        if (dirichletNodes[i][0]) {
-          dirichletFunction.evaluate(relativeTime, iterate[i][0]);
-          break;
-        }
-        if (dirichletNodes[i][1]) {
-          iterate[i][1] = 0;
-          break;
-        }
-        assert(false);
-      default:
-        assert(false);
-    }
+  iterate = this->v_o;
+
+  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
+    for (size_t j = 0; j < dim; ++j)
+      if (this->dirichletNodes[i][j])
+        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
     Vector const &iterate) {
-  postProcessCalled = true;
+  this->postProcessCalled = true;
 
-  v = iterate;
+  this->v = iterate;
 
-  u = u_o;
-  Arithmetic::addProduct(u, tau, v);
-}
+  this->u = this->u_o;
+  Arithmetic::addProduct(this->u, this->tau, this->v);
 
-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->a = this->v;
+  this->a -= this->v_o;
+  this->a /= this->tau;
 }
 
 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!");
-
-  velocity = v;
+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,
@@ -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,85 +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;
-
-  for (size_t i = 0; i < dirichletNodes.size(); ++i)
-    switch (dirichletNodes[i].count()) {
-      case 0:
-        continue;
-      case dim:
-        iterate[i] = 0;
-        dirichletFunction.evaluate(relativeTime, iterate[i][0]);
-        break;
-      case 1:
-        if (dirichletNodes[i][0]) {
-          dirichletFunction.evaluate(relativeTime, iterate[i][0]);
-          break;
-        }
-        if (dirichletNodes[i][1]) {
-          iterate[i][1] = 0;
-          break;
-        }
-        assert(false);
-      default:
-        assert(false);
-    }
+  iterate = this->v_o;
+
+  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
+    for (size_t j = 0; j < dim; ++j)
+      if (this->dirichletNodes[i][j])
+        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void Newmark<Vector, Matrix, Function, dim>::postProcess(
     Vector const &iterate) {
-  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;
+  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>::extractVelocity(Vector &velocity)
-    const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  velocity = v;
+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