From c7707ff71b53c73f81b041e26bfa34517a06fcb8 Mon Sep 17 00:00:00 2001
From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 18 Dec 2013 20:21:55 +0000
Subject: [PATCH] take over recent changes to step-26

git-svn-id: https://svn.dealii.org/branches/releases/Branch-8-1@32055 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-26/doc/builds-on   |   2 +-
 deal.II/examples/step-26/doc/intro.dox   | 321 ++++++++++++++++++++++-
 deal.II/examples/step-26/doc/kind        |   2 +-
 deal.II/examples/step-26/doc/results.dox | 114 ++++++++
 deal.II/examples/step-26/doc/tooltip     |   2 +-
 deal.II/examples/step-26/step-26.cc      | 280 +++++++++++++-------
 6 files changed, 624 insertions(+), 97 deletions(-)

diff --git a/deal.II/examples/step-26/doc/builds-on b/deal.II/examples/step-26/doc/builds-on
index 48a0f73876..17402734c7 100644
--- a/deal.II/examples/step-26/doc/builds-on
+++ b/deal.II/examples/step-26/doc/builds-on
@@ -1 +1 @@
-step-4
+step-6
diff --git a/deal.II/examples/step-26/doc/intro.dox b/deal.II/examples/step-26/doc/intro.dox
index 51dcc23237..7aa1ab6697 100644
--- a/deal.II/examples/step-26/doc/intro.dox
+++ b/deal.II/examples/step-26/doc/intro.dox
@@ -1,8 +1,259 @@
 <a name="Intro"></a>
 <h1>Introduction</h1>
 
+This program implements the heat equation
+@f{align*}
+  \frac{\partial u(\mathbf x, t)}{\partial t}
+  -
+  \Delta u(\mathbf x, t)
+  &=
+  f(\mathbf x, t),
+  \qquad\qquad &&
+  \forall \mathbf x \in \Omega, t\in (0,T), 
+  \\
+  u(\mathbf x, 0) &= u_0(\mathbf x) &&
+  \forall \mathbf x \in \Omega, \\
+  \\
+  u(\mathbf x, t) &= g(\mathbf x,t) &&
+  \forall \mathbf x \in \partial\Omega, t \in (0,T).
+@f}
+In some sense, this equation is simpler than the ones we have discussed in the
+preceding programs step-23, step-24, step-25, namely the wave equation. This
+is due to the fact that the heat equation smoothes out the solution over time,
+and is consequently more forgiving in many regards. For example, when using
+implicit time stepping methods, we can actually take large time steps, we have
+less trouble with the small disturbances we introduce through adapting the
+mesh every few time steps, etc.
+
+Our goal here will be to solve the equations above using the theta-scheme that
+discretizes the equation in time using the following approach, where we would
+like $u^n(\mathbf x)$ to approximate $u(\mathbf x, t_n)$ at some time $t_n$:
+@f{align*}
+  \frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n}
+  -
+  \left[
+  (1-\theta)\Delta u^{n-1}(\mathbf x)
+  +
+  \theta\Delta u^n(\mathbf x)
+  \right]
+  &=
+  \left[
+  (1-\theta)f(\mathbf x, t_{n-1})
+  +
+  \theta f(\mathbf x, t_n)
+  \right].  
+@f}
+Here, $k_n=t_n-t_{n-1}$ is the time step size. The theta-scheme generalizes
+the explicit Euler ($\theta=0$), implicit Euler ($\theta=1$) and
+Crank-Nicolson ($\theta=\frac 12$) time discretizations. Since the latter has
+the highest convergence order, we will choose $\theta=\frac 12$ in the program
+below, but make it so that playing with this parameter remains simple.
+
+Given this time discretization, space discretization happens as it always
+does, by multiplying with test functions, integrating by parts, and then
+restricting everything to a finite dimensional subspace. This yields the
+following set of fully discrete equations after multiplying through with
+$k_n$:
+@f{align*}
+  M U^n-MU^{n-1}
+  +
+  k_n \left[
+  (1-\theta)A U^{n-1}
+  +
+  \theta A U^n
+  \right]
+  &=
+  k_n
+  \left[
+  (1-\theta)F^{n-1}
+  +
+  \theta F^n
+  \right],
+@f}
+where $M$ is the mass matrix and $A$ is the stiffness matrix that results from
+discretizing the Laplacian. Bringing all known quantities to the right hand
+side yields the linear system we have to solve in every step:
+@f{align*}
+  (M
+  +
+  k_n \theta A) U^n
+  &=
+  MU^{n-1}
+  -
+  k_n
+  (1-\theta)A U^{n-1}
+  +
+  k_n
+  \left[
+  (1-\theta)F^{n-1}
+  +
+  \theta F^n
+  \right].
+@f}
+The linear system on the left hand side is symmetric and positive definite, so
+we should have no trouble solving it with the Conjugate Gradient method.
+
+We can start the iteration above if we have the set of nodal coefficients
+$U^0$ at the initial time. Here, we take the ones we get by interpolating the
+initial values $u_0(\mathbf x)$ onto the mesh used for the first time step. We
+will also need to choose a time step; we will here just choose it as fixed,
+but clearly advanced simulators will want to choose it adaptively. We will
+briefly come back to this in the <a href="#Results">results section
+below</a>.
+
+
+<h3> Adapting meshes for time dependent problems </h3>
+
+When solving the wave equation and its variants in the previous few programs,
+we kept the mesh fixed. Just as for stationary equations, one can make a good
+case that this is not the smartest approach and that significant savings can
+be had by adapting the mesh. There are, however, significant difficulties
+compared to the stationary case. Let us go through them in turn:
 
