step-6.h

Go to the documentation of this file.

 2)
 *     , dof_handler(triangulation)
 *   {}
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_6-Step6setup_system"></a> 
 * <h4>Step6::setup_system</h4>
 * 
 
 * 
 * The next function sets up all the variables that describe the linear
 * finite element problem, such as the DoFHandler, matrices, and
 * vectors. The difference to what we did in @ref step_5 "step-5" is only that we now also
 * have to take care of hanging node constraints. These constraints are
 * handled almost exclusively by the library, i.e. you only need to know
 * that they exist and how to get them, but you do not have to know how they
 * are formed or what exactly is done with them.
 * 
 
 * 
 * At the beginning of the function, you find all the things that are the same
 * as in @ref step_5 "step-5": setting up the degrees of freedom (this time we have
 * quadratic elements, but there is no difference from a user code perspective
 * to the linear -- or any other degree, for that matter -- case), generating
 * the sparsity pattern, and initializing the solution and right hand side
 * vectors. Note that the sparsity pattern will have significantly more
 * entries per row now, since there are now 9 degrees of freedom per cell
 * (rather than only four), that can couple with each other.
 * 
 * @code
 *   template <int dim>
 *   void Step6<dim>::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe);
 *   
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 *   
 * @endcode
 * 
 * We may now populate the AffineConstraints object with the hanging node
 * constraints. Since we will call this function in a loop we first clear
 * the current set of constraints from the last system and then compute new
 * ones:
 * 
 * @code
 *     constraints.clear();
 *     DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *   
 *   
 * @endcode
 * 
 * Now we are ready to interpolate the boundary values with indicator 0 (the
 * whole boundary) and store the resulting constraints in our
 * <code>constraints</code> object. Note that we do not to apply the
 * boundary conditions after assembly, like we did in earlier steps: instead
 * we put all constraints on our function space in the AffineConstraints
 * object. We can add constraints to the AffineConstraints object in either
 * order: if two constraints conflict then the constraint matrix either abort
 * or throw an exception via the Assert macro.
 * 
 * @code
 *     VectorTools::interpolate_boundary_values(dof_handler,
 *                                              types::boundary_id(0),
 *                                              Functions::ZeroFunction<dim>(),
 *                                              constraints);
 *   
 * @endcode
 * 
 * After all constraints have been added, they need to be sorted and
 * rearranged to perform some actions more efficiently. This postprocessing
 * is done using the <code>close()</code> function, after which no further
 * constraints may be added any more:
 * 
 * @code
 *     constraints.close();
 *   
 * @endcode
 * 
 * Now we first build our compressed sparsity pattern like we did in the
 * previous examples. Nevertheless, we do not copy it to the final sparsity
 * pattern immediately.  Note that we call a variant of
 * make_sparsity_pattern that takes the AffineConstraints object as the third
 * argument. We are letting the routine know that we will never write into
 * the locations given by <code>constraints</code> by setting the argument
 * <code>keep_constrained_dofs</code> to false (in other words, that we will
 * never write into entries of the matrix that correspond to constrained
 * degrees of freedom). If we were to condense the
 * constraints after assembling, we would have to pass <code>true</code>
 * instead because then we would first write into these locations only to
 * later set them to zero again during condensation.
 * 
 * @code
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler,
 *                                     dsp,
 *                                     constraints,
 *                                     /*keep_constrained_dofs = */ false);
 *   
 * @endcode
 * 
 * Now all non-zero entries of the matrix are known (i.e. those from
 * regularly assembling the matrix and those that were introduced by
 * eliminating constraints). We may copy our intermediate object to the
 * sparsity pattern:
 * 
 * @code
 *     sparsity_pattern.copy_from(dsp);
 *   
 * @endcode
 * 
 * We may now, finally, initialize the sparse matrix:
 * 
 * @code
 *     system_matrix.reinit(sparsity_pattern);
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_6-Step6assemble_system"></a> 
 * <h4>Step6::assemble_system</h4>
 * 
 
 * 
 * Next, we have to assemble the matrix. However, to copy the local matrix and
 * vector on each cell into the global system, we are no longer using a
 * hand-written loop. Instead, we use
 * AffineConstraints::distribute_local_to_global() that internally executes
 * this loop while performing Gaussian elimination on rows and columns
 * corresponding to constrained degrees on freedom.
 * 
 
 * 
 * The rest of the code that forms the local contributions remains
 * unchanged. It is worth noting, however, that under the hood several things
 * are different than before. First, the variable <code>dofs_per_cell</code>
 * and return value of <code>quadrature_formula.size()</code> now are 9 each,
 * where they were 4 before. Introducing such variables as abbreviations is a
 * good strategy to make code work with different elements without having to
 * change too much code. Secondly, the <code>fe_values</code> object of course
 * needs to do other things as well, since the shape functions are now
 * quadratic, rather than linear, in each coordinate variable. Again, however,
 * this is something that is completely handled by the library.
 * 
 * @code
 *   template <int dim>
 *   void Step6<dim>::assemble_system()
 *   {
 *     const QGauss<dim> quadrature_formula(fe.degree + 1);
 *   
 *     FEValues<dim> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               update_quadrature_points | update_JxW_values);
 *   
 *     const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *   
 *     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *     Vector<double>     cell_rhs(dofs_per_cell);
 *   
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         fe_values.reinit(cell);
 *   
 *         cell_matrix = 0;
 *         cell_rhs    = 0;
 *   
 *         for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *           {
 *             const double current_coefficient =
 *               coefficient(fe_values.quadrature_point(q_index));
 *             for (const unsigned int i : fe_values.dof_indices())
 *               {
 *                 for (const unsigned int j : fe_values.dof_indices())
 *                   cell_matrix(i, j) +=
 *                     (current_coefficient *              // a(x_q)
 *                      fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
 *                      fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
 *                      fe_values.JxW(q_index));           // dx
 *   
 *                 cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
 *                                 1.0 *                               // f(x)
 *                                 fe_values.JxW(q_index));            // dx
 *               }
 *           }
 *   
 * @endcode
 * 
 * Finally, transfer the contributions from @p cell_matrix and
 * @p cell_rhs into the global objects.
 * 
 * @code
 *         cell->get_dof_indices(local_dof_indices);
 *         constraints.distribute_local_to_global(
 *           cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
 *       }
 * @endcode
 * 
 * Now we are done assembling the linear system. The constraint matrix took
 * care of applying the boundary conditions and also eliminated hanging node
 * constraints. The constrained nodes are still in the linear system (there
 * is a nonzero entry, chosen in a way that the matrix is well conditioned,
 * on the diagonal of the matrix and all other entries for this line are set
 * to zero) but the computed values are invalid (i.e., the corresponding
 * entries in <code>system_rhs</code> are currently meaningless). We compute
 * the correct values for these nodes at the end of the <code>solve</code>
 * function.
 * 
 * @code
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_6-Step6solve"></a> 
 * <h4>Step6::solve</h4>
 * 
 
 * 
 * We continue with gradual improvements. The function that solves the linear
 * system again uses the SSOR preconditioner, and is again unchanged except
 * that we have to incorporate hanging node constraints. As mentioned above,
 * the degrees of freedom from the AffineConstraints object corresponding to
 * hanging node constraints and boundary values have been removed from the
 * linear system by giving the rows and columns of the matrix a special
 * treatment. This way, the values for these degrees of freedom have wrong,
 * but well-defined values after solving the linear system. What we then have
 * to do is to use the constraints to assign to them the values that they
 * should have. This process, called <code>distributing</code> constraints,
 * computes the values of constrained nodes from the values of the
 * unconstrained ones, and requires only a single additional function call
 * that you find at the end of this function:
 * 
 
 * 
 * 
 * @code
 *   template <int dim>
 *   void Step6<dim>::solve()
 *   {
 *     SolverControl            solver_control(1000, 1e-6 * system_rhs.l2_norm());
 *     SolverCG<Vector<double>> solver(solver_control);
 *   
 *     PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *     preconditioner.initialize(system_matrix, 1.2);
 *   
 *     solver.solve(system_matrix, solution, system_rhs, preconditioner);
 *   
 *     constraints.distribute(solution);
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_6-Step6refine_grid"></a> 
 * <h4>Step6::refine_grid</h4>
 * 
 
 * 
 * We use a sophisticated error estimation scheme to refine the mesh instead
 * of global refinement. We will use the KellyErrorEstimator class which
 * implements an error estimator for the Laplace equation; it can in principle
 * handle variable coefficients, but we will not use these advanced features,
 * but rather use its most simple form since we are not interested in
 * quantitative results but only in a quick way to generate locally refined
 * grids.
 * 
 
 * 
 * Although the error estimator derived by Kelly et al. was originally
 * developed for the Laplace equation, we have found that it is also well
 * suited to quickly generate locally refined grids for a wide class of
 * problems. This error estimator uses the solution gradient's jump at
 * cell faces (which is a measure for the second derivatives) and
 * scales it by the size of the cell. It is therefore a measure for the local
 * smoothness of the solution at the place of each cell and it is thus
 * understandable that it yields reasonable grids also for hyperbolic
 * transport problems or the wave equation as well, although these grids are
 * certainly suboptimal compared to approaches specially tailored to the
 * problem. This error estimator may therefore be understood as a quick way to
 * test an adaptive program.
 * 
 
 * 
 * The way the estimator works is to take a <code>DoFHandler</code> object
 * describing the degrees of freedom and a vector of values for each degree of
 * freedom as input and compute a single indicator value for each active cell
 * of the triangulation (i.e. one value for each of the active cells). To do
 * so, it needs two additional pieces of information: a face quadrature formula,
 * i.e., a quadrature formula on <code>dim-1</code> dimensional objects. We use
 * a 3-point Gauss rule again, a choice that is consistent and appropriate with
 * the bi-quadratic finite element shape functions in this program.
 * (What constitutes a suitable quadrature rule here of course depends on
 * knowledge of the way the error estimator evaluates the solution field. As
 * said above, the jump of the gradient is integrated over each face, which
 * would be a quadratic function on each face for the quadratic elements in
 * use in this example. In fact, however, it is the square of the jump of the
 * gradient, as explained in the documentation of that class, and that is a
 * quartic function, for which a 3 point Gauss formula is sufficient since it
 * integrates polynomials up to order 5 exactly.)
 * 
 
 * 
 * Secondly, the function wants a list of boundary indicators for those
 * boundaries where we have imposed Neumann values of the kind
 * @f$\partial_n u(\mathbf x) = h(\mathbf x)@f$, along with a function @f$h(\mathbf
 * x)@f$ for each such boundary. This information is represented by a map from
 * boundary indicators to function objects describing the Neumann boundary
 * values. In the present example program, we do not use Neumann boundary
 * values, so this map is empty, and in fact constructed using the default
 * constructor of the map in the place where the function call expects the
 * respective function argument.
 * 
 
 * 
 * The output is a vector of values for all active cells. While it may
 * make sense to compute the <b>value</b> of a solution degree of freedom
 * very accurately, it is usually not necessary to compute the <b>error
 * indicator</b> corresponding to the solution on a cell particularly
 * accurately. We therefore typically use a vector of floats instead of a vector
 * of doubles to represent error indicators.
 * 
 * @code
 *   template <int dim>
 *   void Step6<dim>::refine_grid()
 *   {
 *     Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   
 *     KellyErrorEstimator<dim>::estimate(dof_handler,
 *                                        QGauss<dim - 1>(fe.degree + 1),
 *                                        {},
 *                                        solution,
 *                                        estimated_error_per_cell);
 *   
 * @endcode
 * 
 * The above function returned one error indicator value for each cell in
 * the <code>estimated_error_per_cell</code> array. Refinement is now done
 * as follows: refine those 30 per cent of the cells with the highest error
 * values, and coarsen the 3 per cent of cells with the lowest values.
 *   
 
 * 
 * One can easily verify that if the second number were zero, this would
 * approximately result in a doubling of cells in each step in two space
 * dimensions, since for each of the 30 per cent of cells, four new would be
 * replaced, while the remaining 70 per cent of cells remain untouched. In
 * practice, some more cells are usually produced since it is disallowed
 * that a cell is refined twice while the neighbor cell is not refined; in
 * that case, the neighbor cell would be refined as well.
 *   
 
 * 
 * In many applications, the number of cells to be coarsened would be set to
 * something larger than only three per cent. A non-zero value is useful
 * especially if for some reason the initial (coarse) grid is already rather
 * refined. In that case, it might be necessary to refine it in some
 * regions, while coarsening in some other regions is useful. In our case
 * here, the initial grid is very coarse, so coarsening is only necessary in
 * a few regions where over-refinement may have taken place. Thus a small,
 * non-zero value is appropriate here.
 *   
 
 * 
 * The following function now takes these refinement indicators and flags
 * some cells of the triangulation for refinement or coarsening using the
 * method described above. It is from a class that implements several
 * different algorithms to refine a triangulation based on cell-wise error
 * indicators.
 * 
 * @code
 *     GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                     estimated_error_per_cell,
 *                                                     0.3,
 *                                                     0.03);
 *   
 * @endcode
 * 
 * After the previous function has exited, some cells are flagged for
 * refinement, and some other for coarsening. The refinement or coarsening
 * itself is not performed by now, however, since there are cases where
 * further modifications of these flags is useful. Here, we don't want to do
 * any such thing, so we can tell the triangulation to perform the actions
 * for which the cells are flagged:
 * 
 * @code
 *     triangulation.execute_coarsening_and_refinement();
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_6-Step6output_results"></a> 
 * <h4>Step6::output_results</h4>
 * 
 
 * 
 * At the end of computations on each grid, and just before we continue the
 * next cycle with mesh refinement, we want to output the results from this
 * cycle.
 * 
 
 * 
 * We have already seen in @ref step_1 "step-1" how this can be achieved for the
 * mesh itself. Here, we change a few things:
 * <ol>
 * <li>We use two different formats: gnuplot and VTU.</li>
 * <li>We embed the cycle number in the output file name.</li>
 * <li>For gnuplot output, we set up a GridOutFlags::Gnuplot object to
 * provide a few extra visualization arguments so that edges appear
 * curved. This is explained in further detail in @ref step_10 "step-10".</li>
 * </ol>
 * 
 * @code
 *   template <int dim>
 *   void Step6<dim>::output_results(const unsigned int cycle) const
 *   {
 *     {
 *       GridOut               grid_out;
 *       std::ofstream         output("grid-" + std::to_string(cycle) + ".gnuplot");
 *       GridOutFlags::Gnuplot gnuplot_flags(false, 5);
 *       grid_out.set_flags(gnuplot_flags);
 *       const MappingQ<dim> mapping(3);
 *       grid_out.write_gnuplot(triangulation, output, &mapping);
 *     }
 *   
 *     {
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "solution");
 *       data_out.build_patches();
 *   
 *       std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
 *       data_out.write_vtu(output);
 *     }
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_6-Step6run"></a> 
 * <h4>Step6::run</h4>
 * 
 
 * 
 * The final function before <code>main()</code> is again the main driver of
 * the class, <code>run()</code>. It is similar to the one of @ref step_5 "step-5", except
 * that we generate a file in the program again instead of reading it from
 * disk, in that we adaptively instead of globally refine the mesh, and that
 * we output the solution on the final mesh in the present function.
 * 
 
 * 
 * The first block in the main loop of the function deals with mesh generation.
 * If this is the first cycle of the program, instead of reading the grid from
 * a file on disk as in the previous example, we now again create it using a
 * library function. The domain is again a circle with center at the origin and
 * a radius of one (these are the two hidden arguments to the function, which
 * have default values).
 * 
 
 * 
 * You will notice by looking at the coarse grid that it is of inferior
 * quality than the one which we read from the file in the previous example:
 * the cells are less equally formed. However, using the library function this
 * program works in any space dimension, which was not the case before.
 * 
 
 * 
 * In case we find that this is not the first cycle, we want to refine the
 * grid. Unlike the global refinement employed in the last example program, we
 * now use the adaptive procedure described above.
 * 
 
 * 
 * The rest of the loop looks as before:
 * 
 * @code
 *   template <int dim>
 *   void Step6<dim>::run()
 *   {
 *     for (unsigned int cycle = 0; cycle < 8; ++cycle)
 *       {
 *         std::cout << "Cycle " << cycle << ':' << std::endl;
 *   
 *         if (cycle == 0)
 *           {
 *             GridGenerator::hyper_ball(triangulation);
 *             triangulation.refine_global(1);
 *           }
 *         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="step_6-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * The main function is unaltered in its functionality from the previous
 * example, but we have taken a step of additional caution. Sometimes,
 * something goes wrong (such as insufficient disk space upon writing an
 * output file, not enough memory when trying to allocate a vector or a
 * matrix, or if we can't read from or write to a file for whatever reason),
 * and in these cases the library will throw exceptions. Since these are
 * run-time problems, not programming errors that can be fixed once and for
 * all, this kind of exceptions is not switched off in optimized mode, in
 * contrast to the <code>Assert</code> macro which we have used to test
 * against programming errors. If uncaught, these exceptions propagate the
 * call tree up to the <code>main</code> function, and if they are not caught
 * there either, the program is aborted. In many cases, like if there is not
 * enough memory or disk space, we can't do anything but we can at least print
 * some text trying to explain the reason why the program failed. A way to do
 * so is shown in the following. It is certainly useful to write any larger
 * program in this way, and you can do so by more or less copying this
 * function except for the <code>try</code> block that actually encodes the
 * functionality particular to the present application.
 * 
 * @code
 *   int main()
 *   {
 * @endcode
 * 
 * The general idea behind the layout of this function is as follows: let's
 * try to run the program as we did before...
 * 
 * @code
 *     try
 *       {
 *         Step6<2> laplace_problem_2d;
 *         laplace_problem_2d.run();
 *       }
 * @endcode
 * 
 * ...and if this should fail, try to gather as much information as
 * possible. Specifically, if the exception that was thrown is an object of
 * a class that is derived from the C++ standard class
 * <code>exception</code>, then we can use the <code>what</code> member
 * function to get a string which describes the reason why the exception was
 * thrown.
 *   
 
 * 
 * The deal.II exception classes are all derived from the standard class,
 * and in particular, the <code>exc.what()</code> function will return
 * approximately the same string as would be generated if the exception was
 * thrown using the <code>Assert</code> macro. You have seen the output of
 * such an exception in the previous example, and you then know that it
 * contains the file and line number of where the exception occurred, and
 * some other information. This is also what the following statements would
 * print.
 *   
 
 * 
 * Apart from this, there isn't much that we can do except exiting the
 * program with an error code (this is what the <code>return 1;</code>
 * does):
 * 
 * @code
 *     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;
 *       }
 * @endcode
 * 
 * If the exception that was thrown somewhere was not an object of a class
 * derived from the standard <code>exception</code> class, then we can't do
 * anything at all. We then simply print an error message and exit.
 * 
 * @code
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 * @endcode
 * 
 * If we got to this point, there was no exception which propagated up to
 * the main function (there may have been exceptions, but they were caught
 * somewhere in the program or the library). Therefore, the program
 * performed as was expected and we can return without error.
 * 
 * @code
 *     return 0;
 *   }
 * @endcode
<a name="step_6-Results"></a><h1>Results</h1>
 
 
 
The output of the program looks as follows:
@code
Cycle 0:
   Number of active cells:       20
   Number of degrees of freedom: 89
Cycle 1:
   Number of active cells:       38
   Number of degrees of freedom: 183
Cycle 2:
   Number of active cells:       71
   Number of degrees of freedom: 325
Cycle 3:
   Number of active cells:       146
   Number of degrees of freedom: 689
Cycle 4:
   Number of active cells:       284
   Number of degrees of freedom: 1373
Cycle 5:
   Number of active cells:       599
   Number of degrees of freedom: 2769
Cycle 6:
   Number of active cells:       1241
   Number of degrees of freedom: 5833
Cycle 7:
   Number of active cells:       2507
   Number of degrees of freedom: 11916
@endcode
 
 
 
As intended, the number of cells roughly doubles in each cycle. The
number of degrees is slightly more than four times the number of
cells; one would expect a factor of exactly four in two spatial
dimensions on an infinite grid (since the spacing between the degrees
of freedom is half the cell width: one additional degree of freedom on
each edge and one in the middle of each cell), but it is larger than
that factor due to the finite size of the mesh and due to additional
degrees of freedom which are introduced by hanging nodes and local
refinement.
 
 
 
The program outputs the solution and mesh in each cycle of the
refinement loop. The solution looks as follows:
 
<img src="https://www.dealii.org/images/steps/developer/step-6.solution.9.2.png" alt="">
 
It is interesting to follow how the program arrives at the final mesh:
 
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_0.svg"
         alt="Initial grid: the five-cell circle grid with one global refinement."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_1.svg"
         alt="First grid: the five-cell circle grid with two global refinements."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_2.svg"
         alt="Second grid: the five-cell circle grid with one adaptive refinement."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_3.svg"
         alt="Third grid: the five-cell circle grid with two adaptive
         refinements, showing clustering around the inner circle."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_4.svg"
         alt="Fourth grid: the five-cell circle grid with three adaptive
         refinements, showing clustering around the inner circle."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_5.svg"
         alt="Fifth grid: the five-cell circle grid with four adaptive
         refinements, showing clustering around the inner circle."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_6.svg"
         alt="Sixth grid: the five-cell circle grid with five adaptive
         refinements, showing clustering around the inner circle."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_7.svg"
         alt="Last grid: the five-cell circle grid with six adaptive
         refinements, showing that most cells are clustered around the inner circle."
         width="300" height="300">
  </div>
</div>
 
 
It is clearly visible that the region where the solution has a kink,
i.e. the circle at radial distance 0.5 from the center, is
refined most. Furthermore, the central region where the solution is
very smooth and almost flat, is almost not refined at all, but this
results from the fact that we did not take into account that the
coefficient is large there. The region outside is refined rather
arbitrarily, since the second derivative is constant there and refinement
is therefore mostly based on the size of the cells and their deviation
from the optimal square.
 
 
 
<a name="step-6-extensions"></a>
<a name="step_6-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
<a name="step_6-Solversandpreconditioners"></a><h4>Solvers and preconditioners</h4>
 
 
One thing that is always worth playing around with if one solves
problems of appreciable size (much bigger than the one we have here)
is to try different solvers or preconditioners. In the current case,
the linear system is symmetric and positive definite, which makes the
Conjugate Gradient (CG) algorithm pretty much the canonical choice for solving. However,
the SSOR preconditioner we use in the <code>solve()</code> function is
up for grabs.
 
In deal.II, it is relatively simple to change the preconditioner.
Several simple preconditioner choices are accessible
by changing the following two lines
@code
  PreconditionSSOR<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(system_matrix, 1.2);
@endcode
of code in the program. For example, switching this into
@code
  PreconditionSSOR<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(system_matrix, 1.0);
@endcode
allows us to try out a different relaxation parameter for SSOR (1.0 instead
of the 1.2 in the original program). Similarly, by using
@code
  PreconditionJacobi<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(system_matrix);
@endcode
we can use Jacobi as a preconditioner. And by using
@code
  SparseILU<double> preconditioner;
  preconditioner.initialize(system_matrix);
@endcode
we can use a simple incomplete LU decomposition without any thresholding or
strengthening of the diagonal (to use this preconditioner, you have to also
add the header file <code>deal.II/lac/sparse_ilu.h</code> to the include list
at the top of the file).
 
Using these various different preconditioners, we can compare the
number of CG iterations needed (available through the
<code>solver_control.last_step()</code> call, see
@ref step_4 "step-4") as well as CPU time needed (using the Timer class,
discussed, for example, in @ref step_28 "step-28") and get the
following results (left: iterations; right: CPU time):
 
<table width="60%" align="center">
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_iterations.png" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-6.q2.dofs_vs_time.png" alt="">
    </td>
  </tr>
</table>
 
As we can see, all preconditioners behave pretty much the same on this
simple problem, with the number of iterations growing like @f${\cal
O}(N^{1/2})@f$ and because each iteration requires around @f${\cal
O}(N)@f$ operations the total CPU time grows like @f${\cal
O}(N^{3/2})@f$ (for the few smallest meshes, the CPU time is so small
that it doesn't record). Note that even though it is the simplest
method, Jacobi is the fastest for this problem.
 
The situation changes slightly when the finite element is not a
bi-quadratic one (i.e., polynomial degree two) as selected in the
constructor of this program, but a bi-linear one (polynomial degree one).
If one makes this change, the results are as follows:
 
<table width="60%" align="center">
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_iterations.png" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-6.q1.dofs_vs_time.png" alt="">
    </td>
  </tr>
</table>
 
In other words, while the increase in iterations and CPU time is as
before, Jacobi is now the method that requires the most iterations; it
is still the fastest one, however, owing to the simplicity of the
operations it has to perform. This is not to say that Jacobi
is actually a good preconditioner -- for problems of appreciable size, it is
definitely not, and other methods will be substantially better -- but really
only that it is fast because its implementation is so simple that it can
compensate for a larger number of iterations.
 
The message to take away from this is not that simplicity in
preconditioners is always best. While this may be true for the current
problem, it definitely is not once we move to more complicated
problems (elasticity or Stokes, for examples @ref step_8 "step-8" or
@ref step_22 "step-22"). Secondly, all of these preconditioners still
lead to an increase in the number of iterations as the number @f$N@f$ of
degrees of freedom grows, for example @f${\cal O}(N^\alpha)@f$; this, in
turn, leads to a total growth in effort as @f${\cal O}(N^{1+\alpha})@f$
since each iteration takes @f${\cal O}(N)@f$ work. This behavior is
undesirable: we would really like to solve linear systems with @f$N@f$
unknowns in a total of @f${\cal O}(N)@f$ work; there is a class
of preconditioners that can achieve this, namely geometric (@ref step_16 "step-16",
@ref step_37 "step-37", @ref step_39 "step-39")
or algebraic multigrid (@ref step_31 "step-31", @ref step_40 "step-40", and several others)
preconditioners. They are, however, significantly more complex than
the preconditioners outlined above, and so we will leave their use
to these later tutorial programs. The point to make, however, is
that "real" finite element programs do not use the preconditioners
we mention above: These are simply shown for expository purposes.
 
Finally, the last message to take
home is that when the data shown above was generated (in 2018), linear
systems with 100,000 unknowns are
easily solved on a desktop or laptop machine in about a second, making
the solution of relatively simple 2d problems even to very high
accuracy not that big a task as it used to be in the
past. At the same time, the situation for 3d problems continues to be
quite different: A uniform 2d mesh with 100,000 unknowns corresponds to
a grid with about @f$300 \times 300@f$ nodes; the corresponding 3d mesh has
@f$300 \times 300 \times 300@f$ nodes and 30 million unknowns. Because
finite element matrices in 3d have many more nonzero entries than in 2d,
solving these linear systems will not only take 300 times as much CPU
time, but substantially longer. In other words, achieving the
same resolution in 3d *is* quite a large problem, and solving it
within a reasonable amount of time *will* require much more work to
implement better linear solvers. As mentioned above, multigrid methods
and matrix-free methods (see, for example, @ref step_37 "step-37"), along with
parallelization (@ref step_40 "step-40") will be necessary, but are then also able
to comfortably solve such linear systems.
 
 
 
<a name="step_6-Abettermesh"></a><h4>A better mesh</h4>
 
 
If you look at the meshes above, you will see even though the domain is the
unit disk, and the jump in the coefficient lies along a circle, the cells
that make up the mesh do not track this geometry well. The reason, already hinted
at in @ref step_1 "step-1", is that in the absence of other information,
the Triangulation class only sees a bunch of
coarse grid cells but has, of course, no real idea what kind of geometry they
might represent when looked at together. For this reason, we need to tell
the Triangulation what to do when a cell is refined: where should the new
vertices at the edge midpoints and the cell midpoint be located so that the
child cells better represent the desired geometry than the parent cell.
 
To visualize what the triangulation actually knows about the geometry,
it is not enough to just output the location of vertices and draw a
straight line for each edge; instead, we have to output both interior
and boundary lines as multiple segments so that they look
curved. We can do this by making one change to the gnuplot part of
<code>output_results</code>:
@code
{
  GridOut       grid_out;
  std::ofstream output("grid-" + std::to_string(cycle) + ".gnuplot");
  GridOutFlags::Gnuplot gnuplot_flags(false, 5, /*curved_interior_cells*/true);
  grid_out.set_flags(gnuplot_flags);
  MappingQ<dim> mapping(3);
  grid_out.write_gnuplot(triangulation, output, &mapping);
}
@endcode
 
In the code above, we already do this for faces that sit at the boundary: this
happens automatically since we use GridGenerator::hyper_ball, which attaches a
SphericalManifold to the boundary of the domain. To make the mesh
<i>interior</i> also track a circular domain, we need to work a bit harder,
though. First, recall that our coarse mesh consists of a central square
cell and four cells around it. Now first consider what would happen if we
also attached the SphericalManifold object not only to the four exterior faces
but also the four cells at the perimeter as well as all of their faces. We can
do this by adding the following snippet (testing that the center of a cell is
larger than a small multiple, say one tenth, of the cell diameter away from
center of the mesh only fails for the central square of the mesh):
@code
GridGenerator::hyper_ball(triangulation);
// after GridGenerator::hyper_ball is called the Triangulation has
// a SphericalManifold with id 0. We can use it again on the interior.
const Point<dim> mesh_center;
for (const auto &cell : triangulation.active_cell_iterators())
  if (mesh_center.distance (cell->center()) > cell->diameter()/10)
    cell->set_all_manifold_ids(0);
 
triangulation.refine_global(1);
@endcode
 
After a few global refinement steps, this would lead to a mesh of the following
kind:
 
 
  <div class="onecolumn" style="width: 80%">
    <div>
      <img src="https://www.dealii.org/images/steps/developer/step_6_bad_grid_4.svg"
           alt="Grid where some central cells are nearly triangular."
           width="300" height="300">
    </div>
  </div>
 
This is not a good mesh: the central cell has been refined in such a way that
the children located in the four corners of the original central cell
<i>degenerate</i>: they all tend towards triangles as mesh refinement
continues. This means that the Jacobian matrix of the transformation from
reference cell to actual cell degenerates for these cells, and because
all error estimates for finite element solutions contain the norm of the
inverse of the Jacobian matrix, you will get very large errors on these
cells and, in the limit as mesh refinement, a loss of convergence order because
the cells in these corners become worse and worse under mesh refinement.
 
So we need something smarter. To this end, consider the following solution
originally developed by Konstantin Ladutenko. We will use the following code:
@code
GridGenerator::hyper_ball(triangulation);
 
const Point<dim> mesh_center;
const double core_radius  = 1.0/5.0,
             inner_radius = 1.0/3.0;
 
// Step 1: Shrink the inner cell
//
// We cannot get a circle out of the inner cell because of
// the degeneration problem mentioned above. Rather, shrink
// the inner cell to a core radius of 1/5 that stays
// sufficiently far away from the place where the
// coefficient will have a discontinuity and where we want
// to have cell interfaces that actually lie on a circle.
// We do this shrinking by just scaling the location of each
// of the vertices, given that the center of the circle is
// simply the origin of the coordinate system.
for (const auto &cell : triangulation.active_cell_iterators())
  if (mesh_center.distance(cell->center()) < 1e-5)
    {
      for (const auto v : cell->vertex_indices())
        cell->vertex(v) *= core_radius/mesh_center.distance(cell->vertex(v));
    }
 
// Step 2: Refine all cells except the central one
for (const auto &cell : triangulation.active_cell_iterators())
  if (mesh_center.distance(cell->center()) >= 1e-5)
    cell->set_refine_flag();
triangulation.execute_coarsening_and_refinement();
 
// Step 3: Resize the inner children of the outer cells
//
// The previous step replaced each of the four outer cells
// by its four children, but the radial distance at which we
// have intersected is not what we want to later refinement
// steps. Consequently, move the vertices that were just
// created in radial direction to a place where we need
// them.
for (const auto &cell : triangulation.active_cell_iterators())
  for (const auto v : cell->vertex_indices())
    {
      const double dist = mesh_center.distance(cell->vertex(v));
      if (dist > core_radius*1.0001 && dist < 0.9999)
        cell->vertex(v) *= inner_radius/dist;
    }
 
// Step 4: Apply curved manifold description
//
// As discussed above, we can not expect to subdivide the
// inner four cells (or their faces) onto concentric rings,
// but we can do so for all other cells that are located
// outside the inner radius. To this end, we loop over all
// cells and determine whether it is in this zone. If it
// isn't, then we set the manifold description of the cell
// and all of its bounding faces to the one that describes
// the spherical manifold already introduced above and that
// will be used for all further mesh refinement.
for (const auto &cell : triangulation.active_cell_iterators())
  {
    bool is_in_inner_circle = false;
    for (const auto v : cell->vertex_indices())
      if (mesh_center.distance(cell->vertex(v)) < inner_radius)
        {
          is_in_inner_circle = true;
          break;
        }
 
    if (is_in_inner_circle == false)
    // The Triangulation already has a SphericalManifold with
    // manifold id 0 (see the documentation of
    // GridGenerator::hyper_ball) so we just attach it to the outer
    // ring here:
      cell->set_all_manifold_ids(0);
  }
@endcode
 
This code then generates the following, much better sequence of meshes:
 
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_0_ladutenko.svg"
         alt="Initial grid: the Ladutenko grid with one global refinement."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_1_ladutenko.svg"
         alt="First adaptively refined Ladutenko grid."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_2_ladutenko.svg"
         alt="Second adaptively refined Ladutenko grid."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_3_ladutenko.svg"
         alt="Third adaptively refined Ladutenko grid."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_4_ladutenko.svg"
         alt="Fourth adaptively refined Ladutenko grid. The cells are clustered
         along the inner circle."
         width="300" height="300">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_6_grid_5_ladutenko.svg"
         alt="Fifth adaptively refined Ladutenko grid: the cells are clustered
         along the inner circle."
         width="300" height="300">
  </div>
</div>
 
Creating good meshes, and in particular making them fit the geometry you
want, is a complex topic in itself. You can find much more on this in
@ref step_49 "step-49", @ref step_53 "step-53", and @ref step_54 "step-54", among other tutorial programs that cover
the issue. @ref step_65 "step-65" shows another, less manual way to achieve a mesh
well fit to the problem here.
Information on curved domains can also be found in the
documentation topic on @ref manifold "Manifold descriptions".
 
Why does it make sense to choose a mesh that tracks the internal
interface? There are a number of reasons, but the most essential one
comes down to what we actually integrate in our bilinear
form. Conceptually, we want to integrate the term @f$A_{ij}^K=\int_K
a(\mathbf x) \nabla \varphi_i(\mathbf x) \nabla \varphi_j(\mathbf x) ; dx@f$ as the
contribution of cell @f$K@f$ to the matrix entry @f$A_{ij}@f$. We can not
compute it exactly and have to resort to quadrature. We know that
quadrature is accurate if the integrand is smooth. That is because
quadrature in essence computes a polynomial approximation to the
integrand that coincides with the integrand in the quadrature points,
and then computes the volume under this polynomial as an approximation
to the volume under the original integrand. This polynomial
interpolant is accurate if the integrand is smooth on a cell, but it
is usually rather inaccurate if the integrand is discontinuous on a
cell.
 
Consequently, it is worthwhile to align cells in such a way that the
interfaces across which the coefficient is discontinuous are aligned
with cell interfaces. This way, the coefficient is constant on each
cell, following which the integrand will be smooth, and its polynomial
approximation and the quadrature approximation of the integral will
both be accurate. Note that such an alignment is common in many
practical cases, so deal.II provides a number of functions (such as
@ref GlossMaterialId "material_id") to help manage such a scenario.
Refer to @ref step_28 "step-28" and @ref step_46 "step-46" for examples of how material ids can be
applied.
 
Finally, let us consider the case of a coefficient that has a smooth
and non-uniform distribution in space. We can repeat once again all of
the above discussion on the representation of such a function with the
quadrature. So, to simulate it accurately there are a few readily
available options: you could reduce the cell size, increase the order
of the polynomial used in the quadrature formula, select a more
appropriate quadrature formula, or perform a combination of these
steps. The key is that providing the best fit of the coefficient's
spatial dependence with the quadrature polynomial will lead to a more
accurate finite element solution of the PDE.
 
As a final note: The discussion in the previous paragraphs shows, we here
have a very concrete way of stating what we think of a good mesh -- it should
be aligned with the jump in the coefficient. But one could also have asked
this kind of question in a more general setting: Given some equation with
a smooth solution and smooth coefficients, can we say what a good mesh
would look like? This is a question for which the answer is easier to state
in intuitive terms than mathematically: A good mesh has cells that all,
by and large, look like squares (or cubes, in 3d). A bad mesh would contain
cells that are very elongated in some directions or, more generally, for which
there are cells that have both short and long edges. There are many ways
in which one could assign a numerical quality index to each cell that measures
whether the cell is "good" or "bad"; some of these are often chosen because
they are cheap and easy to compute, whereas others are based on what enters
into proofs of convergence. An example of the former would be the ratio of
the longest to the shortest edge of a cell: In the ideal case, that ratio
would be one; bad cells have values much larger than one. An example of the
latter kind would consider the gradient (the "Jacobian") of the mapping
from the reference cell @f$\hat K=[0,1]^d@f$ to the real cell @f$K@f$; this
gradient is a matrix, and a quantity that enters into error estimates
is the maximum over all points on the reference cell of the ratio of the
largest to the smallest eigenvalue of this matrix. It is again not difficult
to see that this ratio is constant if the cell @f$K@f$ is an affine image of
@f$\hat K@f$, and that it is one for squares and cubes.
 
In practice, it might be interesting to visualize such quality measures.
The function GridTools::compute_aspect_ratio_of_cells() provides one
way to get this kind of information. Even better, visualization tools
such as VisIt often allow you to visualize this sort of information
for a variety of measures from within the visualization software
itself; in the case of VisIt, just add a "pseudo-color" plot and select
one of the mesh quality measures instead of the solution field.
 
 
<a name="step_6-Playingwiththeregularityofthesolution"></a><h4>Playing with the regularity of the solution</h4>
 
 
From a mathematical perspective, solutions of the Laplace equation
@f[
  -\Delta u = f
@f]
on smoothly bounded, convex domains are known to be smooth themselves. The exact degree
of smoothness, i.e., the function space in which the solution lives, depends
on how smooth exactly the boundary of the domain is, and how smooth the right
hand side is. Some regularity of the solution may be lost at the boundary, but
one generally has that the solution is twice more differentiable in
compact subsets of the domain than the right hand side.
If, in particular, the right hand side satisfies @f$f\in C^\infty(\Omega)@f$, then
@f$u \in C^\infty(\Omega_i)@f$ where @f$\Omega_i@f$ is any compact subset of @f$\Omega@f$
(@f$\Omega@f$ is an open domain, so a compact subset needs to keep a positive distance
from @f$\partial\Omega@f$).
 
The situation we chose for the current example is different, however: we look
at an equation with a non-constant coefficient @f$a(\mathbf x)@f$:
@f[
  -\nabla \cdot (a \nabla u) = f.
@f]
Here, if @f$a@f$ is not smooth, then the solution will not be smooth either,
regardless of @f$f@f$. In particular, we expect that wherever @f$a@f$ is discontinuous
along a line (or along a plane in 3d),
the solution will have a kink. This is easy to see: if for example @f$f@f$
is continuous, then @f$f=-\nabla \cdot (a \nabla u)@f$ needs to be
continuous. This means that @f$a \nabla u@f$ must be continuously differentiable
(not have a kink). Consequently, if @f$a@f$ has a discontinuity, then @f$\nabla u@f$
must have an opposite discontinuity so that the two exactly cancel and their
product yields a function without a discontinuity. But for @f$\nabla u@f$ to have
a discontinuity, @f$u@f$ must have a kink. This is of course exactly what is
happening in the current example, and easy to observe in the pictures of the
solution.
 
In general, if the coefficient @f$a(\mathbf x)@f$ is discontinuous along a line in 2d,
or a plane in 3d, then the solution may have a kink, but the gradient of the
solution will not go to infinity. That means, that the solution is at least
still in the <a href="https://en.wikipedia.org/wiki/Sobolev_space">Sobolev space</a>
@f$W^{1,\infty}@f$ (i.e., roughly speaking, in the
space of functions whose derivatives are bounded). On the other hand,
we know that in the most
extreme cases -- i.e., where the domain has reentrant corners, the
right hand side only satisfies @f$f\in H^{-1}@f$, or the coefficient @f$a@f$ is only in
@f$L^\infty@f$ -- all we can expect is that @f$u\in H^1@f$ (i.e., the
<a
href="https://en.wikipedia.org/wiki/Sobolev_space#Sobolev_spaces_with_integer_k">Sobolev
space</a> of functions whose derivative is square integrable), a much larger space than
@f$W^{1,\infty}@f$. It is not very difficult to create cases where
the solution is in a space @f$H^{1+s}@f$ where we can get @f$s@f$ to become as small
as we want. Such cases are often used to test adaptive finite element
methods because the mesh will have to resolve the singularity that causes
the solution to not be in @f$W^{1,\infty}@f$ any more.
 
The typical example one uses for this is called the <i>Kellogg problem</i>
(referring to @cite Kel74), which in the commonly used form has a coefficient
@f$a(\mathbf x)@f$ that has different values in the four quadrants of the plane
(or eight different values in the octants of @f${\mathbb R}^3@f$). The exact degree
of regularity (the @f$s@f$ in the index of the Sobolev space above) depends on the
values of @f$a(\mathbf x)@f$ coming together at the origin, and by choosing the
jumps large enough, the regularity of the solution can be made as close as
desired to @f$H^1@f$.
 
To implement something like this, one could replace the coefficient
function by the following (shown here only for the 2d case):
@code
template <int dim>
double coefficient (const Point<dim> &p)
{
  if ((p[0] < 0) && (p[1] < 0))           // lower left quadrant
    return 1;
  else if ((p[0] >= 0) && (p[1] < 0))     // lower right quadrant
    return 10;
  else if ((p[0] < 0) && (p[1] >= 0))     // upper left quadrant
    return 100;
  else if ((p[0] >= 0) && (p[1] >= 0))    // upper right quadrant
    return 1000;
  else
    DEAL_II_ASSERT_UNREACHABLE();
}
@endcode
(Adding the `DEAL_II_ASSERT_UNREACHABLE();` call at the end ensures that either an exception
is thrown or that the program aborts if we ever get to that point
-- which of course we shouldn't,
but this is a good way to insure yourself: We all make mistakes by
sometimes not thinking of all cases, for example by checking
for <code>p[0]</code> to be less than and greater than zero,
rather than greater-or-equal to zero, and thereby forgetting
some cases that would otherwise lead to bugs that are awkward
to find.)
 
By playing with such cases where four or more sectors come
together and on which the coefficient has different values, one can
construct cases where the solution has singularities at the
origin. One can also see how the meshes are refined in such cases.
 *
 *
<a name="step_6-PlainProg"></a>
<h1> The plain program</h1>
@include "step-6.cc"
*/