diff --git a/src/.gitignore b/src/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..ed1309bc89fd12c000fb36a6c97122cc9cecec01
--- /dev/null
+++ b/src/.gitignore
@@ -0,0 +1 @@
+timestepping.cc
diff --git a/src/Makefile.am b/src/Makefile.am
index 2467aafbdfe3640ba86b5db998b0fa8ee7df0d13..0777d6eee9a04467214e8e35c9b31a42303b0205 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -141,3 +141,8 @@ DISTCHECK_CONFIGURE_FLAGS = \
CXX="$(CXX)" CC="$(CC)"
include $(top_srcdir)/am/global-rules
+
+$(srcdir)/timestepping.cc: $(srcdir)/timestepping.org
+ emacs -Q --batch --eval \
+ "(let (vc-handled-backends) \
+ (org-babel-tangle-file \"$<\" nil 'c++))"
diff --git a/src/timestepping.cc b/src/timestepping.cc
deleted file mode 100644
index 0a3e3dd7046c8791f386896dd54446157642499e..0000000000000000000000000000000000000000
--- a/src/timestepping.cc
+++ /dev/null
@@ -1,211 +0,0 @@
-#include <dune/fufem/arithmetic.hh>
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class TimeSteppingScheme {
-public:
- TimeSteppingScheme(VectorType const &_ell, MatrixType const &_A,
- VectorType const &_u_old,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction, double _time,
- double _tau)
- : ell(_ell),
- A(_A),
- u_old(_u_old),
- dirichletNodes(_dirichletNodes),
- dirichletFunction(_dirichletFunction),
- time(_time),
- tau(_tau) {}
-
- virtual ~TimeSteppingScheme() {}
-
- void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
- MatrixType &problem_A) const = 0;
-
- void virtual extractSolution(VectorType const &problem_iterate,
- VectorType &solution) const = 0;
-
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const = 0;
-
-protected:
- VectorType const ℓ
- MatrixType const &A;
- VectorType const &u_old;
- Dune::BitSetVector<dim> const &dirichletNodes;
- FunctionType const &dirichletFunction;
- double const time;
- double const tau;
-};
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class ImplicitEuler
- : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
-public:
- ImplicitEuler(
- VectorType const &_ell, MatrixType const &_A, VectorType const &_u_old,
- VectorType const *_u_old_old_ptr, // FIXME: this should not be necessary
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction, double _time, double _tau)
- : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
- _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
- u_old_old_ptr(_u_old_old_ptr) {}
-
- void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
- MatrixType &problem_A) const {
- problem_rhs = this->ell;
- this->A.mmv(this->u_old, problem_rhs);
-
- problem_A = this->A;
- problem_A *= this->tau;
-
- // Use the old velocity as an initial iterate if possible
- if (u_old_old_ptr) {
- problem_iterate = this->u_old;
- problem_iterate -= *u_old_old_ptr;
- problem_iterate /= this->tau;
- } else
- problem_iterate = 0.0;
-
- for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
- if (this->dirichletNodes[i].count() == dim) {
- problem_iterate[i] = 0;
- this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
- } else if (this->dirichletNodes[i][1])
- problem_iterate[i][1] = 0; // Y direction prescribed
- }
-
- void virtual extractSolution(VectorType const &problem_iterate,
- VectorType &solution) const {
- solution = this->u_old;
- solution.axpy(this->tau, problem_iterate);
- }
-
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const {
- velocity = problem_iterate;
- }
-
-private:
- VectorType const *u_old_old_ptr;
-};
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class ImplicitTwoStep
- : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
-public:
- ImplicitTwoStep(VectorType const &_ell, MatrixType const &_A,
- VectorType const &_u_old, VectorType const *_u_old_old_ptr,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction, double _time,
- double _tau)
- : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
- _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
- u_old_old_ptr(_u_old_old_ptr) {}
-
- void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
- MatrixType &problem_A) const {
- problem_rhs = this->ell;
- this->A.usmv(-4.0 / 3.0, this->u_old, problem_rhs);
- this->A.usmv(+1.0 / 3.0, *u_old_old_ptr, problem_rhs);
-
- problem_A = this->A;
- problem_A *= 2.0 / 3.0 * this->tau;
-
- // Use the old velocity as an initial iterate
- problem_iterate = this->u_old;
- problem_iterate -= *u_old_old_ptr;
- problem_iterate /= this->tau;
-
- for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
- if (this->dirichletNodes[i].count() == dim) {
- problem_iterate[i] = 0;
- this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
- } else if (this->dirichletNodes[i][1])
- problem_iterate[i][1] = 0; // Y direction prescribed
- }
-
- void virtual extractSolution(VectorType const &problem_iterate,
- VectorType &solution) const {
- solution = problem_iterate;
- solution *= this->tau;
- solution.axpy(2, this->u_old);
- solution.axpy(-.5, *u_old_old_ptr);
- solution *= 2.0 / 3.0;
-
- // Check if we split correctly
- {
- VectorType test = problem_iterate;
- test *= this->tau;
- test.axpy(-1.5, solution);
- test.axpy(+2, this->u_old);
- test.axpy(-.5, *u_old_old_ptr);
- assert(test.two_norm() < 1e-10);
- }
- }
-
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const {
- velocity = problem_iterate;
- }
-
-private:
- VectorType const *u_old_old_ptr;
-};
-
-template <class VectorType, class MatrixType, class FunctionType, int dim>
-class Newmark
- : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
-public:
- Newmark(VectorType const &_ell, MatrixType const &_A, MatrixType const &_B,
- VectorType const &_u_old, VectorType const &_ud_old,
- VectorType const &_udd_old,
- Dune::BitSetVector<dim> const &_dirichletNodes,
- FunctionType const &_dirichletFunction, double _time, double _tau)
- : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>(
- _ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
- B(_B),
- ud_old(_ud_old),
- udd_old(_udd_old) {}
-
- void virtual setup(VectorType &problem_rhs, VectorType &problem_iterate,
- MatrixType &problem_A) const {
- problem_rhs = this->ell;
- /* */ B.usmv(2.0 / this->tau, ud_old, problem_rhs);
- /* */ B.usmv(1.0, udd_old, problem_rhs);
- this->A.usmv(-1.0, this->u_old, problem_rhs);
- this->A.usmv(-this->tau / 2.0, ud_old, problem_rhs);
-
- // For fixed tau, we'd only really have to do this once
- problem_A = this->A;
- problem_A *= this->tau / 2.0;
- Arithmetic::addProduct(problem_A, 2.0 / this->tau, B);
-
- // ud_old makes a good initial iterate; we could use anything, though
- problem_iterate = ud_old;
-
- for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
- if (this->dirichletNodes[i].count() == dim) {
- problem_iterate[i] = 0;
- this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
- } else if (this->dirichletNodes[i][1])
- problem_iterate[i][1] = 0; // Y direction prescribed
- }
-
- // u1 = u0 + tau/2 (du1 + du0)
- void virtual extractSolution(VectorType const &problem_iterate,
- VectorType &solution) const {
- solution = this->u_old;
- Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
- Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
- }
-
- void virtual extractVelocity(VectorType const &problem_iterate,
- VectorType &velocity) const {
- velocity = problem_iterate;
- }
-
-private:
- MatrixType const &B;
- VectorType const &ud_old;
- VectorType const &udd_old;
-};
diff --git a/src/timestepping.org b/src/timestepping.org
new file mode 100644
index 0000000000000000000000000000000000000000..51e8b5d893fdc8f5cea65661bea1f2e88015478a
--- /dev/null
+++ b/src/timestepping.org
@@ -0,0 +1,339 @@
+#+name: includes
+#+begin_src c++
+ #include <dune/fufem/arithmetic.hh>
+#+end_src
+
+#+name: preamble
+#+begin_src latex
+ \documentclass{scrartcl}
+ \usepackage{amsmath}
+ \begin{document}
+#+end_src
+
+* Abstract TimeSteppingScheme
+#+name: TimeSteppingScheme
+#+begin_src c++
+ template <class VectorType, class MatrixType, class FunctionType, int dim>
+ class TimeSteppingScheme
+ {
+ public:
+ TimeSteppingScheme(VectorType const &_ell,
+ MatrixType const &_A,
+ VectorType const &_u_old,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction,
+ double _time,
+ double _tau)
+ : ell(_ell),
+ A(_A),
+ u_old(_u_old),
+ dirichletNodes(_dirichletNodes),
+ dirichletFunction(_dirichletFunction),
+ time(_time),
+ tau(_tau)
+ {}
+
+ virtual ~TimeSteppingScheme()
+ {}
+
+ void virtual setup(VectorType &problem_rhs,
+ VectorType &problem_iterate,
+ MatrixType &problem_A) const = 0;
+
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const = 0;
+
+ void virtual extractVelocity(VectorType const &problem_iterate,
+ VectorType &velocity) const = 0;
+
+ protected:
+ VectorType const ℓ
+ MatrixType const &A;
+ VectorType const &u_old;
+ Dune::BitSetVector<dim> const &dirichletNodes;
+ FunctionType const &dirichletFunction;
+ double const time;
+ double const tau;
+ };
+#+end_src
+
+* TimeSteppingScheme: Implicit Euler
+#+name: ImplicitEuler
+#+begin_src c++
+ template <class VectorType, class MatrixType, class FunctionType, int dim>
+ class ImplicitEuler : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
+ {
+ public:
+ ImplicitEuler(VectorType const &_ell,
+ MatrixType const &_A,
+ VectorType const &_u_old,
+ VectorType const *_u_old_old_ptr, // FIXME: this should not be necessary
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction,
+ double _time,
+ double _tau)
+ : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
+ (_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+ u_old_old_ptr(_u_old_old_ptr)
+ {}
+
+ void virtual
+ setup(VectorType &problem_rhs,
+ VectorType &problem_iterate,
+ MatrixType &problem_A) const
+ {
+ problem_rhs = this->ell;
+ this->A.mmv(this->u_old, problem_rhs);
+
+ problem_A = this->A;
+ problem_A *= this->tau;
+
+ // Use the old velocity as an initial iterate if possible
+ if (u_old_old_ptr) {
+ problem_iterate = this->u_old;
+ problem_iterate -= *u_old_old_ptr;
+ problem_iterate /= this->tau;
+ } else
+ problem_iterate = 0.0;
+
+ for (size_t i=0; i < this->dirichletNodes.size(); ++i)
+ if (this->dirichletNodes[i].count() == dim) {
+ problem_iterate[i] = 0;
+ this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
+ } else if (this->dirichletNodes[i][1])
+ problem_iterate[i][1] = 0; // Y direction prescribed
+ }
+
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const
+ {
+ solution = this->u_old;
+ solution.axpy(this->tau, problem_iterate);
+ }
+
+ void virtual extractVelocity(VectorType const &problem_iterate,
+ VectorType &velocity) const
+ {
+ velocity = problem_iterate;
+ }
+
+ private:
+ VectorType const *u_old_old_ptr;
+ };
+#+end_src
+
+* TimeSteppingScheme: Implicit Two-Step
+#+begin_src latex :tangle twostep.tex :noweb yes
+ \documentclass{scrartcl}
+ \usepackage{amsmath}
+ \begin{document}
+ We start out with
+ \begin{align}
+ a(u_1, v - \dot u_1) + j(v) - j(\dot u_1) \ge \ell(v-\dot u_1)
+ \end{align}
+ With the two-step implicit scheme
+ \begin{equation*}
+ \tau \dot u_1 = \frac 32 u_1 - 2 u_0 + \frac 12 u_{-1}
+ \end{equation*}
+ or equivalently
+ \begin{equation*}
+ \frac 23 \left( \tau \dot u_1 + 2 u_0 - \frac 12 u_{-1} \right) = u_1
+ \end{equation*}
+ we obtain
+ \begin{equation*}
+ \frac 23 \tau a(\dot u_1, v - \dot u_1) + j(v) - j(\dot u_1) \ge \ell(v-\dot u_1) - a\left(\frac 43 u_0 - \frac 13 u_{-1}, v - \dot u_1\right)
+ \end{equation*}
+ \end{document}
+#+end_src
+
+#+name: ImplicitTwoStep
+#+begin_src c++
+ template <class VectorType, class MatrixType, class FunctionType, int dim>
+ class ImplicitTwoStep : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
+ {
+ public:
+ ImplicitTwoStep(VectorType const &_ell,
+ MatrixType const &_A,
+ VectorType const &_u_old,
+ VectorType const *_u_old_old_ptr,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction,
+ double _time,
+ double _tau)
+ : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
+ (_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+ u_old_old_ptr(_u_old_old_ptr)
+ {}
+
+ void virtual setup(VectorType &problem_rhs,
+ VectorType &problem_iterate,
+ MatrixType &problem_A) const
+ {
+ problem_rhs = this->ell;
+ this->A.usmv(-4.0/3.0, this->u_old, problem_rhs);
+ this->A.usmv(+1.0/3.0, *u_old_old_ptr, problem_rhs);
+
+ problem_A = this->A;
+ problem_A *= 2.0/3.0 * this->tau;
+
+ // Use the old velocity as an initial iterate
+ problem_iterate = this->u_old;
+ problem_iterate -= *u_old_old_ptr;
+ problem_iterate /= this->tau;
+
+ for (size_t i=0; i < this->dirichletNodes.size(); ++i)
+ if (this->dirichletNodes[i].count() == dim) {
+ problem_iterate[i] = 0;
+ this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
+ } else if (this->dirichletNodes[i][1])
+ problem_iterate[i][1] = 0; // Y direction prescribed
+ }
+
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const
+ {
+ solution = problem_iterate;
+ solution *= this->tau;
+ solution.axpy(2, this->u_old);
+ solution.axpy(-.5, *u_old_old_ptr);
+ solution *= 2.0/3.0;
+
+ // Check if we split correctly
+ {
+ VectorType test = problem_iterate;
+ test *= this->tau;
+ test.axpy(-1.5, solution);
+ test.axpy(+2, this->u_old);
+ test.axpy(-.5, *u_old_old_ptr);
+ assert(test.two_norm() < 1e-10);
+ }
+ }
+
+ void virtual extractVelocity(VectorType const &problem_iterate,
+ VectorType &velocity) const
+ {
+ velocity = problem_iterate;
+ }
+
+ private:
+ VectorType const *u_old_old_ptr;
+ };
+#+end_src
+
+* TimeSteppingScheme: Newmark
+#+begin_src latex :tangle newmark.tex :noweb yes
+ <<preamble>>
+ \noindent The Newmark scheme in its classical form with $\gamma = 1/2$
+ and $\beta = 1/4$ reads
+ \begin{align}
+ \label{eq:1} \dot u_1
+ &= \dot u_0 + \frac \tau 2 (\ddot u_0 + \ddot u_1 )\\
+ \label{eq:2} u_1
+ &= u_0 + \tau \dot u_0 + \frac {\tau^2}4 ( \ddot u_0 + \ddot u_1 )\text.
+ \intertext{We can also write \eqref{eq:1} as}
+ \label{eq:3} \ddot u_1
+ &= \frac 2\tau (\dot u_1 - \dot u_0) - \ddot u_0
+ \intertext{so that it yields}
+ \label{eq:4} u_1
+ &= u_0 + \tau \dot u_0 + \frac {\tau^2}4 \ddot u_0 + \frac {\tau^2}4 \left( \frac 2\tau (\dot u_1 - \dot u_0) - \ddot u_0\right)\\
+ &= u_0 + \tau \dot u_0 + \frac \tau 2 (\dot u_1 - \dot u_0)\nonumber\\
+ &= u_0 + \frac \tau 2 (\dot u_1 + \dot u_0)\nonumber
+ \end{align}
+ in conjunction with \eqref{eq:2}. If we now consider the EVI
+ \begin{align*}
+ b(\ddot u_1, v - \dot u_1) + a(u_1, v - \dot u_1) + j(v) - j(\dot u_1)
+ &\ge \ell (v - \dot u_1)
+ \intertext{with unknowns $u_1$, $\dot u_1$, and $\ddot u_1$, we first derive}
+ \frac 2\tau b(\dot u_1, v - \dot u_1) + a(u_1, v - \dot u_1) + j(v) - j(\dot u_1)
+ &\ge \ell (v - \dot u_1) + b\left(\frac 2\tau \dot u_0 + \ddot u_0, v - \dot u_1\right)
+ \intertext{from \eqref{eq:3} and then}
+ \frac 2\tau b(\dot u_1, v - \dot u_1) + \frac \tau 2 a(\dot u_1, v - \dot u_1) + j(v) - j(\dot u_1)
+ &\ge \ell (v - \dot u_1) + b\left(\frac 2\tau \dot u_0 + \ddot u_0, v - \dot u_1\right)\\
+ &\quad - a\left(u_0 + \frac \tau 2 \dot u_0, v - \dot u_1\right)
+ \end{align*}
+ from \eqref{eq:4}. The only unknown is now $\dot u_1$.
+ \end{document}
+#+end_src
+
+#+name: Newmark
+#+begin_src c++
+ template <class VectorType, class MatrixType, class FunctionType, int dim>
+ class Newmark : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
+ {
+ public:
+ Newmark(VectorType const &_ell,
+ MatrixType const &_A,
+ MatrixType const &_B,
+ VectorType const &_u_old,
+ VectorType const &_ud_old,
+ VectorType const &_udd_old,
+ Dune::BitSetVector<dim> const &_dirichletNodes,
+ FunctionType const &_dirichletFunction,
+ double _time,
+ double _tau)
+ : TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim>
+ (_ell, _A, _u_old, _dirichletNodes, _dirichletFunction, _time, _tau),
+ B(_B),
+ ud_old(_ud_old),
+ udd_old(_udd_old)
+ {}
+
+ void virtual
+ setup(VectorType &problem_rhs,
+ VectorType &problem_iterate,
+ MatrixType &problem_A) const
+ {
+ problem_rhs = this->ell;
+ /* */ B.usmv( 2.0/this->tau, ud_old, problem_rhs);
+ /* */ B.usmv( 1.0, udd_old, problem_rhs);
+ this->A.usmv( -1.0, this->u_old, problem_rhs);
+ this->A.usmv(-this->tau/2.0, ud_old, problem_rhs);
+
+ // For fixed tau, we'd only really have to do this once
+ problem_A = this->A;
+ problem_A *= this->tau/2.0;
+ Arithmetic::addProduct(problem_A, 2.0/this->tau, B);
+
+ // ud_old makes a good initial iterate; we could use anything, though
+ problem_iterate = ud_old;
+
+ for (size_t i=0; i < this->dirichletNodes.size(); ++i)
+ if (this->dirichletNodes[i].count() == dim) {
+ problem_iterate[i] = 0;
+ this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
+ } else if (this->dirichletNodes[i][1])
+ problem_iterate[i][1] = 0; // Y direction prescribed
+ }
+
+ // u1 = u0 + tau/2 (du1 + du0)
+ void virtual extractSolution(VectorType const &problem_iterate,
+ VectorType &solution) const
+ {
+ solution = this->u_old;
+ Arithmetic::addProduct(solution, this->tau / 2.0, problem_iterate); // ud
+ Arithmetic::addProduct(solution, this->tau / 2.0, ud_old);
+ }
+
+ void virtual extractVelocity(VectorType const &problem_iterate,
+ VectorType &velocity) const
+ {
+ velocity = problem_iterate;
+ }
+
+ private:
+ MatrixType const &B;
+ VectorType const &ud_old;
+ VectorType const &udd_old;
+ };
+#+end_src
+
+* C++ code generation
+#+name: timestepping
+#+begin_src c++ :tangle timestepping.cc :noweb yes
+ <<includes>>
+
+ <<TimeSteppingScheme>>
+ <<ImplicitEuler>>
+ <<ImplicitTwoStep>>
+ <<Newmark>>
+#+end_src