From 4e6db7b053d119312309bca1056f2467a96776fc Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 11 Jul 2013 17:17:12 +0200
Subject: [PATCH] [Cleanup] Split up timestepping into multiple headers

---
 src/timestepping.cc               | 291 +-----------------------------
 src/timestepping.hh               | 100 +---------
 src/timestepping/eulerpair.cc     |  95 ++++++++++
 src/timestepping/eulerpair.hh     |  35 ++++
 src/timestepping/impliciteuler.cc |  87 +++++++++
 src/timestepping/impliciteuler.hh |  34 ++++
 src/timestepping/newmark.cc       | 104 +++++++++++
 src/timestepping/newmark.hh       |  38 ++++
 8 files changed, 400 insertions(+), 384 deletions(-)
 create mode 100644 src/timestepping/eulerpair.cc
 create mode 100644 src/timestepping/eulerpair.hh
 create mode 100644 src/timestepping/impliciteuler.cc
 create mode 100644 src/timestepping/impliciteuler.hh
 create mode 100644 src/timestepping/newmark.cc
 create mode 100644 src/timestepping/newmark.hh

diff --git a/src/timestepping.cc b/src/timestepping.cc
index fe6b09e2..42a9b679 100644
--- a/src/timestepping.cc
+++ b/src/timestepping.cc
@@ -4,294 +4,11 @@
 
 #include <dune/istl/matrixindexset.hh>
 #include <dune/fufem/arithmetic.hh>
-#include "timestepping.hh"
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::ImplicitEuler(
-    MatrixType const &_A, VectorType const &_u_initial,
-    VectorType const &_v_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    FunctionType const &_dirichletFunction)
-    : A(_A),
-      u(_u_initial),
-      v(_v_initial),
-      dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
-  v_old = v;
-  u_old = u;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::setup(
-    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
-    VectorType &problem_iterate, MatrixType &problem_AB) {
-  postProcessCalled = false;
-
-  tau = _tau;
-
-  problem_rhs = ell;
-  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
-
-  // For fixed tau, we'd only really have to do this once
-  problem_AB = A;
-  problem_AB *= tau;
-
-  // v_old makes a good initial iterate; we could use anything, though
-  problem_iterate = v_old;
-
-  for (size_t i = 0; i < dirichletNodes.size(); ++i)
-    switch (dirichletNodes[i].count()) {
-      case 0:
-        continue;
-      case dim:
-        problem_iterate[i] = 0;
-        dirichletFunction.evaluate(time, problem_iterate[i][0]);
-        break;
-      case 1:
-        if (dirichletNodes[i][0]) {
-          dirichletFunction.evaluate(time, problem_iterate[i][0]);
-          break;
-        }
-        if (dirichletNodes[i][1]) {
-          problem_iterate[i][1] = 0;
-          break;
-        }
-        assert(false);
-      default:
-        assert(false);
-    }
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::postProcess(
-    VectorType const &problem_iterate) {
-  postProcessCalled = true;
-
-  v = problem_iterate;
-
-  u = u_old;
-  Arithmetic::addProduct(u, tau, v);
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void ImplicitEuler<VectorType, MatrixType, FunctionType,
-                   dim>::extractDisplacement(VectorType &displacement) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  displacement = u;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
-    VectorType &velocity) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  velocity = v;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
-    MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
-    VectorType const &_v_initial, VectorType const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    FunctionType const &_dirichletFunction)
-    : A(_A),
-      B(_B),
-      u(_u_initial),
-      v(_v_initial),
-      a(_a_initial),
-      dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void Newmark<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
-  a_old = a;
-  v_old = v;
-  u_old = u;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
-    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
-    VectorType &problem_iterate, MatrixType &problem_AB) {
-  postProcessCalled = false;
-
-  tau = _tau;
-
-  problem_rhs = ell;
-  Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, v_old);
-  Arithmetic::addProduct(problem_rhs, B, a_old);
-  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
-  Arithmetic::addProduct(problem_rhs, -tau / 2.0, A, v_old);
-
-  // For fixed tau, we'd only really have to do this once
-  Dune::MatrixIndexSet indices(A.N(), A.M());
-  indices.import(A);
-  indices.import(B);
-  indices.exportIdx(problem_AB);
-  problem_AB = 0.0;
-  Arithmetic::addProduct(problem_AB, tau / 2.0, A);
-  Arithmetic::addProduct(problem_AB, 2.0 / tau, B);
-
-  // v_old makes a good initial iterate; we could use anything, though
-  problem_iterate = 0.0;
 