-<h4> Verifying whether the code is correct </h4>
+<ul>
+  <li><i>Time step size and minimal mesh size</i>: For stationary problems, the
+  general approach is "make the mesh as fine as it is necessary". For problems
+  with singularities, this often leads to situations where we get many levels
+  of refinement into corners or along interfaces. The very first tutorial to
+  use adaptive meshes, step-6, is a point in case already.
+
+  However, for time dependent problems, we typically need to choose the time
+  step related to the mesh size. For explicit time discretizations, this is
+  obvious, since we need to respect a CFL condition that ties the time step
+  size to the smallest mesh size. For implicit time discretizations, no such
+  hard restriction exists, but in practice we still want to make the time step
+  smaller if we make the mesh size smaller since we typically have error
+  estimates of the form $\|e\| \le {\cal O}(k^p + h^q)$ where $p,q$ are the
+  convergence orders of the time and space discretization, respectively. We
+  can only make the error small if we decrease both terms. Ideally, an
+  estimate like this would suggest to choose $k \propto h^{q/p}$. Because, at
+  least for problems with non-smooth solutions, the error is typically
+  localized in the cells with the smallest mesh size, we have to indeed choose
+  $k \propto h_{\text{min}}^{q/p}$, using the <i>smallest</i> mesh size.
+
+  The consequence is that refining the mesh further in one place implies not
+  only the moderate additional effort of increasing the number of degrees of
+  freedom slightly, but also the much larger effort of having the solve the
+  <i>global</i> linear system more often because of the smaller time step.
+
+  In practice, one typically deals with this by acknowledging that we can not
+  make the time step arbitrarily small, and consequently can not make the
+  local mesh size arbitrarily small. Rather, we set a maximal level of
+  refinement and when we flag cells for refinement, we simply do not refine
+  those cells whose children would exceed this maximal level of refinement.
+
+  There is a similar problem in that we will choose a right hand side that
+  will switch on in different parts of the domain at different times. To avoid
+  being caught flat footed with too coarse a mesh in areas where we suddenly
+  need a finer mesh, we will also enforce in our program a <i>minimal</i> mesh
+  refinement level.
+
+  <li><i>Test functions from different meshes</i>: Let us consider again the
+  semi-discrete equations we have written down above:
+  @f{align*}
+    \frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n}
+    -
+    \left[
+    (1-\theta)\Delta u^{n-1}(\mathbf x)
+    +
+    \theta\Delta u^n(\mathbf x)
+    \right]
+    &=
+    \left[
+    (1-\theta)f(\mathbf x, t_{n-1})
+    +
+    \theta f(\mathbf x, t_n)
+    \right].  
+  @f}
+  We can here consider $u^{n-1}$ as data since it has presumably been computed
+  before. Now, let us replace 
+  @f{align*}
+    u^n(\mathbf x)\approx u_h^n(\mathbf x) 
+    =
+    \sum_j U^n \varphi_j(\mathbf x),
+  @f}
+  multiply with test functions $\varphi_i(\mathbf x)$ and integrate by parts
+  where necessary. In a process as outlined above, this would yield
+  @f{align*}
+    \sum_j
+    (M
+    +
+    k_n \theta A)_{ij} U^n_j
+    &=
+    (\varphi_i, u_h^{n-1})
+    -
+    k_n
+    (1-\theta)(\nabla \varphi_i, \nabla u_h^{n-1})
+    +
+    k_n
+    \left[
+    (1-\theta)F^{n-1}
+    +
+    \theta F^n
+    \right].
+  @f}
+  Now imagine that we have changed the mesh between time steps $n-1$ and
+  $n$. Then the problem is that the basis functions we use for $u_h^n$ and
+  $u^{n-1}$ are different! This pertains to the terms on the right hand side,
+  the first of which we could more clearly write as (the second follows the
+  same pattern)
+  @f{align*}
+    (\varphi_i, u_h^{n-1})
+    =
+    (\varphi_i^n, u_h^{n-1})
+    =
+    \sum_{j=1}^{N_{n-1}}
+    (\varphi_i^n, \varphi_j^{n-1}) U^{n-1}_j,
+    \qquad\qquad
+    i=1\ldots N_n.
+  @f}
+  If the meshes used in these two time steps are the same, then
+  $(\varphi_i^n, \varphi_j^{n-1})$ forms a square mass matrix
+  $M_{ij}$. However, if the meshes are not the same, then in general the matrix
+  is rectangular. Worse, it is difficult to even compute these integrals
+  because if we loop over the cells of the mesh at time step $n$, then we need
+  to evaluate $\varphi_j^{n-1}$ at the quadrature points of these cells, but
+  they do not necessarily correspond to the cells of the mesh at time step
+  $n-1$ and $\varphi_j^{n-1}$ is not defined via these cells; the same of
+  course applies if we wanted to compute the integrals via integration on the
+  cells of mesh $n-1$.
+
+  In any case, what we have to face is a situation where we need to integrate
+  shape functions defined on two different meshes. This can be done, and is in
+  fact demonstrated in step-28, but the process is at best described by the
+  word "awkward".
+
+  In practice, one does not typically want to do this. Rather, we avoid the
+  whole situation by interpolating the solution from the old to the new mesh
+  every time we adapt the mesh. In other words, rather than solving the
+  equations above, we instead solve the problem
+  @f{align*}
+    \sum_j
+    (M
+    +
+    k_n \theta A)_{ij} U^n_j
+    &=
+    (\varphi_i, I_h^n u_h^{n-1})
+    -
+    k_n
+    (1-\theta)(\nabla \varphi_i, \nabla I_h^n u_h^{n-1})
+    +
+    k_n
+    \left[
+    (1-\theta)F^{n-1}
+    +
+    \theta F^n
+    \right],
+  @f}
+  where $I_h^n$ is the interpolation operator onto the finite element space
+  used in time step $n$. This is not the optimal approach since it introduces
+  an additional error besides time and space discretization, but it is a
+  pragmatic one that makes it feasible to do time adapting meshes.
+</ul>
+
+
+
+<h3> What could possibly go wrong? Verifying whether the code is correct </h3>
 
 There are a number of things one can typically get wrong when implementing a
 finite element code. In particular, for time dependent problems, the following
