step-61.h

Go to the documentation of this file.

,
 *                                     const unsigned int /*component*/) const
 *   {
 *     return 0;
 *   }
 * 
 * 
 * 
 *   template <int dim>
 *   class RightHandSide : public Function<dim>
 *   {
 *   public:
 *     virtual double value(const Point<dim> & p,
 *                          const unsigned int component = 0) const override;
 *   };
 * 
 * 
 * 
 *   template <int dim>
 *   double RightHandSide<dim>::value(const Point<dim> &p,
 *                                    const unsigned int /*component*/) const
 *   {
 *     return (2 * numbers::PI * numbers::PI * std::sin(numbers::PI * p[0]) *
 *             std::sin(numbers::PI * p[1]));
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * The class that implements the exact pressure solution has an
 * oddity in that we implement it as a vector-valued one with two
 * components. (We say that it has two components in the constructor
 * where we call the constructor of the base Function class.) In the
 * `value()` function, we do not test for the value of the
 * `component` argument, which implies that we return the same value
 * for both components of the vector-valued function. We do this
 * because we describe the finite element in use in this program as
 * a vector-valued system that contains the interior and the
 * interface pressures, and when we compute errors, we will want to
 * use the same pressure solution to test both of these components.
 * 
 * @code
 *   template <int dim>
 *   class ExactPressure : public Function<dim>
 *   {
 *   public:
 *     ExactPressure()
 *       : Function<dim>(2)
 *     {}
 * 
 *     virtual double value(const Point<dim> & p,
 *                          const unsigned int component) const override;
 *   };
 * 
 * 
 * 
 *   template <int dim>
 *   double ExactPressure<dim>::value(const Point<dim> &p,
 *                                    const unsigned int /*component*/) const
 *   {
 *     return std::sin(numbers::PI * p[0]) * std::sin(numbers::PI * p[1]);
 *   }
 * 
 * 
 * 
 *   template <int dim>
 *   class ExactVelocity : public TensorFunction<1, dim>
 *   {
 *   public:
 *     ExactVelocity()
 *       : TensorFunction<1, dim>()
 *     {}
 * 
 *     virtual Tensor<1, dim> value(const Point<dim> &p) const override;
 *   };
 * 
 * 
 * 
 *   template <int dim>
 *   Tensor<1, dim> ExactVelocity<dim>::value(const Point<dim> &p) const
 *   {
 *     Tensor<1, dim> return_value;
 *     return_value[0] = -numbers::PI * std::cos(numbers::PI * p[0]) *
 *                       std::sin(numbers::PI * p[1]);
 *     return_value[1] = -numbers::PI * std::sin(numbers::PI * p[0]) *
 *                       std::cos(numbers::PI * p[1]);
 *     return return_value;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationclassimplementation"></a> 
 * <h3>WGDarcyEquation class implementation</h3>
 * 
 
 * 
 * 
 * <a name="WGDarcyEquationWGDarcyEquation"></a> 
 * <h4>WGDarcyEquation::WGDarcyEquation</h4>
 * 
 
 * 
 * In this constructor, we create a finite element space for vector valued
 * functions, which will here include the ones used for the interior and
 * interface pressures, @f$p^\circ@f$ and @f$p^\partial@f$.
 * 
 * @code
 *   template <int dim>
 *   WGDarcyEquation<dim>::WGDarcyEquation(const unsigned int degree)
 *     : fe(FE_DGQ<dim>(degree), 1, FE_FaceQ<dim>(degree), 1)
 *     , dof_handler(triangulation)
 *     , fe_dgrt(degree)
 *     , dof_handler_dgrt(triangulation)
 *   {}
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationmake_grid"></a> 
 * <h4>WGDarcyEquation::make_grid</h4>
 * 
 
 * 
 * We generate a mesh on the unit square domain and refine it.
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::make_grid()
 *   {
 *     GridGenerator::hyper_cube(triangulation, 0, 1);
 *     triangulation.refine_global(5);
 * 
 *     std::cout << "   Number of active cells: " << triangulation.n_active_cells()
 *               << std::endl
 *               << "   Total number of cells: " << triangulation.n_cells()
 *               << std::endl;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationsetup_system"></a> 
 * <h4>WGDarcyEquation::setup_system</h4>
 * 
 
 * 
 * After we have created the mesh above, we distribute degrees of
 * freedom and resize matrices and vectors. The only piece of
 * interest in this function is how we interpolate the boundary
 * values for the pressure. Since the pressure consists of interior
 * and interface components, we need to make sure that we only
 * interpolate onto that component of the vector-valued solution
 * space that corresponds to the interface pressures (as these are
 * the only ones that are defined on the boundary of the domain). We
 * do this via a component mask object for only the interface
 * pressures.
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe);
 *     dof_handler_dgrt.distribute_dofs(fe_dgrt);
 * 
 *     std::cout << "   Number of pressure degrees of freedom: "
 *               << dof_handler.n_dofs() << std::endl;
 * 
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 * 
 * 
 *     {
 *       constraints.clear();
 *       const FEValuesExtractors::Scalar interface_pressure(1);
 *       const ComponentMask              interface_pressure_mask =
 *         fe.component_mask(interface_pressure);
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                BoundaryValues<dim>(),
 *                                                constraints,
 *                                                interface_pressure_mask);
 *       constraints.close();
 *     }
 * 
 * 
 * @endcode
 * 
 * In the bilinear form, there is no integration term over faces
 * between two neighboring cells, so we can just use
 * <code>DoFTools::make_sparsity_pattern</code> to calculate the sparse
 * matrix.
 * 
 * @code
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
 *     sparsity_pattern.copy_from(dsp);
 * 
 *     system_matrix.reinit(sparsity_pattern);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationassemble_system"></a> 
 * <h4>WGDarcyEquation::assemble_system</h4>
 * 
 
 * 
 * This function is more interesting. As detailed in the
 * introduction, the assembly of the linear system requires us to
 * evaluate the weak gradient of the shape functions, which is an
 * element in the Raviart-Thomas space. As a consequence, we need to
 * define a Raviart-Thomas finite element object, and have FEValues
 * objects that evaluate it at quadrature points. We then need to
 * compute the matrix @f$C^K@f$ on every cell @f$K@f$, for which we need the
 * matrices @f$M^K@f$ and @f$G^K@f$ mentioned in the introduction.
 *   
 
 * 
 * A point that may not be obvious is that in all previous tutorial
 * programs, we have always called FEValues::reinit() with a cell
 * iterator from a DoFHandler. This is so that one can call
 * functions such as FEValuesBase::get_function_values() that
 * extract the values of a finite element function (represented by a
 * vector of DoF values) on the quadrature points of a cell. For
 * this operation to work, one needs to know which vector elements
 * correspond to the degrees of freedom on a given cell -- i.e.,
 * exactly the kind of information and operation provided by the
 * DoFHandler class.
 *   
 
 * 
 * We could create a DoFHandler object for the "broken" Raviart-Thomas space
 * (using the FE_DGRT class), but we really don't want to here: At
 * least in the current function, we have no need for any globally defined
 * degrees of freedom associated with this broken space, but really only
 * need to reference the shape functions of such a space on the current
 * cell. As a consequence, we use the fact that one can call
 * FEValues::reinit() also with cell iterators into Triangulation
 * objects (rather than DoFHandler objects). In this case, FEValues
 * can of course only provide us with information that only
 * references information about cells, rather than degrees of freedom
 * enumerated on these cells. So we can't use
 * FEValuesBase::get_function_values(), but we can use
 * FEValues::shape_value() to obtain the values of shape functions
 * at quadrature points on the current cell. It is this kind of
 * functionality we will make use of below. The variable that will
 * give us this information about the Raviart-Thomas functions below
 * is then the `fe_values_rt` (and corresponding `fe_face_values_rt`)
 * object.
 *   
 
 * 
 * Given this introduction, the following declarations should be
 * pretty obvious:
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::assemble_system()
 *   {
 *     const QGauss<dim>     quadrature_formula(fe_dgrt.degree + 1);
 *     const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
 * 
 *     FEValues<dim>     fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_quadrature_points |
 *                               update_JxW_values);
 *     FEFaceValues<dim> fe_face_values(fe,
 *                                      face_quadrature_formula,
 *                                      update_values | update_normal_vectors |
 *                                        update_quadrature_points |
 *                                        update_JxW_values);
 * 
 *     FEValues<dim>     fe_values_dgrt(fe_dgrt,
 *                                  quadrature_formula,
 *                                  update_values | update_gradients |
 *                                    update_quadrature_points |
 *                                    update_JxW_values);
 *     FEFaceValues<dim> fe_face_values_dgrt(fe_dgrt,
 *                                           face_quadrature_formula,
 *                                           update_values |
 *                                             update_normal_vectors |
 *                                             update_quadrature_points |
 *                                             update_JxW_values);
 * 
 *     const unsigned int dofs_per_cell      = fe.n_dofs_per_cell();
 *     const unsigned int dofs_per_cell_dgrt = fe_dgrt.n_dofs_per_cell();
 * 
 *     const unsigned int n_q_points      = fe_values.get_quadrature().size();
 *     const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
 * 
 *     const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
 * 
 *     RightHandSide<dim>  right_hand_side;
 *     std::vector<double> right_hand_side_values(n_q_points);
 * 
 *     const Coefficient<dim>      coefficient;
 *     std::vector<Tensor<2, dim>> coefficient_values(n_q_points);
 * 
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 * 
 * 
 * @endcode
 * 
 * Next, let us declare the various cell matrices discussed in the
 * introduction:
 * 
 * @code
 *     FullMatrix<double> cell_matrix_M(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
 *     FullMatrix<double> cell_matrix_G(dofs_per_cell_dgrt, dofs_per_cell);
 *     FullMatrix<double> cell_matrix_C(dofs_per_cell, dofs_per_cell_dgrt);
 *     FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
 *     Vector<double>     cell_rhs(dofs_per_cell);
 *     Vector<double>     cell_solution(dofs_per_cell);
 * 
 * @endcode
 * 
 * We need <code>FEValuesExtractors</code> to access the @p interior and
 * @p face component of the shape functions.
 * 
 * @code
 *     const FEValuesExtractors::Vector velocities(0);
 *     const FEValuesExtractors::Scalar pressure_interior(0);
 *     const FEValuesExtractors::Scalar pressure_face(1);
 * 
 * @endcode
 * 
 * This finally gets us in position to loop over all cells. On
 * each cell, we will first calculate the various cell matrices
 * used to construct the local matrix -- as they depend on the
 * cell in question, they need to be re-computed on each cell. We
 * need shape functions for the Raviart-Thomas space as well, for
 * which we need to create first an iterator to the cell of the
 * triangulation, which we can obtain by assignment from the cell
 * pointing into the DoFHandler.
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         fe_values.reinit(cell);
 * 
 *         const typename Triangulation<dim>::active_cell_iterator cell_dgrt =
 *           cell;
 *         fe_values_dgrt.reinit(cell_dgrt);
 * 
 *         right_hand_side.value_list(fe_values.get_quadrature_points(),
 *                                    right_hand_side_values);
 *         coefficient.value_list(fe_values.get_quadrature_points(),
 *                                coefficient_values);
 * 
 * @endcode
 * 
 * The first cell matrix we will compute is the mass matrix
 * for the Raviart-Thomas space.  Hence, we need to loop over
 * all the quadrature points for the velocity FEValues object.
 * 
 * @code
 *         cell_matrix_M = 0;
 *         for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
 *           for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
 *             {
 *               const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
 *               for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
 *                 {
 *                   const Tensor<1, dim> v_k =
 *                     fe_values_dgrt[velocities].value(k, q);
 *                   cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
 *                 }
 *             }
 * @endcode
 * 
 * Next we take the inverse of this matrix by using
 * FullMatrix::gauss_jordan(). It will be used to calculate
 * the coefficient matrix @f$C^K@f$ later. It is worth recalling
 * later that `cell_matrix_M` actually contains the *inverse*
 * of @f$M^K@f$ after this call.
 * 
 * @code
 *         cell_matrix_M.gauss_jordan();
 * 
 * @endcode
 * 
 * From the introduction, we know that the right hand side
 * @f$G^K@f$ of the equation that defines @f$C^K@f$ is the difference
 * between a face integral and a cell integral. Here, we
 * approximate the negative of the contribution in the
 * interior. Each component of this matrix is the integral of
 * a product between a basis function of the polynomial space
 * and the divergence of a basis function of the
 * Raviart-Thomas space. These basis functions are defined in
 * the interior.
 * 
 * @code
 *         cell_matrix_G = 0;
 *         for (unsigned int q = 0; q < n_q_points; ++q)
 *           for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
 *             {
 *               const double div_v_i =
 *                 fe_values_dgrt[velocities].divergence(i, q);
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 {
 *                   const double phi_j_interior =
 *                     fe_values[pressure_interior].value(j, q);
 * 
 *                   cell_matrix_G(i, j) -=
 *                     (div_v_i * phi_j_interior * fe_values.JxW(q));
 *                 }
 *             }
 * 
 * 
 * @endcode
 * 
 * Next, we approximate the integral on faces by quadrature.
 * Each component is the integral of a product between a basis function
 * of the polynomial space and the dot product of a basis function of
 * the Raviart-Thomas space and the normal vector. So we loop over all
 * the faces of the element and obtain the normal vector.
 * 
 * @code
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             fe_face_values.reinit(cell, face);
 *             fe_face_values_dgrt.reinit(cell_dgrt, face);
 * 
 *             for (unsigned int q = 0; q < n_face_q_points; ++q)
 *               {
 *                 const Tensor<1, dim> &normal = fe_face_values.normal_vector(q);
 * 
 *                 for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
 *                   {
 *                     const Tensor<1, dim> v_i =
 *                       fe_face_values_dgrt[velocities].value(i, q);
 *                     for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                       {
 *                         const double phi_j_face =
 *                           fe_face_values[pressure_face].value(j, q);
 * 
 *                         cell_matrix_G(i, j) +=
 *                           ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
 *                       }
 *                   }
 *               }
 *           }
 * 
 * @endcode
 * 
 * @p cell_matrix_C is then the matrix product between the
 * transpose of @f$G^K@f$ and the inverse of the mass matrix
 * (where this inverse is stored in @p cell_matrix_M):
 * 
 * @code
 *         cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
 * 
 * @endcode
 * 
 * Finally we can compute the local matrix @f$A^K@f$.  Element
 * @f$A^K_{ij}@f$ is given by @f$\int_{E} \sum_{k,l} C_{ik} C_{jl}
 * (\mathbf{K} \mathbf{v}_k) \cdot \mathbf{v}_l
 * \mathrm{d}x@f$. We have calculated the coefficients @f$C@f$ in
 * the previous step, and so obtain the following after
 * suitably re-arranging the loops:
 * 
 * @code
 *         local_matrix = 0;
 *         for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
 *           {
 *             for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
 *               {
 *                 const Tensor<1, dim> v_k =
 *                   fe_values_dgrt[velocities].value(k, q);
 *                 for (unsigned int l = 0; l < dofs_per_cell_dgrt; ++l)
 *                   {
 *                     const Tensor<1, dim> v_l =
 *                       fe_values_dgrt[velocities].value(l, q);
 * 
 *                     for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                       for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                         local_matrix(i, j) +=
 *                           (coefficient_values[q] * cell_matrix_C[i][k] * v_k) *
 *                           cell_matrix_C[j][l] * v_l * fe_values_dgrt.JxW(q);
 *                   }
 *               }
 *           }
 * 
 * @endcode
 * 
 * Next, we calculate the right hand side, @f$\int_{K} f q \mathrm{d}x@f$:
 * 
 * @code
 *         cell_rhs = 0;
 *         for (unsigned int q = 0; q < n_q_points; ++q)
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               cell_rhs(i) += (fe_values[pressure_interior].value(i, q) *
 *                               right_hand_side_values[q] * fe_values.JxW(q));
 *             }
 * 
 * @endcode
 * 
 * The last step is to distribute components of the local
 * matrix into the system matrix and transfer components of
 * the cell right hand side into the system right hand side:
 * 
 * @code
 *         cell->get_dof_indices(local_dof_indices);
 *         constraints.distribute_local_to_global(
 *           local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
 *       }
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationdimsolve"></a> 
 * <h4>WGDarcyEquation<dim>::solve</h4>
 * 
 
 * 
 * This step is rather trivial and the same as in many previous
 * tutorial programs:
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::solve()
 *   {
 *     SolverControl            solver_control(1000, 1e-8 * system_rhs.l2_norm());
 *     SolverCG<Vector<double>> solver(solver_control);
 *     solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
 *     constraints.distribute(solution);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationdimcompute_postprocessed_velocity"></a> 
 * <h4>WGDarcyEquation<dim>::compute_postprocessed_velocity</h4>
 * 
 
 * 
 * In this function, compute the velocity field from the pressure
 * solution previously computed. The
 * velocity is defined as @f$\mathbf{u}_h = \mathbf{Q}_h \left(
 * -\mathbf{K}\nabla_{w,d}p_h \right)@f$, which requires us to compute
 * many of the same terms as in the assembly of the system matrix.
 * There are also the matrices @f$E^K,D^K@f$ we need to assemble (see
 * the introduction) but they really just follow the same kind of
 * pattern.
 *   
 
 * 
 * Computing the same matrices here as we have already done in the
 * `assemble_system()` function is of course wasteful in terms of
 * CPU time. Likewise, we copy some of the code from there to this
 * function, and this is also generally a poor idea. A better
 * implementation might provide for a function that encapsulates
 * this duplicated code. One could also think of using the classic
 * trade-off between computing efficiency and memory efficiency to
 * only compute the @f$C^K@f$ matrices once per cell during the
 * assembly, storing them somewhere on the side, and re-using them
 * here. (This is what @ref step_51 "step-51" does, for example, where the
 * `assemble_system()` function takes an argument that determines
 * whether the local matrices are recomputed, and a similar approach
 * -- maybe with storing local matrices elsewhere -- could be
 * adapted for the current program.)
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::compute_postprocessed_velocity()
 *   {
 *     darcy_velocity.reinit(dof_handler_dgrt.n_dofs());
 * 
 *     const QGauss<dim>     quadrature_formula(fe_dgrt.degree + 1);
 *     const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
 * 
 *     FEValues<dim> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_quadrature_points |
 *                               update_JxW_values);
 * 
 *     FEFaceValues<dim> fe_face_values(fe,
 *                                      face_quadrature_formula,
 *                                      update_values | update_normal_vectors |
 *                                        update_quadrature_points |
 *                                        update_JxW_values);
 * 
 *     FEValues<dim> fe_values_dgrt(fe_dgrt,
 *                                  quadrature_formula,
 *                                  update_values | update_gradients |
 *                                    update_quadrature_points |
 *                                    update_JxW_values);
 * 
 *     FEFaceValues<dim> fe_face_values_dgrt(fe_dgrt,
 *                                           face_quadrature_formula,
 *                                           update_values |
 *                                             update_normal_vectors |
 *                                             update_quadrature_points |
 *                                             update_JxW_values);
 * 
 *     const unsigned int dofs_per_cell      = fe.n_dofs_per_cell();
 *     const unsigned int dofs_per_cell_dgrt = fe_dgrt.n_dofs_per_cell();
 * 
 *     const unsigned int n_q_points      = fe_values.get_quadrature().size();
 *     const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
 * 
 *     const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
 * 
 * 
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *     std::vector<types::global_dof_index> local_dof_indices_dgrt(
 *       dofs_per_cell_dgrt);
 * 
 *     FullMatrix<double> cell_matrix_M(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
 *     FullMatrix<double> cell_matrix_G(dofs_per_cell_dgrt, dofs_per_cell);
 *     FullMatrix<double> cell_matrix_C(dofs_per_cell, dofs_per_cell_dgrt);
 *     FullMatrix<double> cell_matrix_D(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
 *     FullMatrix<double> cell_matrix_E(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
 * 
 *     Vector<double> cell_solution(dofs_per_cell);
 *     Vector<double> cell_velocity(dofs_per_cell_dgrt);
 * 
 *     const Coefficient<dim>      coefficient;
 *     std::vector<Tensor<2, dim>> coefficient_values(n_q_points_dgrt);
 * 
 *     const FEValuesExtractors::Vector velocities(0);
 *     const FEValuesExtractors::Scalar pressure_interior(0);
 *     const FEValuesExtractors::Scalar pressure_face(1);
 * 
 * @endcode
 * 
 * In the introduction, we explained how to calculate the numerical velocity
 * on the cell. We need the pressure solution values on each cell,
 * coefficients of the Gram matrix and coefficients of the @f$L_2@f$ projection.
 * We have already calculated the global solution, so we will extract the
 * cell solution from the global solution. The coefficients of the Gram
 * matrix have been calculated when we assembled the system matrix for the
 * pressures. We will do the same way here. For the coefficients of the
 * projection, we do matrix multiplication, i.e., the inverse of the Gram
 * matrix times the matrix with @f$(\mathbf{K} \mathbf{w}, \mathbf{w})@f$ as
 * components. Then, we multiply all these coefficients and call them beta.
 * The numerical velocity is the product of beta and the basis functions of
 * the Raviart-Thomas space.
 * 
 * @code
 *     typename DoFHandler<dim>::active_cell_iterator
 *       cell = dof_handler.begin_active(),
 *       endc = dof_handler.end(), cell_dgrt = dof_handler_dgrt.begin_active();
 *     for (; cell != endc; ++cell, ++cell_dgrt)
 *       {
 *         fe_values.reinit(cell);
 *         fe_values_dgrt.reinit(cell_dgrt);
 * 
 *         coefficient.value_list(fe_values_dgrt.get_quadrature_points(),
 *                                coefficient_values);
 * 
 * @endcode
 * 
 * The component of this <code>cell_matrix_E</code> is the integral of
 * @f$(\mathbf{K} \mathbf{w}, \mathbf{w})@f$. <code>cell_matrix_M</code> is
 * the Gram matrix.
 * 
 * @code
 *         cell_matrix_M = 0;
 *         cell_matrix_E = 0;
 *         for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
 *           for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
 *             {
 *               const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
 *               for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
 *                 {
 *                   const Tensor<1, dim> v_k =
 *                     fe_values_dgrt[velocities].value(k, q);
 * 
 *                   cell_matrix_E(i, k) +=
 *                     (coefficient_values[q] * v_i * v_k * fe_values_dgrt.JxW(q));
 * 
 *                   cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
 *                 }
 *             }
 * 
 * @endcode
 * 
 * To compute the matrix @f$D@f$ mentioned in the introduction, we
 * then need to evaluate @f$D=M^{-1}E@f$ as explained in the
 * introduction:
 * 
 * @code
 *         cell_matrix_M.gauss_jordan();
 *         cell_matrix_M.mmult(cell_matrix_D, cell_matrix_E);
 * 
 * @endcode
 * 
 * Then we also need, again, to compute the matrix @f$C@f$ that is
 * used to evaluate the weak discrete gradient. This is the
 * exact same code as used in the assembly of the system
 * matrix, so we just copy it from there:
 * 
 * @code
 *         cell_matrix_G = 0;
 *         for (unsigned int q = 0; q < n_q_points; ++q)
 *           for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
 *             {
 *               const double div_v_i =
 *                 fe_values_dgrt[velocities].divergence(i, q);
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 {
 *                   const double phi_j_interior =
 *                     fe_values[pressure_interior].value(j, q);
 * 
 *                   cell_matrix_G(i, j) -=
 *                     (div_v_i * phi_j_interior * fe_values.JxW(q));
 *                 }
 *             }
 * 
 *         for (const auto &face : cell->face_iterators())
 *           {
 *             fe_face_values.reinit(cell, face);
 *             fe_face_values_dgrt.reinit(cell_dgrt, face);
 * 
 *             for (unsigned int q = 0; q < n_face_q_points; ++q)
 *               {
 *                 const Tensor<1, dim> &normal = fe_face_values.normal_vector(q);
 * 
 *                 for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
 *                   {
 *                     const Tensor<1, dim> v_i =
 *                       fe_face_values_dgrt[velocities].value(i, q);
 *                     for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                       {
 *                         const double phi_j_face =
 *                           fe_face_values[pressure_face].value(j, q);
 * 
 *                         cell_matrix_G(i, j) +=
 *                           ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
 *                       }
 *                   }
 *               }
 *           }
 *         cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
 * 
 * @endcode
 * 
 * Finally, we need to extract the pressure unknowns that
 * correspond to the current cell:
 * 
 * @code
 *         cell->get_dof_values(solution, cell_solution);
 * 
 * @endcode
 * 
 * We are now in a position to compute the local velocity
 * unknowns (with respect to the Raviart-Thomas space we are
 * projecting the term @f$-\mathbf K \nabla_{w,d} p_h@f$ into):
 * 
 * @code
 *         cell_velocity = 0;
 *         for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
 *           for (unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               cell_velocity(k) +=
 *                 -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
 * 
 * @endcode
 * 
 * We compute Darcy velocity.
 * This is same as cell_velocity but used to graph Darcy velocity.
 * 
 * @code
 *         cell_dgrt->get_dof_indices(local_dof_indices_dgrt);
 *         for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
 *           for (unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               darcy_velocity(local_dof_indices_dgrt[k]) +=
 *                 -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
 *       }
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationdimcompute_pressure_error"></a> 
 * <h4>WGDarcyEquation<dim>::compute_pressure_error</h4>
 * 
 
 * 
 * This part is to calculate the @f$L_2@f$ error of the pressure.  We
 * define a vector that holds the norm of the error on each cell.
 * Next, we use VectorTool::integrate_difference() to compute the
 * error in the @f$L_2@f$ norm on each cell. However, we really only
 * care about the error in the interior component of the solution
 * vector (we can't even evaluate the interface pressures at the
 * quadrature points because these are all located in the interior
 * of cells) and consequently have to use a weight function that
 * ensures that the interface component of the solution variable is
 * ignored. This is done by using the ComponentSelectFunction whose
 * arguments indicate which component we want to select (component
 * zero, i.e., the interior pressures) and how many components there
 * are in total (two).
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::compute_pressure_error()
 *   {
 *     Vector<float> difference_per_cell(triangulation.n_active_cells());
 *     const ComponentSelectFunction<dim> select_interior_pressure(0, 2);
 *     VectorTools::integrate_difference(dof_handler,
 *                                       solution,
 *                                       ExactPressure<dim>(),
 *                                       difference_per_cell,
 *                                       QGauss<dim>(fe.degree + 2),
 *                                       VectorTools::L2_norm,
 *                                       &select_interior_pressure);
 * 
 *     const double L2_error = difference_per_cell.l2_norm();
 *     std::cout << "L2_error_pressure " << L2_error << std::endl;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationdimcompute_velocity_error"></a> 
 * <h4>WGDarcyEquation<dim>::compute_velocity_error</h4>
 * 
 
 * 
 * In this function, we evaluate @f$L_2@f$ errors for the velocity on
 * each cell, and @f$L_2@f$ errors for the flux on faces. The function
 * relies on the `compute_postprocessed_velocity()` function having
 * previous computed, which computes the velocity field based on the
 * pressure solution that has previously been computed.
 *   
 
 * 
 * We are going to evaluate velocities on each cell and calculate
 * the difference between numerical and exact velocities.
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::compute_velocity_errors()
 *   {
 *     const QGauss<dim>     quadrature_formula(fe_dgrt.degree + 1);
 *     const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
 * 
 *     FEValues<dim> fe_values_dgrt(fe_dgrt,
 *                                  quadrature_formula,
 *                                  update_values | update_gradients |
 *                                    update_quadrature_points |
 *                                    update_JxW_values);
 * 
 *     FEFaceValues<dim> fe_face_values_dgrt(fe_dgrt,
 *                                           face_quadrature_formula,
 *                                           update_values |
 *                                             update_normal_vectors |
 *                                             update_quadrature_points |
 *                                             update_JxW_values);
 * 
 *     const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
 *     const unsigned int n_face_q_points_dgrt =
 *       fe_face_values_dgrt.get_quadrature().size();
 * 
 *     std::vector<Tensor<1, dim>> velocity_values(n_q_points_dgrt);
 *     std::vector<Tensor<1, dim>> velocity_face_values(n_face_q_points_dgrt);
 * 
 *     const FEValuesExtractors::Vector velocities(0);
 * 
 *     const ExactVelocity<dim> exact_velocity;
 * 
 *     double L2_err_velocity_cell_sqr_global = 0;
 *     double L2_err_flux_sqr                 = 0;
 * 
 * @endcode
 * 
 * Having previously computed the postprocessed velocity, we here
 * only have to extract the corresponding values on each cell and
 * face and compare it to the exact values.
 * 
 * @code
 *     for (const auto &cell_dgrt : dof_handler_dgrt.active_cell_iterators())
 *       {
 *         fe_values_dgrt.reinit(cell_dgrt);
 * 
 * @endcode
 * 
 * First compute the @f$L_2@f$ error between the postprocessed velocity
 * field and the exact one:
 * 
 * @code
 *         fe_values_dgrt[velocities].get_function_values(darcy_velocity,
 *                                                        velocity_values);
 *         double L2_err_velocity_cell_sqr_local = 0;
 *         for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
 *           {
 *             const Tensor<1, dim> velocity = velocity_values[q];
 *             const Tensor<1, dim> true_velocity =
 *               exact_velocity.value(fe_values_dgrt.quadrature_point(q));
 * 
 *             L2_err_velocity_cell_sqr_local +=
 *               ((velocity - true_velocity) * (velocity - true_velocity) *
 *                fe_values_dgrt.JxW(q));
 *           }
 *         L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local;
 * 
 * @endcode
 * 
 * For reconstructing the flux we need the size of cells and
 * faces. Since fluxes are calculated on faces, we have the
 * loop over all four faces of each cell. To calculate the
 * face velocity, we extract values at the quadrature points from the
 * `darcy_velocity` which we have computed previously. Then, we
 * calculate the squared velocity error in normal direction. Finally, we
 * calculate the @f$L_2@f$ flux error on the cell by appropriately scaling
 * with face and cell areas and add it to the global error.
 * 
 * @code
 *         const double cell_area = cell_dgrt->measure();
 *         for (const auto &face_dgrt : cell_dgrt->face_iterators())
 *           {
 *             const double face_length = face_dgrt->measure();
 *             fe_face_values_dgrt.reinit(cell_dgrt, face_dgrt);
 *             fe_face_values_dgrt[velocities].get_function_values(
 *               darcy_velocity, velocity_face_values);
 * 
 *             double L2_err_flux_face_sqr_local = 0;
 *             for (unsigned int q = 0; q < n_face_q_points_dgrt; ++q)
 *               {
 *                 const Tensor<1, dim> velocity = velocity_face_values[q];
 *                 const Tensor<1, dim> true_velocity =
 *                   exact_velocity.value(fe_face_values_dgrt.quadrature_point(q));
 * 
 *                 const Tensor<1, dim> &normal =
 *                   fe_face_values_dgrt.normal_vector(q);
 * 
 *                 L2_err_flux_face_sqr_local +=
 *                   ((velocity * normal - true_velocity * normal) *
 *                    (velocity * normal - true_velocity * normal) *
 *                    fe_face_values_dgrt.JxW(q));
 *               }
 *             const double err_flux_each_face =
 *               L2_err_flux_face_sqr_local / face_length * cell_area;
 *             L2_err_flux_sqr += err_flux_each_face;
 *           }
 *       }
 * 
 * @endcode
 * 
 * After adding up errors over all cells and faces, we take the
 * square root and get the @f$L_2@f$ errors of velocity and
 * flux. These we output to screen.
 * 
 * @code
 *     const double L2_err_velocity_cell =
 *       std::sqrt(L2_err_velocity_cell_sqr_global);
 *     const double L2_err_flux_face = std::sqrt(L2_err_flux_sqr);
 * 
 *     std::cout << "L2_error_vel:  " << L2_err_velocity_cell << std::endl
 *               << "L2_error_flux: " << L2_err_flux_face << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationoutput_results"></a> 
 * <h4>WGDarcyEquation::output_results</h4>
 * 
 
 * 
 * We have two sets of results to output: the interior solution and
 * the skeleton solution. We use <code>DataOut</code> to visualize
 * interior results. The graphical output for the skeleton results
 * is done by using the DataOutFaces class.
 *   
 
 * 
 * In both of the output files, both the interior and the face
 * variables are stored. For the interface output, the output file
 * simply contains the interpolation of the interior pressures onto
 * the faces, but because it is undefined which of the two interior
 * pressure variables you get from the two adjacent cells, it is
 * best to ignore the interior pressure in the interface output
 * file. Conversely, for the cell interior output file, it is of
 * course impossible to show any interface pressures @f$p^\partial@f$,
 * because these are only available on interfaces and not cell
 * interiors. Consequently, you will see them shown as an invalid
 * value (such as an infinity).
 *   
 
 * 
 * For the cell interior output, we also want to output the velocity
 * variables. This is a bit tricky since it lives on the same mesh
 * but uses a different DoFHandler object (the pressure variables live
 * on the `dof_handler` object, the Darcy velocity on the `dof_handler_dgrt`
 * object). Fortunately, there are variations of the
 * DataOut::add_data_vector() function that allow specifying which
 * DoFHandler a vector corresponds to, and consequently we can visualize
 * the data from both DoFHandler objects within the same file.
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::output_results() const
 *   {
 *     {
 *       DataOut<dim> data_out;
 * 
 * @endcode
 * 
 * First attach the pressure solution to the DataOut object:
 * 
 * @code
 *       const std::vector<std::string> solution_names = {"interior_pressure",
 *                                                        "interface_pressure"};
 *       data_out.add_data_vector(dof_handler, solution, solution_names);
 * 
 * @endcode
 * 
 * Then do the same with the Darcy velocity field, and continue
 * with writing everything out into a file.
 * 
 * @code
 *       const std::vector<std::string> velocity_names(dim, "velocity");
 *       const std::vector<
 *         DataComponentInterpretation::DataComponentInterpretation>
 *         velocity_component_interpretation(
 *           dim, DataComponentInterpretation::component_is_part_of_vector);
 *       data_out.add_data_vector(dof_handler_dgrt,
 *                                darcy_velocity,
 *                                velocity_names,
 *                                velocity_component_interpretation);
 * 
 *       data_out.build_patches(fe.degree);
 *       std::ofstream output("solution_interior.vtu");
 *       data_out.write_vtu(output);
 *     }
 * 
 *     {
 *       DataOutFaces<dim> data_out_faces(false);
 *       data_out_faces.attach_dof_handler(dof_handler);
 *       data_out_faces.add_data_vector(solution, "Pressure_Face");
 *       data_out_faces.build_patches(fe.degree);
 *       std::ofstream face_output("solution_interface.vtu");
 *       data_out_faces.write_vtu(face_output);
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="WGDarcyEquationrun"></a> 
 * <h4>WGDarcyEquation::run</h4>
 * 
 
 * 
 * This is the final function of the main class. It calls the other functions
 * of our class.
 * 
 * @code
 *   template <int dim>
 *   void WGDarcyEquation<dim>::run()
 *   {
 *     std::cout << "Solving problem in " << dim << " space dimensions."
 *               << std::endl;
 *     make_grid();
 *     setup_system();
 *     assemble_system();
 *     solve();
 *     compute_postprocessed_velocity();
 *     compute_pressure_error();
 *     compute_velocity_errors();
 *     output_results();
 *   }
 * 
 * } // namespace Step61
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * This is the main function. We can change the dimension here to run in 3d.
 * 
 * @code
 * int main()
 * {
 *   try
 *     {
 *       Step61::WGDarcyEquation<2> wg_darcy(0);
 *       wg_darcy.run();
 *     }
 *   catch (std::exception &exc)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Exception on processing: " << std::endl
 *                 << exc.what() << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 *   catch (...)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Unknown exception!" << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 * 
 *   return 0;
 * }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
We run the program with a right hand side that will produce the
solution @f$p = \sin(\pi x) \sin(\pi y)@f$ and with homogeneous Dirichlet
boundary conditions in the domain @f$\Omega = (0,1)^2@f$. In addition, we
choose the coefficient matrix in the differential operator
@f$\mathbf{K}@f$ as the identity matrix. We test this setup using
@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ and
@f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ element combinations, which one can
select by using the appropriate constructor argument for the
`WGDarcyEquation` object in `main()`. We will then visualize pressure
values in interiors of cells and on faces. We want to see that the
pressure maximum is around 1 and the minimum is around 0. With mesh
refinement, the convergence rates of pressure, velocity and flux
should then be around 1 for @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ , 2 for
@f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$, and 3 for
@f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$.
 
 
<a name="TestresultsoniWGQsub0subQsub0subRTsub0subi"></a><h3>Test results on <i>WG(Q<sub>0</sub>,Q<sub>0</sub>;RT<sub>[0]</sub>)</i></h3>
 
 
The following figures show interior pressures and face pressures using the
@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ element. The mesh is refined 2 times (top)
and 4 times (bottom), respectively. (This number can be adjusted in the
`make_grid()` function.) When the mesh is coarse, one can see
the face pressures @f$p^\partial@f$ neatly between the values of the interior
pressures @f$p^\circ@f$ on the two adjacent cells.
 
<table align="center">
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_2d_2.png" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_3d_2.png" alt=""></td>
  </tr>
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_2d_4.png" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_3d_4.png" alt=""></td>
  </tr>
</table>
 
From the figures, we can see that with the mesh refinement, the maximum and
minimum pressure values are approaching the values we expect.
Since the mesh is a rectangular mesh and numbers of cells in each direction is even, we
have symmetric solutions. From the 3d figures on the right,
we can see that on @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, the pressure is a constant
in the interior of the cell, as expected.
 
<a name="Convergencetableforik0i"></a><h4>Convergence table for <i>k=0</i></h4>
 
 
We run the code with differently refined meshes (chosen in the `make_grid()` function)
and get the following convergence rates of pressure,
velocity, and flux (as defined in the introduction).
 
<table align="center" class="doxtable">
  <tr>
   <th>number of refinements </th><th>  @f$\|p-p_h^\circ\|@f$  </th><th>  @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
  </tr>
  <tr>
   <td>   2                  </td><td>    1.587e-01        </td><td>        5.113e-01               </td><td>   7.062e-01 </td>
  </tr>
  <tr>
   <td>   3                  </td><td>    8.000e-02        </td><td>        2.529e-01               </td><td>   3.554e-01 </td>
  </tr>
  <tr>
   <td>   4                  </td><td>    4.006e-02        </td><td>        1.260e-01               </td><td>   1.780e-01 </td>
  </tr>
  <tr>
   <td>   5                  </td><td>    2.004e-02        </td><td>        6.297e-02               </td><td>   8.902e-02 </td>
  </tr>
  <tr>
   <th>Conv.rate             </th><th>      1.00           </th><th>          1.00                  </th><th>      1.00   </th>
  </tr>
</table>
 
We can see that the convergence rates of @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ are around 1.
This, of course, matches our theoretical expectations.
 
 
<a name="TestresultsoniWGQsub1subQsub1subRTsub1subi"></a><h3>Test results on <i>WG(Q<sub>1</sub>,Q<sub>1</sub>;RT<sub>[1]</sub>)</i></h3>
 
 
We can repeat the experiment from above using the next higher polynomial
degree:
The following figures are interior pressures and face pressures implemented using
@f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$. The mesh is refined 4 times.  Compared to the
previous figures using
@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, on each cell, the solution is no longer constant
on each cell, as we now use bilinear polynomials to do the approximation.
Consequently, there are 4 pressure values in one interior, 2 pressure values on
each face.
 
<table align="center">
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg111_2d_4.png" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg111_3d_4.png" alt=""></td>
  </tr>
</table>
 
Compared to the corresponding image for the @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$
combination, the solution is now substantially more accurate and, in
particular so close to being continuous at the interfaces that we can
no longer distinguish the interface pressures @f$p^\partial@f$ from the
interior pressures @f$p^\circ@f$ on the adjacent cells.
 
<a name="Convergencetableforik1i"></a><h4>Convergence table for <i>k=1</i></h4>
 
 
The following are the convergence rates of pressure, velocity, and flux
we obtain from using the @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ element combination:
 
<table align="center" class="doxtable">
  <tr>
   <th>number of refinements </th><th>  @f$\|p-p_h^\circ\|@f$  </th><th>  @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
  </tr>
  <tr>
    <td>  2           </td><td>           1.613e-02      </td><td>          5.093e-02     </td><td>             7.167e-02   </td>
  </tr>
  <tr>
    <td>  3           </td><td>           4.056e-03      </td><td>          1.276e-02     </td><td>             1.802e-02    </td>
  </tr>
  <tr>
    <td>  4           </td><td>           1.015e-03      </td><td>          3.191e-03     </td><td>             4.512e-03  </td>
  </tr>
  <tr>
    <td>  5           </td><td>           2.540e-04      </td><td>          7.979e-04     </td><td>             1.128e-03  </td>
  </tr>
  <tr>
    <th>Conv.rate     </th><th>              2.00        </th><th>             2.00       </th><th>                 2.00    </th>
  </tr>
</table>
 
The convergence rates of @f$WG(Q_1,Q_1;RT_{[1]})@f$ are around 2, as expected.
 
 
 
<a name="TestresultsoniWGQsub2subQsub2subRTsub2subi"></a><h3>Test results on <i>WG(Q<sub>2</sub>,Q<sub>2</sub>;RT<sub>[2]</sub>)</i></h3>
 
 
Let us go one polynomial degree higher.
The following are interior pressures and face pressures implemented using
@f$WG(Q_2,Q_2;RT_{[2]})@f$, with mesh size @f$h = 1/32@f$ (i.e., 5 global mesh
refinement steps). In the program, we use
`data_out_face.build_patches(fe.degree)` when generating graphical output
(see the documentation of DataOut::build_patches()), which here implies that
we divide each 2d cell interior into 4 subcells in order to provide a better
visualization of the quadratic polynomials.
<table align="center">
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg222_2d_5.png" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg222_3d_5.png" alt=""></td>
  </tr>
</table>
 
 
<a name="Convergencetableforik2i"></a><h4>Convergence table for <i>k=2</i></h4>
 
 
As before, we can generate convergence data for the
@f$L_2@f$ errors of pressure, velocity, and flux
using the @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ combination:
 
<table align="center" class="doxtable">
  <tr>
   <th>number of refinements </th><th>  @f$\|p-p_h^\circ\|@f$  </th><th>  @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
  </tr>
  <tr>
     <td>  2               </td><td>       1.072e-03       </td><td>         3.375e-03       </td><td>           4.762e-03   </td>
  </tr>
  <tr>
    <td>   3               </td><td>       1.347e-04       </td><td>         4.233e-04       </td><td>           5.982e-04    </td>
  </tr>
  <tr>
    <td>   4               </td><td>       1.685e-05      </td><td>          5.295e-05       </td><td>           7.487e-05  </td>
  </tr>
  <tr>
    <td>   5               </td><td>       2.107e-06      </td><td>          6.620e-06       </td><td>           9.362e-06  </td>
  </tr>
  <tr>
    <th>Conv.rate          </th><th>         3.00         </th><th>            3.00          </th><th>              3.00    </th>
  </tr>
</table>
 
Once more, the convergence rates of @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ is
as expected, with values around 3.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-61.cc"
*/