step-4.h

Go to the documentation of this file.

) const
 *   {
 *     double return_value = 0.0;
 *     for (unsigned int i = 0; i < dim; ++i)
 *       return_value += 4.0 * std::pow(p[i], 4.0);
 *   
 *     return return_value;
 *   }
 *   
 *   
 * @endcode
 * 
 * As boundary values, we choose @f$x^2+y^2@f$ in 2d, and @f$x^2+y^2+z^2@f$ in 3d. This
 * happens to be equal to the square of the vector from the origin to the
 * point at which we would like to evaluate the function, irrespective of the
 * dimension. So that is what we return:
 * 
 * @code
 *   template <int dim>
 *   double BoundaryValues<dim>::value(const Point<dim> &p,
 *                                     const unsigned int /*component*/) const
 *   {
 *     return p.square();
 *   }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-ImplementationofthecodeStep4codeclass"></a> 
 * <h3>Implementation of the <code>Step4</code> class</h3>
 * 
 
 * 
 * Next for the implementation of the class template that makes use of the
 * functions above. As before, we will write everything as templates that have
 * a formal parameter <code>dim</code> that we assume unknown at the time we
 * define the template functions. Only later, the compiler will find a
 * declaration of <code>Step4@<2@></code> (in the <code>main</code> function,
 * actually) and compile the entire class with <code>dim</code> replaced by 2,
 * a process referred to as `instantiation of a template'. When doing so, it
 * will also replace instances of <code>RightHandSide@<dim@></code> by
 * <code>RightHandSide@<2@></code> and instantiate the latter class from the
 * class template.
 * 
 
 * 
 * In fact, the compiler will also find a declaration <code>Step4@<3@></code>
 * in <code>main()</code>. This will cause it to again go back to the general
 * <code>Step4@<dim@></code> template, replace all occurrences of
 * <code>dim</code>, this time by 3, and compile the class a second time. Note
 * that the two instantiations <code>Step4@<2@></code> and
 * <code>Step4@<3@></code> are completely independent classes; their only
 * common feature is that they are both instantiated from the same general
 * template, but they are not convertible into each other, for example, and
 * share no code (both instantiations are compiled completely independently).
 * 
 
 * 
 * 
 
 * 
 * 
 * <a name="step_4-Step4Step4"></a> 
 * <h4>Step4::Step4</h4>
 * 
 
 * 
 * After this introduction, here is the constructor of the <code>Step4</code>
 * class. It specifies the desired polynomial degree of the finite elements
 * and associates the DoFHandler to the triangulation just as in the previous
 * example program, @ref step_3 "step-3":
 * 
 * @code
 *   template <int dim>
 *   Step4<dim>::Step4()
 *     : fe(/* polynomial degree = */ 1)
 *     , dof_handler(triangulation)
 *   {}
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Step4make_grid"></a> 
 * <h4>Step4::make_grid</h4>
 * 
 
 * 
 * Grid creation is something inherently dimension dependent. However, as long
 * as the domains are sufficiently similar in 2d or 3d, the library can
 * abstract for you. In our case, we would like to again solve on the square
 * @f$[-1,1]\times [-1,1]@f$ in 2d, or on the cube @f$[-1,1] \times [-1,1] \times
 * [-1,1]@f$ in 3d; both can be termed GridGenerator::hyper_cube(), so we may
 * use the same function in whatever dimension we are. Of course, the
 * functions that create a hypercube in two and three dimensions are very much
 * different, but that is something you need not care about. Let the library
 * handle the difficult things.
 * 
 * @code
 *   template <int dim>
 *   void Step4<dim>::make_grid()
 *   {
 *     GridGenerator::hyper_cube(triangulation, -1, 1);
 *     triangulation.refine_global(4);
 *   
 *     std::cout << "   Number of active cells: " << triangulation.n_active_cells()
 *               << std::endl
 *               << "   Total number of cells: " << triangulation.n_cells()
 *               << std::endl;
 *   }
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Step4setup_system"></a> 
 * <h4>Step4::setup_system</h4>
 * 
 
 * 
 * This function looks exactly like in the previous example, although it
 * performs actions that in their details are quite different if
 * <code>dim</code> happens to be 3. The only significant difference from a
 * user's perspective is the number of cells resulting, which is much higher
 * in three than in two space dimensions!
 * 
 * @code
 *   template <int dim>
 *   void Step4<dim>::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe);
 *   
 *     std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *               << std::endl;
 *   
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *     sparsity_pattern.copy_from(dsp);
 *   
 *     system_matrix.reinit(sparsity_pattern);
 *   
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Step4assemble_system"></a> 
 * <h4>Step4::assemble_system</h4>
 * 
 
 * 
 * Unlike in the previous example, we would now like to use a non-constant
 * right hand side function and non-zero boundary values. Both are tasks that
 * are readily achieved with only a few new lines of code in the assemblage of
 * the matrix and right hand side.
 * 
 
 * 
 * More interesting, though, is the way we assemble matrix and right hand side
 * vector dimension independently: there is simply no difference to the
 * two-dimensional case. Since the important objects used in this function
 * (quadrature formula, FEValues) depend on the dimension by way of a template
 * parameter as well, they can take care of setting up properly everything for
 * the dimension for which this function is compiled. By declaring all classes
 * which might depend on the dimension using a template parameter, the library
 * can make nearly all work for you and you don't have to care about most
 * things.
 * 
 * @code
 *   template <int dim>
 *   void Step4<dim>::assemble_system()
 *   {
 *     const QGauss<dim> quadrature_formula(fe.degree + 1);
 *   
 * @endcode
 * 
 * We wanted to have a non-constant right hand side, so we use an object of
 * the class declared above to generate the necessary data. Since this right
 * hand side object is only used locally in the present function, we declare
 * it here as a local variable:
 * 
 * @code
 *     RightHandSide<dim> right_hand_side;
 *   
 * @endcode
 * 
 * Compared to the previous example, in order to evaluate the non-constant
 * right hand side function we now also need the quadrature points on the
 * cell we are presently on (previously, we only required values and
 * gradients of the shape function from the FEValues object, as well as the
 * quadrature weights, FEValues::JxW() ). We can tell the FEValues object to
 * do for us by also giving it the #update_quadrature_points flag:
 * 
 * @code
 *     FEValues<dim> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               update_quadrature_points | update_JxW_values);
 *   
 * @endcode
 * 
 * We then again define the same abbreviation as in the previous program.
 * The value of this variable of course depends on the dimension which we
 * are presently using, but the FiniteElement class does all the necessary
 * work for you and you don't have to care about the dimension dependent
 * parts:
 * 
 * @code
 *     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);
 *   
 * @endcode
 * 
 * Next, we again have to loop over all cells and assemble local
 * contributions.  Note, that a cell is a quadrilateral in two space
 * dimensions, but a hexahedron in 3d. In fact, the
 * <code>active_cell_iterator</code> data type is something different,
 * depending on the dimension we are in, but to the outside world they look
 * alike and you will probably never see a difference. In any case, the real
 * type is hidden by using `auto`:
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         fe_values.reinit(cell);
 *   
 *         cell_matrix = 0;
 *         cell_rhs    = 0;
 *   
 * @endcode
 * 
 * Now we have to assemble the local matrix and right hand side. This is
 * done exactly like in the previous example, but now we revert the
 * order of the loops (which we can safely do since they are independent
 * of each other) and merge the loops for the local matrix and the local
 * vector as far as possible to make things a bit faster.
 *       
 
 * 
 * Assembling the right hand side presents the only significant
 * difference to how we did things in @ref step_3 "step-3": Instead of using a
 * constant right hand side with value 1, we use the object representing
 * the right hand side and evaluate it at the quadrature points:
 * 
 * @code
 *         for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *           for (const unsigned int i : fe_values.dof_indices())
 *             {
 *               for (const unsigned int j : fe_values.dof_indices())
 *                 cell_matrix(i, j) +=
 *                   (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
 *   
 *               const auto &x_q = fe_values.quadrature_point(q_index);
 *               cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
 *                               right_hand_side.value(x_q) *        // f(x_q)
 *                               fe_values.JxW(q_index));            // dx
 *             }
 * @endcode
 * 
 * As a final remark to these loops: when we assemble the local
 * contributions into <code>cell_matrix(i,j)</code>, we have to multiply
 * the gradients of shape functions @f$i@f$ and @f$j@f$ at point number
 * q_index and
 * multiply it with the scalar weights JxW. This is what actually
 * happens: <code>fe_values.shape_grad(i,q_index)</code> returns a
 * <code>dim</code> dimensional vector, represented by a
 * <code>Tensor@<1,dim@></code> object, and the operator* that
 * multiplies it with the result of
 * <code>fe_values.shape_grad(j,q_index)</code> makes sure that the
 * <code>dim</code> components of the two vectors are properly
 * contracted, and the result is a scalar floating point number that
 * then is multiplied with the weights. Internally, this operator* makes
 * sure that this happens correctly for all <code>dim</code> components
 * of the vectors, whether <code>dim</code> be 2, 3, or any other space
 * dimension; from a user's perspective, this is not something worth
 * bothering with, however, making things a lot simpler if one wants to
 * write code dimension independently.
 * 
 
 * 
 * With the local systems assembled, the transfer into the global matrix
 * and right hand side is done exactly as before, but here we have again
 * merged some loops for efficiency:
 * 
 * @code
 *         cell->get_dof_indices(local_dof_indices);
 *         for (const unsigned int i : fe_values.dof_indices())
 *           {
 *             for (const unsigned int j : fe_values.dof_indices())
 *               system_matrix.add(local_dof_indices[i],
 *                                 local_dof_indices[j],
 *                                 cell_matrix(i, j));
 *   
 *             system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *           }
 *       }
 *   
 * @endcode
 * 
 * As the final step in this function, we wanted to have non-homogeneous
 * boundary values in this example, unlike the one before. This is a simple
 * task, we only have to replace the Functions::ZeroFunction used there by an
 * object of the class which describes the boundary values we would like to
 * use (i.e. the <code>BoundaryValues</code> class declared above):
 *   
 
 * 
 * The function VectorTools::interpolate_boundary_values() will only work
 * on faces that have been marked with boundary indicator 0 (because that's
 * what we say the function should work on with the second argument below).
 * If there are faces with boundary id other than 0, then the function
 * interpolate_boundary_values() will do nothing on these faces. For
 * the Laplace equation doing nothing is equivalent to assuming that
 * on those parts of the boundary a zero Neumann boundary condition holds.
 * 
 * @code
 *     std::map<types::global_dof_index, double> boundary_values;
 *     VectorTools::interpolate_boundary_values(dof_handler,
 *                                              types::boundary_id(0),
 *                                              BoundaryValues<dim>(),
 *                                              boundary_values);
 *     MatrixTools::apply_boundary_values(boundary_values,
 *                                        system_matrix,
 *                                        solution,
 *                                        system_rhs);
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Step4solve"></a> 
 * <h4>Step4::solve</h4>
 * 
 
 * 
 * Solving the linear system of equations is something that looks almost
 * identical in most programs. In particular, it is dimension independent, so
 * this function is copied verbatim from the previous example.
 * 
 * @code
 *   template <int dim>
 *   void Step4<dim>::solve()
 *   {
 *     SolverControl            solver_control(1000, 1e-6 * system_rhs.l2_norm());
 *     SolverCG<Vector<double>> solver(solver_control);
 *     solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
 *   
 *     std::cout << "   " << solver_control.last_step()
 *               << " CG iterations needed to obtain convergence." << std::endl;
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Step4output_results"></a> 
 * <h4>Step4::output_results</h4>
 * 
 
 * 
 * This function also does what the respective one did in @ref step_3 "step-3". No changes
 * here for dimension independence either.
 * 
 
 * 
 * Since the program will run both 2d and 3d versions of the Laplace solver,
 * we use the dimension in the filename to generate distinct filenames for
 * each run (in a better program, one would check whether <code>dim</code> can
 * have other values than 2 or 3, but we neglect this here for the sake of
 * brevity).
 * 
 * @code
 *   template <int dim>
 *   void Step4<dim>::output_results() const
 *   {
 *     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(dim == 2 ? "solution-2d.vtk" : "solution-3d.vtk");
 *     data_out.write_vtk(output);
 *   }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Step4run"></a> 
 * <h4>Step4::run</h4>
 * 
 
 * 
 * This is the function which has the top-level control over everything. Apart
 * from one line of additional output, it is the same as for the previous
 * example.
 * 
 * @code
 *   template <int dim>
 *   void Step4<dim>::run()
 *   {
 *     std::cout << "Solving problem in " << dim << " space dimensions."
 *               << std::endl;
 *   
 *     make_grid();
 *     setup_system();
 *     assemble_system();
 *     solve();
 *     output_results();
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_4-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * And this is the main function. It also looks mostly like in @ref step_3 "step-3", but if
 * you look at the code below, note how we first create a variable of type
 * <code>Step4@<2@></code> (forcing the compiler to compile the class template
 * with <code>dim</code> replaced by <code>2</code>) and run a 2d simulation,
 * and then we do the whole thing over in 3d.
 * 
 
 * 
 * In practice, this is probably not what you would do very frequently (you
 * probably either want to solve a 2d problem, or one in 3d, but not both at
 * the same time). However, it demonstrates the mechanism by which we can
 * simply change which dimension we want in a single place, and thereby force
 * the compiler to recompile the dimension independent class templates for the
 * dimension we request. The emphasis here lies on the fact that we only need
 * to change a single place. This makes it rather trivial to debug the program
 * in 2d where computations are fast, and then switch a single place to a 3 to
 * run the much more computing intensive program in 3d for `real'
 * computations.
 * 
 
 * 
 * Each of the two blocks is enclosed in braces to make sure that the
 * <code>laplace_problem_2d</code> variable goes out of scope (and releases
 * the memory it holds) before we move on to allocate memory for the 3d
 * case. Without the additional braces, the <code>laplace_problem_2d</code>
 * variable would only be destroyed at the end of the function, i.e. after
 * running the 3d problem, and would needlessly hog memory while the 3d run
 * could actually use it.
 * 
 * @code
 *   int main()
 *   {
 *     {
 *       Step4<2> laplace_problem_2d;
 *       laplace_problem_2d.run();
 *     }
 *   
 *     {
 *       Step4<3> laplace_problem_3d;
 *       laplace_problem_3d.run();
 *     }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="step_4-Results"></a><h1>Results</h1>
 
 
 
The output of the program looks as follows (the number of iterations
may vary by one or two, depending on your computer, since this is
often dependent on the round-off accuracy of floating point
operations, which differs between processors):
@code
Solving problem in 2 space dimensions.
   Number of active cells: 256
   Total number of cells: 341
   Number of degrees of freedom: 289
   19 CG iterations needed to obtain convergence.
Solving problem in 3 space dimensions.
   Number of active cells: 4096
   Total number of cells: 4681
   Number of degrees of freedom: 4913
   20 CG iterations needed to obtain convergence.
@endcode
It is obvious that in three spatial dimensions the number of cells and
therefore also the number of degrees of freedom is
much higher. What cannot be seen here, is that besides this higher
number of rows and columns in the matrix, there are also significantly
more entries per row of the matrix in three space
dimensions. Together, this leads to a much higher numerical effort for
solving the system of equation, which you can feel in the run time of the two
solution steps when you actually run the program.
 
 
 
The program produces two files: <code>solution-2d.vtk</code> and
<code>solution-3d.vtk</code>, which can be viewed using the programs
VisIt or Paraview (in case you do not have these programs, you can easily
change the
output format in the program to something which you can view more
easily). Visualizing solutions is a bit of an art, but it can also be fun, so
you should play around with your favorite visualization tool to get familiar
with its functionality. Here's what I have come up with for the 2d solution:
 
<p align="center">
  <img src="https://www.dealii.org/images/steps/developer/step-4.solution-2d.png" alt="">
</p>
 
(See also <a href="https://www.math.colostate.edu/~bangerth/videos.676.11.html">video lecture 11</a>, <a href="https://www.math.colostate.edu/~bangerth/videos.676.32.html">video lecture 32</a>.)
The picture shows the solution of the problem under consideration as
a 3D plot. As can be seen, the solution is almost flat in the interior
of the domain and has a higher curvature near the boundary. This, of
course, is due to the fact that for Laplace's equation the curvature
of the solution is equal to the right hand side and that was chosen as
a quartic polynomial which is nearly zero in the interior and is only
rising sharply when approaching the boundaries of the domain; the
maximal values of the right hand side function are at the corners of
the domain, where also the solution is moving most rapidly.
It is also nice to see that the solution follows the desired quadratic
boundary values along the boundaries of the domain.
It can also be useful to verify a computed solution against an analytical
solution. For an explanation of this technique, see @ref step_7 "step-7".
 
On the other hand, even though the picture does not show the mesh lines
explicitly, you can see them as little kinks in the solution. This clearly
indicates that the solution hasn't been computed to very high accuracy and
that to get a better solution, we may have to compute on a finer mesh.
 
In three spatial dimensions, visualization is a bit more difficult. The left
picture shows the solution and the mesh it was computed on on the surface of
the domain. This is nice, but it has the drawback that it completely hides
what is happening on the inside. The picture on the right is an attempt at
visualizing the interior as well, by showing surfaces where the solution has
constant values (as indicated by the legend at the top left). Isosurface
pictures look best if one makes the individual surfaces slightly transparent
so that it is possible to see through them and see what's behind.
 
<table width="60%" align="center">
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-4.solution-3d.png" alt="">
    </td>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-4.contours-3d.png" alt="">
    </td>
  </tr>
</table>
 
@note
A final remark on visualization: the idea of visualization is to give insight,
which is not the same as displaying information. In particular, it is easy to
overload a picture with information, but while it shows more information it
makes it also more difficult to glean insight. As an example, the program I
used to generate these pictures, VisIt, by default puts tick marks on every
axis, puts a big fat label "X Axis" on the @f$x@f$ axis and similar for the other
axes, shows the file name from which the data was taken in the top left and
the name of the user doing so and the time and date on the bottom right. None
of this is important
here: the axes are equally easy to make out because the tripod at the bottom
left is still visible, and we know from the program that the domain is
@f$[-1,1]^3@f$, so there is no need for tick marks. As a consequence, I have
switched off all the extraneous stuff in the picture: the art of visualization
is to reduce the picture to those parts that are important to see what one
wants to see, but no more.
 
 
 
<a name="step-4-postprocessing"></a>
<a name="step_4-PostprocessingWhattodowiththesolution"></a><h3> Postprocessing: What to do with the solution? </h3>
 
 
This tutorial -- like most of the other programs -- principally only shows how
to numerically approximate the solution of a partial differential equation, and
then how to visualize this solution graphically. But
solving a PDE is of course not the goal in most practical applications (unless
you are a numerical methods developer and the *method* is the goal): We generally
want to solve a PDE because we want to *extract information* from it. Examples
for what people are interested in from solutions include the following:
- Let's say you solve the equations of elasticity (which we will do in @ref step_8 "step-8"),
  then that's presumably because you want to know about the deformation of an
  elastic object under a given load. From an engineering perspective, what you
  then presumably want to learn is the degree of deformation of the object,
  say at a specific point; or you may want to know the maximum
  [stress](https://en.wikipedia.org/wiki/Stress_(mechanics)) in order to
  determine whether the applied load exceeds the safe maximal stress the
  material can withstand.
- If you are solving fluid flow problems (such as in @ref step_22 "step-22", @ref step_57 "step-57", @ref step_67 "step-67",
  and @ref step_69 "step-69"), then you might be interested in the fluid velocity at specific
  points, and oftentimes the forces the fluid exerts on the boundary of the
  fluid domain. The latter is important in many applications: If the fluid
  in question is the air flowing around an airplane, then we are typically
  interested in the drag and lift forces on the fuselage and wings. If the
  fluid is water flowing around a ship, then we typically care about
  the drag force on the ship.
- If you are solving the Maxwell equations of electromagnetics, you are
  typically interested in how much energy is radiated in certain directions
  (say, in order to know the range of information transmission via an
  antenna, or to determine the
  [radar cross section](https://en.wikipedia.org/wiki/Radar_cross_section)
  of planes or ships).
 
The point here is that from an engineering perspective, *solving* the
PDE is only the first step. The second step is to *evaluate* the computed
solution in order to extract relevant numbers that allow us to
either *optimize a design*, or to *make decisions*. This second step is often
called "postprocessing the solution".
 
This program does not solve a solid or fluid mechanics problem, so we
should try to illustrate postprocessing with something that makes sense
in the context of the equation we solve here. The Poisson equation in two
space dimensions is a model for the vertical deformation of a membrane
that is clamped at the boundary and is subject to a vertical force. For this
kind of situation, it makes sense to evaluate the *average vertical
displacement*,
@f[
  \bar u_h = \frac{\int_\Omega u_h(\mathbf x) \, dx}{|\Omega|},
@f]
where @f$|\Omega| = \int_\Omega 1 \, dx@f$ is the area of the domain. To compute
@f$\bar u_h@f$, as usual we replace integrals over the domain by a sum of integrals
over cells,
@f[
  \int_\Omega u_h(\mathbf x) \, dx
  =
  \sum_K \int_K u_h(\mathbf x) \, dx,
@f]
and then integrals over cells are approximated by quadrature:
@f{align*}{
  \int_\Omega u_h(\mathbf x) \, dx
  &=
  \sum_K \int_K u_h(\mathbf x) \, dx,
  \\
  &=
  \sum_K \sum_q u_h(\mathbf x_q^K) w_q^K,
@f}
where @f$w_q^K@f$ is the weight of the @f$q@f$th quadrature point evaluated on
cell @f$K@f$. All of this is as always provided by the FEValues class -- the
entry point for all integrals in deal.II.
 
The actual implementation of this is straightforward once you know how
to get the values of the solution @f$u@f$ at the quadrature points of a cell.
This functionality is provided by FEValues::get_function_values(), a
function that takes a global vector of nodal values as input and returns
a vector of function values at the quadrature points of the current cell.
Using this function, to see how it all works together you can
place the following code snippet anywhere in the program after the
solution has been computed (the `output_results()` function seems like a good
place to also do postprocessing, for example):
@code
  QGauss<dim>   quadrature_formula(fe.degree + 1);
  FEValues<dim> fe_values(fe,
                          quadrature_formula,
                          update_values | update_JxW_values);
 
  std::vector<double> solution_values(quadrature_formula.size());
  double              integral_of_u   = 0;
  double              volume_of_omega = 0;
 
  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      fe_values.reinit(cell);
      fe_values.get_function_values(solution, solution_values);
 
      for (const unsigned int q_point : fe_values.quadrature_point_indices())
        {
          integral_of_u += solution_values[q_point] * fe_values.JxW(q_point);
          volume_of_omega += 1 * fe_values.JxW(q_point);
        }
    }
  std::cout << "   Mean value of u=" << integral_of_u / volume_of_omega
            << std::endl;
@endcode
In this code snippet, we also compute the volume (or, since we are currently
thinking about a two-dimensional situation: the area) @f$|\Omega|@f$ by computing
the integral @f$|\Omega| = \int_\Omega 1 \, dx@f$ in exactly the same way, via
quadrature. (We could avoid having to compute @f$|\Omega|@f$ by hand here, using the
fact that deal.II has a function for this: GridTools::volume(). That said,
it is efficient to compute the two integrals
concurrently in the same loop, and so that's what we do.)
 
This program of course also solves the same Poisson equation in three space
dimensions. In this situation, the Poisson equation is often used as a model
for diffusion of either a physical species (say, of ink in a tank of water,
or a pollutant in the air) or of energy (specifically, of thermal energy in
a solid body). In that context, the quantity
@f[
  \Phi_h = \int_{\partial\Omega} \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx
@f]
is the *flux* of this species or energy across the boundary. (In actual
physical models, one would also have to multiply the right hand side by
a diffusivity or conductivity constant, but let us ignore this here.) In
much the same way as before, we compute such integrals by splitting
it over integrals of *faces* of cells, and then applying quadrature:
@f{align*}{
  \Phi_h
  &=
  \int_{\partial\Omega} \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx
  \\
  &=
  \sum_K
  \sum_{f \in \text{faces of @f$K@f$}, f\subset\partial\Omega}
  \int_f \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx
  \\
  &=
  \sum_K
  \sum_{f \in \text{faces of @f$K@f$}, f\subset\partial\Omega}
  \sum_q \nabla u_h(\mathbf x_q^f) \cdot \mathbf n(\mathbf x_q^f) w_q^f,
@f}
where now @f$\mathbf x_q^f@f$ are the quadrature points located on face @f$f@f$,
and @f$w_q^f@f$ are the weights associated with these faces. The second
of the sum symbols loops over all faces of cell @f$K@f$, but restricted to
those that are actually at the boundary.
 
This all is easily implemented by the following code that replaces the use of the
FEValues class (which is used for integrating over cells -- i.e., domain integrals)
by the FEFaceValues class (which is used for integrating over faces -- i.e.,
boundary integrals):
@code
  QGauss<dim - 1>   face_quadrature_formula(fe.degree + 1);
  FEFaceValues<dim> fe_face_values(fe,
                                   face_quadrature_formula,
                                   update_gradients | update_normal_vectors |
                                     update_JxW_values);
 
  std::vector<Tensor<1, dim>> solution_gradients(face_quadrature_formula.size());
  double                      flux = 0;
 
  for (const auto &cell : dof_handler.active_cell_iterators())
    for (const auto &face : cell->face_iterators())
      if (face->at_boundary())
        {
          fe_face_values.reinit(cell, face);
          fe_face_values.get_function_gradients(solution, solution_gradients);
 
          for (const unsigned int q_point :
               fe_face_values.quadrature_point_indices())
            {
              flux += solution_gradients[q_point] *
                      fe_face_values.normal_vector(q_point) *
                      fe_face_values.JxW(q_point);
            }
        }
  std::cout << "   Flux=" << flux << std::endl;
@endcode
 
If you add these two code snippets to the code, you will get output like the
following when you run the program:
@code
Solving problem in 2 space dimensions.
   Number of active cells: 256
   Total number of cells: 341
   Number of degrees of freedom: 289
   26 CG iterations needed to obtain convergence.
   Mean value of u=1.33303
   Flux=-3.68956
Solving problem in 3 space dimensions.
   Number of active cells: 4096
   Total number of cells: 4681
   Number of degrees of freedom: 4913
   30 CG iterations needed to obtain convergence.
   Mean value of u=1.58058
   Flux=-8.29435
@endcode
 
This makes some sense: If you look, for example, at the 2d output above,
the solution varies between values of 1 and 2, but with a larger part of the
solution closer to one than two; so an average value of 1.33 for the mean value
is reasonable. For the flux, recall that @f$\nabla u \cdot \mathbf n@f$ is the
directional derivative in the normal direction -- in other words, how the
solution changes as we move from the interior of the domain towards the
boundary. If you look at the 2d solution, you will realize that for most parts
of the boundary, the solution *decreases* as we approach the boundary, so the
normal derivative is negative -- so if we integrate along the boundary, we
should expect (and obtain!) a negative value.
 
 
 
<a name="step-4-extensions"></a>
<a name="step_4-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
There are many ways with which one can play with this program. The simpler
ones include essentially all the possibilities already discussed in the
<a href="step_3.html#extensions" target="body">Possibilities for extensions in the documentation of step 3</a>,
except that you will have to think about whether something now also applies
to the 3d case discussed in the current program.
 
It is also worthwhile considering the postprocessing options discussed
above. The documentation states two numbers (the mean value and the
normal flux) for both the 2d and 3d cases. Can we trust these
numbers? We have convinced ourselves that at least the mean value
is reasonable, and that the sign of the flux is probably correct.
But are these numbers accurate?
 
A general rule is that we should never trust a number unless we have
verified it in some way. From the theory of finite element methods,
we know that as we make the mesh finer and finer, the numerical
solution @f$u_h@f$ we compute here must converge to the exact solution
@f$u@f$. As a consequence, we also expect that @f$\bar u_h \rightarrow \bar u@f$
and @f$\Phi_h \rightarrow \Phi@f$, but that does not mean that for any
given mesh @f$\bar u_h@f$ or @f$\Phi_h@f$ are particularly accurate approximations.
 
To test this kind of thing, we have already considered the convergence of
a point value in @ref step_3 "step-3". We can do the same here by selecting how many
times the mesh is globally refined in the `make_grid()` function of this
program. For the mean value of the solution, we then get the following
numbers:
  <table align="center" class="doxtable">
    <tr> <th># of refinements</th>
         <th>@f$\bar u_h@f$ in 2d</th>
         <th>@f$\bar u_h@f$ in 3d</th>
    </tr>
    <tr> <td>4</td> <td>1.33303</td> <td>1.58058</td> </tr>
    <tr> <td>5</td> <td>1.33276</td> <td>1.57947</td> </tr>
    <tr> <td>6</td> <td>1.3327</td>  <td>1.5792</td> </tr>
    <tr> <td>7</td> <td>1.33269</td> <td>1.57914</td> </tr>
    <tr> <td>8</td> <td>1.33268</td> <td></td> </tr>
    <tr> <td>9</td> <td>1.33268</td> <td></td> </tr>
  </table>
I did not have the patience to run the last two values for the 3d case --
one needs quite a fine mesh for this, with correspondingly long run times.
But we can be reasonably assured that values around 1.33 (for the 2d case)
and 1.58 (for the 3d case) are about right -- and at least for engineering
applications, three digits of accuracy are good enough.
 
The situation looks very different for the flux. Here, we get results
such as the following:
  <table align="center" class="doxtable">
    <tr> <th># of refinements</th>
         <th>@f$\Phi_h@f$ in 2d</th>
         <th>@f$\Phi_h@f$ in 3d</th>
    </tr>
    <tr> <td>4</td>   <td>-3.68956</td>   <td>-8.29435</td>  </tr>
    <tr> <td>5</td>   <td>-4.90147</td>   <td>-13.0691</td>  </tr>
    <tr> <td>6</td>   <td>-5.60745</td>   <td>-15.9171</td>  </tr>
    <tr> <td>7</td>   <td>-5.99111</td>   <td>-17.4918</td>  </tr>
    <tr> <td>8</td>   <td>-6.19196</td>   <td></td>          </tr>
    <tr> <td>9</td>   <td>-6.29497</td>   <td></td>          </tr>
    <tr> <td>10</td>   <td>-6.34721</td>  <td></td>          </tr>
    <tr> <td>11</td>   <td>-6.37353</td>  <td></td>          </tr>
  </table>
So this is not great. For the 2d case, we might infer that perhaps
a value around -6.4 might be right if we just refine the mesh enough --
though 11 refinements already leads to some 4,194,304 cells. In any
case, the first number (the one shown in the beginning where we
discussed postprocessing) was off by almost a factor of 2!
 
For the 3d case, the last number shown was on a mesh with 2,097,152
cells; the next one would have had 8 times as many cells. In any case, the
numbers mean that we can't even be sure
that the first digit of that last number is correct! In other words,
it was worth checking, or we would have just believed all of these
numbers. In fact, that last column isn't even doing a particularly
good job convincing us that the code might be correctly implemented.
 
If you keep reading through the other tutorial programs, you will find many ways
to make these sorts of computations more accurate and to come to
believe that the flux actually does converge to its correct value.
For example, we can dramatically increase the accuracy of the computation
by using adaptive mesh refinement (@ref step_6 "step-6") near the boundary, and
in particular by using higher polynomial degree finite elements (also
@ref step_6 "step-6", but also @ref step_7 "step-7"). Using the latter, using cubic elements
(polynomial degree 3), we can actually compute the flux pretty
accurately even in 3d: @f$\Phi_h=-19.0148@f$ with 4 global refinement steps,
and @f$\Phi_h=-19.1533@f$ with 5 refinement steps. These numbers are already
pretty close together and give us a reasonable idea of the first
two correct digits of the "true" answer.
 
@note We would be remiss to not also comment on the fact that there
  are good theoretical reasons why computing the flux accurately
  appears to be so much more difficult than the average value.
  This has to do with the fact that finite element theory
  provides us with the estimate
  @f$\|u-u_h\|_{L_2(\Omega)} \le C h^2 \|\nabla^2u\|_{L_2(\Omega)}@f$
  when using the linear elements this program uses -- that is, for
  every global mesh refinement, @f$h@f$ is reduced by a factor of two
  and the error goes down by a factor of 4. Now, the @f$L_2@f$ error is
  not equivalent to the error in the mean value, but the two are
  related: They are both integrals over the domain, using the *value*
  of the solution. We expect the mean value to converge no worse than
  the @f$L_2@f$ norm of the error. At the same time, theory also provides
  us with this estimate:
  @f$\|\nabla (u-u_h)\|_{L_2(\partial\Omega)} \le
    C h^{1/2} \|\nabla^2u\|_{L_2(\Omega)}@f$. The move from values to
  gradients reduces the convergence rates by one order, and the move
  from domain to boundary by another half order. Here, then, each
  refinement step reduces the error not by a factor of 4 any more,
  by only by a factor of @f$\sqrt{2} \approx 1.4@f$. It takes a lot
  of global refinement steps to reduce the error by, say, a factor
  ten or hundred, and this is reflected in the very slow convergence
  evidenced by the table. On the other hand, for cubic elements (i.e.,
  polynomial degree 3), we would get
  @f$\|u-u_h\|_{L_2(\Omega)} \le C h^4 \|\nabla^4u\|_{L_2(\Omega)}@f$
  and after reduction by 1.5 orders, we would still have
  @f$\|\nabla (u-u_h)\|_{L_2(\partial\Omega)} \le
    C h^{2+1/2} \|\nabla^4u\|_{L_2(\Omega)}@f$. This rate,
  @f${\cal O}(h^{2.5})@f$ is still quite rapid, and it is perhaps not
  surprising that we get much better answers with these higher
  order elements. This also illustrates that when trying to
  approximate anything that relates to a gradient of the solution,
  using linear elements (polynomial degree one) is really not a
  good choice at all.
 
@note In this very specific case, it turns out that we can actually
  compute the exact value of @f$\Phi@f$. This is because for the Poisson
  equation we compute the solution of here, @f$-\Delta u = f@f$, we can
  integrate over the domain, @f$-\int_\Omega \Delta u = \int_\Omega f@f$,
  and then use that @f$\Delta = \text{div}\;\text{grad}@f$; this allows
  us to use the
  divergence theorem followed by multiplying by minus one to find
  @f$\int_{\partial\Omega} \nabla u \cdot n = -\int_\Omega f@f$. The
  left hand side happens to be @f$\Phi@f$. For the specific right
  hand side @f$f(x_1,x_2)=4(x_1^4+x_2^4)@f$ we use in 2d, we then
  get @f$-\int_\Omega f = -\int_{-1}^{1} \int_{-1}^{1} 4(x_1^4+x_2^4) \; dx_2\; dx_1
  = -16 \left[\int_{-1}^{1} x^4 \; dx\right] = -16\times\frac 25@f$,
  which has a numerical value of exactly -6.4 -- right on with our
  guess above. In 3d, we can do the same and get that the exact
  value is
  @f$-\int_\Omega f =
   -\int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} 4(x_1^4+x_2^4+x_3^4) \; dx_3 \; dx_2\; dx_1
   = -48\times\frac 25=-19.2@f$. What we found with cubic elements
  is then quite close to this exact value. Of course, in practice
  we almost never have exact values to compare with: If we could
  compute something on a piece of paper, we wouldn't have to solve
  the PDE numerically. But these sorts of situations make for excellent
  test cases that help us verify that our numerical solver works
  correctly. In many other cases, the literature contains
  numbers where others have already computed an answer accurately
  using their own software, and these are also often useful to
  compare against in verifying the correctness of our codes.
 *
 *
<a name="step_4-PlainProg"></a>
<h1> The plain program</h1>
@include "step-4.cc"
*/