-  for (size_t i = 0; i < dirichletNodes.size(); ++i)
-    switch (dirichletNodes[i].count()) {
-      case 0:
-        continue;
-      case dim:
-        problem_iterate[i] = 0;
-        dirichletFunction.evaluate(time, problem_iterate[i][0]);
-        break;
-      case 1:
-        if (dirichletNodes[i][0]) {
-          dirichletFunction.evaluate(time, problem_iterate[i][0]);
-          break;
-        }
-        if (dirichletNodes[i][1]) {
-          problem_iterate[i][1] = 0;
-          break;
-        }
-        assert(false);
-      default:
-        assert(false);
-    }
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void Newmark<VectorType, MatrixType, FunctionType, dim>::postProcess(
-    VectorType const &problem_iterate) {
-  postProcessCalled = true;
-
-  v = problem_iterate;
-
-  u = u_old;
-  Arithmetic::addProduct(u, tau / 2.0, v);
-  Arithmetic::addProduct(u, tau / 2.0, v_old);
-
-  a = 0;
-  Arithmetic::addProduct(a, 2.0 / tau, v);
-  Arithmetic::addProduct(a, -2.0 / tau, v_old);
-  Arithmetic::addProduct(a, -1.0, a_old);
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void Newmark<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
-    VectorType &displacement) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  displacement = u;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void Newmark<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
-    VectorType &velocity) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  velocity = v;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-EulerPair<VectorType, MatrixType, FunctionType, dim>::EulerPair(
-    MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
-    VectorType const &_v_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    FunctionType const &_dirichletFunction)
-    : A(_A),
-      B(_B),
-      u(_u_initial),
-      v(_v_initial),
-      dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
-  v_old = v;
-  u_old = u;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
-    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
-    VectorType &problem_iterate, MatrixType &problem_AB) {
-  postProcessCalled = false;
-
-  tau = _tau;
-
-  problem_rhs = ell;
-  Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, v_old);
-  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
-
-  // For fixed tau, we'd only really have to do this once
-  Dune::MatrixIndexSet indices(A.N(), A.M());
-  indices.import(A);
-  indices.import(B);
-  indices.exportIdx(problem_AB);
-  problem_AB = 0.0;
-  Arithmetic::addProduct(problem_AB, tau, A);
-  Arithmetic::addProduct(problem_AB, 1.0 / tau, B);
-
-  // v_old makes a good initial iterate; we could use anything, though
-  problem_iterate = v_old;
-
-  for (size_t i = 0; i < dirichletNodes.size(); ++i)
-    switch (dirichletNodes[i].count()) {
-      case 0:
-        continue;
-      case dim:
-        problem_iterate[i] = 0;
-        dirichletFunction.evaluate(time, problem_iterate[i][0]);
-        break;
-      case 1:
-        if (dirichletNodes[i][0]) {
-          dirichletFunction.evaluate(time, problem_iterate[i][0]);
-          break;
-        }
-        if (dirichletNodes[i][1]) {
-          problem_iterate[i][1] = 0;
-          break;
-        }
-        assert(false);
-      default:
-        assert(false);
-    }
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void EulerPair<VectorType, MatrixType, FunctionType, dim>::postProcess(
-    VectorType const &problem_iterate) {
-  postProcessCalled = true;
-
-  v = problem_iterate;
-
-  u = u_old;
-  Arithmetic::addProduct(u, tau, v);
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
-    VectorType &displacement) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-
-  displacement = u;
-}
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
-    VectorType &velocity) const {
-  if (!postProcessCalled)
-    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+#include "timestepping.hh"
 
-  velocity = v;
-}
+#include "timestepping/eulerpair.cc"
+#include "timestepping/impliciteuler.cc"
+#include "timestepping/newmark.cc"
 
 #include "timestepping_tmpl.cc"
