step-9.h

Go to the documentation of this file.

false);
 *     sparsity_pattern.copy_from(dsp);
 * 
 *     system_matrix.reinit(sparsity_pattern);
 * 
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * In the following function, the matrix and right hand side are
 * assembled. As stated in the documentation of the main class above, it
 * does not do this itself, but rather delegates to the function following
 * next, utilizing the WorkStream concept discussed in @ref threads .
 *   
 
 * 
 * If you have looked through the @ref threads module, you will have
 * seen that assembling in parallel does not take an incredible
 * amount of extra code as long as you diligently describe what the
 * scratch and copy data objects are, and if you define suitable
 * functions for the local assembly and the copy operation from local
 * contributions to global objects. This done, the following will do
 * all the heavy lifting to get these operations done on multiple
 * threads on as many cores as you have in your system:
 * 
 * @code
 *   template <int dim>
 *   void AdvectionProblem<dim>::assemble_system()
 *   {
 *     WorkStream::run(dof_handler.begin_active(),
 *                     dof_handler.end(),
 *                     *this,
 *                     &AdvectionProblem::local_assemble_system,
 *                     &AdvectionProblem::copy_local_to_global,
 *                     AssemblyScratchData(fe),
 *                     AssemblyCopyData());
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * As already mentioned above, we need to have scratch objects for
 * the parallel computation of local contributions. These objects
 * contain FEValues and FEFaceValues objects (as well as some arrays), and so
 * we will need to have constructors and copy constructors that allow us to
 * create them. For the cell terms we need the values
 * and gradients of the shape functions, the quadrature points in
 * order to determine the source density and the advection field at
 * a given point, and the weights of the quadrature points times the
 * determinant of the Jacobian at these points. In contrast, for the
 * boundary integrals, we don't need the gradients, but rather the
 * normal vectors to the cells. This determines which update flags
 * we will have to pass to the constructors of the members of the
 * class:
 * 
 * @code
 *   template <int dim>
 *   AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
 *     const FiniteElement<dim> &fe)
 *     : fe_values(fe,
 *                 QGauss<dim>(fe.degree + 1),
 *                 update_values | update_gradients | update_quadrature_points |
 *                   update_JxW_values)
 *     , fe_face_values(fe,
 *                      QGauss<dim - 1>(fe.degree + 1),
 *                      update_values | update_quadrature_points |
 *                        update_JxW_values | update_normal_vectors)
 *     , rhs_values(fe_values.get_quadrature().size())
 *     , advection_directions(fe_values.get_quadrature().size())
 *     , face_boundary_values(fe_face_values.get_quadrature().size())
 *     , face_advection_directions(fe_face_values.get_quadrature().size())
 *   {}
 * 
 * 
 * 
 *   template <int dim>
 *   AdvectionProblem<dim>::AssemblyScratchData::AssemblyScratchData(
 *     const AssemblyScratchData &scratch_data)
 *     : fe_values(scratch_data.fe_values.get_fe(),
 *                 scratch_data.fe_values.get_quadrature(),
 *                 update_values | update_gradients | update_quadrature_points |
 *                   update_JxW_values)
 *     , fe_face_values(scratch_data.fe_face_values.get_fe(),
 *                      scratch_data.fe_face_values.get_quadrature(),
 *                      update_values | update_quadrature_points |
 *                        update_JxW_values | update_normal_vectors)
 *     , rhs_values(scratch_data.rhs_values.size())
 *     , advection_directions(scratch_data.advection_directions.size())
 *     , face_boundary_values(scratch_data.face_boundary_values.size())
 *     , face_advection_directions(scratch_data.face_advection_directions.size())
 *   {}
 * 
 * 
 * 
 * @endcode
 * 
 * Now, this is the function that does the actual work. It is not very
 * different from the <code>assemble_system</code> functions of previous
 * example programs, so we will again only comment on the differences. The
 * mathematical stuff closely follows what we have said in the introduction.
 *   
 
 * 
 * There are a number of points worth mentioning here, though. The
 * first one is that we have moved the FEValues and FEFaceValues
 * objects into the ScratchData object. We have done so because the
 * alternative would have been to simply create one every time we
 * get into this function -- i.e., on every cell. It now turns out
 * that the FEValues classes were written with the explicit goal of
 * moving everything that remains the same from cell to cell into
 * the construction of the object, and only do as little work as
 * possible in FEValues::reinit() whenever we move to a new
 * cell. What this means is that it would be very expensive to
 * create a new object of this kind in this function as we would
 * have to do it for every cell -- exactly the thing we wanted to
 * avoid with the FEValues class. Instead, what we do is create it
 * only once (or a small number of times) in the scratch objects and
 * then re-use it as often as we can.
 *   
 
 * 
 * This begs the question of whether there are other objects we
 * create in this function whose creation is expensive compared to
 * its use. Indeed, at the top of the function, we declare all sorts
 * of objects. The <code>AdvectionField</code>,
 * <code>RightHandSide</code> and <code>BoundaryValues</code> do not
 * cost much to create, so there is no harm here. However,
 * allocating memory in creating the <code>rhs_values</code> and
 * similar variables below typically costs a significant amount of
 * time, compared to just accessing the (temporary) values we store
 * in them. Consequently, these would be candidates for moving into
 * the <code>AssemblyScratchData</code> class. We will leave this as
 * an exercise.
 * 
 * @code
 *   template <int dim>
 *   void AdvectionProblem<dim>::local_assemble_system(
 *     const typename DoFHandler<dim>::active_cell_iterator &cell,
 *     AssemblyScratchData &                                 scratch_data,
 *     AssemblyCopyData &                                    copy_data)
 *   {
 * @endcode
 * 
 * We define some abbreviations to avoid unnecessarily long lines:
 * 
 * @code
 *     const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *     const unsigned int n_q_points =
 *       scratch_data.fe_values.get_quadrature().size();
 *     const unsigned int n_face_q_points =
 *       scratch_data.fe_face_values.get_quadrature().size();
 * 
 * @endcode
 * 
 * We declare cell matrix and cell right hand side...
 * 
 * @code
 *     copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
 *     copy_data.cell_rhs.reinit(dofs_per_cell);
 * 
 * @endcode
 * 
 * ... an array to hold the global indices of the degrees of freedom of
 * the cell on which we are presently working...
 * 
 * @code
 *     copy_data.local_dof_indices.resize(dofs_per_cell);
 * 
 * @endcode
 * 
 * ... then initialize the <code>FEValues</code> object...
 * 
 * @code
 *     scratch_data.fe_values.reinit(cell);
 * 
 * @endcode
 * 
 * ... obtain the values of right hand side and advection directions
 * at the quadrature points...
 * 
 * @code
 *     scratch_data.advection_field.value_list(
 *       scratch_data.fe_values.get_quadrature_points(),
 *       scratch_data.advection_directions);
 *     scratch_data.right_hand_side.value_list(
 *       scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
 * 
 * @endcode
 * 
 * ... set the value of the streamline diffusion parameter as
 * described in the introduction...
 * 
 * @code
 *     const double delta = 0.1 * cell->diameter();
 * 
 * @endcode
 * 
 * ... and assemble the local contributions to the system matrix and
 * right hand side as also discussed above:
 * 
 * @code
 *     for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *       for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *         {
 * @endcode
 * 
 * Alias the AssemblyScratchData object to keep the lines from
 * getting too long:
 * 
 * @code
 *           const auto &sd = scratch_data;
 *           for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *             copy_data.cell_matrix(i, j) +=
 *               ((sd.fe_values.shape_value(i, q_point) +           // (phi_i +
 *                 delta * (sd.advection_directions[q_point] *      // delta beta
 *                          sd.fe_values.shape_grad(i, q_point))) * // grad phi_i)
 *                sd.advection_directions[q_point] *                // beta
 *                sd.fe_values.shape_grad(j, q_point)) *            // grad phi_j
 *               sd.fe_values.JxW(q_point);                         // dx
 * 
 *           copy_data.cell_rhs(i) +=
 *             (sd.fe_values.shape_value(i, q_point) +           // (phi_i +
 *              delta * (sd.advection_directions[q_point] *      // delta beta
 *                       sd.fe_values.shape_grad(i, q_point))) * // grad phi_i)
 *             sd.rhs_values[q_point] *                          // f
 *             sd.fe_values.JxW(q_point);                        // dx
 *         }
 * 
 * @endcode
 * 
 * Besides the cell terms which we have built up now, the bilinear
 * form of the present problem also contains terms on the boundary of
 * the domain. Therefore, we have to check whether any of the faces of
 * this cell are on the boundary of the domain, and if so assemble the
 * contributions of this face as well. Of course, the bilinear form
 * only contains contributions from the <code>inflow</code> part of
 * the boundary, but to find out whether a certain part of a face of
 * the present cell is part of the inflow boundary, we have to have
 * information on the exact location of the quadrature points and on
 * the direction of flow at this point; we obtain this information
 * using the FEFaceValues object and only decide within the main loop
 * whether a quadrature point is on the inflow boundary.
 * 
 * @code
 *     for (const auto &face : cell->face_iterators())
 *       if (face->at_boundary())
 *         {
 * @endcode
 * 
 * Ok, this face of the present cell is on the boundary of the
 * domain. Just as for the usual FEValues object which we have
 * used in previous examples and also above, we have to
 * reinitialize the FEFaceValues object for the present face:
 * 
 * @code
 *           scratch_data.fe_face_values.reinit(cell, face);
 * 
 * @endcode
 * 
 * For the quadrature points at hand, we ask for the values of
 * the inflow function and for the direction of flow:
 * 
 * @code
 *           scratch_data.boundary_values.value_list(
 *             scratch_data.fe_face_values.get_quadrature_points(),
 *             scratch_data.face_boundary_values);
 *           scratch_data.advection_field.value_list(
 *             scratch_data.fe_face_values.get_quadrature_points(),
 *             scratch_data.face_advection_directions);
 * 
 * @endcode
 * 
 * Now loop over all quadrature points and see whether this face is on
 * the inflow or outflow part of the boundary. The normal
 * vector points out of the cell: since the face is at
 * the boundary, the normal vector points out of the domain,
 * so if the advection direction points into the domain, its
 * scalar product with the normal vector must be negative (to see why
 * this is true, consider the scalar product definition that uses a
 * cosine):
 * 
 * @code
 *           for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
 *             if (scratch_data.fe_face_values.normal_vector(q_point) *
 *                   scratch_data.face_advection_directions[q_point] <
 *                 0.)
 * @endcode
 * 
 * If the face is part of the inflow boundary, then compute the
 * contributions of this face to the global matrix and right
 * hand side, using the values obtained from the
 * FEFaceValues object and the formulae discussed in the
 * introduction:
 * 
 * @code
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                 {
 *                   for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                     copy_data.cell_matrix(i, j) -=
 *                       (scratch_data.face_advection_directions[q_point] *
 *                        scratch_data.fe_face_values.normal_vector(q_point) *
 *                        scratch_data.fe_face_values.shape_value(i, q_point) *
 *                        scratch_data.fe_face_values.shape_value(j, q_point) *
 *                        scratch_data.fe_face_values.JxW(q_point));
 * 
 *                   copy_data.cell_rhs(i) -=
 *                     (scratch_data.face_advection_directions[q_point] *
 *                      scratch_data.fe_face_values.normal_vector(q_point) *
 *                      scratch_data.face_boundary_values[q_point] *
 *                      scratch_data.fe_face_values.shape_value(i, q_point) *
 *                      scratch_data.fe_face_values.JxW(q_point));
 *                 }
 *         }
 * 
 * @endcode
 * 
 * The final piece of information the copy routine needs is the global
 * indices of the degrees of freedom on this cell, so we end by writing
 * them to the local array:
 * 
 * @code
 *     cell->get_dof_indices(copy_data.local_dof_indices);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * The second function we needed to write was the one that copies
 * the local contributions the previous function computed (and
 * put into the AssemblyCopyData object) into the global matrix and right
 * hand side vector objects. This is essentially what we always had
 * as the last block of code when assembling something on every
 * cell. The following should therefore be pretty obvious:
 * 
 * @code
 *   template <int dim>
 *   void
 *   AdvectionProblem<dim>::copy_local_to_global(const AssemblyCopyData &copy_data)
 *   {
 *     hanging_node_constraints.distribute_local_to_global(
 *       copy_data.cell_matrix,
 *       copy_data.cell_rhs,
 *       copy_data.local_dof_indices,
 *       system_matrix,
 *       system_rhs);
 *   }
 * 
 * @endcode
 * 
 * Here comes the linear solver routine. As the system is no longer
 * symmetric positive definite as in all the previous examples, we cannot
 * use the Conjugate Gradient method anymore. Rather, we use a solver that
 * is more general and does not rely on any special properties of the
 * matrix: the GMRES method. GMRES, like the conjugate gradient method,
 * requires a decent preconditioner: we use a Jacobi preconditioner here,
 * which works well enough for this problem.
 * 
 * @code
 *   template <int dim>
 *   void AdvectionProblem<dim>::solve()
 *   {
 *     SolverControl               solver_control(std::max<std::size_t>(1000,
 *                                                        system_rhs.size() / 10),
 *                                  1e-10 * system_rhs.l2_norm());
 *     SolverGMRES<Vector<double>> solver(solver_control);
 *     PreconditionJacobi<SparseMatrix<double>> preconditioner;
 *     preconditioner.initialize(system_matrix, 1.0);
 *     solver.solve(system_matrix, solution, system_rhs, preconditioner);
 * 
 *     Vector<double> residual(dof_handler.n_dofs());
 * 
 *     system_matrix.vmult(residual, solution);
 *     residual -= system_rhs;
 *     std::cout << "   Iterations required for convergence: "
 *               << solver_control.last_step() << '\n'
 *               << "   Max norm of residual:                "
 *               << residual.linfty_norm() << '\n';
 * 
 *     hanging_node_constraints.distribute(solution);
 *   }
 * 
 * @endcode
 * 
 * The following function refines the grid according to the quantity
 * described in the introduction. The respective computations are made in
 * the class <code>GradientEstimation</code>.
 * 
 * @code
 *   template <int dim>
 *   void AdvectionProblem<dim>::refine_grid()
 *   {
 *     Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 * 
 *     GradientEstimation::estimate(dof_handler,
 *                                  solution,
 *                                  estimated_error_per_cell);
 * 
 *     GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                     estimated_error_per_cell,
 *                                                     0.3,
 *                                                     0.03);
 * 
 *     triangulation.execute_coarsening_and_refinement();
 *   }
 * 
 * @endcode
 * 
 * This function is similar to the one in step 6, but since we use a higher
 * degree finite element we save the solution in a different
 * way. Visualization programs like VisIt and Paraview typically only
 * understand data that is associated with nodes: they cannot plot
 * fifth-degree basis functions, which results in a very inaccurate picture
 * of the solution we computed. To get around this we save multiple
 * <em>patches</em> per cell: in 2D we save 64 bilinear `cells' to the VTU
 * file for each cell, and in 3D we save 512. The end result is that the
 * visualization program will use a piecewise linear interpolation of the
 * cubic basis functions: this captures the solution detail and, with most
 * screen resolutions, looks smooth. We save the grid in a separate step
 * with no extra patches so that we have a visual representation of the cell
 * faces.
 *   
 
 * 
 * Version 9.1 of deal.II gained the ability to write higher degree
 * polynomials (i.e., write piecewise bicubic visualization data for our
 * piecewise bicubic solution) VTK and VTU output: however, not all recent
 * versions of ParaView and VisIt (as of 2018) can read this format, so we
 * use the older, more general (but less efficient) approach here.
 * 
 * @code
 *   template <int dim>
 *   void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
 *   {
 *     {
 *       GridOut       grid_out;
 *       std::ofstream output("grid-" + std::to_string(cycle) + ".vtu");
 *       grid_out.write_vtu(triangulation, output);
 *     }
 * 
 *     {
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "solution");
 *       data_out.build_patches(8);
 * 
 * @endcode
 * 
 * VTU output can be expensive, both to compute and to write to
 * disk. Here we ask ZLib, a compression library, to compress the data
 * in a way that maximizes throughput.
 * 
 * @code
 *       DataOutBase::VtkFlags vtk_flags;
 *       vtk_flags.compression_level =
 *         DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
 *       data_out.set_flags(vtk_flags);
 * 
 *       std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
 *       data_out.write_vtu(output);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * ... as is the main loop (setup -- solve -- refine), aside from the number
 * of cycles and the initial grid:
 * 
 * @code
 *   template <int dim>
 *   void AdvectionProblem<dim>::run()
 *   {
 *     for (unsigned int cycle = 0; cycle < 10; ++cycle)
 *       {
 *         std::cout << "Cycle " << cycle << ':' << std::endl;
 * 
 *         if (cycle == 0)
 *           {
 *             GridGenerator::hyper_cube(triangulation, -1, 1);
 *             triangulation.refine_global(3);
 *           }
 *         else
 *           {
 *             refine_grid();
 *           }
 * 
 * 
 *         std::cout << "   Number of active cells:              "
 *                   << triangulation.n_active_cells() << std::endl;
 * 
 *         setup_system();
 * 
 *         std::cout << "   Number of degrees of freedom:        "
 *                   << dof_handler.n_dofs() << std::endl;
 * 
 *         assemble_system();
 *         solve();
 *         output_results(cycle);
 *       }
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="GradientEstimationclassimplementation"></a> 
 * <h3>GradientEstimation class implementation</h3>
 * 
 
 * 
 * Now for the implementation of the <code>GradientEstimation</code> class.
 * Let us start by defining constructors for the
 * <code>EstimateScratchData</code> class used by the
 * <code>estimate_cell()</code> function:
 * 
 * @code
 *   template <int dim>
 *   GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
 *     const FiniteElement<dim> &fe,
 *     const Vector<double> &    solution,
 *     Vector<float> &           error_per_cell)
 *     : fe_midpoint_value(fe,
 *                         QMidpoint<dim>(),
 *                         update_values | update_quadrature_points)
 *     , solution(solution)
 *     , error_per_cell(error_per_cell)
 *     , cell_midpoint_value(1)
 *     , neighbor_midpoint_value(1)
 *   {
 * @endcode
 * 
 * We allocate a vector to hold iterators to all active neighbors of
 * a cell. We reserve the maximal number of active neighbors in order to
 * avoid later reallocations. Note how this maximal number of active
 * neighbors is computed here.
 * 
 * @code
 *     active_neighbors.reserve(GeometryInfo<dim>::faces_per_cell *
 *                              GeometryInfo<dim>::max_children_per_face);
 *   }
 * 
 * 
 *   template <int dim>
 *   GradientEstimation::EstimateScratchData<dim>::EstimateScratchData(
 *     const EstimateScratchData &scratch_data)
 *     : fe_midpoint_value(scratch_data.fe_midpoint_value.get_fe(),
 *                         scratch_data.fe_midpoint_value.get_quadrature(),
 *                         update_values | update_quadrature_points)
 *     , solution(scratch_data.solution)
 *     , error_per_cell(scratch_data.error_per_cell)
 *     , cell_midpoint_value(1)
 *     , neighbor_midpoint_value(1)
 *   {}
 * 
 * 
 * @endcode
 * 
 * Next comes the implementation of the <code>GradientEstimation</code>
 * class. The first function does not much except for delegating work to the
 * other function, but there is a bit of setup at the top.
 *   
 
 * 
 * Before starting with the work, we check that the vector into
 * which the results are written has the right size. Programming
 * mistakes in which one forgets to size arguments correctly at the
 * calling site are quite common. Because the resulting damage from
 * not catching such errors is often subtle (e.g., corruption of
 * data somewhere in memory, or non-reproducible results), it is
 * well worth the effort to check for such things.
 * 
 * @code
 *   template <int dim>
 *   void GradientEstimation::estimate(const DoFHandler<dim> &dof_handler,
 *                                     const Vector<double> & solution,
 *                                     Vector<float> &        error_per_cell)
 *   {
 *     Assert(
 *       error_per_cell.size() == dof_handler.get_triangulation().n_active_cells(),
 *       ExcInvalidVectorLength(error_per_cell.size(),
 *                              dof_handler.get_triangulation().n_active_cells()));
 * 
 *     WorkStream::run(dof_handler.begin_active(),
 *                     dof_handler.end(),
 *                     &GradientEstimation::template estimate_cell<dim>,
 *                     std::function<void(const EstimateCopyData &)>(),
 *                     EstimateScratchData<dim>(dof_handler.get_fe(),
 *                                              solution,
 *                                              error_per_cell),
 *                     EstimateCopyData());
 *   }
 * 
 * 
 * @endcode
 * 
 * Here comes the function that estimates the local error by computing the
 * finite difference approximation of the gradient. The function first
 * computes the list of active neighbors of the present cell and then
 * computes the quantities described in the introduction for each of
 * the neighbors. The reason for this order is that it is not a one-liner
 * to find a given neighbor with locally refined meshes. In principle, an
 * optimized implementation would find neighbors and the quantities
 * depending on them in one step, rather than first building a list of
 * neighbors and in a second step their contributions but we will gladly
 * leave this as an exercise. As discussed before, the worker function
 * passed to WorkStream::run works on "scratch" objects that keep all
 * temporary objects. This way, we do not need to create and initialize
 * objects that are expensive to initialize within the function that does
 * the work every time it is called for a given cell. Such an argument is
 * passed as the second argument. The third argument would be a "copy-data"
 * object (see @ref threads for more information) but we do not actually use
 * any of these here. Since WorkStream::run() insists on passing three
 * arguments, we declare this function with three arguments, but simply
 * ignore the last one.
 *   
 
 * 
 * (This is unsatisfactory from an aesthetic perspective. It can be avoided
 * by using an anonymous (lambda) function. If you allow, let us here show
 * how. First, assume that we had declared this function to only take two
 * arguments by omitting the unused last one. Now, WorkStream::run still
 * wants to call this function with three arguments, so we need to find a
 * way to "forget" the third argument in the call. Simply passing
 * WorkStream::run the pointer to the function as we do above will not do
 * this -- the compiler will complain that a function declared to have two
 * arguments is called with three arguments. However, we can do this by
 * passing the following as the third argument to WorkStream::run():
 * <div class=CodeFragmentInTutorialComment>
 * @code
 * [](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *    EstimateScratchData<dim> &                            scratch_data,
 *    EstimateCopyData &)
 * {
 *   GradientEstimation::estimate_cell<dim>(cell, scratch_data);
 * }
 * @endcode
 * </div>
 * This is not much better than the solution implemented below: either the
 * routine itself must take three arguments or it must be wrapped by
 * something that takes three arguments. We don't use this since adding the
 * unused argument at the beginning is simpler.
 *   
 
 * 
 * Now for the details:
 * 
 * @code
 *   template <int dim>
 *   void GradientEstimation::estimate_cell(
 *     const typename DoFHandler<dim>::active_cell_iterator &cell,
 *     EstimateScratchData<dim> &                            scratch_data,
 *     const EstimateCopyData &)
 *   {
 * @endcode
 * 
 * We need space for the tensor <code>Y</code>, which is the sum of
 * outer products of the y-vectors.
 * 
 * @code
 *     Tensor<2, dim> Y;
 * 
 * @endcode
 * 
 * First initialize the <code>FEValues</code> object, as well as the
 * <code>Y</code> tensor:
 * 
 * @code
 *     scratch_data.fe_midpoint_value.reinit(cell);
 * 
 * @endcode
 * 
 * Now, before we go on, we first compute a list of all active neighbors
 * of the present cell. We do so by first looping over all faces and see
 * whether the neighbor there is active, which would be the case if it
 * is on the same level as the present cell or one level coarser (note
 * that a neighbor can only be once coarser than the present cell, as
 * we only allow a maximal difference of one refinement over a face in
 * deal.II). Alternatively, the neighbor could be on the same level
 * and be further refined; then we have to find which of its children
 * are next to the present cell and select these (note that if a child
 * of a neighbor of an active cell that is next to this active cell,
 * needs necessarily be active itself, due to the one-refinement rule
 * cited above).
 *     
 
 * 
 * Things are slightly different in one space dimension, as there the
 * one-refinement rule does not exist: neighboring active cells may
 * differ in as many refinement levels as they like. In this case, the
 * computation becomes a little more difficult, but we will explain
 * this below.
 *     
 
 * 
 * Before starting the loop over all neighbors of the present cell, we
 * have to clear the array storing the iterators to the active
 * neighbors, of course.
 * 
 * @code
 *     scratch_data.active_neighbors.clear();
 *     for (const auto face_n : cell->face_indices())
 *       if (!cell->at_boundary(face_n))
 *         {
 * @endcode
 * 
 * First define an abbreviation for the iterator to the face and
 * the neighbor
 * 
 * @code
 *           const auto face     = cell->face(face_n);
 *           const auto neighbor = cell->neighbor(face_n);
 * 
 * @endcode
 * 
 * Then check whether the neighbor is active. If it is, then it
 * is on the same level or one level coarser (if we are not in
 * 1D), and we are interested in it in any case.
 * 
 * @code
 *           if (neighbor->is_active())
 *             scratch_data.active_neighbors.push_back(neighbor);
 *           else
 *             {
 * @endcode
 * 
 * If the neighbor is not active, then check its children.
 * 
 * @code
 *               if (dim == 1)
 *                 {
 * @endcode
 * 
 * To find the child of the neighbor which bounds to the
 * present cell, successively go to its right child if
 * we are left of the present cell (n==0), or go to the
 * left child if we are on the right (n==1), until we
 * find an active cell.
 * 
 * @code
 *                   auto neighbor_child = neighbor;
 *                   while (neighbor_child->has_children())
 *                     neighbor_child = neighbor_child->child(face_n == 0 ? 1 : 0);
 * 
 * @endcode
 * 
 * As this used some non-trivial geometrical intuition,
 * we might want to check whether we did it right,
 * i.e., check whether the neighbor of the cell we found
 * is indeed the cell we are presently working
 * on. Checks like this are often useful and have
 * frequently uncovered errors both in algorithms like
 * the line above (where it is simple to involuntarily
 * exchange <code>n==1</code> for <code>n==0</code> or
 * the like) and in the library (the assumptions
 * underlying the algorithm above could either be wrong,
 * wrongly documented, or are violated due to an error
 * in the library). One could in principle remove such
 * checks after the program works for some time, but it
 * might be a good things to leave it in anyway to check
 * for changes in the library or in the algorithm above.
 *                   
 
 * 
 * Note that if this check fails, then this is certainly
 * an error that is irrecoverable and probably qualifies
 * as an internal error. We therefore use a predefined
 * exception class to throw here.
 * 
 * @code
 *                   Assert(neighbor_child->neighbor(face_n == 0 ? 1 : 0) == cell,
 *                          ExcInternalError());
 * 
 * @endcode
 * 
 * If the check succeeded, we push the active neighbor
 * we just found to the stack we keep:
 * 
 * @code
 *                   scratch_data.active_neighbors.push_back(neighbor_child);
 *                 }
 *               else
 * @endcode
 * 
 * If we are not in 1d, we collect all neighbor children
 * `behind' the subfaces of the current face and move on:
 * 
 * @code
 *                 for (unsigned int subface_n = 0; subface_n < face->n_children();
 *                      ++subface_n)
 *                   scratch_data.active_neighbors.push_back(
 *                     cell->neighbor_child_on_subface(face_n, subface_n));
 *             }
 *         }
 * 
 * @endcode
 * 
 * OK, now that we have all the neighbors, lets start the computation
 * on each of them. First we do some preliminaries: find out about the
 * center of the present cell and the solution at this point. The
 * latter is obtained as a vector of function values at the quadrature
 * points, of which there are only one, of course. Likewise, the
 * position of the center is the position of the first (and only)
 * quadrature point in real space.
 * 
 * @code
 *     const Point<dim> this_center =
 *       scratch_data.fe_midpoint_value.quadrature_point(0);
 * 
 *     scratch_data.fe_midpoint_value.get_function_values(
 *       scratch_data.solution, scratch_data.cell_midpoint_value);
 * 
 * @endcode
 * 
 * Now loop over all active neighbors and collect the data we
 * need.
 * 
 * @code
 *     Tensor<1, dim> projected_gradient;
 *     for (const auto &neighbor : scratch_data.active_neighbors)
 *       {
 * @endcode
 * 
 * Then get the center of the neighbor cell and the value of the
 * finite element function at that point. Note that for this
 * information we have to reinitialize the <code>FEValues</code>
 * object for the neighbor cell.
 * 
 * @code
 *         scratch_data.fe_midpoint_value.reinit(neighbor);
 *         const Point<dim> neighbor_center =
 *           scratch_data.fe_midpoint_value.quadrature_point(0);
 * 
 *         scratch_data.fe_midpoint_value.get_function_values(
 *           scratch_data.solution, scratch_data.neighbor_midpoint_value);
 * 
 * @endcode
 * 
 * Compute the vector <code>y</code> connecting the centers of the
 * two cells. Note that as opposed to the introduction, we denote
 * by <code>y</code> the normalized difference vector, as this is
 * the quantity used everywhere in the computations.
 * 
 * @code
 *         Tensor<1, dim> y        = neighbor_center - this_center;
 *         const double   distance = y.norm();
 *         y /= distance;
 * 
 * @endcode
 * 
 * Then add up the contribution of this cell to the Y matrix...
 * 
 * @code
 *         for (unsigned int i = 0; i < dim; ++i)
 *           for (unsigned int j = 0; j < dim; ++j)
 *             Y[i][j] += y[i] * y[j];
 * 
 * @endcode
 * 
 * ... and update the sum of difference quotients:
 * 
 * @code
 *         projected_gradient += (scratch_data.neighbor_midpoint_value[0] -
 *                                scratch_data.cell_midpoint_value[0]) /
 *                               distance * y;
 *       }
 * 
 * @endcode
 * 
 * If now, after collecting all the information from the neighbors, we
 * can determine an approximation of the gradient for the present
 * cell, then we need to have passed over vectors <code>y</code> which
 * span the whole space, otherwise we would not have all components of
 * the gradient. This is indicated by the invertibility of the matrix.
 *     
 
 * 
 * If the matrix is not invertible, then the present
 * cell had an insufficient number of active neighbors. In contrast to
 * all previous cases (where we raised exceptions) this is, however,
 * not a programming error: it is a runtime error that can happen in
 * optimized mode even if it ran well in debug mode, so it is
 * reasonable to try to catch this error also in optimized mode. For
 * this case, there is the <code>AssertThrow</code> macro: it checks
 * the condition like the <code>Assert</code> macro, but not only in
 * debug mode; it then outputs an error message, but instead of
 * aborting the program as in the case of the <code>Assert</code>
 * macro, the exception is thrown using the <code>throw</code> command
 * of C++. This way, one has the possibility to catch this error and
 * take reasonable counter actions. One such measure would be to
 * refine the grid globally, as the case of insufficient directions
 * can not occur if every cell of the initial grid has been refined at
 * least once.
 * 
 * @code
 *     AssertThrow(determinant(Y) != 0, ExcInsufficientDirections());
 * 
 * @endcode
 * 
 * If, on the other hand, the matrix is invertible, then invert it,
 * multiply the other quantity with it, and compute the estimated error
 * using this quantity and the correct powers of the mesh width:
 * 
 * @code
 *     const Tensor<2, dim> Y_inverse = invert(Y);
 * 
 *     const Tensor<1, dim> gradient = Y_inverse * projected_gradient;
 * 
 * @endcode
 * 
 * The last part of this function is the one where we write into
 * the element of the output vector what we have just
 * computed. The address of this vector has been stored in the
 * scratch data object, and all we have to do is know how to get
 * at the correct element inside this vector -- but we can ask the
 * cell we're on the how-manyth active cell it is for this:
 * 
 * @code
 *     scratch_data.error_per_cell(cell->active_cell_index()) =
 *       (std::pow(cell->diameter(), 1 + 1.0 * dim / 2) * gradient.norm());
 *   }
 * } // namespace Step9
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Mainfunction"></a> 
 * <h3>Main function</h3>
 * 
 
 * 
 * The <code>main</code> function is similar to the previous examples. The
 * primary difference is that we use MultithreadInfo to set the maximum
 * number of threads (see the documentation module @ref threads
 * "Parallel computing with multiple processors accessing shared memory"
 * for more information). The number of threads used is the minimum of the
 * environment variable DEAL_II_NUM_THREADS and the parameter of
 * <code>set_thread_limit</code>. If no value is given to
 * <code>set_thread_limit</code>, the default value from the Intel Threading
 * Building Blocks (TBB) library is used. If the call to
 * <code>set_thread_limit</code> is omitted, the number of threads will be
 * chosen by TBB independently of DEAL_II_NUM_THREADS.
 * 
 * @code
 * int main()
 * {
 *   using namespace dealii;
 *   try
 *     {
 *       MultithreadInfo::set_thread_limit();
 * 
 *       Step9::AdvectionProblem<2> advection_problem_2d;
 *       advection_problem_2d.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>
 
 
 
The results of this program are not particularly spectacular. They
consist of the console output, some grid files, and the solution on
each of these grids. First for the console output:
@code
Cycle 0:
   Number of active cells:              64
   Number of degrees of freedom:        1681
   Iterations required for convergence: 298
   Max norm of residual:                3.60316e-12
Cycle 1:
   Number of active cells:              124
   Number of degrees of freedom:        3537
   Iterations required for convergence: 415
   Max norm of residual:                3.70682e-12
Cycle 2:
   Number of active cells:              247
   Number of degrees of freedom:        6734
   Iterations required for convergence: 543
   Max norm of residual:                7.19716e-13
Cycle 3:
   Number of active cells:              502
   Number of degrees of freedom:        14105
   Iterations required for convergence: 666
   Max norm of residual:                3.45628e-13
Cycle 4:
   Number of active cells:              1003
   Number of degrees of freedom:        27462
   Iterations required for convergence: 1064
   Max norm of residual:                1.86495e-13
Cycle 5:
   Number of active cells:              1993
   Number of degrees of freedom:        55044
   Iterations required for convergence: 1251
   Max norm of residual:                1.28765e-13
Cycle 6:
   Number of active cells:              3985
   Number of degrees of freedom:        108492
   Iterations required for convergence: 2035
   Max norm of residual:                6.78085e-14
Cycle 7:
   Number of active cells:              7747
   Number of degrees of freedom:        210612
   Iterations required for convergence: 2187
   Max norm of residual:                2.61457e-14
Cycle 8:
   Number of active cells:              15067
   Number of degrees of freedom:        406907
   Iterations required for convergence: 3079
   Max norm of residual:                2.9932e-14
Cycle 9:
   Number of active cells:              29341
   Number of degrees of freedom:        780591
   Iterations required for convergence: 3913
   Max norm of residual:                8.15689e-15
@endcode
 
Quite a number of cells are used on the finest level to resolve the features of
the solution. Here are the fourth and tenth grids:
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-9-grid-3.png"
         alt="Fourth grid in the refinement cycle, showing some adaptivity to features."
         width="400" height="400">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-9-grid-9.png"
         alt="Tenth grid in the refinement cycle, showing that the waves are fully captured."
         width="400" height="400">
  </div>
</div>
and the fourth and tenth solutions:
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-9-solution-3.png"
         alt="Fourth solution, showing that we resolve most features but some
         are sill unresolved and appear blury."
         width="400" height="400">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-9-solution-9.png"
         alt="Tenth solution, showing a fully resolved flow."
         width="400" height="400">
  </div>
</div>
and both the grid and solution zoomed in:
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-9-solution-3-zoom.png"
         alt="Detail of the fourth solution, showing that we resolve most
         features but some are sill unresolved and appear blury. In particular,
         the larger cells need to be refined."
         width="400" height="400">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-9-solution-9-zoom.png"
         alt="Detail of the tenth solution, showing that we needed a lot more
         cells than were present in the fourth solution."
         width="400" height="400">
  </div>
</div>
 
The solution is created by that part that is transported along the wiggly
advection field from the left and lower boundaries to the top right, and the
part that is created by the source in the lower left corner, and the results of
which are also transported along. The grid shown above is well-adapted to
resolve these features. The comparison between plots shows that, even though we
are using a high-order approximation, we still need adaptive mesh refinement to
fully resolve the wiggles.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-9.cc"
*/