From b8476178917fa8f5d71152728f7a23812d3cbbec Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 17 Sep 2012 15:51:30 +0200
Subject: [PATCH] Compile time stepping separately

---
 src/.gitignore         |   2 +
 src/Makefile.am        |   3 +-
 src/one-body-sample.cc |   2 +-
 src/timestepping.org   | 653 ++++++++++++++++++++++++-----------------
 4 files changed, 381 insertions(+), 279 deletions(-)

diff --git a/src/.gitignore b/src/.gitignore
index ed1309bc..c52dff94 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -1 +1,3 @@
 timestepping.cc
+timestepping.hh
+timestepping_tmpl.cc
diff --git a/src/Makefile.am b/src/Makefile.am
index b49ae312..9b222c90 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -10,6 +10,7 @@ SOURCES = \
 	compute_state_ruina.cc \
 	mysolver.cc \
 	one-body-sample.cc \
+	timestepping.cc \
 	vtk.cc
 
 ## 2D
@@ -101,7 +102,7 @@ DISTCHECK_CONFIGURE_FLAGS = \
 
 include $(top_srcdir)/am/global-rules
 
-BUILT_SOURCES = timestepping.cc
+BUILT_SOURCES = timestepping.hh timestepping.cc
 $(srcdir)/timestepping.cc: $(srcdir)/timestepping.org
 	emacs -Q --batch --eval \
 	  "(let (vc-handled-backends) \
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index dafab2dd..3fc2e190 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -75,7 +75,7 @@
 #include "enum_state_model.cc"
 #include "enum_scheme.cc"
 
-#include "timestepping.cc"
+#include "timestepping.hh"
 
 int const dim = DIM;
 
diff --git a/src/timestepping.org b/src/timestepping.org
index 911af73e..ced45917 100644
--- a/src/timestepping.org
+++ b/src/timestepping.org
@@ -1,17 +1,33 @@
-#+name: includes
+* Preamble
+#+name: includes.hh
 #+begin_src c++
+  #include <dune/common/bitsetvector.hh>
+#+end_src
+#+name: preamble.cc
+#+begin_src c++
+  #ifdef HAVE_CONFIG_H
+  #include "config.h"
+  #endif
+  
   #include <dune/fufem/arithmetic.hh>
+  #include "timestepping.hh"
 #+end_src
-
-#+name: preamble
+#+name: preamble.tex
 #+begin_src latex
   \documentclass{scrartcl}
   \usepackage{amsmath}
   \begin{document}
 #+end_src
-
+#+name: virtual_function_declaration
+#+begin_src c++
+  void virtual setup(VectorType const &, double, double, VectorType &,
+                     VectorType &, MatrixType &);
+  void virtual postProcess(VectorType const &);
+  void virtual extractDisplacement(VectorType &) const;
+  void virtual extractVelocity(VectorType &) const;
+#+end_src
 * Abstract TimeSteppingScheme
-#+name: TimeSteppingScheme
+#+name: TimeSteppingScheme.hh
 #+begin_src c++
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   class TimeSteppingScheme
@@ -25,20 +41,14 @@
           VectorType &problem_iterate,
           MatrixType &problem_A) = 0;
   
-    void virtual
-    postProcess(VectorType const &problem_iterate) = 0;
-  
-    void virtual
-    extractDisplacement(VectorType &displacement) const = 0;
-  
-    void virtual
-    extractVelocity(VectorType &velocity) const = 0;
+    void virtual postProcess(VectorType const &problem_iterate) = 0;
+    void virtual extractDisplacement(VectorType &displacement) const = 0;
+    void virtual extractVelocity(VectorType &velocity) const = 0;
   };
 #+end_src
-
 * TimeSteppingScheme: Implicit Euler