diff --git a/src/timestepping.hh b/src/timestepping.hh
index efbcb416..b09e6b30 100644
--- a/src/timestepping.hh
+++ b/src/timestepping.hh
@@ -16,102 +16,8 @@ class TimeSteppingScheme {
   void virtual extractVelocity(VectorType &velocity) const = 0;
 };
 
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class ImplicitEuler
-    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
-public:
-  ImplicitEuler(MatrixType const &_A, VectorType const &_u_initial,
-                VectorType const &_v_initial,
-                Dune::BitSetVector<dim> const &_dirichletNodes,
-                FunctionType const &_dirichletFunction);
-
-  void virtual nextTimeStep() override;
-  void virtual setup(VectorType const &, double, double, VectorType &,
-                     VectorType &, MatrixType &) override;
-  void virtual postProcess(VectorType const &) override;
-  void virtual extractDisplacement(VectorType &) const override;
-  void virtual extractVelocity(VectorType &) const override;
-
-private:
-  MatrixType const &A;
-  VectorType u;
-  VectorType v;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  FunctionType const &dirichletFunction;
-
-  VectorType u_old;
-  VectorType v_old;
-
-  double tau;
-
-  bool postProcessCalled = false;
-};
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class Newmark
-    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
-public:
-  Newmark(MatrixType const &_A, MatrixType const &_B,
-          VectorType const &_u_initial, VectorType const &_v_initial,
-          VectorType const &_a_initial,
-          Dune::BitSetVector<dim> const &_dirichletNodes,
-          FunctionType const &_dirichletFunction);
-
-  void virtual nextTimeStep() override;
-  void virtual setup(VectorType const &, double, double, VectorType &,
-                     VectorType &, MatrixType &) override;
-  void virtual postProcess(VectorType const &) override;
-  void virtual extractDisplacement(VectorType &) const override;
-  void virtual extractVelocity(VectorType &) const override;
-
-private:
-  MatrixType const &A;
-  MatrixType const &B;
-  VectorType u;
-  VectorType v;
-  VectorType a;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  FunctionType const &dirichletFunction;
-
-  VectorType u_old;
-  VectorType v_old;
-  VectorType a_old;
-
-  double tau;
-
-  bool postProcessCalled = false;
-};
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class EulerPair
-    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
-public:
-  EulerPair(MatrixType const &_A, MatrixType const &_B,
-            VectorType const &_u_initial, VectorType const &_v_initial,
-            Dune::BitSetVector<dim> const &_dirichletNodes,
-            FunctionType const &_dirichletFunction);
-
-  void virtual nextTimeStep() override;
-  void virtual setup(VectorType const &, double, double, VectorType &,
-                     VectorType &, MatrixType &) override;
-  void virtual postProcess(VectorType const &) override;
-  void virtual extractDisplacement(VectorType &) const override;
-  void virtual extractVelocity(VectorType &) const override;
-
-private:
-  MatrixType const &A;
-  MatrixType const &B;
-  VectorType u;
-  VectorType v;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  FunctionType const &dirichletFunction;
-
-  VectorType u_old;
-  VectorType v_old;
-
-  double tau;
-
-  bool postProcessCalled = false;
-};
+#include "timestepping/impliciteuler.hh"
+#include "timestepping/newmark.hh"
+#include "timestepping/eulerpair.hh"
 
 #endif
