step-3.h

Go to the documentation of this file.

 1)
 *     , dof_handler(triangulation)
 *   {}
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_3-Step3make_grid"></a> 
 * <h4>Step3::make_grid</h4>
 * 
 
 * 
 * Now, the first thing we've got to do is to generate the triangulation on
 * which we would like to do our computation and number each vertex with a
 * degree of freedom. We have seen these two steps in @ref step_1 "step-1" and @ref step_2 "step-2"
 * before, respectively.
 * 
 
 * 
 * This function does the first part, creating the mesh.  We create the grid
 * and refine all cells five times. Since the initial grid (which is the
 * square @f$[-1,1] \times [-1,1]@f$) consists of only one cell, the final grid
 * has 32 times 32 cells, for a total of 1024.
 * 
 
 * 
 * Unsure that 1024 is the correct number? We can check that by outputting the
 * number of cells using the <code>n_active_cells()</code> function on the
 * triangulation.
 * 
 * @code
 *   void Step3::make_grid()
 *   {
 *     GridGenerator::hyper_cube(triangulation, -1, 1);
 *     triangulation.refine_global(5);
 *   
 *     std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *               << std::endl;
 *   }
 *   
 * @endcode
 * 
 * @note We call the Triangulation::n_active_cells() function, rather than
 * Triangulation::n_cells(). Here, <i>active</i> means the cells that aren't
 * refined any further. We stress the adjective "active" since there are more
 * cells, namely the parent cells of the finest cells, their parents, etc, up
 * to the one cell which made up the initial grid. Of course, on the next
 * coarser level, the number of cells is one quarter that of the cells on the
 * finest level, i.e. 256, then 64, 16, 4, and 1. If you called
 * <code>triangulation.n_cells()</code> instead in the code above, you would
 * consequently get a value of 1365 instead. On the other hand, the number of
 * cells (as opposed to the number of active cells) is not typically of much
 * interest, so there is no good reason to print it.
 * 
 
 * 
 * 
 
 * 
 * 
 * <a name="step_3-Step3setup_system"></a> 
 * <h4>Step3::setup_system</h4>
 * 
 
 * 
 * Next we enumerate all the degrees of freedom and set up matrix and vector
 * objects to hold the system data. Enumerating is done by using
 * DoFHandler::distribute_dofs(), as we have seen in the @ref step_2 "step-2" example. Since
 * we use the FE_Q class and have set the polynomial degree to 1 in the
 * constructor, i.e. bilinear elements, this associates one degree of freedom
 * with each vertex. While we're at generating output, let us also take a look
 * at how many degrees of freedom are generated:
 * 
 * @code
 *   void Step3::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe);
 *     std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *               << std::endl;
 * @endcode
 * 
 * There should be one DoF for each vertex. Since we have a 32 times 32
 * grid, the number of DoFs should be 33 times 33, or 1089.
 * 
 
 * 
 * As we have seen in the previous example, we set up a sparsity pattern by
 * first creating a temporary structure, tagging those entries that might be
 * nonzero, and then copying the data over to the SparsityPattern object
 * that can then be used by the system matrix.
 * 
 * @code
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *     sparsity_pattern.copy_from(dsp);
 *   
 * @endcode
 * 
 * Note that the SparsityPattern object does not hold the values of the
 * matrix, it only stores the places where entries are. The entries
 * themselves are stored in objects of type SparseMatrix, of which our
 * variable system_matrix is one.
 *   
 
 * 
 * The distinction between sparsity pattern and matrix was made to allow
 * several matrices to use the same sparsity pattern. This may not seem
 * relevant here, but when you consider the size which matrices can have,
 * and that it may take some time to build the sparsity pattern, this
 * becomes important in large-scale problems if you have to store several
 * matrices in your program.
 * 
 * @code
 *     system_matrix.reinit(sparsity_pattern);
 *   
 * @endcode
 * 
 * The last thing to do in this function is to set the sizes of the right
 * hand side vector and the solution vector to the right values:
 * 
 * @code
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 *   }
 *   
 * @endcode
 * 
 * 
 * <a name="step_3-Step3assemble_system"></a> 
 * <h4>Step3::assemble_system</h4>
 * 
 
 * 
 * 
 
 * 
 * The next step is to compute the entries of the matrix and right hand side
 * that form the linear system from which we compute the solution. This is the
 * central function of each finite element program and we have discussed the
 * primary steps in the introduction already.
 * 
 
 * 
 * The general approach to assemble matrices and vectors is to loop over all
 * cells, and on each cell compute the contribution of that cell to the global
 * matrix and right hand side by quadrature. The point to realize now is that
 * we need the values of the shape functions at the locations of quadrature
 * points on the real cell. However, both the finite element shape functions
 * as well as the quadrature points are only defined on the reference
 * cell. They are therefore of little help to us, and we will in fact hardly
 * ever query information about finite element shape functions or quadrature
 * points from these objects directly.
 * 
 
 * 
 * Rather, what is required is a way to map this data from the reference cell
 * to the real cell. Classes that can do that are derived from the Mapping
 * class, though one again often does not have to deal with them directly:
 * many functions in the library can take a mapping object as argument, but
 * when it is omitted they simply resort to the standard bilinear Q1
 * mapping. We will go this route, and not bother with it for the moment (we
 * come back to this in @ref step_10 "step-10", @ref step_11 "step-11", and @ref step_12 "step-12").
 * 
 
 * 
 * So what we now have is a collection of three classes to deal with: finite
 * element, quadrature, and mapping objects. That's too much, so there is one
 * type of class that orchestrates information exchange between these three:
 * the FEValues class. If given one instance of each three of these objects
 * (or two, and an implicit linear mapping), it will be able to provide you
 * with information about values and gradients of shape functions at
 * quadrature points on a real cell.
 * 
 
 * 
 * Using all this, we will assemble the linear system for this problem in the
 * following function:
 * 
 * @code
 *   void Step3::assemble_system()
 *   {
 * @endcode
 * 
 * Ok, let's start: we need a quadrature formula for the evaluation of the
 * integrals on each cell. Let's take a Gauss formula with two quadrature
 * points in each direction, i.e. a total of four points since we are in
 * 2d. This quadrature formula integrates polynomials of degrees up to three
 * exactly (in 1d). It is easy to check that this is sufficient for the
 * present problem:
 * 
 * @code
 *     const QGauss<2> quadrature_formula(fe.degree + 1);
 * @endcode
 * 
 * And we initialize the object which we have briefly talked about above. It
 * needs to be told which finite element we want to use, and the quadrature
 * points and their weights (jointly described by a Quadrature object). As
 * mentioned, we use the implied Q1 mapping, rather than specifying one
 * ourselves explicitly. Finally, we have to tell it what we want it to
 * compute on each cell: we need the values of the shape functions at the
 * quadrature points (for the right hand side @f$(\varphi_i,f)@f$), their
 * gradients (for the matrix entries @f$(\nabla \varphi_i, \nabla
 * \varphi_j)@f$), and also the weights of the quadrature points and the
 * determinants of the Jacobian transformations from the reference cell to
 * the real cells.
 *   
 
 * 
 * This list of what kind of information we actually need is given as a
 * collection of flags as the third argument to the constructor of
 * FEValues. Since these values have to be recomputed, or updated, every
 * time we go to a new cell, all of these flags start with the prefix
 * <code>update_</code> and then indicate what it actually is that we want
 * updated. The flag to give if we want the values of the shape functions
 * computed is #update_values; for the gradients it is
 * #update_gradients. The determinants of the Jacobians and the quadrature
 * weights are always used together, so only the products (Jacobians times
 * weights, or short <code>JxW</code>) are computed; since we need them, we
 * have to list #update_JxW_values as well:
 * 
 * @code
 *     FEValues<2> fe_values(fe,
 *                           quadrature_formula,
 *                           update_values | update_gradients | update_JxW_values);
 * @endcode
 * 
 * The advantage of this approach is that we can specify what kind of
 * information we actually need on each cell. It is easily understandable
 * that this approach can significantly speed up finite element computations,
 * compared to approaches where everything, including second derivatives,
 * normal vectors to cells, etc are computed on each cell, regardless of
 * whether they are needed or not.
 *   
 
 * 
 * @note The syntax <code>update_values | update_gradients |
 * update_JxW_values</code> is not immediately obvious to anyone not
 * used to programming bit operations in C for years already. First,
 * <code>operator|</code> is the <i>bitwise or operator</i>, i.e.,
 * it takes two integer arguments that are interpreted as bit
 * patterns and returns an integer in which every bit is set for
 * which the corresponding bit is set in at least one of the two
 * arguments. For example, consider the operation
 * <code>9|10</code>. In binary, <code>9=0b1001</code> (where the
 * prefix <code>0b</code> indicates that the number is to be
 * interpreted as a binary number) and <code>10=0b1010</code>. Going
 * through each bit and seeing whether it is set in one of the
 * argument, we arrive at <code>0b1001|0b1010=0b1011</code> or, in
 * decimal notation, <code>9|10=11</code>. The second piece of
 * information you need to know is that the various
 * <code>update_*</code> flags are all integers that have <i>exactly
 * one bit set</i>. For example, assume that
 * <code>update_values=0b00001=1</code>,
 * <code>update_gradients=0b00010=2</code>,
 * <code>update_JxW_values=0b10000=16</code>. Then
 * <code>update_values | update_gradients | update_JxW_values =
 * 0b10011 = 19</code>. In other words, we obtain a number that
 * <i>encodes a binary mask representing all of the operations you
 * want to happen</i>, where each operation corresponds to exactly
 * one bit in the integer that, if equal to one, means that a
 * particular piece should be updated on each cell and, if it is
 * zero, means that we need not compute it. In other words, even
 * though <code>operator|</code> is the <i>bitwise OR operation</i>,
 * what it really represents is <i>I want this AND that AND the
 * other</i>. Such binary masks are quite common in C programming,
 * but maybe not so in higher level languages like C++, but serve
 * the current purpose quite well.
 * 
 
 * 
 * For use further down below, we define a shortcut for a value that will
 * be used very frequently. Namely, an abbreviation for the number of degrees
 * of freedom on each cell (since we are in 2d and degrees of freedom are
 * associated with vertices only, this number is four, but we rather want to
 * write the definition of this variable in a way that does not preclude us
 * from later choosing a different finite element that has a different
 * number of degrees of freedom per cell, or work in a different space
 * dimension).
 *   
 
 * 
 * In general, it is a good idea to use a symbolic name instead of
 * hard-coding these numbers even if you know them, since for example,
 * you may want to change the finite element at some time. Changing the
 * element would have to be done in a different function and it is easy
 * to forget to make a corresponding change in another part of the program.
 * It is better to not rely on your own calculations, but instead ask
 * the right object for the information: Here, we ask the finite element
 * to tell us about the number of degrees of freedom per cell and we
 * will get the correct number regardless of the space dimension or
 * polynomial degree we may have chosen elsewhere in the program.
 *   
 
 * 
 * The shortcut here, defined primarily to discuss the basic concept
 * and not because it saves a lot of typing, will then make the following
 * loops a bit more readable. You will see such shortcuts in many places in
 * larger programs, and `dofs_per_cell` is one that is more or less the
 * conventional name for this kind of object.
 * 
 * @code
 *     const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *   
 * @endcode
 * 
 * Now, we said that we wanted to assemble the global matrix and vector
 * cell-by-cell. We could write the results directly into the global matrix,
 * but this is not very efficient since access to the elements of a sparse
 * matrix is slow. Rather, we first compute the contribution of each cell in
 * a small matrix with the degrees of freedom on the present cell, and only
 * transfer them to the global matrix when the computations are finished for
 * this cell. We do the same for the right hand side vector. So let's first
 * allocate these objects (these being local objects, all degrees of freedom
 * are coupling with all others, and we should use a full matrix object
 * rather than a sparse one for the local operations; everything will be
 * transferred to a global sparse matrix later on):
 * 
 * @code
 *     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *     Vector<double>     cell_rhs(dofs_per_cell);
 *   
 * @endcode
 * 
 * When assembling the contributions of each cell, we do this with the local
 * numbering of the degrees of freedom (i.e. the number running from zero
 * through dofs_per_cell-1). However, when we transfer the result into the
 * global matrix, we have to know the global numbers of the degrees of
 * freedom. When we query them, we need a scratch (temporary) array for
 * these numbers (see the discussion at the end of the introduction for
 * the type, types::global_dof_index, used here):
 * 
 * @code
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 * @endcode
 * 
 * Now for the loop over all cells. We have seen before how this works for a
 * triangulation. A DoFHandler has cell iterators that are exactly analogous
 * to those of a Triangulation, but with extra information about the degrees
 * of freedom for the finite element you're using. Looping over the active
 * cells of a degree-of-freedom handler works the same as for a triangulation.
 *   
 
 * 
 * Note that we declare the type of the cell as `const auto &` instead of
 * `auto` this time around. In step 1, we were modifying the cells of the
 * triangulation by flagging them with refinement indicators. Here we're only
 * examining the cells without modifying them, so it's good practice to
 * declare `cell` as `const` in order to enforce this invariant.
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 * @endcode
 * 
 * We are now sitting on one cell, and we would like the values and
 * gradients of the shape functions be computed, as well as the
 * determinants of the Jacobian matrices of the mapping between
 * reference cell and true cell, at the quadrature points. Since all
 * these values depend on the geometry of the cell, we have to have the
 * FEValues object re-compute them on each cell:
 * 
 * @code
 *         fe_values.reinit(cell);
 *   
 * @endcode
 * 
 * Next, reset the local cell's contributions to global matrix and
 * global right hand side to zero, before we fill them:
 * 
 * @code
 *         cell_matrix = 0;
 *         cell_rhs    = 0;
 *   
 * @endcode
 * 
 * Now it is time to start integration over the cell, which we
 * do by looping over all quadrature points, which we will
 * number by q_index.
 * 
 * @code
 *         for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *           {
 * @endcode
 * 
 * First assemble the matrix: For the Laplace problem, the
 * matrix on each cell is the integral over the gradients of
 * shape function i and j. Since we do not integrate, but
 * rather use quadrature, this is the sum over all
 * quadrature points of the integrands times the determinant
 * of the Jacobian matrix at the quadrature point times the
 * weight of this quadrature point. You can get the gradient
 * of shape function @f$i@f$ at quadrature point with number q_index by
 * using <code>fe_values.shape_grad(i,q_index)</code>; this
 * gradient is a 2-dimensional vector (in fact it is of type
 * Tensor@<1,dim@>, with here dim=2) and the product of two
 * such vectors is the scalar product, i.e. the product of
 * the two shape_grad function calls is the dot
 * product. This is in turn multiplied by the Jacobian
 * determinant and the quadrature point weight (that one
 * gets together by the call to FEValues::JxW() ). Finally,
 * this is repeated for all shape functions @f$i@f$ and @f$j@f$:
 * 
 * @code
 *             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
 *   
 * @endcode
 * 
 * We then do the same thing for the right hand side. Here,
 * the integral is over the shape function i times the right
 * hand side function, which we choose to be the function
 * with constant value one (more interesting examples will
 * be considered in the following programs).
 * 
 * @code
 *             for (const unsigned int i : fe_values.dof_indices())
 *               cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
 *                               1. *                                // f(x_q)
 *                               fe_values.JxW(q_index));            // dx
 *           }
 * @endcode
 * 
 * Now that we have the contribution of this cell, we have to transfer
 * it to the global matrix and right hand side. To this end, we first
 * have to find out which global numbers the degrees of freedom on this
 * cell have. Let's simply ask the cell for that information:
 * 
 * @code
 *         cell->get_dof_indices(local_dof_indices);
 *   
 * @endcode
 * 
 * Then again loop over all shape functions i and j and transfer the
 * local elements to the global matrix. The global numbers can be
 * obtained using local_dof_indices[i]:
 * 
 * @code
 *         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));
 *   
 * @endcode
 * 
 * And again, we do the same thing for the right hand side vector.
 * 
 * @code
 *         for (const unsigned int i : fe_values.dof_indices())
 *           system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *       }
 *   
 *   
 * @endcode
 * 
 * Now almost everything is set up for the solution of the discrete
 * system. However, we have not yet taken care of boundary values (in fact,
 * Laplace's equation without Dirichlet boundary values is not even uniquely
 * solvable, since you can add an arbitrary constant to the discrete
 * solution). We therefore have to do something about the situation.
 *   
 
 * 
 * For this, we first obtain a list of the degrees of freedom on the
 * boundary and the value the shape function shall have there. For
 * simplicity, we only interpolate the boundary value function, rather than
 * projecting it onto the boundary. There is a function in the library which
 * does exactly this: VectorTools::interpolate_boundary_values(). Its
 * parameters are (omitting parameters for which default values exist and
 * that we don't care about): the DoFHandler object to get the global
 * numbers of the degrees of freedom on the boundary; the component of the
 * boundary where the boundary values shall be interpolated; the boundary
 * value function itself; and the output object.
 *   
 
 * 
 * The component of the boundary is meant as follows: in many cases, you may
 * want to impose certain boundary values only on parts of the boundary. For
 * example, you may have inflow and outflow boundaries in fluid dynamics, or
 * clamped and free parts of bodies in deformation computations of
 * bodies. Then you will want to denote these different parts of the
 * boundary by indicators, and tell the interpolate_boundary_values
 * function to only compute the boundary values on a certain part of the
 * boundary (e.g. the clamped part, or the inflow boundary). By default,
 * all boundaries have a 0 boundary indicator, unless otherwise specified.
 * (For example, many functions in namespace GridGenerator specify otherwise.)
 * If sections of the boundary have different boundary conditions, you have to
 * number those parts with different boundary indicators. The function call
 * below will then only determine boundary values for those parts of the
 * boundary for which the boundary indicator is in fact the zero specified as
 * the second argument.
 *   
 
 * 
 * The function describing the boundary values is an object of type Function
 * or of a derived class. One of the derived classes is
 * Functions::ZeroFunction, which describes (not unexpectedly) a function
 * which is zero everywhere. We create such an object in-place and pass it to
 * the VectorTools::interpolate_boundary_values() function.
 *   
 
 * 
 * Finally, the output object is a list of pairs of global degree of freedom
 * numbers (i.e. the number of the degrees of freedom on the boundary) and
 * their boundary values (which are zero here for all entries). This mapping
 * of DoF numbers to boundary values is done by the <code>std::map</code>
 * class.
 * 
 * @code
 *     std::map<types::global_dof_index, double> boundary_values;
 *     VectorTools::interpolate_boundary_values(dof_handler,
 *                                              types::boundary_id(0),
 *                                              Functions::ZeroFunction<2>(),
 *                                              boundary_values);
 * @endcode
 * 
 * Now that we got the list of boundary DoFs and their respective boundary
 * values, let's use them to modify the system of equations
 * accordingly. This is done by the following function call:
 * 
 * @code
 *     MatrixTools::apply_boundary_values(boundary_values,
 *                                        system_matrix,
 *                                        solution,
 *                                        system_rhs);
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_3-Step3solve"></a> 
 * <h4>Step3::solve</h4>
 * 
 
 * 
 * The following function solves the discretized equation. As discussed in
 * the introduction, we want to use an iterative solver to do this,
 * specifically the Conjugate Gradient (CG) method.
 * 
 
 * 
 * The way to do this in deal.II is a three-step process:
 * - First, we need to have an object that knows how to tell the CG algorithm
 * when to stop. This is done by using a SolverControl object, and as
 * stopping criterion we say: stop after a maximum of 1000 iterations (which
 * is far more than is needed for 1089 variables; see the results section to
 * find out how many were really used), and stop if the norm of the residual
 * is below @f$\tau=10^{-6}\|\mathbf b\|@f$ where @f$\mathbf b@f$ is the right hand
 * side vector. In practice, this latter criterion will be the one
 * which stops the iteration.
 * - Then we need the solver itself. The template parameter to the SolverCG
 * class is the type of the vectors we are using.
 * - The last step is to actually solve the system of equations. The CG solver
 * takes as arguments the components of the linear system @f$Ax=b@f$ (in the
 * order in which they appear in this equation), and a preconditioner
 * as the fourth argument. We don't feel ready to delve into preconditioners
 * yet, so we tell it to use the identity operation as preconditioner. Later
 * tutorial programs will spend significant amount of time and space on
 * constructing better preconditioners.
 * 
 
 * 
 * At the end of this process, the `solution` variable contains the
 * nodal values of the solution function. At the end of the function, we
 * output how many Conjugate Gradients iterations it took to solve the
 * linear system.
 * 
 * @code
 *   void Step3::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_3-Step3output_results"></a> 
 * <h4>Step3::output_results</h4>
 * 
 
 * 
 * The last part of a typical finite element program is to output the results
 * and maybe do some postprocessing (for example compute the maximal stress
 * values at the boundary, or the average flux across the outflow, etc). We
 * have no such postprocessing here, but we would like to write the solution
 * to a file.
 * 
 * @code
 *   void Step3::output_results() const
 *   {
 * @endcode
 * 
 * To write the output to a file, we need an object which knows about output
 * formats and the like. This is the DataOut class, and we need an object of
 * that type:
 * 
 * @code
 *     DataOut<2> data_out;
 * @endcode
 * 
 * Now we have to tell it where to take the values from which it shall
 * write. We tell it which DoFHandler object to use, and the solution vector
 * (and the name by which the solution variable shall appear in the output
 * file). If we had more than one vector which we would like to look at in
 * the output (for example right hand sides, errors per cell, etc) we would
 * add them as well:
 * 
 * @code
 *     data_out.attach_dof_handler(dof_handler);
 *     data_out.add_data_vector(solution, "solution");
 * @endcode
 * 
 * After the DataOut object knows which data it is to work on, we have to
 * tell it to process them into something the back ends can handle. The
 * reason is that we have separated the frontend (which knows about how to
 * treat DoFHandler objects and data vectors) from the back end (which knows
 * many different output formats) and use an intermediate data format to
 * transfer data from the front- to the backend. The data is transformed
 * into this intermediate format by the following function:
 * 
 * @code
 *     data_out.build_patches();
 *   
 * @endcode
 * 
 * Now we have everything in place for the actual output. Just open a file
 * and write the data into it, using VTK format (there are many other
 * functions in the DataOut class we are using here that can write the
 * data in postscript, AVS, GMV, Gnuplot, or some other file
 * formats):
 * 
 * @code
 *     const std::string filename = "solution.vtk";
 *     std::ofstream     output(filename);
 *     data_out.write_vtk(output);
 *     std::cout << "Output written to " << filename << std::endl;
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_3-Step3run"></a> 
 * <h4>Step3::run</h4>
 * 
 
 * 
 * Finally, the last function of this class is the main function which calls
 * all the other functions of the <code>Step3</code> class. The order in which
 * this is done resembles the order in which most finite element programs
 * work. Since the names are mostly self-explanatory, there is not much to
 * comment about:
 * 
 * @code
 *   void Step3::run()
 *   {
 *     make_grid();
 *     setup_system();
 *     assemble_system();
 *     solve();
 *     output_results();
 *   }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_3-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * This is the main function of the program. Since the concept of a
 * main function is mostly a remnant from the pre-object oriented era
 * before C++ programming, it often does not do much more than
 * creating an object of the top-level class and calling its principle
 * function.
 * 
 * @code
 *   int main()
 *   {
 *     Step3 laplace_problem;
 *     laplace_problem.run();
 *   
 *     return 0;
 *   }
 * @endcode
<a name="step_3-Results"></a><h1>Results</h1>
 
 
The output of the program looks as follows:
@code
Number of active cells: 1024
Number of degrees of freedom: 1089
36 CG iterations needed to obtain convergence.
Output written to solution.vtk
@endcode
 
The last line is the output we generated at the bottom of the
`output_results()` function: The program generated the file
<code>solution.vtk</code>, which is in the VTK format that is widely
used by many visualization programs today -- including the two
heavy-weights <a href="https://www.llnl.gov/visit">VisIt</a> and
<a href="https://www.paraview.org">Paraview</a> that are the most
commonly used programs for this purpose today.
 
Using VisIt, it is not very difficult to generate a picture of the
solution like this:
<table width="60%" align="center">
  <tr>
    <td align="center">
      <img src="https://dealii.org/images/steps/developer/step-3.solution-3.png" alt="Visualization of the solution of step-3">
    </td>
  </tr>
</table>
It shows both the solution and the mesh, elevated above the @f$x@f$-@f$y@f$ plane
based on the value of the solution at each point. Of course the solution
here is not particularly exciting, but that is a result of both what the
Laplace equation represents and the right hand side @f$f(\mathbf x)=1@f$ we
have chosen for this program: The Laplace equation describes (among many
other uses) the vertical deformation of a membrane subject to an external,
vertical force. In the current example, the membrane's borders
are clamped to a square frame with no vertical variation; a constant
force density will therefore intuitively lead to a membrane that
simply bulges upward -- like the one shown above.
 
VisIt and Paraview both allow playing with various kinds of visualizations
of the solution. Several video lectures show how to use these programs.
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>.
 
 
 
@anchor step_3-Extensions
<a name="step_3-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
If you want to play around a little bit with this program, here are a few
suggestions:
</p>
 
<ul>
  <li>
  Change the geometry and mesh: In the program, we have generated a square
  domain and mesh by using the GridGenerator::hyper_cube()
  function. However, the <code>GridGenerator</code> has a good number of other
  functions as well. Try an L-shaped domain, a ring, or other domains you find
  there.
  </li>
 
  <li>
  Change the boundary condition: The code uses the Functions::ZeroFunction
  function to generate zero boundary conditions. However, you may want to try
  non-zero constant boundary values using
  <code>Functions::ConstantFunction&lt;2&gt;(1)</code> instead of
  <code>Functions::ZeroFunction&lt;2&gt;()</code> to have unit Dirichlet boundary
  values. More exotic functions are described in the documentation of the
  Functions namespace, and you may pick one to describe your particular boundary
  values.
  </li>
 
  <li> Modify the type of boundary condition: Presently, what happens
  is that we use Dirichlet boundary values all around, since the
  default is that all boundary parts have boundary indicator zero, and
  then we tell the
  VectorTools::interpolate_boundary_values() function to
  interpolate boundary values to zero on all boundary components with
  indicator zero.  <p> We can change this behavior if we assign parts
  of the boundary different indicators. For example, try this
  immediately after calling GridGenerator::hyper_cube():
  @code
  triangulation.begin_active()->face(0)->set_boundary_id(1);
  @endcode
 
  What this does is it first asks the triangulation to
  return an iterator that points to the first active cell. Of course,
  this being the coarse mesh for the triangulation of a square, the
  triangulation has only a single cell at this moment, and it is
  active. Next, we ask the cell to return an iterator to its first
  face, and then we ask the face to reset the boundary indicator of
  that face to 1. What then follows is this: When the mesh is refined,
  faces of child cells inherit the boundary indicator of their
  parents, i.e. even on the finest mesh, the faces on one side of the
  square have boundary indicator 1. Later, when we get to
  interpolating boundary conditions, the
  VectorTools::interpolate_boundary_values() call will only produce boundary
  values for those faces that have zero boundary indicator, and leave
  those faces alone that have a different boundary indicator. What
  this then does is to impose Dirichlet boundary conditions on the
  former, and homogeneous Neumann conditions on the latter (i.e. zero
  normal derivative of the solution, unless one adds additional terms
  to the right hand side of the variational equality that deal with
  potentially non-zero Neumann conditions). You will see this if you
  run the program.
 
  An alternative way to change the boundary indicator is to label
  the boundaries based on the Cartesian coordinates of the face centers.
  For example, we can label all of the cells along the top and
  bottom boundaries with a boundary indicator 1 by checking to
  see if the cell centers' y-coordinates are within a tolerance
  (here `1e-12`) of -1 and 1. Try this immediately after calling
  GridGenerator::hyper_cube(), as before:
  @code
  for (auto &face : triangulation.active_face_iterators())
    if (face->at_boundary())
      if (std::fabs(face->center()(1) - (-1.0)) < 1e-12 ||
          std::fabs(face->center()(1) - (1.0)) < 1e-12)
        face->set_boundary_id(1);
  @endcode
  Although this code is a bit longer than before, it is useful for
  complex geometries, as it does not require knowledge of face labels.
 
  <li>
  A slight variation of the last point would be to set different boundary
  values as above, but then use a different boundary value function for
  boundary indicator one. In practice, what you have to do is to add a second
  call to <code>interpolate_boundary_values</code> for boundary indicator one:
  @code
  VectorTools::interpolate_boundary_values(dof_handler,
                                           types::boundary_id(1),
                                           Functions::ConstantFunction<2>(1.),
                                           boundary_values);
  @endcode
  If you have this call immediately after the first one to this function, then
  it will interpolate boundary values on faces with boundary indicator 1 to the
  unit value, and merge these interpolated values with those previously
  computed for boundary indicator 0. The result will be that we will get
  discontinuous boundary values, zero on three sides of the square, and one on
  the fourth.
 
  <li>
  Use triangles: As mentioned in the results section of @ref step_1 "step-1", for
  historical reasons, almost all tutorial programs for deal.II are
  written using quadrilateral or hexahedral meshes. But deal.II also
  supports triangular and tetrahedral meshes. So a good experiment would
  be to replace the mesh used here by a triangular mesh.
 
  This is *almost* trivial. First, as discussed in @ref step_1 "step-1", we may want
  to start with the quadrilateral mesh we are already creating, and
  then convert it into a triangular one. You can do that by replacing
  the first line of `Step3::make_grid()` by the following code:
  @code
  Triangulation<2> triangulation_quad;
  GridGenerator::hyper_cube(triangulation_quad, -1, 1);
  GridGenerator::convert_hypercube_to_simplex_mesh (triangulation_quad,
                                                    triangulation);
  @endcode
  The GridGenerator::convert_hypercube_to_simplex_mesh() replaces each
  quadrilateral by eight triangles with half the diameter of the original
  quadrilateral; as a consequence, the resulting mesh is substantially
  finer and one might expect that the solution is consequently more
  accurate (but also has many more degrees of freedom). That is a question
  you can explore with the techniques discussed in the "Results" section
  of @ref step_4 "step-4", but that goes beyond what we want to demonstrate here.
 
  If you run this program, you will run into an error message that
  will look something like this:
  @code
--------------------------------------------------------
An error occurred in line <2633> of file </home/bangerth/p/deal.II/1/dealii/include/deal.II/dofs/dof_accessor.templates.h> in function
    const ::FiniteElement<dimension_, space_dimension_>& ::DoFCellAccessor<dim, spacedim, lda>::get_fe() const [with int dimension_ = 2; int space_dimension_ = 2; bool level_dof_access = false]
The violated condition was:
    this->reference_cell() == fe.reference_cell()
Additional information:
    The reference-cell type used on this cell (Tri) does not match the
    reference-cell type of the finite element associated with this cell
    (Quad). Did you accidentally use simplex elements on hypercube meshes
    (or the other way around), or are you using a mixed mesh and assigned
    a simplex element to a hypercube cell (or the other way around) via
    the active_fe_index?
  @endcode
  It is worth carefully reading the error message. It doesn't just
  state that there is an error, but also how it may have
  arisen. Specifically, it asks whether we are using a finite element
  for simplex meshes (in 2d simplices are triangles) with a hypercube
  mesh (in 2d hypercubes are quadrilaterals) or the other way around?
 
  Of course, this is exactly what we are doing, though this may
  perhaps not be clear to you. But if you look up the documentation,
  you will find that the FE_Q element we use in the main class can
  only be used on hypercube meshes; what we *want* to use instead now
  that we are using a simplex mesh is the FE_SimplexP class that is the
  equivalent to FE_Q for simplex cells. (To do this, you will also
  have to add `#include <deal.II/fe/fe_simplex_p.h>` at the top of the file.)
 
  The last thing you need to change (which at the time of writing is
  unfortunately not prompted by getting an error message) is that when
  we integrate, we need to use a quadrature formula that is
  appropriate for triangles. This is done by changing QGauss by
  QGaussSimplex in the code.
 
  With all of these steps, you then get the following solution:
  <img src="https://dealii.org/images/steps/developer/step-3.solution-triangles.png" alt="Visualization of the solution of step-3 using triangles">
 
  <li>
  Verify correctness and observe convergence: When you compute something,
  your next step should be to make sure it is actually correct. That's
  probably good life advice anyway, but it is a particularly useful skill in scientific
  computing. In the current context, it is made more difficult by the fact
  that we do not actually know what the exact solution to the Laplace
  equation with the given right hand is. In the @ref step_7 "step-7" tutorial program,
  we will come back to this issue and come up with ways to *choose* the
  right hand side in a way so that we know what the exact solution is; we can
  then compare the computed (approximate) solution against this known (exact) solution.
 
  Here, we can at least check one aspect of a correct implementation:
  If we make the mesh finer and finer, finite element theory tells us
  the solution needs to *converge* to the exact solution. Because we
  don't know what the exact solution is, so for the moment all we can
  say is that the computed solutions need to converge to *something*
  as the mesh is refined. This is true for the solution as a function
  of @f$x,y@f$, and it is specifically also true for the solution at a
  particular point. This is easy to check here: We can evaluate the
  value of the solution at a
  single point and compare the value for different %numbers of global
  refinement (the number of global refinement steps is set in
  <code>LaplaceProblem::make_grid()</code> above). To evaluate the
  solution at a point, say at @f$(\frac 13, \frac 13)@f$, we could add the
  following code to the <code>LaplaceProblem::output_results</code> function:
  @code
    std::cout << "Solution at (1/3,1/3): "
              << VectorTools::point_value(dof_handler, solution,
                                          Point<2>(1./3, 1./3))
              << std::endl;
  @endcode
  For 1 through 9 global refinement steps, we then get the following sequence
  of point values:
 
  <table align="center" class="doxtable">
    <tr> <th># of refinements</th> <th>@f$u_h(\frac 13,\frac13)@f$</th> </tr>
    <tr> <td>1</td> <td>0.166667</td> </tr>
    <tr> <td>2</td> <td>0.227381</td> </tr>
    <tr> <td>3</td> <td>0.237375</td> </tr>
    <tr> <td>4</td> <td>0.240435</td> </tr>
    <tr> <td>5</td> <td>0.241140</td> </tr>
    <tr> <td>6</td> <td>0.241324</td> </tr>
    <tr> <td>7</td> <td>0.241369</td> </tr>
    <tr> <td>8</td> <td>0.241380</td> </tr>
    <tr> <td>9</td> <td>0.241383</td> </tr>
  </table>
 
  By noticing that the difference between each two consecutive values reduces
  by about a factor of 4, we can conjecture that the "correct" value may be
  @f$u(\frac 13, \frac 13)\approx 0.241384@f$. In fact, if we assumed this to be
  the correct value, we could show that the sequence above indeed shows @f${\cal
  O}(h^2)@f$ convergence &mdash; theoretically, the convergence order should be
  @f${\cal O}(h^2 |\log h|)@f$ (because that's the theoretical error bound in the
  @f$L_\infty@f$ norm) but the symmetry of the domain and the mesh may lead
  to the better convergence order observed.
 
  A slight variant of this would be to repeat the test with quadratic
  elements. All you need to do is to set the polynomial degree of the finite
  element to two in the constructor
  <code>LaplaceProblem::LaplaceProblem</code>.
 
  <li>Convergence of the mean: A different way to see that the solution
  actually converges (to something) is to compute the mean of the solution.
  To this end, add
  the following code to <code>LaplaceProblem::output_results</code>:
  @code
    std::cout << "Mean value: "
              << VectorTools::compute_mean_value (dof_handler,
                                                  QGauss<2>(fe.degree + 1),
                                                  solution,
                                                  0)
              << std::endl;
  @endcode
  The documentation of the VectorTools::compute_mean_value() function
  explains what the second and fourth
  parameters mean, while the first and third should be obvious. Doing the same
  study again where we change the number of global refinement steps, we get
  the following result:
 
  <table align="center" class="doxtable">
    <tr> <th># of refinements</th> <th>@f$\int_\Omega u_h(x)\; dx@f$</th> </tr>
    <tr> <td>0</td> <td>0.09375000</td> </tr>
    <tr> <td>1</td> <td>0.12790179</td> </tr>
    <tr> <td>2</td> <td>0.13733440</td> </tr>
    <tr> <td>3</td> <td>0.13976069</td> </tr>
    <tr> <td>4</td> <td>0.14037251</td> </tr>
    <tr> <td>5</td> <td>0.14052586</td> </tr>
    <tr> <td>6</td> <td>0.14056422</td> </tr>
    <tr> <td>7</td> <td>0.14057382</td> </tr>
    <tr> <td>8</td> <td>0.14057622</td> </tr>
  </table>
 
  Again, the difference between two adjacent values goes down by about a
  factor of four, indicating convergence as @f${\cal O}(h^2)@f$.
</ul>
 
 
 
<a name="step_3-UsingHDF5tooutputthesolutionandadditionaldata"></a><h3>Using %HDF5 to output the solution and additional data</h3>
 
 
%HDF5 is a commonly used format that can be read by many scripting
languages (e.g. R or Python). It is not difficult to get deal.II to
produce some %HDF5 files that can then be used in external scripts to
postprocess some of the data generated by this program. Here are some
ideas on what is possible.
 
 
<a name="step_3-Changingtheoutputtoh5"></a><h4> Changing the output to .h5</h4>
 
 
To fully make use of the automation we first need to introduce a private variable for the number of
global refinement steps <code>unsigned int n_refinement_steps </code>, which will be used for the output filename.
In <code>make_grid()</code> we then replace <code>triangulation.refine_global(5);</code> with
@code
n_refinement_steps = 5;
triangulation.refine_global(n_refinement_steps);
@endcode
The deal.II library has two different %HDF5 bindings, one in the HDF5
namespace (for interfacing to general-purpose data files)
and another one in DataOut (specifically for writing files for the
visualization of solutions).
Although the HDF5 deal.II binding supports both serial and MPI, the %HDF5 DataOut binding
only supports parallel output.
For this reason we need to initialize an MPI
communicator with only one processor. This is done by adding the following code.
@code
int main(int argc, char* argv[])
{
  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
  ...
}
@endcode
Next we change the `Step3::output_results()` output routine as
described in the DataOutBase namespace documentation:
@code
const std::string filename_h5 = "solution_" + std::to_string(n_refinement_steps) + ".h5";
DataOutBase::DataOutFilterFlags flags(true, true);
DataOutBase::DataOutFilter data_filter(flags);
data_out.write_filtered_data(data_filter);
data_out.write_hdf5_parallel(data_filter, filename_h5, MPI_COMM_WORLD);
@endcode
The resulting file can then be visualized just like the VTK file that
the original version of the tutorial produces; but, since %HDF5 is a
more general file format, it can also easily be processed in scripting
languages for other purposes.
 
 
<a name="step_3-Addingthepointvalueandthemeanseeextensionaboveintotheh5file"></a><h4> Adding the point value and the mean (see extension above) into the .h5 file</h4>
 
 
After outputting the solution, the file can be opened again to include
more datasets.  This allows us to keep all the necessary information
of our experiment in a single result file, which can then be read and
processed by some postprocessing script.
(Have a look at HDF5::Group::write_dataset() for further
information on the possible output options.)
 
To make this happen, we first include the necessary header into our file :
@code
#include <deal.II/base/hdf5.h>
@endcode
Adding the following lines to the end
of our output routine adds the information about the value of the
solution at a particular point, as well as the mean value of the
solution, to our %HDF5 file :
@code
HDF5::File data_file(filename_h5, HDF5::File::FileAccessMode::open, MPI_COMM_WORLD);
Vector<double> point_value(1);
point_value[0] = VectorTools::point_value(dof_handler, solution,
                                          Point<2>(1./3, 1./3));
data_file.write_dataset("point_value", point_value);
Vector<double> mean_value(1);
mean_value[0] = VectorTools::compute_mean_value(dof_handler,
                                                QGauss<2>(fe.degree + 1),
                                                solution, 0);
data_file.write_dataset("mean_value",mean_value);
@endcode
 
 
 
<a name="step_3-UsingRandggplot2togenerateplots"></a><h3> Using R and ggplot2 to generate plots</h3>
 
@note Alternatively, one could use the Python code in the next subsection.
 
The data put into %HDF5 files above can then be used from scripting
languages for further postprocessing. In the following, let us show
how this can, in particular, be done with the
<a href="https://en.wikipedia.org/wiki/R_(programming_language)">R
programming language</a>, a widely used language in statistical data
analysis. (Similar things can also be done in Python, for example.)
If you are unfamiliar with R and ggplot2 you could check out the data carpentry course on R
<a href="https://datacarpentry.org/R-ecology-lesson/index.html">here</a>.
Furthermore, since most search engines struggle with searches of the form "R + topic",
we recommend using the specializes service <a
href="http://rseek.org">RSeek </a> instead.
 
The most prominent difference between R and other languages is that
the assignment operator (`a = 5`) is typically written as
`a <- 5`. As the latter is considered standard we will use it in our examples as well.
To open the `.h5` file in R you have to install the <a href="https://bioconductor.org/packages/release/bioc/html/rhdf5.html">rhdf5</a> package, which is a part of the Bioconductor package.
 
First we will include all necessary packages and have a look at how the data is structured in our file.
@code{.r}
library(rhdf5)     # library for handling HDF5 files
library(ggplot2)   # main plotting library
library(grDevices) # needed for output to PDF
library(viridis)   # contains good colormaps for sequential data
 
refinement <- 5
h5f <- H5Fopen(paste("solution_",refinement,".h5",sep=""))
print(h5f)
@endcode
This gives the following output
@code{.unparsed}
HDF5 FILE
   name /
filename
 
    name       otype  dclass     dim
0 cells       H5I_DATASET INTEGER  x 1024
1 mean_value  H5I_DATASET FLOAT   1
2 nodes       H5I_DATASET FLOAT    x 1089
3 point_value H5I_DATASET FLOAT   1
4 solution    H5I_DATASET FLOAT    x 1089
@endcode
The datasets can be accessed by <code>h5f\$name</code>. The function
<code>dim(h5f\$cells)</code> gives us the dimensions of the matrix
that is used to store our cells.
We can see the following three matrices, as well as the two
additional data points we added.
<ul>
<li> <code>cells</code>: a 4x1024 matrix that stores the  (C++) vertex indices for each cell
<li> <code>nodes</code>: a 2x1089 matrix storing the position values (x,y) for our cell vertices
<li> <code>solution</code>: a 1x1089 matrix storing the values of our solution at each vertex
</ul>
Now we can use this data to generate various plots. Plotting with ggplot2 usually splits into two steps.
At first the data needs to be manipulated and added to a <code>data.frame</code>.
After that, a <code>ggplot</code> object is constructed and manipulated by adding plot elements to it.
 
<code>nodes</code> and <code>cells</code> contain all the information we need to plot our grid.
The following code wraps all the data into one dataframe for plotting our grid:
@code{.r}
# Counting in R starts at 1 instead of 0, so we need to increment all
# vertex indices by one:
cell_ids <- h5f@f$cells+1
 
# Store the x and y positions of each vertex in one big vector in a
# cell by cell fashion (every 4 entries belong to one cell):
cells_x <- h5f@f$nodes[1,][cell_ids]
cells_y <- h5f@f$nodes[2,][cell_ids]
 
# Construct a vector that stores the matching cell by cell grouping
# (1,1,1,1,2,2,2,2,...):
groups <- rep(1:ncol(cell_ids),each=4)
 
# Finally put everything into one dataframe:
meshdata <- data.frame(x = cells_x, y = cells_y, id = groups)
@endcode
 
With the finished dataframe we have everything we need to plot our grid:
@code{.r}
pdf (paste("grid_",refinement,".pdf",sep=""),width = 5,height = 5) # Open new PDF file
plt <- ggplot(meshdata,aes(x=x,y=y,group=id))                      # Construction of our plot
                                                                   # object, at first only data
 
plt <- plt + geom_polygon(fill="white",colour="black")             # Actual plotting of the grid as polygons
plt <- plt + ggtitle(paste("grid at refinement level #",refinement))
 
print(plt)                                                         # Show the current state of the plot/add it to the pdf
dev.off()                                                          # Close PDF file
@endcode
 
The contents of this file then look as follows (not very exciting, but
you get the idea):
<table width="60%" align="center">
  <tr>
   <td align="center">
     <img src="https://dealii.org/images/steps/developer/step-3.extensions.grid_5.png" alt="Grid after 5 refinement steps of step-3">
   </td>
  </tr>
</table>
 
We can also visualize the solution itself, and this is going to look
more interesting.
To make a 2D pseudocolor plot of our solution we will use <code>geom_raster</code>.
This function needs a structured grid, i.e. uniform in x and y directions.
Luckily our data at this point is structured in the right way.
The following code plots a pseudocolor representation of our surface into a new PDF:
@code{.r}
pdf (paste("pseudocolor_",refinement,".pdf",sep=""),width = 5,height = 4.2) # Open new PDF file
colordata <- data.frame(x = h5f@f$nodes[1,],y = h5f@f$nodes[2,] , solution = h5f@f$solution[1,])
plt <- ggplot(colordata,aes(x=x,y=y,fill=solution))
plt <- plt + geom_raster(interpolate=TRUE)
plt <- plt + scale_fill_viridis()
plt <- plt + ggtitle(paste("solution at refinement level #",refinement))
 
print(plt)
dev.off()
H5Fclose(h5f) # Close the HDF5 file
@endcode
This is now going to look as follows:
<table width="60%" align="center">
 <tr>
   <td align="center">
     <img src="https://dealii.org/images/steps/developer/step-3.extensions.pseudocolor_5.png" alt="Solution after 5 refinement steps of step-3">
   </td>
 </tr>
</table>
 
For plotting the convergence curves we need to re-run the C++ code multiple times with different values for <code>n_refinement_steps</code>
starting from 1.
Since every file only contains a single data point we need to loop over them and concatenate the results into a single vector.
@code{.r}
n_ref <- 8   # Maximum refinement level for which results are existing
 
# First we initiate all vectors with the results of the first level
h5f   <- H5Fopen("solution_1.h5")
dofs  <- dim(h5f@f$solution)[2]
mean  <- h5f@f$mean_value
point <- h5f@f$point_value
H5Fclose(h5f)
 
for (reflevel in 2:n_ref)
{
   h5f   <- H5Fopen(paste("solution_",reflevel,".h5",sep=""))
   dofs  <- c(dofs,dim(h5f\$solution)[2])
   mean  <- c(mean,h5f\$mean_value)
   point <- c(point,h5f\$point_value)
   H5Fclose(h5f)
}
@endcode
As we are not interested in the values themselves but rather in the error compared to a "exact" solution we will
assume our highest refinement level to be that solution and omit it from the data.
@code{.r}
# Calculate the error w.r.t. our maximum refinement step
mean_error  <- abs(mean[1:n_ref-1]-mean[n_ref])
point_error <- abs(point[1:n_ref-1]-point[n_ref])
 
# Remove the highest value from our DoF data
dofs     <- dofs[1:n_ref-1]
convdata <- data.frame(dofs = dofs, mean_value= mean_error, point_value = point_error)
@endcode
Now we have all the data available to generate our plots.
It is often useful to plot errors on a log-log scale, which is
accomplished in the following code:
@code
pdf (paste("convergence.pdf",sep=""),width = 5,height = 4.2)
plt <- ggplot(convdata,mapping=aes(x = dofs, y = mean_value))
plt <- plt+geom_line()
plt <- plt+labs(x="#DoFs",y = "mean value error")
plt <- plt+scale_x_log10()+scale_y_log10()
print(plt)
 
plt <- ggplot(convdata,mapping=aes(x = dofs, y = point_value))
plt <- plt+geom_line()
plt <- plt+labs(x="#DoFs",y = "point value error")
plt <- plt+scale_x_log10()+scale_y_log10()
print(plt)
 
dev.off()
@endcode
This results in the following plot that shows how the errors in the
mean value and the solution value at the chosen point nicely converge
to zero:
<table style="width:50%" align="center">
  <tr>
    <td><img src="https://dealii.org/images/steps/developer/step-3.extensions.convergence_mean.png" alt=""></td>
    <td><img src="https://dealii.org/images/steps/developer/step-3.extensions.convergence_point.png" alt=""></td>
  </tr>
</table>
 
<a name="step_3-UsingPythontogenerateplots"></a><h3> Using Python to generate plots</h3>
 
 
In this section we discuss the postprocessing of the data stored in %HDF5 files using the "python" programming language.
The necessary packages to import are
@code{.py}
import numpy as np                        # to work with multidimensional arrays
import h5py                               # to work with %HDF5 files
 
import pandas as pd                       # for data frames
import matplotlib.pyplot as plt           # plotting
from matplotlib.patches import Polygon
 
from scipy.interpolate import griddata    # interpolation function
from matplotlib import cm                 # for colormaps
 
@endcode
The %HDF5 solution file corresponding to `refinement = 5` can be opened as
@code{.py}
refinement = 5
filename = "solution_%d.h5" % refinement
file = h5py.File(filename, "r")
@endcode
The following prints out the tables in the %HDF5 file
@code{.py}
for key, value in file.items():
    print(key, " : ", value)
@endcode
which prints out
@code{.unparsed}
cells  :  <HDF5 dataset "cells": shape (1024, 4), type "<u4">
mean_value  :  <HDF5 dataset "mean_value": shape (1,), type "<f8">
nodes  :  <HDF5 dataset "nodes": shape (1089, 2), type "<f8">
point_value  :  <HDF5 dataset "point_value": shape (1,), type "<f8">
solution  :  <HDF5 dataset "solution": shape (1089, 1), type "<f8">
@endcode
There are @f$(32+1)\times(32+1) = 1089@f$ nodes.
The coordinates of these nodes @f$(x,y)@f$ are stored in the table `nodes` in the %HDF5 file.
There are a total of @f$32\times 32 = 1024@f$ cells. The nodes which make up each cell are
marked in the table `cells` in the %HDF5 file.
 
Next, we extract the data into multidimensional arrays
@code{.py}
nodes = np.array(file["/nodes"])
cells = np.array(file["/cells"])
solution = np.array(file["/solution"])
 
x, y = nodes.T
@endcode
 
The following stores the @f$x@f$ and @f$y@f$ coordinates of each node of each cell in one flat array.
@code{.py}
cell_x = x[cells.flatten()]
cell_y = y[cells.flatten()]
@endcode
The following tags the cell ids. Each four entries correspond to one cell.
Then we collect the coordinates and ids into a data frame
@code{.py}
n_cells = cells.shape[0]
cell_ids = np.repeat(np.arange(n_cells), 4)
meshdata = pd.DataFrame({"x": cell_x, "y": cell_y, "ids": cell_ids})
@endcode
The data frame looks
@code{.unparsed}
print(meshdata)
 
      x       y       ids
0       0.00000 0.00000 0
1       0.03125 0.00000 0
2       0.03125 0.03125 0
3       0.00000 0.03125 0
4       0.03125 0.00000 1
...     ...     ...     ...
4091    0.93750 1.00000 1022
4092    0.96875 0.96875 1023
4093    1.00000 0.96875 1023
4094    1.00000 1.00000 1023
4095    0.96875 1.00000 1023
 
4096 rows × 3 columns
@endcode
 
To plot the mesh, we loop over all cells and connect the first and last node to complete the polygon
@code{.py}
fig, ax = plt.subplots()
ax.set_aspect("equal", "box")
ax.set_title("grid at refinement level #%d" % refinement)
 
for cell_id, cell in meshdata.groupby(["ids"]):
    cell = pd.concat([cell, cell.head(1)])
    plt.plot(cell["x"], cell["y"], c="k")
@endcode
Alternatively one could use the matplotlib built-in Polygon class
@code{.py}
fig, ax = plt.subplots()
ax.set_aspect("equal", "box")
ax.set_title("grid at refinement level #%d" % refinement)
for cell_id, cell in meshdata.groupby(["ids"]):
    patch = Polygon(cell[["x", "y"]], facecolor="w", edgecolor="k")
    ax.add_patch(patch)
@endcode
 
To plot the solution, we first create a finer grid and then interpolate the solution values
into the grid and then plot it.
@code{.py}
nx = int(np.sqrt(n_cells)) + 1
nx *= 10
xg = np.linspace(x.min(), x.max(), nx)
yg = np.linspace(y.min(), y.max(), nx)
 
xgrid, ygrid = np.meshgrid(xg, yg)
solution_grid = griddata((x, y), solution.flatten(), (xgrid, ygrid), method="linear")
 
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_title("solution at refinement level #%d" % refinement)
c = ax.pcolor(xgrid, ygrid, solution_grid, cmap=cm.viridis)
fig.colorbar(c, ax=ax)
 
plt.show()
@endcode
 
To check the convergence of `mean_value` and `point_value`
we loop over data of all refinements and store into vectors <code>mean_values</code> and <code>mean_values</code>
@code{.py}
mean_values = np.zeros((8,))
point_values = np.zeros((8,))
dofs = np.zeros((8,))
 
for refinement in range(1, 9):
    filename = "solution_%d.h5" % refinement
    file = h5py.File(filename, "r")
    mean_values[refinement - 1] = np.array(file["/mean_value"])[0]
    point_values[refinement - 1] = np.array(file["/point_value"])[0]
    dofs[refinement - 1] = np.array(file["/solution"]).shape[0]
@endcode
 
Following is the plot of <code>mean_values</code> on `log-log` scale
@code{.py}
mean_error = np.abs(mean_values[1:] - mean_values[:-1])
plt.loglog(dofs[:-1], mean_error)
plt.grid()
plt.xlabel("#DoFs")
plt.ylabel("mean value error")
plt.show()
@endcode
 
Following is the plot of <code>point_values</code> on `log-log` scale
@code{.py}
point_error = np.abs(point_values[1:] - point_values[:-1])
plt.loglog(dofs[:-1], point_error)
plt.grid()
plt.xlabel("#DoFs")
plt.ylabel("point value error")
plt.show()
@endcode
 
A Python package which mimics the `R` ggplot2 (which is based on specifying the grammar of the graphics) is
<a href="https://plotnine.org/">plotnine</a>.
@code{.py}
We need to import the following from the <code>plotnine</code> package
from plotnine import (
    ggplot,
    aes,
    geom_raster,
    geom_polygon,
    geom_line,
    labs,
    scale_x_log10,
    scale_y_log10,
    ggtitle,
)
@endcode
Then plot the mesh <code>meshdata</code> dataframe
@code{.py}
plot = (
    ggplot(meshdata, aes(x="x", y="y", group="ids"))
    + geom_polygon(fill="white", colour="black")
    + ggtitle("grid at refinement level #%d" % refinement)
)
 
print(plot)
@endcode
Collect the solution into a dataframe
@code{.py}
colordata = pd.DataFrame({"x": x, "y": y, "solution": solution.flatten()})
@endcode
Plot of the solution
@code{.py}
plot = (
    ggplot(colordata, aes(x="x", y="y", fill="solution"))
    + geom_raster(interpolate=True)
    + ggtitle("solution at refinement level #%d" % refinement)
)
 
print(plot)
@endcode
 
Collect the convergence data into a dataframe
@code{.py}
convdata = pd.DataFrame(
    {"dofs": dofs[:-1], "mean_value": mean_error, "point_value": point_error}
)
 
@endcode
Following is the plot of <code>mean_values</code> on `log-log` scale
@code{.py}
plot = (
    ggplot(convdata, mapping=aes(x="dofs", y="mean_value"))
    + geom_line()
    + labs(x="#DoFs", y="mean value error")
    + scale_x_log10()
    + scale_y_log10()
)
 
plot.save("mean_error.pdf", dpi=60)
print(plot)
@endcode
 
Following is the plot of <code>point_values</code> on `log-log` scale
@code{.py}
plot = (
    ggplot(convdata, mapping=aes(x="dofs", y="point_value"))
    + geom_line()
    + labs(x="#DoFs", y="point value error")
    + scale_x_log10()
    + scale_y_log10()
)
 
plot.save("point_error.pdf", dpi=60)
print(plot)
@endcode
 *
 *
<a name="step_3-PlainProg"></a>
<h1> The plain program</h1>
@include "step-3.cc"
*/