-#+name: ImplicitEuler
-#+begin_src c++
+#+name: ImplicitEuler.hh
+#+begin_src c++ noweb: yes
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   class ImplicitEuler : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
   {
@@ -47,90 +57,109 @@
                   VectorType const &_u_initial,
                   VectorType const &_ud_initial,
                   Dune::BitSetVector<dim> const &_dirichletNodes,
-                  FunctionType const &_dirichletFunction)
-      : A(_A),
-        u(_u_initial),
-        ud(_ud_initial),
-        dirichletNodes(_dirichletNodes),
-        dirichletFunction(_dirichletFunction)
-    {}
+                  FunctionType const &_dirichletFunction);
+
+  <<virtual_function_declaration>>
   
-    void virtual
-    setup(VectorType const &ell,
-          double _tau,
-          double time,
-          VectorType &problem_rhs,
-          VectorType &problem_iterate,
-          MatrixType &problem_A)
-    {
-      postProcessCalled = false;
+  private:
+    MatrixType const &A;
+    VectorType u;
+    VectorType ud;
+    Dune::BitSetVector<dim> const &dirichletNodes;
+    FunctionType const &dirichletFunction;
   
-      tau = _tau;
+    VectorType u_old;
+    VectorType ud_old;
   
-      // This function is only called once for every timestep
-      ud_old = ud;
-      u_old = u;
+    double tau;
   
-      problem_rhs = ell;
-      A.mmv(u_old, problem_rhs);
+    bool postProcessCalled = false;
+  };
+#+end_src
+#+name: ImplicitEuler.cc
+#+begin_src c++
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
+  ImplicitEuler(MatrixType const &_A,
+                VectorType const &_u_initial,
+                VectorType const &_ud_initial,
+                Dune::BitSetVector<dim> const &_dirichletNodes,
+                FunctionType const &_dirichletFunction)
+    : A(_A),
+      u(_u_initial),
+      ud(_ud_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction)
+  {}
   
-      // For fixed tau, we'd only really have to do this once
-      problem_A = A;
-      problem_A *= tau;
+  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_A)
+  {
+    postProcessCalled = false;
   
-      // ud_old makes a good initial iterate; we could use anything, though
-      problem_iterate = ud_old;
+    tau = _tau;
   
-      for (size_t i=0; i < dirichletNodes.size(); ++i)
-        if (dirichletNodes[i].count() == dim) {
-          problem_iterate[i] = 0;
-          dirichletFunction.evaluate(time, problem_iterate[i][0]);
-        } else if (dirichletNodes[i][1])
-          problem_iterate[i][1] = 0; // Y direction prescribed
-    }
+    // This function is only called once for every timestep
+    ud_old = ud;
+    u_old = u;
   
-    void virtual postProcess(VectorType const &problem_iterate)
-    {
-      postProcessCalled = true;
+    problem_rhs = ell;
+    A.mmv(u_old, problem_rhs);
   
-      ud = problem_iterate;
+    // For fixed tau, we'd only really have to do this once
+    problem_A = A;
+    problem_A *= tau;
   
-      u = u_old;
-      Arithmetic::addProduct(u, tau, ud);
-    }
+    // ud_old makes a good initial iterate; we could use anything, though
+    problem_iterate = ud_old;
   
-    void virtual extractDisplacement(VectorType &displacement) const
-    {
-      if (!postProcessCalled)
-        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+    for (size_t i=0; i < dirichletNodes.size(); ++i)
+      if (dirichletNodes[i].count() == dim) {
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+      } else if (dirichletNodes[i][1])
+        problem_iterate[i][1] = 0; // Y direction prescribed
+  }
   
-      displacement = u;
-    }
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::
+  postProcess(VectorType const &problem_iterate)
+  {
+    postProcessCalled = true;
   
-    void virtual extractVelocity(VectorType &velocity) const
-    {
-      if (!postProcessCalled)
-        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+    ud = problem_iterate;
   
-      velocity = ud;
-    }
+    u = u_old;
+    Arithmetic::addProduct(u, tau, ud);
+  }
   
