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 hp::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.
 *   

 * 
 * Finally, we initialize FESeries::Fourier object which will be used to
 * calculate coefficient in Fourier series as described in the introduction.
 * In addition to the hp::FECollection, we need to provide quadrature rules
 * hp::QCollection for integration on the reference cell.
 *   

 * 
 * In order to resize fourier_coefficients Table, we use the following
 * auxiliary function
 * 
 * @code
 *   template <int dim, typename T>
 *   void resize(Table<dim, T> &coeff, const unsigned int N)
 *   {
 *     TableIndices<dim> size;
 *     for (unsigned int d = 0; d < dim; d++)
 *       size[d] = N;
 *     coeff.reinit(size);
 *   }
 * 
 *   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
 * 
 * As described in the introduction, we define the Fourier vectors @f${\bf
 * k}@f$ for which we want to compute Fourier coefficients of the solution
 * on each cell as follows. In 2d, we will need coefficients corresponding
 * to vectors @f${\bf k}=(2 \pi i, 2\pi j)^T@f$ for which @f$\sqrt{i^2+j^2}\le N@f$,
 * with @f$i,j@f$ integers and @f$N@f$ being the maximal polynomial degree we use
 * for the finite elements in this program. The FESeries::Fourier class'
 * constructor first parameter @f$N@f$ defines the number of coefficients in 1D
 * with the total number of coefficients being @f$N^{dim}@f$. Although we will
 * not use coefficients corresponding to
 * @f$\sqrt{i^2+j^2}> N@f$ and @f$i+j==0@f$, the overhead of their calculation is
 * minimal. The transformation matrices for each FiniteElement will be
 * calculated only once the first time they are required in the course of
 * hp-adaptive refinement. Because we work on the unit cell, we can do all
 * this work without a mapping from reference to real cell and consequently
 * can precalculate these matrices. The calculation of expansion
 * coefficients for a particular set of local degrees of freedom on a given
 * cell then follows as a simple matrix-vector product.
 * The 3d case is handled analogously.
 * 
 * @code
 *     const unsigned int N = max_degree;
 * 
 * @endcode
 * 
 * We will need to assemble the matrices that do the Fourier transforms
 * for each of the finite elements we deal with, i.e. the matrices @f${\cal
 * F}_{{\bf k},j}@f$ defined in the introduction. We have to do that for
 * each of the finite elements in use. To that end we need a quadrature
 * rule. In this example we use the same quadrature formula for each
 * finite element, namely that is obtained by iterating a
 * 2-point Gauss formula as many times as the maximal exponent we use for
 * the term @f$e^{i{\bf k}\cdot{\bf x}}@f$:
 * 
 * @code
 *     QGauss<1>      base_quadrature(2);
 *     QIterated<dim> quadrature(base_quadrature, N);
 *     for (unsigned int i = 0; i < fe_collection.size(); i++)
 *       fourier_q_collection.push_back(quadrature);
 * 
 * @endcode
 * 
 * Now we are ready to set-up the FESeries::Fourier object
 * 
 * @code
 *     const std::vector<unsigned int> n_coefficients_per_direction(
 *       fe_collection.size(), N);
 *     fourier = std_cxx14::make_unique<FESeries::Fourier<dim>>(
 *       n_coefficients_per_direction, fe_collection, fourier_q_collection);
 * 
 * @endcode
 * 
 * We need to resize the matrix of fourier coefficients according to the
 * number of modes N.
 * 
 * @code
 *     resize(fourier_coefficients, N);
 *   }
 * 
 * 
 * @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 an hp object.
 * 
 * @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 implement 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().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. Estimating the smoothness is done in the respective
 * function of this class; that function is discussed further down below:
 * 
 * @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);
 * 
 * 
 *     Vector<float> smoothness_indicators(triangulation.n_active_cells());
 *     estimate_smoothness(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 it 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. Note that
 * the DataOut class has a second template argument (which defaults to
 * DoFHandler@<dim@>, which is why we have never seen it in previous
 * tutorial programs) that indicates the type of DoF handler to be
 * used. Here, we have to use the hp::DoFHandler class:
 * 
 * @code
 *       DataOut<dim, hp::DoFHandler<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
 * threshold. For this, we first have to determine the maximal and
 * minimal values of the smoothness indicators of all flagged cells,
 * which we do using a loop over all cells and comparing current minimal
 * and maximal values. (We start with the minimal and maximal values of
 * <i>all</i> cells, a range within which the minimal and maximal values
 * on cells flagged for refinement must surely lie.) Absent any better
 * strategies, we will then set the threshold above which will increase
 * @f$p@f$ instead of reducing @f$h@f$ as the mean value between minimal and
 * maximal smoothness indicators on cells flagged for refinement:
 * 
 * @code
 *       float max_smoothness = *std::min_element(smoothness_indicators.begin(),
 *                                                smoothness_indicators.end()),
 *             min_smoothness = *std::max_element(smoothness_indicators.begin(),
 *                                                smoothness_indicators.end());
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         if (cell->refine_flag_set())
 *           {
 *             max_smoothness =
 *               std::max(max_smoothness,
 *                        smoothness_indicators(cell->active_cell_index()));
 *             min_smoothness =
 *               std::min(min_smoothness,
 *                        smoothness_indicators(cell->active_cell_index()));
 *           }
 *       const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;
 * 
 * @endcode
 * 
 * With this, we can go back, loop over all cells again, and for those
 * cells 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 then increase the polynomial
 * degree and in return remove the flag indicating that the cell should
 * undergo bisection. For all other cells, the refinement flags remain
 * untouched:
 * 
 * @code
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         if (cell->refine_flag_set() &&
 *             (smoothness_indicators(cell->active_cell_index()) >
 *              threshold_smoothness) &&
 *             (cell->active_fe_index() + 1 < fe_collection.size()))
 *           {
 *             cell->clear_refine_flag();
 *             cell->set_active_fe_index(cell->active_fe_index() + 1);
 *           }
 * 
 * @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:
 * 
 * @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. It is a
 * specialization for the 2d case, i.e. a corresponding function needs to be
 * implemented if the program is run in anything other then 2d. The function
 * is actually stolen from @ref step_14 "step-14" and generates the same mesh used already
 * there, i.e. the square domain with the square hole in the middle. The
 * meaning of the different parts of this function are explained in the
 * documentation of @ref step_14 "step-14":
 * 
 * @code
 *   template <>
 *   void LaplaceProblem<2>::create_coarse_grid()
 *   {
 *     const unsigned int dim = 2;
 * 
 *     const std::vector<Point<2>> vertices = {
 *       {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0}, 
 *       {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5}, 
 *       {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0},               
 *       {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5}, 
 *       {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}};
 * 
 *     const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
 *       cell_vertices = {{{0, 1, 5, 6}},
 *                        {{1, 2, 6, 7}},
 *                        {{2, 3, 7, 8}},
 *                        {{3, 4, 8, 9}},
 *                        {{5, 6, 10, 11}},
 *                        {{8, 9, 12, 13}},
 *                        {{10, 11, 14, 15}},
 *                        {{12, 13, 17, 18}},
 *                        {{14, 15, 19, 20}},
 *                        {{15, 16, 20, 21}},
 *                        {{16, 17, 21, 22}},
 *                        {{17, 18, 22, 23}}};
 * 
 *     const unsigned int n_cells = cell_vertices.size();
 * 
 *     std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
 *     for (unsigned int i = 0; i < n_cells; ++i)
 *       {
 *         for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell; ++j)
 *           cells[i].vertices[j] = cell_vertices[i][j];
 *         cells[i].material_id = 0;
 *       }
 * 
 *     triangulation.create_triangulation(vertices, cells, SubCellData());
 *     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);
 *       }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemestimate_smoothness"></a> 
 * <h4>LaplaceProblem::estimate_smoothness</h4>
 * 

 * 
 * As described in the introduction, we will need to take the maximum
 * absolute value of fourier coefficients which correspond to @f$k@f$-vector
 * @f$|{\bf k}|= const@f$. To filter the coefficients Table we
 * will use the FESeries::process_coefficients() which requires a predicate
 * to be specified. The predicate should operate on TableIndices and return
 * a pair of <code>bool</code> and <code>unsigned int</code>. The latter
 * is the value of the map from TableIndicies to unsigned int.  It is
 * used to define subsets of coefficients from which we search for the one
 * with highest absolute value, i.e. @f$l^\infty@f$-norm. The <code>bool</code>
 * parameter defines which indices should be used in processing. In the
 * current case we are interested in coefficients which correspond to
 * @f$0 < i*i+j*j < N*N@f$ and @f$0 < i*i+j*j+k*k < N*N@f$ in 2D and 3D, respectively.
 * 
 * @code
 *   template <int dim>
 *   std::pair<bool, unsigned int>
 *   LaplaceProblem<dim>::predicate(const TableIndices<dim> &ind)
 *   {
 *     unsigned int v = 0;
 *     for (unsigned int i = 0; i < dim; i++)
 *       v += ind[i] * ind[i];
 *     if (v > 0 && v < max_degree * max_degree)
 *       return std::make_pair(true, v);
 *     else
 *       return std::make_pair(false, v);
 *   }
 * 
 * @endcode
 * 
 * This last function of significance implements the algorithm to estimate
 * the smoothness exponent using the algorithms explained in detail in the
 * introduction. We will therefore only comment on those points that are of
 * implementational importance.
 * 
 * @code
 *   template <int dim>
 *   void
 *   LaplaceProblem<dim>::estimate_smoothness(Vector<float> &smoothness_indicators)
 *   {
 * @endcode
 * 
 * Since most of the hard work is done for us in FESeries::Fourier and
 * we set up the object of this class in the constructor, what we are left
 * to do here is apply this class to calculate coefficients and then
 * perform linear regression to fit their decay slope.
 * 

 * 
 * 

 * 
 * First thing to do is to loop over all cells and do our work there, i.e.
 * to locally do the Fourier transform and estimate the decay coefficient.
 * We will use the following array as a scratch array in the loop to store
 * local DoF values:
 * 
 * @code
 *     Vector<double> local_dof_values;
 * 
 * @endcode
 * 
 * Then here is the loop:
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 * @endcode
 * 
 * Inside the loop, we first need to get the values of the local
 * degrees of freedom (which we put into the
 * <code>local_dof_values</code> array after setting it to the right
 * size) and then need to compute the Fourier transform by multiplying
 * this vector with the matrix @f${\cal F}@f$ corresponding to this finite
 * element. This is done by calling FESeries::Fourier::calculate(),
 * that has to be provided with the <code>local_dof_values</code>,
 * <code>cell->active_fe_index()</code> and a Table to store
 * coefficients.
 * 
 * @code
 *         local_dof_values.reinit(cell->get_fe().dofs_per_cell);
 *         cell->get_dof_values(solution, local_dof_values);
 * 
 *         fourier->calculate(local_dof_values,
 *                            cell->active_fe_index(),
 *                            fourier_coefficients);
 * 
 * @endcode
 * 
 * The next thing, as explained in the introduction, is that we wanted
 * to only fit our exponential decay of Fourier coefficients to the
 * largest coefficients for each possible value of @f$|{\bf k}|@f$. To
 * this end, we use FESeries::process_coefficients() to rework
 * coefficients into the desired format. We'll only take those Fourier
 * coefficients with the largest magnitude for a given value of @f$|{\bf
 * k}|@f$ and thereby need to use VectorTools::Linfty_norm:
 * 
 * @code
 *         std::pair<std::vector<unsigned int>, std::vector<double>> res =
 *           FESeries::process_coefficients<dim>(
 *             fourier_coefficients,
 *             [this](const TableIndices<dim> &indices) {
 *               return this->predicate(indices);
 *             },
 *             VectorTools::Linfty_norm);
 * 
 *         Assert(res.first.size() == res.second.size(), ExcInternalError());
 * 
 * @endcode
 * 
 * The first vector in the <code>std::pair</code> will store values of
 * the predicate, that is @f$i*i+j*j= const@f$ or @f$i*i+j*j+k*k = const@f$ in
 * 2D or 3D respectively. This vector will be the same for all the cells
 * so we can calculate logarithms of the corresponding Fourier vectors
 * @f$|{\bf k}|@f$ only once in the whole hp-refinement cycle:
 * 
 * @code
 *         if (ln_k.size() == 0)
 *           {
 *             ln_k.resize(res.first.size(), 0);
 *             for (unsigned int f = 0; f < ln_k.size(); f++)
 *               ln_k[f] =
 *                 std::log(2.0 * numbers::PI * std::sqrt(1. * res.first[f]));
 *           }
 * 
 * @endcode
 * 
 * We have to calculate the logarithms of absolute values of
 * coefficients and use it in a linear regression fit to obtain @f$\mu@f$.
 * 
 * @code
 *         for (double &residual_element : res.second)
 *           residual_element = std::log(residual_element);
 * 
 *         std::pair<double, double> fit =
 *           FESeries::linear_regression(ln_k, res.second);
 * 
 * @endcode
 * 
 * The final step is to compute the Sobolev index @f$s=\mu-\frac d2@f$ and
 * store it in the vector of estimated values for each cell:
 * 
 * @code
 *         smoothness_indicators(cell->active_cell_index()) =
 *           -fit.first - 1. * dim / 2;
 *       }
 *   }
 * } // 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 .

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      : 966
   Number of degrees of freedom: 5245
   Number of constraints       : 936
Cycle 2:
   Number of active cells      : 1143
   Number of degrees of freedom: 8441
   Number of constraints       : 1929
Cycle 3:
   Number of active cells      : 1356
   Number of degrees of freedom: 12349
   Number of constraints       : 3046
Cycle 4:
   Number of active cells      : 1644
   Number of degrees of freedom: 18178
   Number of constraints       : 4713
Cycle 5:
   Number of active cells      : 1728
   Number of degrees of freedom: 22591
   Number of constraints       : 6095
@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-0.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-1.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-2.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-3.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-4.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-5.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
grey corresponds to degree two and pink corresponds to degree seven:

<div class="threecolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-27-cell-degree-0.png"
         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-cell-degree-1.png"
         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-cell-degree-2.png"
         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-cell-degree-3.png"
         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-cell-degree-4.png"
         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-cell-degree-5.png"
         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-0.png"
         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-1.png"
         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-2.png"
         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-3.png"
         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-4.png"
         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-5.png"
         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="PlainProg"></a>
<h1> The plain program</h1>
@include "step-27.cc"
*/