@@ -113,3 +364,71 @@ $\frac{1}{(n_x^2+n_y^2)\pi^2}$.
 Once we have verified that the time integration and right hand side handling
 are correct using this scheme, we can go on to verifying that we have the
 boundary values correct, using a very similar approach.
+
+
+
+<h3> The testcase </h3>
+
+Solving the heat equation on a simple domain with a simple right hand side
+almost always leads to solutions that are exceedingly boring, since they
+become very smooth very quickly and then do not move very much any
+more. Rather, we here solve the equation on the L-shaped domain with zero
+Dirichlet boundary values and zero initial conditions, but as right hand side
+we choose
+@f{align*}
+  f(\mathbf x, t)
+  =
+  \left\{
+  \begin{array}{ll}
+    \chi_1(\mathbf x)
+    & \text{if \(0\le t \le 0.2\tau\) or \(\tau\le t \le 1.2\tau\) or \(2\tau\le t
+    \le 2.2\tau\), etc}
+    \\
+    \chi_2(\mathbf x)
+    & \text{if \(0.5\le t \le 0.7\tau\) or \(1.5\tau\le t \le 1.7\tau\) or \(2.5\tau\le t
+    \le 2.7\tau\), etc}
+    \\
+    0
+    & \text{otherwise}
+  \end{array}
+  \right.
+@f}
+Here,
+@f{align*}
+  \chi_1(\mathbf x) &= 
+  \left\{
+  \begin{array}{ll}
+    1
+    & \text{if \(x>0.5\) and \(y>-0.5\)}
+    \\
+    0
+    & \text{otherwise}
+  \end{array}
+  \right.
+  \\
+  \chi_2(\mathbf x) &= 
+  \left\{
+  \begin{array}{ll}
+    1
+    & \text{if \(x>-0.5\) and \(y>0.5\)}
+    \\
+    0
+    & \text{otherwise}
+  \end{array}
+  \right.
+@f}
+In other words, in every period of lenght $\tau$, the right hand side first
+flashes on in domain 1, then off completely, then on in domain 2, then off
+completely again. This pattern is probably best observed via the little
+animation of the solution shown in the <a href="#Results">results
+section</a>.
+
+If you interpret the heat equation as finding the spatially and temporally
+variable temperature distribution of a conducting solid, then the test case
+above corresponds to an L-shaped body where we keep the boundary at zero
+temperature, and heat alternatingly in two parts of the domain. While heating
+is in effect, the temperature rises in these places, after which it diffuses
+and diminishes again. The point of these initial conditions is that they
+provide us with a solution that has singularities both in time (when sources
+switch on and off) as well as time (at the reentrant corner as well as at the
+edges and corners of the regions where the source acts).
diff --git a/deal.II/examples/step-26/doc/kind b/deal.II/examples/step-26/doc/kind
index c1d9154931..86a44aa1ef 100644
--- a/deal.II/examples/step-26/doc/kind
+++ b/deal.II/examples/step-26/doc/kind
@@ -1 +1 @@
-techniques
+time dependent
diff --git a/deal.II/examples/step-26/doc/results.dox b/deal.II/examples/step-26/doc/results.dox
index f4c6feefb5..798cc8f298 100644
--- a/deal.II/examples/step-26/doc/results.dox
+++ b/deal.II/examples/step-26/doc/results.dox
@@ -1 +1,115 @@
 <h1>Results</h1>
+
+As in many of the tutorials, the actual output of the program matters less
+than how we arrived there. Nonetheless, here it is:
+@code
+===========================================
+Number of active cells: 48
+Number of degrees of freedom: 65
+
+Time step 1 at t=0.002
+     7 CG iterations.
+
+===========================================
+Number of active cells: 60
+Number of degrees of freedom: 81
+
+
+Time step 1 at t=0.002
+     7 CG iterations.
+
+===========================================
+Number of active cells: 105
+Number of degrees of freedom: 136
+
+
+Time step 1 at t=0.002
+     7 CG iterations.
+
+[...]
+
+Time step 249 at t=0.498
+     13 CG iterations.
+Time step 250 at t=0.5
+     14 CG iterations.
+
+===========================================
+Number of active cells: 1803
+Number of degrees of freedom: 2109
+@endcode
+
+Maybe of more interest is a visualization of the solution and the mesh on which
+it was computed:
+
+<img
+src="http://www.dealii.org/images/steps/developer/step-26.movie.gif"
+alt="Animation of the solution of step 26.">
+
+The movie shows how the two sources switch on and off and how the mesh reacts
+to this. It is quite obvious that the mesh as is is probably not the best we
+could come up with. We'll get back to this in the next section.
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+There are at least two areas where one can improve this program significantly:
+adaptive time stepping and a better choice of the mesh.
+
+<h4>Adaptive time stepping</h4>
+
+Having chosen an implicit time stepping scheme, we are not bound by any
+CFL-like condition on the time step. Furthermore, because the time scales on
+which change happens on a given cell in the heat equation are not bound to the
+cells diameter (unlike the case with the wave equation, where we had a fixed
+speed of information transport that couples the temporal and spatial scales),
+we can choose the time step as we please. Or, better, choose it as we deem
+necessary for accuracy.
+
+Looking at the solution, it is clear that the action does not happen uniformly
+over time: a lot is changing around the time we switch on a source, things
+become less dramatic once a source is on for a little while, and we enter a
+long phase of decline when both sources are off. During these times, we could
+surely get away with a larger time step than before without sacrificing too
+much accuracy.
+
+The literature has many suggestions on how to choose the time step size
+adaptively. Much can be learned, for example, from the way ODE solvers choose
+their time steps. One can also be inspired by a posteriori error estimators
+that can, ideally, be written in a way that the consist of a temporal and a
+spatial contribution to the overall error. If the temporal one is too large,
+we should choose a smaller time step. Ideas in this direction can be found,
+for example, in the PhD thesis of a former principal developer of deal.II,
+Ralf Hartmann, published by the University of Heidelberg, Germany, in 2002.
+
+
+<h4>Better refinement criteria</h4>
+
+If you look at the meshes in the movie above, it is clear that they are not
+particularly well suited to the task at hand. In fact, they look rather
+random. 
+
+There are two factors at play. First, there are some islands where cells
+have been refined but that are surrounded by non-refined cells (and there
+are probably also a few occasional coarsened islands). These are not terrible,
+as they most of the time do not affect the approximation quality of the mesh,
+but they also don't help because so many of their additional degrees of
+freedom are in fact constrained by hanging node constraints. That said,
+this is easy to fix: the Triangulation class takes an argument to its
+constructor indicating a level of "mesh smoothing". Passing one of many 
+possible flags, this instructs the triangulation to refine some additional 
+cells, or not to refine some cells, so that the resulting mesh does not have 
+these artifacts.
+
+The second problem is more severe: the mesh appears to lag the solution. 
+The underlying reason is that we only adapt the mesh once every fifth
+time step, and only allow for a single refinement in these cases. Whenever a
+source switches on, the solution had been very smooth in this area before and
+the mesh was consequently rather coarse. This implies that the next time step
+when we refine the mesh, we will get one refinement level more in this area,
+and five time steps later another level, etc. But this is not enough: first,
+we should refine immediately when a source switches on (after all, in the
+current context we at least know what the right hand side is), and we should
+allow for more than one refinement level. Of course, all of this can be done
+using deal.II, it just requires a bit of algorithmic thinking in how to make
+this work!
diff --git a/deal.II/examples/step-26/doc/tooltip b/deal.II/examples/step-26/doc/tooltip
index 7b3e8b447a..9bd689e03e 100644
--- a/deal.II/examples/step-26/doc/tooltip
+++ b/deal.II/examples/step-26/doc/tooltip
@@ -1 +1 @@
-Complicated boundary descriptions.
+The heat equation. Time dependent meshes.
diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc
index d047fc7b5c..72adfa1b85 100644
--- a/deal.II/examples/step-26/step-26.cc
+++ b/deal.II/examples/step-26/step-26.cc
@@ -19,6 +19,8 @@
  */
 
 
+// The program starts with the usual include files, all of which you should
+// have seen before by now:
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/function.h>
@@ -51,12 +53,30 @@
 #include <iostream>
 
 
+// Then the usual placing of all content of this program into a namespace and
+// the importation of the deal.II namespace into the one we will work in:
 namespace Step26
 {
   using namespace dealii;
 
 
-
+  // @sect3{The <code>HeatEquation</code> class}
+  //
+  // The next piece is the declaration of the main class of this program. It
+  // follows the well trodden path of previous examples. If you have looked at
+  // step-6, for example, the only thing worth noting here is that we need to
+  // build two matrices (the mass and Laplace matrix) and keep the current and
+  // previous time step's solution. We then also need to store the current
+  // time, the size of the time step, and the number of the current time
+  // step. The last of the member variables denotes the theta parameter
+  // discussed in the introduction that allows us to treat the explicit and
+  // implicit Euler methods as well as the Crank-Nicolson method and other
+  // generalizations all in one program.
+  //
+  // As far as member functions are concerned, the only possible surprise is
+  // that the <code>refine_mesh</code> function takes arguments for the
+  // minimal and maximal mesh refinement level. The purpose of this is
+  // discussed in the introduction.
   template<int dim>
   class HeatEquation
   {
@@ -95,18 +115,26 @@ namespace Step26
 
 
 
+  // @sect3{Equation data}
+
+  // In the following classes and functions, we implement the various pieces
+  // of data that define this problem (right hand side and boundary values)
+  // that are used in this program and for which we need function objects. The
+  // right hand side is chosen as discussed at the end of the
+  // introduction. For boundary values, we choose zero values, but this is
+  // easily changed below.
   template<int dim>
-  class RightHandSide: public Function<dim>
+  class RightHandSide : public Function<dim>
   {
   public:
-    RightHandSide()
+    RightHandSide ()
       :
       Function<dim>(),
       period (0.2)
     {}
 
-    virtual double value(const Point<dim> &p,
-                         const unsigned int component = 0) const;
+    virtual double value (const Point<dim> &p,
+			  const unsigned int component = 0) const;
 
   private:
     const double period;
@@ -115,8 +143,8 @@ namespace Step26
 
 
   template<int dim>
-  double RightHandSide<dim>::value(const Point<dim> &p,
-                                   const unsigned int component) const
+  double RightHandSide<dim>::value (const Point<dim> &p,
+				    const unsigned int component) const
   {
     Assert (component == 0, ExcInternalError());
     Assert (dim == 2, ExcNotImplemented());
@@ -145,22 +173,18 @@ namespace Step26
 
 
   template<int dim>
-  class BoundaryValues: public Function<dim>
+  class BoundaryValues : public Function<dim>
   {
   public:
-    BoundaryValues()
-      :
-      Function<dim>()
-    {
-    }
-
-    virtual double value(const Point<dim> &p,
-                         const unsigned int component = 0) const;
+    virtual double value (const Point<dim>  &p,
+			  const unsigned int component = 0) const;
   };
 
+
+  
   template<int dim>
-  double BoundaryValues<dim>::value(const Point<dim> &/*p*/,
-                                    const unsigned int component) const
+  double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
+				     const unsigned int component) const
   {
     Assert(component == 0, ExcInternalError());
     return 0;
@@ -168,18 +192,45 @@ namespace Step26
 
 
 
+  // @sect3{The <code>HeatEquation</code> implementation}
+  //
+  // It is time now for the implementation of the main class. Let's
+  // start with the constructor which selects a linear element, a time
+  // step constant at 1/500 (remember that one period of the source
+  // on the right hand side was set to 0.2 above, so we resolve each
+  // period with 100 time steps) and chooses the Crank Nicolson method
+  // by setting $\theta=1/2$.
   template<int dim>
-  HeatEquation<dim>::HeatEquation()
+  HeatEquation<dim>::HeatEquation ()
     :
     fe(1),
     dof_handler(triangulation),
     time_step(1. / 500),
     theta(0.5)
-  {
-  }
+  {}
 
 
 
+  // @sect4{<code>HeatEquation::setup_system</code>}
+  //
+  // The next function is the one that sets up the DoFHandler object,
+  // computes the constraints, and sets the linear algebra objects
+  // to their correct sizes. We also compute the mass and Laplace
+  // matrix here by simply calling two functions in the library.
+  //
+  // Note that we compute these matrices taking into account already the
+  // constraints due to hanging nodes. These are all homogenous, i.e.,
+  // they only consist of constraints of the form $U_i = \alpha_{ij} U_j
+  // + \alpha_{ik} U_k$ (whereas inhomogenous constraints would also
+  // have a term not proportional to $U$, i.e., $U_i = \alpha_{ij} U_j
+  // + \alpha_{ik} U_k + c_i$). For this kind of constraint, we can
+  // eliminate hanging nodes independently in the matrix and the
+  // right hand side vectors, but this is not the case for inhomogenous
+  // constraints for which we can eliminate constrained degrees of freedom
+  // only by looking at both the system matrix and corresponding right
+  // right hand side at the same time. This may become a problem when
+  // dealing with non-zero Dirichlet boundary conditions, though we
+  // do not do this here in the current program.
   template<int dim>
   void HeatEquation<dim>::setup_system()
   {
@@ -227,7 +278,10 @@ namespace Step26
   }
 
 
-
+  // @sect4{<code>HeatEquation::solve_time_step</code>}
+  //
+  // The next function is the one that solves the actual linear system
+  // for a single time step. There is nothing surprising here:
   template<int dim>
   void HeatEquation<dim>::solve_time_step()
   {
@@ -237,7 +291,8 @@ namespace Step26
     PreconditionSSOR<> preconditioner;
     preconditioner.initialize(system_matrix, 1.0);
 
-    cg.solve(system_matrix, solution, system_rhs, preconditioner);
+    cg.solve(system_matrix, solution, system_rhs,
+             preconditioner);
 
     constraints.distribute(solution);
 
@@ -247,6 +302,9 @@ namespace Step26
 
 
 
+  // @sect4{<code>HeatEquation::output_results</code>}
+  //
+  // Neither is there anything new in generating graphical output:
   template<int dim>
   void HeatEquation<dim>::output_results() const
   {
@@ -265,40 +323,29 @@ namespace Step26
   }
 
 
-  // @sect4{BoussinesqFlowProblem::refine_mesh}
+  // @sect4{<code>HeatEquation::refine_mesh</code>}
   //
-  // This function takes care of the adaptive mesh refinement. The three tasks
+  // This function is the interesting part of the program. It takes care of
+  // the adaptive mesh refinement. The three tasks
   // this function performs is to first find out which cells to
   // refine/coarsen, then to actually do the refinement and eventually
   // transfer the solution vectors between the two different grids. The first
   // task is simply achieved by using the well-established Kelly error
-  // estimator on the temperature (it is the temperature we're mainly
-  // interested in for this program, and we need to be accurate in regions of
-  // high temperature gradients, also to not have too much numerical
-  // diffusion). The second task is to actually do the remeshing. That
-  // involves only basic functions as well, such as the
+  // estimator on the solution. The second task is to actually do the
+  // remeshing. That involves only basic functions as well, such as the
   // <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
-  // with the largest estimated error that together make up 80 per cent of the
+  // with the largest estimated error that together make up 60 per cent of the
   // error, and coarsens those cells with the smallest error that make up for
-  // a combined 10 per cent of the error.
+  // a combined 40 per cent of the error. Note that for problems such as the
+  // current one where the areas where something is going on are shifting
+  // around, we want to aggressively coarsen so that we can move cells
+  // around to where it is necessary.
   //
-  // If implemented like this, we would get a program that will not make much
-  // progress: Remember that we expect temperature fields that are nearly
-  // discontinuous (the diffusivity $\kappa$ is very small after all) and
-  // consequently we can expect that a freely adapted mesh will refine further
-  // and further into the areas of large gradients. This decrease in mesh size
-  // will then be accompanied by a decrease in time step, requiring an
-  // exceedingly large number of time steps to solve to a given final time. It
-  // will also lead to meshes that are much better at resolving
-  // discontinuities after several mesh refinement cycles than in the
-  // beginning.
-  //
-  // In particular to prevent the decrease in time step size and the
-  // correspondingly large number of time steps, we limit the maximal
-  // refinement depth of the mesh. To this end, after the refinement indicator
-  // has been applied to the cells, we simply loop over all cells on the
-  // finest level and unselect them from refinement if they would result in
-  // too high a mesh level.
+  // As already discussed in the introduction, too small a mesh leads to
+  // too small a time step, whereas too large a mesh leads to too little
+  // resolution. Consequently, after the first two steps, we have two
+  // loops that limit refinement and coarsening to an allowable range of
+  // cells:
   template <int dim>
   void HeatEquation<dim>::refine_mesh (const unsigned int min_grid_level,
                                        const unsigned int max_grid_level)
@@ -314,6 +361,7 @@ namespace Step26
     GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
                                                        estimated_error_per_cell,
                                                        0.6, 0.4);
+
     if (triangulation.n_levels() > max_grid_level)
       for (typename Triangulation<dim>::active_cell_iterator
            cell = triangulation.begin_active(max_grid_level);
@@ -330,58 +378,53 @@ namespace Step26
     // SolutionTransfer class and we have to prepare the solution vectors that
     // should be transferred to the new grid (we will lose the old grid once
     // we have done the refinement so the transfer has to happen concurrently
-    // with refinement). What we definitely need are the current and the old
-    // temperature (BDF-2 time stepping requires two old solutions). Since the
-    // SolutionTransfer objects only support to transfer one object per dof
-    // handler, we need to collect the two temperature solutions in one data
-    // structure. Moreover, we choose to transfer the Stokes solution, too,
-    // since we need the velocity at two previous time steps, of which only
-    // one is calculated on the fly.
+    // with refinement). At the point where we call this function, we will
+    // have just computed the solution, so we no longer need the old_solution
+    // variable (it will be overwritten by the solution just after the mesh
+    // may have been refined, i.e., at the end of the time step; see below).
+    // In other words, we only need the one solution vector, and we copy it
+    // to a temporary object where it is safe from being reset when we further
+    // down below call <code>setup_system()</code>.
     //
-    // Consequently, we initialize two SolutionTransfer objects for the Stokes
-    // and temperature DoFHandler objects, by attaching them to the old dof
-    // handlers. With this at place, we can prepare the triangulation and the
-    // data vectors for refinement (in this order).
-    std::vector<Vector<double> > x_solution (2);
-    x_solution[0] = solution;
-    x_solution[1] = old_solution;
-
+    // Consequently, we initialize a SolutionTransfer object by attaching
+    // it to the old DoF handler. We then prepare the triangulation and the
+    // data vector for refinement (in this order).
     SolutionTransfer<dim> solution_trans(dof_handler);
 
+    Vector<double> previous_solution;
+    previous_solution = solution;
     triangulation.prepare_coarsening_and_refinement();
-    solution_trans.prepare_for_coarsening_and_refinement(x_solution);
+    solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
 
     // Now everything is ready, so do the refinement and recreate the dof
     // structure on the new grid, and initialize the matrix structures and the
-    // new vectors in the <code>setup_dofs</code> function. Next, we actually
-    // perform the interpolation of the solutions between the grids. We create
-    // another copy of temporary vectors for temperature (now corresponding to
-    // the new grid), and let the interpolate function do the job. Then, the
-    // resulting array of vectors is written into the respective vector member
-    // variables. For the Stokes vector, everything is just the same &ndash;
-    // except that we do not need another temporary vector since we just
-    // interpolate a single vector. In the end, we have to tell the program
-    // that the matrices and preconditioners need to be regenerated, since the
-    // mesh has changed.
+    // new vectors in the <code>setup_system</code> function. Next, we actually
+    // perform the interpolation of the solution from old to new grid.
     triangulation.execute_coarsening_and_refinement ();
     setup_system ();
 
-    std::vector<Vector<double> > tmp (2);
-    tmp[0].reinit (solution);
-    tmp[1].reinit (solution);
-    solution_trans.interpolate(x_solution, tmp);
-
-    solution = tmp[0];
-    old_solution = tmp[1];
+    solution_trans.interpolate(previous_solution, solution);
   }
 
 
 
+  // @sect4{<code>HeatEquation::run</code>}
+  //
+  // This is the main driver of the program, where we loop over all
+  // time steps. At the top of the function, we set the number of
+  // initial global mesh refinements and the number of initial cycles of
+  // adaptive mesh refinement by repeating the first time step a few
+  // times. Then we create a mesh, initialize the various objects we will
+  // work with, set a label for where we should start when re-running
+  // the first time step, and interpolate the initial solution onto
+  // out mesh (we choose the zero function here, which of course we could
+  // do in a simpler way by just setting the solution vector to zero). We
+  // also output the initial time step once.
   template<int dim>
   void HeatEquation<dim>::run()
   {
-    const unsigned int initial_global_refinement = (dim == 2 ? 1 : 2);
-    const unsigned int n_adaptive_pre_refinement_steps = 1;
+    const unsigned int initial_global_refinement = 2;
+    const unsigned int n_adaptive_pre_refinement_steps = 4;
 
     GridGenerator::hyper_L (triangulation);
     triangulation.refine_global (initial_global_refinement);
@@ -395,6 +438,10 @@ namespace Step26
 
 start_time_iteration:
 
+    tmp.reinit (solution.size());
+    forcing_terms.reinit (solution.size());
+
+
     VectorTools::interpolate(dof_handler,
                              ZeroFunction<dim>(),
                              old_solution);
@@ -405,6 +452,12 @@ start_time_iteration:
 
     output_results();
 
+    // Then we start the main loop until the computed time exceeds our
+    // end time of 0.5. The first task is to build the right hand
+    // side of the linear system we need to solve in each time step.
+    // Recall that it contains the term $MU^{n-1}-(1-\theta)k_n AU^{n-1}$.
+    // We put these terms into the variable system_rhs, with the
+    // help of a temporary vector:
     while (time <= 0.5)
       {
         time += time_step;
@@ -413,14 +466,18 @@ start_time_iteration:
         std::cout << "Time step " << timestep_number << " at t=" << time
                   << std::endl;
 
-        tmp.reinit (solution.size());
-        forcing_terms.reinit (solution.size());
-
         mass_matrix.vmult(system_rhs, old_solution);
 
         laplace_matrix.vmult(tmp, old_solution);
         system_rhs.add(-(1 - theta) * time_step, tmp);
 
+        // The second piece is to compute the contributions of the source
+        // terms. This corresponds to the term $k_n
+        // \left[ (1-\theta)F^{n-1} + \theta F^n \right]$. The following
+        // code calls VectorTools::create_right_hand_side to compute the
+        // vectors $F$, where we set the time of the right hand side
+        // (source) function before we evaluate it. The result of this
+        // all ends up in the forcing_terms variable:
         RightHandSide<dim> rhs_function;
         rhs_function.set_time(time);
         VectorTools::create_right_hand_side(dof_handler,
@@ -438,8 +495,25 @@ start_time_iteration:
 
         forcing_terms.add(time_step * (1 - theta), tmp);
 
+        // Next, we add the forcing terms to the ones that
+        // come from the time stepping, and also build the matrix
+        // $M+k_n\theta A$ that we have to invert in each time step.
+        // The final piece of these operations is to eliminate
+        // hanging node constrained degrees of freedom from the
+        // linear system:
         system_rhs += forcing_terms;
 
+        system_matrix.copy_from(mass_matrix);
+        system_matrix.add(theta * time_step, laplace_matrix);
+
+        constraints.condense (system_matrix, system_rhs);
+
+        // There is one more operation we need to do before we
+        // can solve it: boundary values. To this end, we create
+        // a boundary value object, set the proper time to the one
+        // of the current time step, and evaluate it as we have
+        // done many times before. The result is used to also
+        // set the correct boundary values in the linear system:
         {
           BoundaryValues<dim> boundary_values_function;
           boundary_values_function.set_time(time);
@@ -450,20 +524,27 @@ start_time_iteration:
                                                    boundary_values_function,
                                                    boundary_values);
 
-          system_matrix.copy_from(mass_matrix);
-          system_matrix.add(theta * time_step, laplace_matrix);
           MatrixTools::apply_boundary_values(boundary_values,
                                              system_matrix,
                                              solution,
                                              system_rhs);
         }
 
-        constraints.condense (system_rhs);
-
+        // With this out of the way, all we have to do is solve the
+        // system, generate graphical data, and...
         solve_time_step();
 
         output_results();
 
+        // ...take care of mesh refinement. Here, what we want to do is
+        // (i) refine the requested number of times at the very beginning
+        // of the solution procedure, after which we jump to the top to
+        // restart the time iteration, (ii) refine every fifth time
+        // step after that.
+        //
+        // The time loop and, indeed, the main part of the program ends
+        // with starting into the next time step by setting old_solution
+        // to the solution we have just computed.
         if ((timestep_number == 1) &&
             (pre_refinement_step < n_adaptive_pre_refinement_steps))
           {
@@ -471,19 +552,32 @@ start_time_iteration:
                          initial_global_refinement + n_adaptive_pre_refinement_steps);
             ++pre_refinement_step;
 
+            tmp.reinit (solution.size());
+            forcing_terms.reinit (solution.size());
+
             std::cout << std::endl;
 
             goto start_time_iteration;
           }
         else if ((timestep_number > 0) && (timestep_number % 5 == 0))
-          refine_mesh (initial_global_refinement,
-                       initial_global_refinement + n_adaptive_pre_refinement_steps);
+          {
+            refine_mesh (initial_global_refinement,
+                         initial_global_refinement + n_adaptive_pre_refinement_steps);
+            tmp.reinit (solution.size());
+            forcing_terms.reinit (solution.size());
+          }
 
         old_solution = solution;
       }
   }
 }
 
+
+// @sect3{The <code>main</code> function}
+//
+// Having made it this far,  there is, again, nothing
+// much to discuss for the main function of this
+// program: it looks like all such functions since step-6.
 int main()
 {
   try
-- 
2.39.5