-  private:
-    MatrixType const &A;
-    VectorType u;
-    VectorType ud;
-    Dune::BitSetVector<dim> const &dirichletNodes;
-    FunctionType const &dirichletFunction;
+  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!");
   
-    VectorType u_old;
-    VectorType ud_old;
+    displacement = u;
+  }
   
-    double tau;
+  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 = ud;
+  }
   
-    bool postProcessCalled = false;
-  };
 #+end_src
-
 * TimeSteppingScheme: Implicit Two-Step
 #+begin_src latex :tangle twostep.tex :noweb yes
   \documentclass{scrartcl}
@@ -154,9 +183,8 @@
   \end{equation*}
   \end{document}
 #+end_src
-
-#+name: ImplicitTwoStep
-#+begin_src c++
+#+name: ImplicitTwoStep.hh
+#+begin_src c++ noweb: yes
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   class ImplicitTwoStep : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
   {
@@ -165,108 +193,9 @@
                     VectorType const &_u_initial,
                     VectorType const &_ud_initial,
                     Dune::BitSetVector<dim> const &_dirichletNodes,
-                    FunctionType const &_dirichletFunction)
-      : A(_A),
-        u(_u_initial),
-        ud(_ud_initial),
-        dirichletNodes(_dirichletNodes),
-        dirichletFunction(_dirichletFunction)
-    {}
-  
-    // FIXME: handle case run == 1
-    void virtual
-    setup(VectorType const &ell,
-          double _tau,
-          double time,
-          VectorType &problem_rhs,
-          VectorType &problem_iterate,
-          MatrixType &problem_A)
-    {
-      postProcessCalled = false;
-  
-      tau = _tau;
-  
-      // This function is only called once for every timestep
-      u_old_old = u_old;
-      ud_old = ud;
-      u_old = u;
-  
-      switch (state) {
-        // Perform an implicit Euler step since we lack information
-      case NO_SETUP:
-        state = FIRST_SETUP;
-  
-        problem_rhs = ell;
-        A.mmv(u_old, problem_rhs);
-  
-        problem_A = A;
-        problem_A *= tau;
-        break;
-      case FIRST_SETUP:
-        state = SECOND_SETUP;
-        // FALLTHROUGH
-      case SECOND_SETUP:
-        problem_rhs = ell;
-        A.usmv(-4.0/3.0, u_old, problem_rhs);
-        A.usmv(+1.0/3.0, u_old_old, problem_rhs);
-  
-        // For fixed tau, we'd only really have to do this once
-        problem_A = A;
-        problem_A *= 2.0/3.0 * tau;
-        break;
-      default:
-        assert(false);
-      }
-  
-      // ud_old makes a good initial iterate; we could use anything, though
-      problem_iterate = ud_old;
-  
-      for (size_t i=0; i < dirichletNodes.size(); ++i)
-        if (dirichletNodes[i].count() == dim) {
-          problem_iterate[i] = 0;
-          dirichletFunction.evaluate(time, problem_iterate[i][0]);
-        } else if (dirichletNodes[i][1])
-          problem_iterate[i][1] = 0; // Y direction prescribed
-    }
-  
-    void virtual postProcess(VectorType const &problem_iterate)
-    {
-      postProcessCalled = true;
-  
-      ud = problem_iterate;
-  
-      switch (state) {
-      case FIRST_SETUP:
-        u = u_old;
-        Arithmetic::addProduct(u, tau, ud);
-        break;
-      case SECOND_SETUP:
-        u = 0.0;
-        Arithmetic::addProduct(u, tau, ud);
-        Arithmetic::addProduct(u, 2.0, u_old);
-        Arithmetic::addProduct(u, -.5, u_old_old);
-        u *= 2.0/3.0;
-        break;
-      default:
-        assert(false);
-      }
-    }
-  
-    void virtual extractDisplacement(VectorType &displacement) const
-    {
-      if (!postProcessCalled)
-        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-  
-      displacement = u;
-    }
-  
-    void virtual extractVelocity(VectorType &velocity) const
-    {
-      if (!postProcessCalled)
-        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-  
-      velocity = ud;
-    }
+                    FunctionType const &_dirichletFunction);
+
+  <<virtual_function_declaration>>
   
   private:
     MatrixType const &A;
