step-27.h

Go to the documentation of this file.

) const
 *   {
 *     double product = 1;
 *     for (unsigned int d = 0; d < dim; ++d)
 *       product *= (p[d] + 1);
 *     return product;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Implementationofthemainclass"></a> 
 * <h3>Implementation of the main class</h3>
 * 
 
 * 
 * 
 * <a name="LaplaceProblemLaplaceProblemconstructor"></a> 
 * <h4>LaplaceProblem::LaplaceProblem constructor</h4>
 * 
 
 * 
 * The constructor of this class is fairly straightforward. It associates
 * the DoFHandler object with the triangulation, and then sets the
 * maximal polynomial degree to 7 (in 1d and 2d) or 5 (in 3d and higher). We
 * do so because using higher order polynomial degrees becomes prohibitively
 * expensive, especially in higher space dimensions.
 *   
 
 * 
 * Following this, we fill the collections of finite element, and cell and
 * face quadrature objects. We start with quadratic elements, and each
 * quadrature formula is chosen so that it is appropriate for the matching
 * finite element in the hp::FECollection object.
 * 
 * @code
 *   template <int dim>
 *   LaplaceProblem<dim>::LaplaceProblem()
 *     : dof_handler(triangulation)
 *     , max_degree(dim <= 2 ? 7 : 5)
 *   {
 *     for (unsigned int degree = 2; degree <= max_degree; ++degree)
 *       {
 *         fe_collection.push_back(FE_Q<dim>(degree));
 *         quadrature_collection.push_back(QGauss<dim>(degree + 1));
 *         face_quadrature_collection.push_back(QGauss<dim - 1>(degree + 1));
 *       }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemLaplaceProblemdestructor"></a> 
 * <h4>LaplaceProblem::~LaplaceProblem destructor</h4>
 * 
 
 * 
 * The destructor is unchanged from what we already did in @ref step_6 "step-6":
 * 
 * @code
 *   template <int dim>
 *   LaplaceProblem<dim>::~LaplaceProblem()
 *   {
 *     dof_handler.clear();
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemsetup_system"></a> 
 * <h4>LaplaceProblem::setup_system</h4>
 *   
 
 * 
 * This function is again a verbatim copy of what we already did in
 * @ref step_6 "step-6". Despite function calls with exactly the same names and arguments,
 * the algorithms used internally are different in some aspect since the
 * dof_handler variable here is in @f$hp@f$-mode.
 * 
 * @code
 *   template <int dim>
 *   void LaplaceProblem<dim>::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe_collection);
 * 
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 * 
 *     constraints.clear();
 *     DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *     VectorTools::interpolate_boundary_values(dof_handler,
 *                                              0,
 *                                              Functions::ZeroFunction<dim>(),
 *                                              constraints);
 *     constraints.close();
 * 
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
 *     sparsity_pattern.copy_from(dsp);
 * 
 *     system_matrix.reinit(sparsity_pattern);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemassemble_system"></a> 
 * <h4>LaplaceProblem::assemble_system</h4>
 * 
 
 * 
 * This is the function that assembles the global matrix and right hand side
 * vector from the local contributions of each cell. Its main working is as
 * has been described in many of the tutorial programs before. The
 * significant deviations are the ones necessary for <i>hp</i> finite
 * element methods. In particular, that we need to use a collection of
 * FEValues object (implemented through the hp::FEValues class), and that we
 * have to eliminate constrained degrees of freedom already when copying
 * local contributions into global objects. Both of these are explained in
 * detail in the introduction of this program.
 *   
 
 * 
 * One other slight complication is the fact that because we use different
 * polynomial degrees on different cells, the matrices and vectors holding
 * local contributions do not have the same size on all cells. At the
 * beginning of the loop over all cells, we therefore each time have to
 * resize them to the correct size (given by <code>dofs_per_cell</code>).
 * Because these classes are implemented in such a way that reducing the size
 * of a matrix or vector does not release the currently allocated memory
 * (unless the new size is zero), the process of resizing at the beginning of
 * the loop will only require re-allocation of memory during the first few
 * iterations. Once we have found in a cell with the maximal finite element
 * degree, no more re-allocations will happen because all subsequent
 * <code>reinit</code> calls will only set the size to something that fits the
 * currently allocated memory. This is important since allocating memory is
 * expensive, and doing so every time we visit a new cell would take
 * significant compute time.
 * 
 * @code
 *   template <int dim>
 *   void LaplaceProblem<dim>::assemble_system()
 *   {
 *     hp::FEValues<dim> hp_fe_values(fe_collection,
 *                                    quadrature_collection,
 *                                    update_values | update_gradients |
 *                                      update_quadrature_points |
 *                                      update_JxW_values);
 * 
 *     RightHandSide<dim> rhs_function;
 * 
 *     FullMatrix<double> cell_matrix;
 *     Vector<double>     cell_rhs;
 * 
 *     std::vector<types::global_dof_index> local_dof_indices;
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
 * 
 *         cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
 *         cell_matrix = 0;
 * 
 *         cell_rhs.reinit(dofs_per_cell);
 *         cell_rhs = 0;
 * 
 *         hp_fe_values.reinit(cell);
 * 
 *         const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
 * 
 *         std::vector<double> rhs_values(fe_values.n_quadrature_points);
 *         rhs_function.value_list(fe_values.get_quadrature_points(), rhs_values);
 * 
 *         for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
 *              ++q_point)
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 cell_matrix(i, j) +=
 *                   (fe_values.shape_grad(i, q_point) * // grad phi_i(x_q)
 *                    fe_values.shape_grad(j, q_point) * // grad phi_j(x_q)
 *                    fe_values.JxW(q_point));           // dx
 * 
 *               cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
 *                               rhs_values[q_point] *               // f(x_q)
 *                               fe_values.JxW(q_point));            // dx
 *             }
 * 
 *         local_dof_indices.resize(dofs_per_cell);
 *         cell->get_dof_indices(local_dof_indices);
 * 
 *         constraints.distribute_local_to_global(
 *           cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
 *       }
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemsolve"></a> 
 * <h4>LaplaceProblem::solve</h4>
 * 
 
 * 
 * The function solving the linear system is entirely unchanged from
 * previous examples. We simply try to reduce the initial residual (which
 * equals the @f$l_2@f$ norm of the right hand side) by a certain factor:
 * 
 * @code
 *   template <int dim>
 *   void LaplaceProblem<dim>::solve()
 *   {
 *     SolverControl            solver_control(system_rhs.size(),
 *                                  1e-12 * system_rhs.l2_norm());
 *     SolverCG<Vector<double>> cg(solver_control);
 * 
 *     PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *     preconditioner.initialize(system_matrix, 1.2);
 * 
 *     cg.solve(system_matrix, solution, system_rhs, preconditioner);
 * 
 *     constraints.distribute(solution);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblempostprocess"></a> 
 * <h4>LaplaceProblem::postprocess</h4>
 * 
 
 * 
 * After solving the linear system, we will want to postprocess the
 * solution. Here, all we do is to estimate the error, estimate the local
 * smoothness of the solution as described in the introduction, then write
 * graphical output, and finally refine the mesh in both @f$h@f$ and @f$p@f$
 * according to the indicators computed before. We do all this in the same
 * function because we want the estimated error and smoothness indicators
 * not only for refinement, but also include them in the graphical output.
 * 
 * @code
 *   template <int dim>
 *   void LaplaceProblem<dim>::postprocess(const unsigned int cycle)
 *   {
 * @endcode
 * 
 * Let us start with computing estimated error and smoothness indicators,
 * which each are one number for each active cell of our
 * triangulation. For the error indicator, we use the KellyErrorEstimator
 * class as always.
 * 
 * @code
 *     Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *     KellyErrorEstimator<dim>::estimate(
 *       dof_handler,
 *       face_quadrature_collection,
 *       std::map<types::boundary_id, const Function<dim> *>(),
 *       solution,
 *       estimated_error_per_cell);
 * 
 * @endcode
 * 
 * Estimating the smoothness is performed with the method of decaying
 * expansion coefficients as outlined in the introduction. We will first
 * need to create an object capable of transforming the finite element
 * solution on every single cell into a sequence of Fourier series
 * coefficients. The SmoothnessEstimator namespace offers a factory function
 * for such a FESeries::Fourier object that is optimized for the process of
 * estimating smoothness. The actual determination of the decay of Fourier
 * coefficients on every individual cell then happens in the last function.
 * 
 * @code
 *     Vector<float> smoothness_indicators(triangulation.n_active_cells());
 *     FESeries::Fourier<dim> fourier =
 *       SmoothnessEstimator::Fourier::default_fe_series(fe_collection);
 *     SmoothnessEstimator::Fourier::coefficient_decay(fourier,
 *                                                     dof_handler,
 *                                                     solution,
 *                                                     smoothness_indicators);
 * 
 * @endcode
 * 
 * Next we want to generate graphical output. In addition to the two
 * estimated quantities derived above, we would also like to output the
 * polynomial degree of the finite elements used on each of the elements
 * on the mesh.
 *     
 
 * 
 * The way to do that requires that we loop over all cells and poll the
 * active finite element index of them using
 * <code>cell-@>active_fe_index()</code>. We then use the result of this
 * operation and query the finite element collection for the finite
 * element with that index, and finally determine the polynomial degree of
 * that element. The result we put into a vector with one element per
 * cell. The DataOut class requires this to be a vector of
 * <code>float</code> or <code>double</code>, even though our values are
 * all integers, so that is what we use:
 * 
 * @code
 *     {
 *       Vector<float> fe_degrees(triangulation.n_active_cells());
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         fe_degrees(cell->active_cell_index()) =
 *           fe_collection[cell->active_fe_index()].degree;
 * 
 * @endcode
 * 
 * With now all data vectors available -- solution, estimated errors and
 * smoothness indicators, and finite element degrees --, we create a
 * DataOut object for graphical output and attach all data:
 * 
 * @code
 *       DataOut<dim> data_out;
 * 
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "solution");
 *       data_out.add_data_vector(estimated_error_per_cell, "error");
 *       data_out.add_data_vector(smoothness_indicators, "smoothness");
 *       data_out.add_data_vector(fe_degrees, "fe_degree");
 *       data_out.build_patches();
 * 
 * @endcode
 * 
 * The final step in generating output is to determine a file name, open
 * the file, and write the data into it (here, we use VTK format):
 * 
 * @code
 *       const std::string filename =
 *         "solution-" + Utilities::int_to_string(cycle, 2) + ".vtk";
 *       std::ofstream output(filename);
 *       data_out.write_vtk(output);
 *     }
 * 
 * @endcode
 * 
 * After this, we would like to actually refine the mesh, in both @f$h@f$ and
 * @f$p@f$. The way we are going to do this is as follows: first, we use the
 * estimated error to flag those cells for refinement that have the
 * largest error. This is what we have always done:
 * 
 * @code
 *     {
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       estimated_error_per_cell,
 *                                                       0.3,
 *                                                       0.03);
 * 
 * @endcode
 * 
 * Next we would like to figure out which of the cells that have been
 * flagged for refinement should actually have @f$p@f$ increased instead of
 * @f$h@f$ decreased. The strategy we choose here is that we look at the
 * smoothness indicators of those cells that are flagged for refinement,
 * and increase @f$p@f$ for those with a smoothness larger than a certain
 * relative threshold. In other words, for every cell for which (i) the
 * refinement flag is set, (ii) the smoothness indicator is larger than
 * the threshold, and (iii) we still have a finite element with a
 * polynomial degree higher than the current one in the finite element
 * collection, we will assign a future FE index that corresponds to a
 * polynomial with degree one higher than it currently is. The following
 * function is capable of doing exactly this. Absent any better
 * strategies, we will set the threshold via interpolation between the
 * minimal and maximal smoothness indicators on cells flagged for
 * refinement. Since the corner singularities are strongly localized, we
 * will favor @f$p@f$- over @f$h@f$-refinement quantitatively. We achieve this
 * with a low threshold by setting a small interpolation factor of 0.2. In
 * the same way, we deal with cells that are going to be coarsened and
 * decrease their polynomial degree when their smoothness indicator is
 * below the corresponding threshold determined on cells to be coarsened.
 * 
 * @code
 *       hp::Refinement::p_adaptivity_from_relative_threshold(
 *         dof_handler, smoothness_indicators, 0.2, 0.2);
 * 
 * @endcode
 * 
 * The above function only determines whether the polynomial degree will
 * change via future FE indices, but does not manipulate the
 * @f$h@f$-refinement flags. So for cells that are flagged for both refinement
 * categories, we prefer @f$p@f$- over @f$h@f$-refinement. The following function
 * call ensures that only one of @f$p@f$- or @f$h@f$-refinement is imposed, and
 * not both at once.
 * 
 * @code
 *       hp::Refinement::choose_p_over_h(dof_handler);
 * 
 * @endcode
 * 
 * For grid adaptive refinement, we ensure a 2:1 mesh balance by limiting
 * the difference of refinement levels of neighboring cells to one by
 * calling Triangulation::prepare_coarsening_and_refinement(). We would
 * like to achieve something similar for the p-levels of neighboring
 * cells: levels of future finite elements are not allowed to differ by
 * more than a specified difference. With its default parameters, a call
 * of hp::Refinement::limit_p_level_difference() ensures that their level
 * difference is limited to one. This will not necessarily decrease the
 * number of hanging nodes in the domain, but makes sure that high order
 * polynomials are not constrained to much lower polynomials on faces,
 * e.g., fifth order to second order polynomials.
 * 
 * @code
 *       triangulation.prepare_coarsening_and_refinement();
 *       hp::Refinement::limit_p_level_difference(dof_handler);
 * 
 * @endcode
 * 
 * At the end of this procedure, we then refine the mesh. During this
 * process, children of cells undergoing bisection inherit their mother
 * cell's finite element index. Further, future finite element indices
 * will turn into active ones, so that the new finite elements will be
 * assigned to cells after the next call of DoFHandler::distribute_dofs().
 * 
 * @code
 *       triangulation.execute_coarsening_and_refinement();
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemcreate_coarse_grid"></a> 
 * <h4>LaplaceProblem::create_coarse_grid</h4>
 * 
 
 * 
 * The following function is used when creating the initial grid. The grid we
 * would like to create is actually similar to the one from @ref step_14 "step-14", i.e., the
 * square domain with the square hole in the middle. It can be generated by
 * exactly the same function. However, since its implementation is only a
 * specialization of the 2d case, we will present a different way of creating
 * this domain which is dimension independent.
 *   
 
 * 
 * We first create a hypercube triangulation with enough cells so that it
 * already holds our desired domain @f$[-1,1]^d@f$, subdivided into @f$4^d@f$ cells.
 * We then remove those cells in the center of the domain by testing the
 * coordinate values of the vertices on each cell. In the end, we refine the
 * so created grid globally as usual.
 * 
 * @code
 *   template <int dim>
 *   void LaplaceProblem<dim>::create_coarse_grid()
 *   {
 *     Triangulation<dim> cube;
 *     GridGenerator::subdivided_hyper_cube(cube, 4, -1., 1.);
 * 
 *     std::set<typename Triangulation<dim>::active_cell_iterator> cells_to_remove;
 *     for (const auto &cell : cube.active_cell_iterators())
 *       for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
 *         if (cell->vertex(v).square() < .1)
 *           cells_to_remove.insert(cell);
 * 
 *     GridGenerator::create_triangulation_with_removed_cells(cube,
 *                                                            cells_to_remove,
 *                                                            triangulation);
 * 
 *     triangulation.refine_global(3);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemrun"></a> 
 * <h4>LaplaceProblem::run</h4>
 * 
 
 * 
 * This function implements the logic of the program, as did the respective
 * function in most of the previous programs already, see for example @ref step_6 "step-6".
 *   
 
 * 
 * Basically, it contains the adaptive loop: in the first iteration create a
 * coarse grid, and then set up the linear system, assemble it, solve, and
 * postprocess the solution including mesh refinement. Then start over
 * again. In the meantime, also output some information for those staring at
 * the screen trying to figure out what the program does:
 * 
 * @code
 *   template <int dim>
 *   void LaplaceProblem<dim>::run()
 *   {
 *     for (unsigned int cycle = 0; cycle < 6; ++cycle)
 *       {
 *         std::cout << "Cycle " << cycle << ':' << std::endl;
 * 
 *         if (cycle == 0)
 *           create_coarse_grid();
 * 
 *         setup_system();
 * 
 *         std::cout << "   Number of active cells      : "
 *                   << triangulation.n_active_cells() << std::endl
 *                   << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *                   << std::endl
 *                   << "   Number of constraints       : "
 *                   << constraints.n_constraints() << std::endl;
 * 
 *         assemble_system();
 *         solve();
 *         postprocess(cycle);
 *       }
 *   }
 * } // namespace Step27
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Themainfunction"></a> 
 * <h3>The main function</h3>
 * 
 
 * 
 * The main function is again verbatim what we had before: wrap creating and
 * running an object of the main class into a <code>try</code> block and catch
 * whatever exceptions are thrown, thereby producing meaningful output if
 * anything should go wrong:
 * 
 * @code
 * int main()
 * {
 *   try
 *     {
 *       using namespace Step27;
 * 
 *       LaplaceProblem<2> laplace_problem;
 *       laplace_problem.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>
 
 
In this section, we discuss a few results produced from running the
current tutorial program. More results, in particular the extension to
3d calculations and determining how much compute time the individual
components of the program take, are given in the @ref hp_paper "hp-paper".
 
When run, this is what the program produces:
 
@code
> make run
[ 66%] Built target @ref step_27 "step-27"
[100%] Run @ref step_27 "step-27" with Release configuration
Cycle 0:
   Number of active cells      : 768
   Number of degrees of freedom: 3264
   Number of constraints       : 384
Cycle 1:
   Number of active cells      : 807
   Number of degrees of freedom: 4764
   Number of constraints       : 756
Cycle 2:
   Number of active cells      : 927
   Number of degrees of freedom: 8226
   Number of constraints       : 1856
Cycle 3:
   Number of active cells      : 978
   Number of degrees of freedom: 12146
   Number of constraints       : 2944
Cycle 4:
   Number of active cells      : 1104
   Number of degrees of freedom: 16892
   Number of constraints       : 3998
Cycle 5:
   Number of active cells      : 1149
   Number of degrees of freedom: 22078
   Number of constraints       : 5230
@endcode
 
The first thing we learn from this is that the number of constrained degrees
of freedom is on the order of 20-25% of the total number of degrees of
freedom, at least on the later grids when we have elements of relatively
high order (in 3d, the fraction of constrained degrees of freedom can be up
to 30%). This is, in fact, on the same order of magnitude as for
non-@f$hp@f$-discretizations. For example, in the last step of the @ref step_6 "step-6"
program, we have 18353 degrees of freedom, 4432 of which are
constrained. The difference is that in the latter program, each constrained
hanging node is constrained against only the two adjacent degrees of
freedom, whereas in the @f$hp@f$-case, constrained nodes are constrained against
many more degrees of freedom. Note also that the current program also
includes nodes subject to Dirichlet boundary conditions in the list of
constraints. In cycle 0, all the constraints are actually because of
boundary conditions.
 
Of maybe more interest is to look at the graphical output. First, here is the
solution of the problem:
 
<img src="https://www.dealii.org/images/steps/developer/step-27-solution.png"
     alt="Elevation plot of the solution, showing the lack of regularity near
          the interior (reentrant) corners."
     width="200" height="200">
 
Secondly, let us look at the sequence of meshes generated:
 
<div class="threecolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-00.svg"
         alt="Triangulation containing reentrant corners without adaptive refinement."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-01.svg"
         alt="Triangulation containing reentrant corners with one level of
         refinement. New cells are placed near the corners."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-02.svg"
         alt="Triangulation containing reentrant corners with two levels of
         refinement. New cells are placed near the corners."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-03.svg"
         alt="Triangulation containing reentrant corners with three levels of
         refinement. New cells are placed near the corners."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-04.svg"
         alt="Triangulation containing reentrant corners with four levels of
         refinement. New cells are placed near the corners."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.mesh-05.svg"
         alt="Triangulation containing reentrant corners with five levels of
         refinement. New cells are placed near the corners."
         width="200" height="200">
  </div>
</div>
 
It is clearly visible how the mesh is refined near the corner singularities,
as one would expect it. More interestingly, we should be curious to see the
distribution of finite element polynomial degrees to these mesh cells, where
the lightest color corresponds to degree two and the darkest one corresponds
to degree seven:
 
<div class="threecolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-00.svg"
         alt="Initial grid where all cells contain just biquadratic functions."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-01.svg"
         alt="Depiction of local approximation degrees after one refinement."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-02.svg"
         alt="Depiction of local approximation degrees after two refinements."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-03.svg"
         alt="Depiction of local approximation degrees after three refinements."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-04.svg"
         alt="Depiction of local approximation degrees after four refinements."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.fe_degree-05.svg"
         alt="Depiction of local approximation degrees after five refinements."
         width="200" height="200">
  </div>
</div>
 
While this is certainly not a perfect arrangement, it does make some sense: we
use low order elements close to boundaries and corners where regularity is
low. On the other hand, higher order elements are used where (i) the error was
at one point fairly large, i.e., mainly in the general area around the corner
singularities and in the top right corner where the solution is large, and
(ii) where the solution is smooth, i.e., far away from the boundary.
 
This arrangement of polynomial degrees of course follows from our smoothness
estimator. Here is the estimated smoothness of the solution, with darker colors
indicating least smoothness and lighter indicating the smoothest areas:
 
<div class="threecolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-00.svg"
         alt="Estimated regularity per cell on the initial grid."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-01.svg"
         alt="Depiction of the estimated regularity per cell after one refinement."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-02.svg"
         alt="Depiction of the estimated regularity per cell after two refinements."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-03.svg"
         alt="Depiction of the estimated regularity per cell after three refinements."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-04.svg"
         alt="Depiction of the estimated regularity per cell after four refinements."
         width="200" height="200">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27.smoothness-05.svg"
         alt="Depiction of the estimated regularity per cell after five refinements."
         width="200" height="200">
  </div>
</div>
 
The primary conclusion one can draw from this is that the loss of regularity at
the internal corners is a highly localized phenomenon; it only seems to impact
the cells adjacent to the corner itself, so when we refine the mesh the black
coloring is no longer visible. Besides the corners, this sequence of plots
implies that the smoothness estimates are somewhat independent of the mesh
refinement, particularly when we are far away from boundaries.
It is also obvious that the smoothness estimates are independent of the actual
size of the solution (see the picture of the solution above), as it should be.
A point of larger concern, however, is that one realizes on closer inspection
that the estimator we have overestimates the smoothness of the solution on
cells with hanging nodes. This in turn leads to higher polynomial degrees in
these areas, skewing the allocation of finite elements onto cells.
 
We have no good explanation for this effect at the moment. One theory is that
the numerical solution on cells with hanging nodes is, of course, constrained
and therefore not entirely free to explore the function space to get close to
the exact solution. This lack of degrees of freedom may manifest itself by
yielding numerical solutions on these cells with suppressed oscillation,
meaning a higher degree of smoothness. The estimator picks this signal up and
the estimated smoothness overestimates the actual value. However, a definite
answer to what is going on currently eludes the authors of this program.
 
The bigger question is, of course, how to avoid this problem. Possibilities
include estimating the smoothness not on single cells, but cell assemblies or
patches surrounding each cell. It may also be possible to find simple
correction factors for each cell depending on the number of constrained
degrees of freedom it has. In either case, there are ample opportunities for
further research on finding good @f$hp@f$-refinement criteria. On the other hand,
the main point of the current program was to demonstrate using the
@f$hp@f$-technology in deal.II, which is unaffected by our use of a possible
sub-optimal refinement criterion.
 
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
<a name="Differenthpdecisionstrategies"></a><h4>Different hp-decision strategies</h4>
 
 
This tutorial demonstrates only one particular strategy to decide between @f$h@f$- and
@f$p@f$-adaptation. In fact, there are many more ways to automatically decide on the
adaptation type, of which a few are already implemented in deal.II:
<ul>
  <li><i>Fourier coefficient decay:</i> This is the strategy currently
  implemented in this tutorial. For more information on this strategy, see
  the general documentation of the SmoothnessEstimator::Fourier namespace.</li>
 
  <li><i>Legendre coefficient decay:</i> This strategy is quite similar
  to the current one, but uses Legendre series expansion rather than the
  Fourier one: instead of sinusoids as basis functions, this strategy uses
  Legendre polynomials. Of course, since we approximate the solution using a
  finite-dimensional polynomial on each cell, the expansion of the solution in
  Legendre polynomials is also finite and, consequently, when we talk about the
  "decay" of this expansion, we can only consider the finitely many nonzero
  coefficients of this expansion, rather than think about it in asymptotic terms.
  But, if we have enough of these coefficients, we can certainly think of the
  decay of these coefficients as characteristic of the decay of the coefficients
  of the exact solution (which is, in general, not polynomial and so will have an
  infinite Legendre expansion), and considering the coefficients we have should
  reveal something about the properties of the exact solution.
 
  The transition from the Fourier strategy to the Legendre one is quite simple:
  You just need to change the series expansion class and the corresponding
  smoothness estimation function to be part of the proper namespaces
  FESeries::Legendre and SmoothnessEstimator::Legendre. This strategy is used
  in @ref step_75 "step-75". For the theoretical background of this strategy, consult the
  general documentation of the SmoothnessEstimator::Legendre namespace, as well
  as @cite mavriplis1994hp , @cite eibner2007hp and @cite davydov2017hp .</li>
 
  <li><i>Refinement history:</i> The last strategy is quite different
  from the other two. In theory, we know how the error will converge
  after changing the discretization of the function space. With
  @f$h@f$-refinement the solution converges algebraically as already pointed
  out in @ref step_7 "step-7". If the solution is sufficiently smooth, though, we
  expect that the solution will converge exponentially with increasing
  polynomial degree of the finite element. We can compare a proper
  prediction of the error with the actual error in the following step to
  see if our choice of adaptation type was justified.
 
  The transition to this strategy is a bit more complicated. For this, we need
  an initialization step with pure @f$h@f$- or @f$p@f$-refinement and we need to
  transfer the predicted errors over adapted meshes. The extensive
  documentation of the hp::Refinement::predict_error() function describes not
  only the theoretical details of this approach, but also presents a blueprint
  on how to implement this strategy in your code. For more information, see
  @cite melenk2001hp .
 
  Note that with this particular function you cannot predict the error for
  the next time step in time-dependent problems. Therefore, this strategy
  cannot be applied to this type of problem without further ado. Alternatively,
  the following approach could be used, which works for all the other
  strategies as well: start each time step with a coarse mesh, keep refining
  until happy with the result, and only then move on to the next time step.</li>
</ul>
 
Try implementing one of these strategies into this tutorial and observe the
subtle changes to the results. You will notice that all strategies are
capable of identifying the singularities near the reentrant corners and
will perform @f$h@f$-refinement in these regions, while preferring @f$p@f$-refinement
in the bulk domain. A detailed comparison of these strategies is presented
in @cite fehling2020 .
 
 
<a name="Parallelhpadaptivefiniteelements"></a><h4>Parallel hp-adaptive finite elements</h4>
 
 
All functionality presented in this tutorial already works for both
sequential and parallel applications. It is possible without too much
effort to change to either the parallel::shared::Triangulation or the
parallel::distributed::Triangulation classes. If you feel eager to try
it, we recommend reading @ref step_18 "step-18" for the former and @ref step_40 "step-40" for the
latter case first for further background information on the topic, and
then come back to this tutorial to try out your newly acquired skills.
 
We go one step further in @ref step_75 "step-75": Here, we combine hp-adapative and
MatrixFree methods in combination with
parallel::distributed::Triangulation objects.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-27.cc"
*/