diff --git a/src/timestepping/eulerpair.cc b/src/timestepping/eulerpair.cc
new file mode 100644
index 00000000..85cee621
--- /dev/null
+++ b/src/timestepping/eulerpair.cc
@@ -0,0 +1,95 @@
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+EulerPair<VectorType, MatrixType, FunctionType, dim>::EulerPair(
+    MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+    VectorType const &_v_initial,
+    Dune::BitSetVector<dim> const &_dirichletNodes,
+    FunctionType const &_dirichletFunction)
+    : A(_A),
+      B(_B),
+      u(_u_initial),
+      v(_v_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+  v_old = v;
+  u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
+    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+    VectorType &problem_iterate, MatrixType &problem_AB) {
+  postProcessCalled = false;
+
+  tau = _tau;
+
+  problem_rhs = ell;
+  Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, v_old);
+  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
+
+  // For fixed tau, we'd only really have to do this once
+  Dune::MatrixIndexSet indices(A.N(), A.M());
+  indices.import(A);
+  indices.import(B);
+  indices.exportIdx(problem_AB);
+  problem_AB = 0.0;
+  Arithmetic::addProduct(problem_AB, tau, A);
+  Arithmetic::addProduct(problem_AB, 1.0 / tau, B);
+
+  // v_old makes a good initial iterate; we could use anything, though
+  problem_iterate = v_old;
+
+  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+    switch (dirichletNodes[i].count()) {
+      case 0:
+        continue;
+      case dim:
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+        break;
+      case 1:
+        if (dirichletNodes[i][0]) {
+          dirichletFunction.evaluate(time, problem_iterate[i][0]);
+          break;
+        }
+        if (dirichletNodes[i][1]) {
+          problem_iterate[i][1] = 0;
+          break;
+        }
+        assert(false);
+      default:
+        assert(false);
+    }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::postProcess(
+    VectorType const &problem_iterate) {
+  postProcessCalled = true;
+
+  v = problem_iterate;
+
+  u = u_old;
+  Arithmetic::addProduct(u, tau, v);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
+    VectorType &displacement) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void EulerPair<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+    VectorType &velocity) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  velocity = v;
+}
diff --git a/src/timestepping/eulerpair.hh b/src/timestepping/eulerpair.hh
new file mode 100644
index 00000000..34e684a5
--- /dev/null
+++ b/src/timestepping/eulerpair.hh
@@ -0,0 +1,35 @@
+#ifndef DUNE_TECTONIC_TIMESTEPPING_EULERPAIR_HH
+#define DUNE_TECTONIC_TIMESTEPPING_EULERPAIR_HH
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class EulerPair
+    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+  EulerPair(MatrixType const &_A, MatrixType const &_B,
+            VectorType const &_u_initial, VectorType const &_v_initial,
+            Dune::BitSetVector<dim> const &_dirichletNodes,
+            FunctionType const &_dirichletFunction);
+
+  void virtual nextTimeStep() override;
+  void virtual setup(VectorType const &, double, double, VectorType &,
+                     VectorType &, MatrixType &) override;
+  void virtual postProcess(VectorType const &) override;
+  void virtual extractDisplacement(VectorType &) const override;
+  void virtual extractVelocity(VectorType &) const override;
+
+private:
+  MatrixType const &A;
+  MatrixType const &B;
+  VectorType u;
+  VectorType v;
+  Dune::BitSetVector<dim> const &dirichletNodes;
+  FunctionType const &dirichletFunction;
+
+  VectorType u_old;
+  VectorType v_old;
+
+  double tau;
+
+  bool postProcessCalled = false;
+};
+#endif
diff --git a/src/timestepping/impliciteuler.cc b/src/timestepping/impliciteuler.cc
new file mode 100644
index 00000000..f2dd8be8
--- /dev/null
+++ b/src/timestepping/impliciteuler.cc
@@ -0,0 +1,87 @@
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::ImplicitEuler(
+    MatrixType const &_A, VectorType const &_u_initial,
+    VectorType const &_v_initial,
+    Dune::BitSetVector<dim> const &_dirichletNodes,
+    FunctionType const &_dirichletFunction)
+    : A(_A),
+      u(_u_initial),
+      v(_v_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+  v_old = v;
+  u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::setup(
+    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+    VectorType &problem_iterate, MatrixType &problem_AB) {
+  postProcessCalled = false;
+
+  tau = _tau;
+
+  problem_rhs = ell;
+  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
+
+  // For fixed tau, we'd only really have to do this once
+  problem_AB = A;
+  problem_AB *= tau;
+
+  // v_old makes a good initial iterate; we could use anything, though
+  problem_iterate = v_old;
+
+  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+    switch (dirichletNodes[i].count()) {
+      case 0:
+        continue;
+      case dim:
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+        break;
+      case 1:
+        if (dirichletNodes[i][0]) {
+          dirichletFunction.evaluate(time, problem_iterate[i][0]);
+          break;
+        }
+        if (dirichletNodes[i][1]) {
+          problem_iterate[i][1] = 0;
+          break;
+        }
+        assert(false);
+      default:
+        assert(false);
+    }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::postProcess(
+    VectorType const &problem_iterate) {
+  postProcessCalled = true;
+
+  v = problem_iterate;
+
+  u = u_old;
+  Arithmetic::addProduct(u, tau, v);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType,
+                   dim>::extractDisplacement(VectorType &displacement) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+    VectorType &velocity) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  velocity = v;
+}
diff --git a/src/timestepping/impliciteuler.hh b/src/timestepping/impliciteuler.hh
new file mode 100644
index 00000000..93f2d3dc
--- /dev/null
+++ b/src/timestepping/impliciteuler.hh
@@ -0,0 +1,34 @@
+#ifndef DUNE_TECTONIC_TIMESTEPPING_IMPLICITEULER_HH
+#define DUNE_TECTONIC_TIMESTEPPING_IMPLICITEULER_HH
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class ImplicitEuler
+    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+  ImplicitEuler(MatrixType const &_A, VectorType const &_u_initial,
+                VectorType const &_v_initial,
+                Dune::BitSetVector<dim> const &_dirichletNodes,
+                FunctionType const &_dirichletFunction);
+
+  void virtual nextTimeStep() override;
+  void virtual setup(VectorType const &, double, double, VectorType &,
+                     VectorType &, MatrixType &) override;
+  void virtual postProcess(VectorType const &) override;
+  void virtual extractDisplacement(VectorType &) const override;
+  void virtual extractVelocity(VectorType &) const override;
+
+private:
+  MatrixType const &A;
+  VectorType u;
+  VectorType v;
+  Dune::BitSetVector<dim> const &dirichletNodes;
+  FunctionType const &dirichletFunction;
+
+  VectorType u_old;
+  VectorType v_old;
+
+  double tau;
+
+  bool postProcessCalled = false;
+};
+#endif
diff --git a/src/timestepping/newmark.cc b/src/timestepping/newmark.cc
new file mode 100644
index 00000000..a479b08f
--- /dev/null
+++ b/src/timestepping/newmark.cc
@@ -0,0 +1,104 @@
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+Newmark<VectorType, MatrixType, FunctionType, dim>::Newmark(
+    MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+    VectorType const &_v_initial, VectorType const &_a_initial,
+    Dune::BitSetVector<dim> const &_dirichletNodes,
+    FunctionType const &_dirichletFunction)
+    : A(_A),
+      B(_B),
+      u(_u_initial),
+      v(_v_initial),
+      a(_a_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction) {}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
+  a_old = a;
+  v_old = v;
+  u_old = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::setup(
+    VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
+    VectorType &problem_iterate, MatrixType &problem_AB) {
+  postProcessCalled = false;
+
+  tau = _tau;
+
+  problem_rhs = ell;
+  Arithmetic::addProduct(problem_rhs, 2.0 / tau, B, v_old);
+  Arithmetic::addProduct(problem_rhs, B, a_old);
+  Arithmetic::addProduct(problem_rhs, -1.0, A, u_old);
+  Arithmetic::addProduct(problem_rhs, -tau / 2.0, A, v_old);
+
+  // For fixed tau, we'd only really have to do this once
+  Dune::MatrixIndexSet indices(A.N(), A.M());
+  indices.import(A);
+  indices.import(B);
+  indices.exportIdx(problem_AB);
+  problem_AB = 0.0;
+  Arithmetic::addProduct(problem_AB, tau / 2.0, A);
+  Arithmetic::addProduct(problem_AB, 2.0 / tau, B);
+
+  // v_old makes a good initial iterate; we could use anything, though
+  problem_iterate = 0.0;
+
+  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+    switch (dirichletNodes[i].count()) {
+      case 0:
+        continue;
+      case dim:
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+        break;
+      case 1:
+        if (dirichletNodes[i][0]) {
+          dirichletFunction.evaluate(time, problem_iterate[i][0]);
+          break;
+        }
+        if (dirichletNodes[i][1]) {
+          problem_iterate[i][1] = 0;
+          break;
+        }
+        assert(false);
+      default:
+        assert(false);
+    }
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::postProcess(
+    VectorType const &problem_iterate) {
+  postProcessCalled = true;
+
+  v = problem_iterate;
+
+  u = u_old;
+  Arithmetic::addProduct(u, tau / 2.0, v);
+  Arithmetic::addProduct(u, tau / 2.0, v_old);
+
+  a = 0;
+  Arithmetic::addProduct(a, 2.0 / tau, v);
+  Arithmetic::addProduct(a, -2.0 / tau, v_old);
+  Arithmetic::addProduct(a, -1.0, a_old);
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::extractDisplacement(
+    VectorType &displacement) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  displacement = u;
+}
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void Newmark<VectorType, MatrixType, FunctionType, dim>::extractVelocity(
+    VectorType &velocity) const {
+  if (!postProcessCalled)
+    DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+  velocity = v;
+}
diff --git a/src/timestepping/newmark.hh b/src/timestepping/newmark.hh
new file mode 100644
index 00000000..b9c1e91d
--- /dev/null
+++ b/src/timestepping/newmark.hh
@@ -0,0 +1,38 @@
+#ifndef DUNE_TECTONIC_TIMESTEPPING_NEWMARK_HH
+#define DUNE_TECTONIC_TIMESTEPPING_NEWMARK_HH
+
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+class Newmark
+    : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
+public:
+  Newmark(MatrixType const &_A, MatrixType const &_B,
+          VectorType const &_u_initial, VectorType const &_v_initial,
+          VectorType const &_a_initial,
+          Dune::BitSetVector<dim> const &_dirichletNodes,
+          FunctionType const &_dirichletFunction);
+
+  void virtual nextTimeStep() override;
+  void virtual setup(VectorType const &, double, double, VectorType &,
+                     VectorType &, MatrixType &) override;
+  void virtual postProcess(VectorType const &) override;
+  void virtual extractDisplacement(VectorType &) const override;
+  void virtual extractVelocity(VectorType &) const override;
+
+private:
+  MatrixType const &A;
+  MatrixType const &B;
+  VectorType u;
+  VectorType v;
+  VectorType a;
+  Dune::BitSetVector<dim> const &dirichletNodes;
+  FunctionType const &dirichletFunction;
+
+  VectorType u_old;
+  VectorType v_old;
+  VectorType a_old;
+
+  double tau;
+
+  bool postProcessCalled = false;
+};
+#endif
-- 
GitLab