@@ -288,10 +217,126 @@
     state_type state = NO_SETUP;
   };
 #+end_src
-
+#+name: ImplicitTwoStep.cc
+#+begin_src c++
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
+  ImplicitTwoStep(MatrixType const &_A,
+                  VectorType const &_u_initial,
+                  VectorType const &_ud_initial,
+                  Dune::BitSetVector<dim> const &_dirichletNodes,
+                  FunctionType const &_dirichletFunction)
+    : A(_A),
+      u(_u_initial),
+      ud(_ud_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction)
+  {}
+  
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
+  setup(VectorType const &ell,
+        double _tau,
+        double time,
+        VectorType &problem_rhs,
+        VectorType &problem_iterate,
+        MatrixType &problem_A)
+  {
+    postProcessCalled = false;
+  
+    tau = _tau;
+  
+    // This function is only called once for every timestep
+    u_old_old = u_old;
+    ud_old = ud;
+    u_old = u;
+  
+    switch (state) {
+      // Perform an implicit Euler step since we lack information
+    case NO_SETUP:
+      state = FIRST_SETUP;
+  
+      problem_rhs = ell;
+      A.mmv(u_old, problem_rhs);
+  
+      problem_A = A;
+      problem_A *= tau;
+      break;
+    case FIRST_SETUP:
+      state = SECOND_SETUP;
+      // FALLTHROUGH
+    case SECOND_SETUP:
+      problem_rhs = ell;
+      A.usmv(-4.0/3.0, u_old, problem_rhs);
+      A.usmv(+1.0/3.0, u_old_old, problem_rhs);
+  
+      // For fixed tau, we'd only really have to do this once
+      problem_A = A;
+      problem_A *= 2.0/3.0 * tau;
+      break;
+    default:
+      assert(false);
+    }
+  
+    // ud_old makes a good initial iterate; we could use anything, though
+    problem_iterate = ud_old;
+  
+    for (size_t i=0; i < dirichletNodes.size(); ++i)
+      if (dirichletNodes[i].count() == dim) {
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+      } else if (dirichletNodes[i][1])
+        problem_iterate[i][1] = 0; // Y direction prescribed
+  }
+  
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
+  postProcess(VectorType const &problem_iterate)
+  {
+    postProcessCalled = true;
+  
+    ud = problem_iterate;
+  
+    switch (state) {
+    case FIRST_SETUP:
+      u = u_old;
+      Arithmetic::addProduct(u, tau, ud);
+      break;
+    case SECOND_SETUP:
+      u = 0.0;
+      Arithmetic::addProduct(u, tau, ud);
+      Arithmetic::addProduct(u, 2.0, u_old);
+      Arithmetic::addProduct(u, -.5, u_old_old);
+      u *= 2.0/3.0;
+      break;
+    default:
+      assert(false);
+    }
+  }
+  
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void ImplicitTwoStep<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 ImplicitTwoStep<VectorType, MatrixType, FunctionType, dim>::
+  extractVelocity(VectorType &velocity) const
+  {
+    if (!postProcessCalled)
+      DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+  
+    velocity = ud;
+  }
+#+end_src
 * TimeSteppingScheme: Newmark
 #+begin_src latex :tangle newmark.tex :noweb yes
-  <<preamble>>
+  <<preamble.tex>>
   \noindent The Newmark scheme in its classical form with $\gamma = 1/2$
   and $\beta = 1/4$ reads
   \begin{align}
@@ -323,9 +368,8 @@
   from \eqref{eq:4}. The only unknown is now $\dot u_1$.
   \end{document}
 #+end_src
