step-23.h

Go to the documentation of this file.

,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         Assert(component == 0, ExcIndexRange(component, 0, 1));
 *         return 0;
 *       }
 *     };
 *   
 *   
 *   
 *     template <int dim>
 *     class InitialValuesV : public Function<dim>
 *     {
 *     public:
 *       virtual double value(const Point<dim> & /*p*/,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         Assert(component == 0, ExcIndexRange(component, 0, 1));
 *         return 0;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * Secondly, we have the right hand side forcing term. Boring as we are, we
 * choose zero here as well:
 * 
 * @code
 *     template <int dim>
 *     class RightHandSide : public Function<dim>
 *     {
 *     public:
 *       virtual double value(const Point<dim> & /*p*/,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         Assert(component == 0, ExcIndexRange(component, 0, 1));
 *         return 0;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * Finally, we have boundary values for @f$u@f$ and @f$v@f$. They are as described
 * in the introduction, one being the time derivative of the other:
 * 
 * @code
 *     template <int dim>
 *     class BoundaryValuesU : public Function<dim>
 *     {
 *     public:
 *       virtual double value(const Point<dim> & p,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         Assert(component == 0, ExcIndexRange(component, 0, 1));
 *   
 *         if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
 *             (p[1] > -1. / 3))
 *           return std::sin(this->get_time() * 4 * numbers::PI);
 *         else
 *           return 0;
 *       }
 *     };
 *   
 *   
 *   
 *     template <int dim>
 *     class BoundaryValuesV : public Function<dim>
 *     {
 *     public:
 *       virtual double value(const Point<dim> & p,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         Assert(component == 0, ExcIndexRange(component, 0, 1));
 *   
 *         if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
 *             (p[1] > -1. / 3))
 *           return (std::cos(this->get_time() * 4 * numbers::PI) * 4 * numbers::PI);
 *         else
 *           return 0;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ImplementationofthecodeWaveEquationcodeclass"></a> 
 * <h3>Implementation of the <code>WaveEquation</code> class</h3>
 * 
 
 * 
 * The implementation of the actual logic is actually fairly short, since we
 * relegate things like assembling the matrices and right hand side vectors
 * to the library. The rest boils down to not much more than 130 lines of
 * actual code, a significant fraction of which is boilerplate code that can
 * be taken from previous example programs (e.g. the functions that solve
 * linear systems, or that generate output).
 *   
 
 * 
 * Let's start with the constructor (for an explanation of the choice of
 * time step, see the section on Courant, Friedrichs, and Lewy in the
 * introduction):
 * 
 * @code
 *     template <int dim>
 *     WaveEquation<dim>::WaveEquation()
 *       : fe(1)
 *       , dof_handler(triangulation)
 *       , time_step(1. / 64)
 *       , time(time_step)
 *       , timestep_number(1)
 *       , theta(0.5)
 *     {}
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="WaveEquationsetup_system"></a> 
 * <h4>WaveEquation::setup_system</h4>
 * 
 
 * 
 * The next function is the one that sets up the mesh, DoFHandler, and
 * matrices and vectors at the beginning of the program, i.e. before the
 * first time step. The first few lines are pretty much standard if you've
 * read through the tutorial programs at least up to @ref step_6 "step-6":
 * 
 * @code
 *     template <int dim>
 *     void WaveEquation<dim>::setup_system()
 *     {
 *       GridGenerator::hyper_cube(triangulation, -1, 1);
 *       triangulation.refine_global(7);
 *   
 *       std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl;
 *   
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl
 *                 << std::endl;
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *       sparsity_pattern.copy_from(dsp);
 *   
 * @endcode
 * 
 * Then comes a block where we have to initialize the 3 matrices we need
 * in the course of the program: the mass matrix, the Laplace matrix, and
 * the matrix @f$M+k^2\theta^2A@f$ used when solving for @f$U^n@f$ in each time
 * step.
 *     
 
 * 
 * When setting up these matrices, note that they all make use of the same
 * sparsity pattern object. Finally, the reason why matrices and sparsity
 * patterns are separate objects in deal.II (unlike in many other finite
 * element or linear algebra classes) becomes clear: in a significant
 * fraction of applications, one has to hold several matrices that happen
 * to have the same sparsity pattern, and there is no reason for them not
 * to share this information, rather than re-building and wasting memory
 * on it several times.
 *     
 
 * 
 * After initializing all of these matrices, we call library functions
 * that build the Laplace and mass matrices. All they need is a DoFHandler
 * object and a quadrature formula object that is to be used for numerical
 * integration. Note that in many respects these functions are better than
 * what we would usually do in application programs, for example because
 * they automatically parallelize building the matrices if multiple
 * processors are available in a machine: for more information see the
 * documentation of WorkStream or the
 * @ref threads "Parallel computing with multiple processors"
 * module. The matrices for solving linear systems will be filled in the
 * run() method because we need to re-apply boundary conditions every time
 * step.
 * 
 * @code
 *       mass_matrix.reinit(sparsity_pattern);
 *       laplace_matrix.reinit(sparsity_pattern);
 *       matrix_u.reinit(sparsity_pattern);
 *       matrix_v.reinit(sparsity_pattern);
 *   
 *       MatrixCreator::create_mass_matrix(dof_handler,
 *                                         QGauss<dim>(fe.degree + 1),
 *                                         mass_matrix);
 *       MatrixCreator::create_laplace_matrix(dof_handler,
 *                                            QGauss<dim>(fe.degree + 1),
 *                                            laplace_matrix);
 *   
 * @endcode
 * 
 * The rest of the function is spent on setting vector sizes to the
 * correct value. The final line closes the hanging node constraints
 * object. Since we work on a uniformly refined mesh, no constraints exist
 * or have been computed (i.e. there was no need to call
 * DoFTools::make_hanging_node_constraints as in other programs), but we
 * need a constraints object in one place further down below anyway.
 * 
 * @code
 *       solution_u.reinit(dof_handler.n_dofs());
 *       solution_v.reinit(dof_handler.n_dofs());
 *       old_solution_u.reinit(dof_handler.n_dofs());
 *       old_solution_v.reinit(dof_handler.n_dofs());
 *       system_rhs.reinit(dof_handler.n_dofs());
 *   
 *       constraints.close();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="WaveEquationsolve_uandWaveEquationsolve_v"></a> 
 * <h4>WaveEquation::solve_u and WaveEquation::solve_v</h4>
 * 
 
 * 
 * The next two functions deal with solving the linear systems associated
 * with the equations for @f$U^n@f$ and @f$V^n@f$. Both are not particularly
 * interesting as they pretty much follow the scheme used in all the
 * previous tutorial programs.
 *   
 
 * 
 * One can make little experiments with preconditioners for the two matrices
 * we have to invert. As it turns out, however, for the matrices at hand
 * here, using Jacobi or SSOR preconditioners reduces the number of
 * iterations necessary to solve the linear system slightly, but due to the
 * cost of applying the preconditioner it is no win in terms of run-time. It
 * is not much of a loss either, but let's keep it simple and just do
 * without:
 * 
 * @code
 *     template <int dim>
 *     void WaveEquation<dim>::solve_u()
 *     {
 *       SolverControl            solver_control(1000, 1e-8 * system_rhs.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       cg.solve(matrix_u, solution_u, system_rhs, PreconditionIdentity());
 *   
 *       std::cout << "   u-equation: " << solver_control.last_step()
 *                 << " CG iterations." << std::endl;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void WaveEquation<dim>::solve_v()
 *     {
 *       SolverControl            solver_control(1000, 1e-8 * system_rhs.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       cg.solve(matrix_v, solution_v, system_rhs, PreconditionIdentity());
 *   
 *       std::cout << "   v-equation: " << solver_control.last_step()
 *                 << " CG iterations." << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="WaveEquationoutput_results"></a> 
 * <h4>WaveEquation::output_results</h4>
 * 
 
 * 
 * Likewise, the following function is pretty much what we've done
 * before. The only thing worth mentioning is how here we generate a string
 * representation of the time step number padded with leading zeros to 3
 * character length using the Utilities::int_to_string function's second
 * argument.
 * 
 * @code
 *     template <int dim>
 *     void WaveEquation<dim>::output_results() const
 *     {
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution_u, "U");
 *       data_out.add_data_vector(solution_v, "V");
 *   
 *       data_out.build_patches();
 *   
 *       const std::string filename =
 *         "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
 * @endcode
 * 
 * Like @ref step_15 "step-15", since we write output at every time step (and the system
 * we have to solve is relatively easy), we instruct DataOut to use the
 * zlib compression algorithm that is optimized for speed instead of disk
 * usage since otherwise plotting the output becomes a bottleneck:
 * 
 * @code
 *       DataOutBase::VtkFlags vtk_flags;
 *       vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
 *       data_out.set_flags(vtk_flags);
 *       std::ofstream output(filename);
 *       data_out.write_vtu(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="WaveEquationrun"></a> 
 * <h4>WaveEquation::run</h4>
 * 
 
 * 
 * The following is really the only interesting function of the program. It
 * contains the loop over all time steps, but before we get to that we have
 * to set up the grid, DoFHandler, and matrices. In addition, we have to
 * somehow get started with initial values. To this end, we use the
 * VectorTools::project function that takes an object that describes a
 * continuous function and computes the @f$L^2@f$ projection of this function
 * onto the finite element space described by the DoFHandler object. Can't
 * be any simpler than that:
 * 
 * @code
 *     template <int dim>
 *     void WaveEquation<dim>::run()
 *     {
 *       setup_system();
 *   
 *       VectorTools::project(dof_handler,
 *                            constraints,
 *                            QGauss<dim>(fe.degree + 1),
 *                            InitialValuesU<dim>(),
 *                            old_solution_u);
 *       VectorTools::project(dof_handler,
 *                            constraints,
 *                            QGauss<dim>(fe.degree + 1),
 *                            InitialValuesV<dim>(),
 *                            old_solution_v);
 *   
 * @endcode
 * 
 * The next thing is to loop over all the time steps until we reach the
 * end time (@f$T=5@f$ in this case). In each time step, we first have to
 * solve for @f$U^n@f$, using the equation @f$(M^n + k^2\theta^2 A^n)U^n =@f$
 * @f$(M^{n,n-1} - k^2\theta(1-\theta) A^{n,n-1})U^{n-1} + kM^{n,n-1}V^{n-1}
 * +@f$ @f$k\theta \left[k \theta F^n + k(1-\theta) F^{n-1} \right]@f$. Note
 * that we use the same mesh for all time steps, so that @f$M^n=M^{n,n-1}=M@f$
 * and @f$A^n=A^{n,n-1}=A@f$. What we therefore have to do first is to add up
 * @f$MU^{n-1} - k^2\theta(1-\theta) AU^{n-1} + kMV^{n-1}@f$ and the forcing
 * terms, and put the result into the <code>system_rhs</code> vector. (For
 * these additions, we need a temporary vector that we declare before the
 * loop to avoid repeated memory allocations in each time step.)
 *     
 
 * 
 * The one thing to realize here is how we communicate the time variable
 * to the object describing the right hand side: each object derived from
 * the Function class has a time field that can be set using the
 * Function::set_time and read by Function::get_time. In essence, using
 * this mechanism, all functions of space and time are therefore
 * considered functions of space evaluated at a particular time. This
 * matches well what we typically need in finite element programs, where
 * we almost always work on a single time step at a time, and where it
 * never happens that, for example, one would like to evaluate a
 * space-time function for all times at any given spatial location.
 * 
 * @code
 *       Vector<double> tmp(solution_u.size());
 *       Vector<double> forcing_terms(solution_u.size());
 *   
 *       for (; time <= 5; time += time_step, ++timestep_number)
 *         {
 *           std::cout << "Time step " << timestep_number << " at t=" << time
 *                     << std::endl;
 *   
 *           mass_matrix.vmult(system_rhs, old_solution_u);
 *   
 *           mass_matrix.vmult(tmp, old_solution_v);
 *           system_rhs.add(time_step, tmp);
 *   
 *           laplace_matrix.vmult(tmp, old_solution_u);
 *           system_rhs.add(-theta * (1 - theta) * time_step * time_step, tmp);
 *   
 *           RightHandSide<dim> rhs_function;
 *           rhs_function.set_time(time);
 *           VectorTools::create_right_hand_side(dof_handler,
 *                                               QGauss<dim>(fe.degree + 1),
 *                                               rhs_function,
 *                                               tmp);
 *           forcing_terms = tmp;
 *           forcing_terms *= theta * time_step;
 *   
 *           rhs_function.set_time(time - time_step);
 *           VectorTools::create_right_hand_side(dof_handler,
 *                                               QGauss<dim>(fe.degree + 1),
 *                                               rhs_function,
 *                                               tmp);
 *   
 *           forcing_terms.add((1 - theta) * time_step, tmp);
 *   
 *           system_rhs.add(theta * time_step, forcing_terms);
 *   
 * @endcode
 * 
 * After so constructing the right hand side vector of the first
 * equation, all we have to do is apply the correct boundary
 * values. As for the right hand side, this is a space-time function
 * evaluated at a particular time, which we interpolate at boundary
 * nodes and then use the result to apply boundary values as we
 * usually do. The result is then handed off to the solve_u()
 * function:
 * 
 * @code
 *           {
 *             BoundaryValuesU<dim> boundary_values_u_function;
 *             boundary_values_u_function.set_time(time);
 *   
 *             std::map<types::global_dof_index, double> boundary_values;
 *             VectorTools::interpolate_boundary_values(dof_handler,
 *                                                      0,
 *                                                      boundary_values_u_function,
 *                                                      boundary_values);
 *   
 * @endcode
 * 
 * The matrix for solve_u() is the same in every time steps, so one
 * could think that it is enough to do this only once at the
 * beginning of the simulation. However, since we need to apply
 * boundary values to the linear system (which eliminate some matrix
 * rows and columns and give contributions to the right hand side),
 * we have to refill the matrix in every time steps before we
 * actually apply boundary data. The actual content is very simple:
 * it is the sum of the mass matrix and a weighted Laplace matrix:
 * 
 * @code
 *             matrix_u.copy_from(mass_matrix);
 *             matrix_u.add(theta * theta * time_step * time_step, laplace_matrix);
 *             MatrixTools::apply_boundary_values(boundary_values,
 *                                                matrix_u,
 *                                                solution_u,
 *                                                system_rhs);
 *           }
 *           solve_u();
 *   
 *   
 * @endcode
 * 
 * The second step, i.e. solving for @f$V^n@f$, works similarly, except
 * that this time the matrix on the left is the mass matrix (which we
 * copy again in order to be able to apply boundary conditions, and
 * the right hand side is @f$MV^{n-1} - k\left[ \theta A U^n +
 * (1-\theta) AU^{n-1}\right]@f$ plus forcing terms. Boundary values
 * are applied in the same way as before, except that now we have to
 * use the BoundaryValuesV class:
 * 
 * @code
 *           laplace_matrix.vmult(system_rhs, solution_u);
 *           system_rhs *= -theta * time_step;
 *   
 *           mass_matrix.vmult(tmp, old_solution_v);
 *           system_rhs += tmp;
 *   
 *           laplace_matrix.vmult(tmp, old_solution_u);
 *           system_rhs.add(-time_step * (1 - theta), tmp);
 *   
 *           system_rhs += forcing_terms;
 *   
 *           {
 *             BoundaryValuesV<dim> boundary_values_v_function;
 *             boundary_values_v_function.set_time(time);
 *   
 *             std::map<types::global_dof_index, double> boundary_values;
 *             VectorTools::interpolate_boundary_values(dof_handler,
 *                                                      0,
 *                                                      boundary_values_v_function,
 *                                                      boundary_values);
 *             matrix_v.copy_from(mass_matrix);
 *             MatrixTools::apply_boundary_values(boundary_values,
 *                                                matrix_v,
 *                                                solution_v,
 *                                                system_rhs);
 *           }
 *           solve_v();
 *   
 * @endcode
 * 
 * Finally, after both solution components have been computed, we
 * output the result, compute the energy in the solution, and go on to
 * the next time step after shifting the present solution into the
 * vectors that hold the solution at the previous time step. Note the
 * function SparseMatrix::matrix_norm_square that can compute
 * @f$\left<V^n,MV^n\right>@f$ and @f$\left<U^n,AU^n\right>@f$ in one step,
 * saving us the expense of a temporary vector and several lines of
 * code:
 * 
 * @code
 *           output_results();
 *   
 *           std::cout << "   Total energy: "
 *                     << (mass_matrix.matrix_norm_square(solution_v) +
 *                         laplace_matrix.matrix_norm_square(solution_u)) /
 *                          2
 *                     << std::endl;
 *   
 *           old_solution_u = solution_u;
 *           old_solution_v = solution_v;
 *         }
 *     }
 *   } // namespace Step23
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * What remains is the main function of the program. There is nothing here
 * that hasn't been shown in several of the previous programs:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step23;
 *   
 *         WaveEquation<2> wave_equation_solver;
 *         wave_equation_solver.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
When the program is run, it produces the following output:
@code
Number of active cells: 16384
Number of degrees of freedom: 16641
 
Time step 1 at t=0.015625
   u-equation: 8 CG iterations.
   v-equation: 22 CG iterations.
   Total energy: 1.17887
Time step 2 at t=0.03125
   u-equation: 8 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 2.9655
Time step 3 at t=0.046875
   u-equation: 8 CG iterations.
   v-equation: 21 CG iterations.
   Total energy: 4.33761
Time step 4 at t=0.0625
   u-equation: 7 CG iterations.
   v-equation: 21 CG iterations.
   Total energy: 5.35499
Time step 5 at t=0.078125
   u-equation: 7 CG iterations.
   v-equation: 21 CG iterations.
   Total energy: 6.18652
Time step 6 at t=0.09375
   u-equation: 7 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 6.6799
 
...
 
Time step 31 at t=0.484375
   u-equation: 7 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 21.9068
Time step 32 at t=0.5
   u-equation: 7 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 23.3394
Time step 33 at t=0.515625
   u-equation: 7 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 23.1019
 
...
 
Time step 319 at t=4.98438
   u-equation: 7 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 23.1019
Time step 320 at t=5
   u-equation: 7 CG iterations.
   v-equation: 20 CG iterations.
   Total energy: 23.1019
@endcode
 
What we see immediately is that the energy is a constant at least after
@f$t=\frac 12@f$ (until which the boundary source term @f$g@f$ is nonzero, injecting
energy into the system).
 
In addition to the screen output, the program writes the solution of each time
step to an output file. If we process them adequately and paste them into a
movie, we get the following:
 
<img src="https://www.dealii.org/images/steps/developer/step-23.movie.gif" alt="Animation of the solution of step 23.">
 
The movie shows the generated wave nice traveling through the domain and back,
being reflected at the clamped boundary. Some numerical noise is trailing the
wave, an artifact of a too-large mesh size that can be reduced by reducing the
mesh width and the time step.
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
If you want to explore a bit, try out some of the following things:
<ul>
  <li>Varying @f$\theta@f$. This gives different time stepping schemes, some of
  which are stable while others are not. Take a look at how the energy
  evolves.
 
  <li>Different initial and boundary conditions, right hand sides.
 
  <li>More complicated domains or more refined meshes. Remember that the time
  step needs to be bounded by the mesh width, so changing the mesh should
  always involve also changing the time step. We will come back to this issue
  in @ref step_24 "step-24".
 
  <li>Variable coefficients: In real media, the wave speed is often
  variable. In particular, the "real" wave equation in realistic media would
  read
  @f[
     \rho(x) \frac{\partial^2 u}{\partial t^2}
     -
     \nabla \cdot
     a(x) \nabla u = f,
  @f]
  where @f$\rho(x)@f$ is the density of the material, and @f$a(x)@f$ is related to the
  stiffness coefficient. The wave speed is then @f$c=\sqrt{a/\rho}@f$.
 
  To make such a change, we would have to compute the mass and Laplace
  matrices with a variable coefficient. Fortunately, this isn't too hard: the
  functions MatrixCreator::create_laplace_matrix and
  MatrixCreator::create_mass_matrix have additional default parameters that can
  be used to pass non-constant coefficient functions to them. The required
  changes are therefore relatively small. On the other hand, care must be
  taken again to make sure the time step is within the allowed range.
 
  <li>In the in-code comments, we discussed the fact that the matrices for
  solving for @f$U^n@f$ and @f$V^n@f$ need to be reset in every time because of
  boundary conditions, even though the actual content does not change. It is
  possible to avoid copying by not eliminating columns in the linear systems,
  which is implemented by appending a @p false argument to the call:
  @code
    MatrixTools::apply_boundary_values(boundary_values,
                                       matrix_u,
                                       solution_u,
                                       system_rhs,
                                       false);
  @endcode
 
  <li>deal.II being a library that supports adaptive meshes it would of course be
  nice if this program supported change the mesh every few time steps. Given the
  structure of the solution &mdash; a wave that travels through the domain &mdash;
  it would seem appropriate if we only refined the mesh where the wave currently is,
  and not simply everywhere. It is intuitively clear that we should be able to
  save a significant amount of cells this way. (Though upon further thought one
  realizes that this is really only the case in the initial stages of the simulation.
  After some time, for wave phenomena, the domain is filled with reflections of
  the initial wave going in every direction and filling every corner of the domain.
  At this point, there is in general little one can gain using local mesh
  refinement.)
 
  To make adaptively changing meshes possible, there are basically two routes.
  The "correct" way would be to go back to the weak form we get using Rothe's
  method. For example, the first of the two equations to be solved in each time
  step looked like this:
  \f{eqnarray*}
  (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
  (u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla u^{n-1},\nabla \varphi)
  +
  k(v^{n-1},\varphi)
  + k^2\theta
  \left[
  \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
  \right].
  \f}
  Now, note that we solve for @f$u^n@f$ on mesh @f${\mathbb T}^n@f$, and
  consequently the test functions @f$\varphi@f$ have to be from the space
  @f$V_h^n@f$ as well. As discussed in the introduction, terms like
  @f$(u^{n-1},\varphi)@f$ then require us to integrate the solution of the
  previous step (which may have been computed on a different mesh
  @f${\mathbb T}^{n-1}@f$) against the test functions of the current mesh,
  leading to a matrix @f$M^{n,n-1}@f$. This process of integrating shape
  functions from different meshes is, at best, awkward. It can be done
  but because it is difficult to ensure that @f${\mathbb T}^{n-1}@f$ and
  @f${\mathbb T}^{n}@f$ differ by at most one level of refinement, one
  has to recursively match cells from both meshes. It is feasible to
  do this, but it leads to lengthy and not entirely obvious code.
 
  The second approach is the following: whenever we change the mesh,
  we simply interpolate the solution from the last time step on the old
  mesh to the new mesh, using the SolutionTransfer class. In other words,
  instead of the equation above, we would solve
  \f{eqnarray*}
  (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
  (I^n u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla I^n u^{n-1},\nabla \varphi)
  +
  k(I^n v^{n-1},\varphi)
  + k^2\theta
  \left[
  \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
  \right],
  \f}
  where @f$I^n@f$ interpolates a given function onto mesh @f${\mathbb T}^n@f$.
  This is a much simpler approach because, in each time step, we no
  longer have to worry whether @f$u^{n-1},v^{n-1}@f$ were computed on the
  same mesh as we are using now or on a different mesh. Consequently,
  the only changes to the code necessary are the addition of a function
  that computes the error, marks cells for refinement, sets up a
  SolutionTransfer object, transfers the solution to the new mesh, and
  rebuilds matrices and right hand side vectors on the new mesh. Neither
  the functions building the matrices and right hand sides, nor the
  solvers need to be changed.
 
  While this second approach is, strictly speaking,
  not quite correct in the Rothe framework (it introduces an addition source
  of error, namely the interpolation), it is nevertheless what
  almost everyone solving time dependent equations does. We will use this
  method in @ref step_31 "step-31", for example.
</ul>
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-23.cc"
*/