-
-#+name: Newmark
-#+begin_src c++
+#+name: Newmark.hh
+#+begin_src c++ noweb: yes
   template <class VectorType, class MatrixType, class FunctionType, int dim>
   class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
   {
@@ -336,86 +380,9 @@
             VectorType const &_ud_initial,
             VectorType const &_udd_initial,
             Dune::BitSetVector<dim> const &_dirichletNodes,
-            FunctionType const &_dirichletFunction)
-      : A(_A),
-        B(_B),
-        u(_u_initial),
-        ud(_ud_initial),
-        udd(_udd_initial),
-        dirichletNodes(_dirichletNodes),
-        dirichletFunction(_dirichletFunction)
-    {}
+            FunctionType const &_dirichletFunction);
   
-    void virtual
-    setup(VectorType const &ell,
-          double _tau,
-          double time,
-          VectorType &problem_rhs,
-          VectorType &problem_iterate,
-          MatrixType &problem_A)
-    {
-      postProcessCalled = false;
-  
-      tau = _tau;
-  
-      // This function is only called once for every timestep
-      udd_old = udd;
-      ud_old = ud;
-      u_old = u;
-  
-      problem_rhs = ell;
-      B.usmv( 2.0/tau,  ud_old, problem_rhs);
-      B.usmv(     1.0, udd_old, problem_rhs);
-      A.usmv(    -1.0,   u_old, problem_rhs);
-      A.usmv(-tau/2.0,  ud_old, problem_rhs);
-  
-      // For fixed tau, we'd only really have to do this once
-      problem_A  = A;
-      problem_A *= tau/2.0;
-      Arithmetic::addProduct(problem_A, 2.0/tau, B);
-  
-      // ud_old makes a good initial iterate; we could use anything, though
-      problem_iterate = ud_old;
-  
-      for (size_t i=0; i < dirichletNodes.size(); ++i)
-        if (dirichletNodes[i].count() == dim) {
-          problem_iterate[i] = 0;
-          dirichletFunction.evaluate(time, problem_iterate[i][0]);
-        } else if (dirichletNodes[i][1])
-          problem_iterate[i][1] = 0; // Y direction prescribed
-    }
-  
-    void virtual postProcess(VectorType const &problem_iterate)
-    {
-      postProcessCalled = true;
-  
-      ud = problem_iterate;
-  
-      u = u_old;
-      Arithmetic::addProduct(u, tau/2.0, ud);
-      Arithmetic::addProduct(u, tau/2.0, ud_old);
-  
-      udd = 0;
-      Arithmetic::addProduct(udd,  2.0/tau, ud);
-      Arithmetic::addProduct(udd, -2.0/tau, ud_old);
-      Arithmetic::addProduct(udd, -1.0,     udd_old);
-    }
-  
-    void virtual extractDisplacement(VectorType &displacement) const
-    {
-      if (!postProcessCalled)
-        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-  
-      displacement = u;
-    }
-  
-    void virtual extractVelocity(VectorType &velocity) const
-    {
-      if (!postProcessCalled)
-        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
-  
-      velocity = ud;
-    }
+  <<virtual_function_declaration>>
   
   private:
     MatrixType const &A;
@@ -434,16 +401,148 @@
   
     bool postProcessCalled = false;
   };
+#+end_src
+#+name: Newmark.cc
+#+begin_src c++
+  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 &_ud_initial,
+          VectorType const &_udd_initial,
+          Dune::BitSetVector<dim> const &_dirichletNodes,
+          FunctionType const &_dirichletFunction)
+    : A(_A),
+      B(_B),
+      u(_u_initial),
+      ud(_ud_initial),
+      udd(_udd_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunction(_dirichletFunction)
+  {}
+  
+  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_A)
+  {
+    postProcessCalled = false;
+  
+    tau = _tau;
   
+    // This function is only called once for every timestep
+    udd_old = udd;
+    ud_old = ud;
+    u_old = u;
+  
+    problem_rhs = ell;
+    B.usmv( 2.0/tau,  ud_old, problem_rhs);
+    B.usmv(     1.0, udd_old, problem_rhs);
+    A.usmv(    -1.0,   u_old, problem_rhs);
+    A.usmv(-tau/2.0,  ud_old, problem_rhs);
+  
+    // For fixed tau, we'd only really have to do this once
+    problem_A  = A;
+    problem_A *= tau/2.0;
+    Arithmetic::addProduct(problem_A, 2.0/tau, B);
+  
+    // ud_old makes a good initial iterate; we could use anything, though
+    problem_iterate = ud_old;
+  
+    for (size_t i=0; i < dirichletNodes.size(); ++i)
+      if (dirichletNodes[i].count() == dim) {
+        problem_iterate[i] = 0;
+        dirichletFunction.evaluate(time, problem_iterate[i][0]);
+      } else if (dirichletNodes[i][1])
+        problem_iterate[i][1] = 0; // Y direction prescribed
+  }
+  
+  template <class VectorType, class MatrixType, class FunctionType, int dim>
+  void Newmark<VectorType, MatrixType, FunctionType, dim>::
+  postProcess(VectorType const &problem_iterate)
+  {
+    postProcessCalled = true;
+  
+    ud = problem_iterate;
+  
+    u = u_old;
+    Arithmetic::addProduct(u, tau/2.0, ud);
+    Arithmetic::addProduct(u, tau/2.0, ud_old);
+  
+    udd = 0;
+    Arithmetic::addProduct(udd,  2.0/tau, ud);
+    Arithmetic::addProduct(udd, -2.0/tau, ud_old);
+    Arithmetic::addProduct(udd, -1.0,     udd_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 = ud;
+  }
 #+end_src
-
 * C++ code generation
-#+name: timestepping
+#+begin_src c++ :tangle timestepping.hh :noweb yes
+  // GENERATED -- DO NOT MODIFY
+  #ifndef DUNE_TECTONIC_TIMESTEPPING_HH
+  #define DUNE_TECTONIC_TIMESTEPPING_HH
+  <<includes.hh>>
+  
+  <<TimeSteppingScheme.hh>>
+  <<ImplicitEuler.hh>>
+  <<ImplicitTwoStep.hh>>
+  <<Newmark.hh>>
+  #endif
+#+end_src
 #+begin_src c++ :tangle timestepping.cc :noweb yes
-  <<includes>>
+  // GENERATED -- DO NOT MODIFY
+  <<preamble.cc>>
   
-  <<TimeSteppingScheme>>
-  <<ImplicitEuler>>
-  <<ImplicitTwoStep>>
-  <<Newmark>>
+  <<ImplicitEuler.cc>>
+  <<ImplicitTwoStep.cc>>
+  <<Newmark.cc>>
+
+  #include "timestepping_tmpl.cc"
+#+end_src
+
+#+begin_src c++ :tangle timestepping_tmpl.cc :noweb yes
+  // GENERATED -- DO NOT MODIFY
+  #ifndef DIM
+  #error DIM unset
+  #endif
+  
+  #include <dune/common/fmatrix.hh>
+  #include <dune/common/function.hh>
+  #include <dune/common/fvector.hh>
+  #include <dune/istl/bcrsmatrix.hh>
+  #include <dune/istl/bvector.hh>
+  
+  typedef Dune::FieldVector<double, DIM> SmallVector;
+  typedef Dune::FieldMatrix<double, DIM, DIM> SmallMatrix;
+  typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
+  typedef Dune::BlockVector<SmallVector> VectorType;
+  typedef Dune::VirtualFunction<double, double> FunctionType;
+  
+  template class ImplicitEuler<VectorType, MatrixType, FunctionType, DIM>;
+  template class ImplicitTwoStep<VectorType, MatrixType, FunctionType, DIM>;
+  template class Newmark<VectorType, MatrixType, FunctionType, DIM>;
 #+end_src
-- 
GitLab