step-20.h

Go to the documentation of this file.

,
 *                                      const unsigned int /*component*/) const
 *     {
 *       return 0;
 *     }
 * 
 * 
 * 
 *     template <int dim>
 *     double
 *     PressureBoundaryValues<dim>::value(const Point<dim> &p,
 *                                        const unsigned int /*component*/) const
 *     {
 *       return -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
 *                alpha * p[0] * p[0] * p[0] / 6);
 *     }
 * 
 * 
 * 
 *     template <int dim>
 *     void ExactSolution<dim>::vector_value(const Point<dim> &p,
 *                                           Vector<double> &  values) const
 *     {
 *       AssertDimension(values.size(), dim + 1);
 * 
 *       values(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2;
 *       values(1) = alpha * p[0] * p[1];
 *       values(2) = -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
 *                     alpha * p[0] * p[0] * p[0] / 6);
 *     }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Theinversepermeabilitytensor"></a> 
 * <h3>The inverse permeability tensor</h3>
 * 
 
 * 
 * In addition to the other equation data, we also want to use a
 * permeability tensor, or better -- because this is all that appears in the
 * weak form -- the inverse of the permeability tensor,
 * <code>KInverse</code>. For the purpose of verifying the exactness of the
 * solution and determining convergence orders, this tensor is more in the
 * way than helpful. We will therefore simply set it to the identity matrix.
 *     
 
 * 
 * However, a spatially varying permeability tensor is indispensable in
 * real-life porous media flow simulations, and we would like to use the
 * opportunity to demonstrate the technique to use tensor valued functions.
 *     
 
 * 
 * Possibly unsurprisingly, deal.II also has a base class not only for
 * scalar and generally vector-valued functions (the <code>Function</code>
 * base class) but also for functions that return tensors of fixed dimension
 * and rank, the <code>TensorFunction</code> template. Here, the function
 * under consideration returns a dim-by-dim matrix, i.e. a tensor of rank 2
 * and dimension <code>dim</code>. We then choose the template arguments of
 * the base class appropriately.
 *     
 
 * 
 * The interface that the <code>TensorFunction</code> class provides is
 * essentially equivalent to the <code>Function</code> class. In particular,
 * there exists a <code>value_list</code> function that takes a list of
 * points at which to evaluate the function, and returns the values of the
 * function in the second argument, a list of tensors:
 * 
 * @code
 *     template <int dim>
 *     class KInverse : public TensorFunction<2, dim>
 *     {
 *     public:
 *       KInverse()
 *         : TensorFunction<2, dim>()
 *       {}
 * 
 *       virtual void
 *       value_list(const std::vector<Point<dim>> &points,
 *                  std::vector<Tensor<2, dim>> &  values) const override;
 *     };
 * 
 * 
 * @endcode
 * 
 * The implementation is less interesting. As in previous examples, we add a
 * check to the beginning of the class to make sure that the sizes of input
 * and output parameters are the same (see @ref step_5 "step-5" for a discussion of this
 * technique). Then we loop over all evaluation points, and for each one
 * set the output tensor to the identity matrix.
 *     
 
 * 
 * There is an oddity at the top of the function (the
 * `(void)points;` statement) that is worth discussing. The values
 * we put into the output `values` array does not actually depend
 * on the `points` arrays of coordinates at which the function is
 * evaluated. In other words, the `points` argument is in fact
 * unused, and we could have just not given it a name if we had
 * wanted. But we want to use the `points` object for checking
 * that the `values` object has the correct size. The problem is
 * that in release mode, `AssertDimension` is defined as a macro
 * that expands to nothing; the compiler will then complain that
 * the `points` object is unused. The idiomatic approach to
 * silencing this warning is to have a statement that evaluates
 * (reads) variable but doesn't actually do anything: That's what
 * `(void)points;` does: It reads from `points`, and then casts
 * the result of the read to `void`, i.e., nothing. This statement
 * is, in other words, completely pointless and implies no actual
 * action except to explain to the compiler that yes, this
 * variable is in fact used even in release mode. (In debug mode,
 * the `AssertDimension` macro expands to something that reads
 * from the variable, and so the funny statement would not be
 * necessary in debug mode.)
 * 
 * @code
 *     template <int dim>
 *     void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
 *                                    std::vector<Tensor<2, dim>> &  values) const
 *     {
 *       (void)points;
 *       AssertDimension(points.size(), values.size());
 * 
 *       for (auto &value : values)
 *         value = unit_symmetric_tensor<dim>();
 *     }
 *   } // namespace PrescribedSolution
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="MixedLaplaceProblemclassimplementation"></a> 
 * <h3>MixedLaplaceProblem class implementation</h3>
 * 
 
 * 
 * 
 * <a name="MixedLaplaceProblemMixedLaplaceProblem"></a> 
 * <h4>MixedLaplaceProblem::MixedLaplaceProblem</h4>
 * 
 
 * 
 * In the constructor of this class, we first store the value that was
 * passed in concerning the degree of the finite elements we shall use (a
 * degree of zero, for example, means to use RT(0) and DG(0)), and then
 * construct the vector valued element belonging to the space @f$X_h@f$ described
 * in the introduction. The rest of the constructor is as in the early
 * tutorial programs.
 *   
 
 * 
 * The only thing worth describing here is the constructor call of the
 * <code>fe</code> variable. The <code>FESystem</code> class to which this
 * variable belongs has a number of different constructors that all refer to
 * binding simpler elements together into one larger element. In the present
 * case, we want to couple a single RT(degree) element with a single
 * DQ(degree) element. The constructor to <code>FESystem</code> that does
 * this requires us to specify first the first base element (the
 * <code>FE_RaviartThomas</code> object of given degree) and then the number
 * of copies for this base element, and then similarly the kind and number
 * of <code>FE_DGQ</code> elements. Note that the Raviart-Thomas element
 * already has <code>dim</code> vector components, so that the coupled
 * element will have <code>dim+1</code> vector components, the first
 * <code>dim</code> of which correspond to the velocity variable whereas the
 * last one corresponds to the pressure.
 *   
 
 * 
 * It is also worth comparing the way we constructed this element from its
 * base elements, with the way we have done so in @ref step_8 "step-8": there, we have
 * built it as <code>fe (FE_Q@<dim@>(1), dim)</code>, i.e. we have simply
 * used <code>dim</code> copies of the <code>FE_Q(1)</code> element, one
 * copy for the displacement in each coordinate direction.
 * 
 * @code
 *   template <int dim>
 *   MixedLaplaceProblem<dim>::MixedLaplaceProblem(const unsigned int degree)
 *     : degree(degree)
 *     , fe(FE_RaviartThomas<dim>(degree), 1, FE_DGQ<dim>(degree), 1)
 *     , dof_handler(triangulation)
 *   {}
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="MixedLaplaceProblemmake_grid_and_dofs"></a> 
 * <h4>MixedLaplaceProblem::make_grid_and_dofs</h4>
 * 
 
 * 
 * This next function starts out with well-known functions calls that create
 * and refine a mesh, and then associate degrees of freedom with it:
 * 
 * @code
 *   template <int dim>
 *   void MixedLaplaceProblem<dim>::make_grid_and_dofs()
 *   {
 *     GridGenerator::hyper_cube(triangulation, -1, 1);
 *     triangulation.refine_global(5);
 * 
 *     dof_handler.distribute_dofs(fe);
 * 
 * @endcode
 * 
 * However, then things become different. As mentioned in the
 * introduction, we want to subdivide the matrix into blocks corresponding
 * to the two different kinds of variables, velocity and pressure. To this
 * end, we first have to make sure that the indices corresponding to
 * velocities and pressures are not intermingled: First all velocity
 * degrees of freedom, then all pressure DoFs. This way, the global matrix
 * separates nicely into a @f$2 \times 2@f$ system. To achieve this, we have to
 * renumber degrees of freedom based on their vector component, an
 * operation that conveniently is already implemented:
 * 
 * @code
 *     DoFRenumbering::component_wise(dof_handler);
 * 
 * @endcode
 * 
 * The next thing is that we want to figure out the sizes of these blocks
 * so that we can allocate an appropriate amount of space. To this end, we
 * call the DoFTools::count_dofs_per_fe_component() function that
 * counts how many shape functions are non-zero for a particular vector
 * component. We have <code>dim+1</code> vector components, and
 * DoFTools::count_dofs_per_fe_component() will count how many shape
 * functions belong to each of these components.
 *     
 
 * 
 * There is one problem here. As described in the documentation of that
 * function, it <i>wants</i> to put the number of @f$x@f$-velocity shape
 * functions into <code>dofs_per_component[0]</code>, the number of
 * @f$y@f$-velocity shape functions into <code>dofs_per_component[1]</code>
 * (and similar in 3d), and the number of pressure shape functions into
 * <code>dofs_per_component[dim]</code>. But, the Raviart-Thomas element
 * is special in that it is non-@ref GlossPrimitive "primitive", i.e.,
 * for Raviart-Thomas elements all velocity shape functions
 * are nonzero in all components. In other words, the function cannot
 * distinguish between @f$x@f$ and @f$y@f$ velocity functions because there
 * <i>is</i> no such distinction. It therefore puts the overall number
 * of velocity into each of <code>dofs_per_component[c]</code>,
 * @f$0\le c\le \text{dim}@f$. On the other hand, the number
 * of pressure variables equals the number of shape functions that are
 * nonzero in the dim-th component.
 *     
 
 * 
 * Using this knowledge, we can get the number of velocity shape
 * functions from any of the first <code>dim</code> elements of
 * <code>dofs_per_component</code>, and then use this below to initialize
 * the vector and matrix block sizes, as well as create output.
 *     
 
 * 
 * @note If you find this concept difficult to understand, you may
 * want to consider using the function DoFTools::count_dofs_per_fe_block()
 * instead, as we do in the corresponding piece of code in @ref step_22 "step-22".
 * You might also want to read up on the difference between
 * @ref GlossBlock "blocks" and @ref GlossComponent "components"
 * in the glossary.
 * 
 * @code
 *     const std::vector<types::global_dof_index> dofs_per_component =
 *       DoFTools::count_dofs_per_fe_component(dof_handler);
 *     const unsigned int n_u = dofs_per_component[0],
 *                        n_p = dofs_per_component[dim];
 * 
 *     std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *               << std::endl
 *               << "Total number of cells: " << triangulation.n_cells()
 *               << std::endl
 *               << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *               << " (" << n_u << '+' << n_p << ')' << std::endl;
 * 
 * @endcode
 * 
 * The next task is to allocate a sparsity pattern for the matrix that we
 * will create. We use a compressed sparsity pattern like in the previous
 * steps, but as <code>system_matrix</code> is a block matrix we use the
 * class <code>BlockDynamicSparsityPattern</code> instead of just
 * <code>DynamicSparsityPattern</code>. This block sparsity pattern has
 * four blocks in a @f$2 \times 2@f$ pattern. The blocks' sizes depend on
 * <code>n_u</code> and <code>n_p</code>, which hold the number of velocity
 * and pressure variables.
 * 
 * @code
 *     const std::vector<types::global_dof_index> block_sizes = {n_u, n_p};
 *     BlockDynamicSparsityPattern                dsp(block_sizes, block_sizes);
 *     DoFTools::make_sparsity_pattern(dof_handler, dsp);
 * 
 * @endcode
 * 
 * We use the compressed block sparsity pattern in the same way as the
 * non-block version to create the sparsity pattern and then the system
 * matrix:
 * 
 * @code
 *     sparsity_pattern.copy_from(dsp);
 *     system_matrix.reinit(sparsity_pattern);
 * 
 * @endcode
 * 
 * Then we have to resize the solution and right hand side vectors in
 * exactly the same way as the block compressed sparsity pattern:
 * 
 * @code
 *     solution.reinit(block_sizes);
 *     system_rhs.reinit(block_sizes);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="MixedLaplaceProblemassemble_system"></a> 
 * <h4>MixedLaplaceProblem::assemble_system</h4>
 * 
 
 * 
 * Similarly, the function that assembles the linear system has mostly been
 * discussed already in the introduction to this example. At its top, what
 * happens are all the usual steps, with the addition that we do not only
 * allocate quadrature and <code>FEValues</code> objects for the cell terms,
 * but also for face terms. After that, we define the usual abbreviations
 * for variables, and the allocate space for the local matrix and right hand
 * side contributions, and the array that holds the global numbers of the
 * degrees of freedom local to the present cell.
 * 
 * @code
 *   template <int dim>
 *   void MixedLaplaceProblem<dim>::assemble_system()
 *   {
 *     QGauss<dim>     quadrature_formula(degree + 2);
 *     QGauss<dim - 1> face_quadrature_formula(degree + 2);
 * 
 *     FEValues<dim>     fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               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);
 * 
 *     const unsigned int dofs_per_cell   = fe.n_dofs_per_cell();
 *     const unsigned int n_q_points      = quadrature_formula.size();
 *     const unsigned int n_face_q_points = face_quadrature_formula.size();
 * 
 *     FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
 *     Vector<double>     local_rhs(dofs_per_cell);
 * 
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 * 
 * @endcode
 * 
 * The next step is to declare objects that represent the source term,
 * pressure boundary value, and coefficient in the equation. In addition
 * to these objects that represent continuous functions, we also need
 * arrays to hold their values at the quadrature points of individual
 * cells (or faces, for the boundary values). Note that in the case of the
 * coefficient, the array has to be one of matrices.
 * 
 * @code
 *     const PrescribedSolution::RightHandSide<dim> right_hand_side;
 *     const PrescribedSolution::PressureBoundaryValues<dim>
 *                                             pressure_boundary_values;
 *     const PrescribedSolution::KInverse<dim> k_inverse;
 * 
 *     std::vector<double>         rhs_values(n_q_points);
 *     std::vector<double>         boundary_values(n_face_q_points);
 *     std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
 * 
 * @endcode
 * 
 * Finally, we need a couple of extractors that we will use to get at the
 * velocity and pressure components of vector-valued shape
 * functions. Their function and use is described in detail in the @ref
 * vector_valued report. Essentially, we will use them as subscripts on
 * the FEValues objects below: the FEValues object describes all vector
 * components of shape functions, while after subscription, it will only
 * refer to the velocities (a set of <code>dim</code> components starting
 * at component zero) or the pressure (a scalar component located at
 * position <code>dim</code>):
 * 
 * @code
 *     const FEValuesExtractors::Vector velocities(0);
 *     const FEValuesExtractors::Scalar pressure(dim);
 * 
 * @endcode
 * 
 * With all this in place, we can go on with the loop over all cells. The
 * body of this loop has been discussed in the introduction, and will not
 * be commented any further here:
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         fe_values.reinit(cell);
 *         local_matrix = 0;
 *         local_rhs    = 0;
 * 
 *         right_hand_side.value_list(fe_values.get_quadrature_points(),
 *                                    rhs_values);
 *         k_inverse.value_list(fe_values.get_quadrature_points(),
 *                              k_inverse_values);
 * 
 *         for (unsigned int q = 0; q < n_q_points; ++q)
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
 *               const double div_phi_i_u = fe_values[velocities].divergence(i, q);
 *               const double phi_i_p     = fe_values[pressure].value(i, q);
 * 
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 {
 *                   const Tensor<1, dim> phi_j_u =
 *                     fe_values[velocities].value(j, q);
 *                   const double div_phi_j_u =
 *                     fe_values[velocities].divergence(j, q);
 *                   const double phi_j_p = fe_values[pressure].value(j, q);
 * 
 *                   local_matrix(i, j) +=
 *                     (phi_i_u * k_inverse_values[q] * phi_j_u 
 *                      - phi_i_p * div_phi_j_u                 
 *                      - div_phi_i_u * phi_j_p)                
 *                     * fe_values.JxW(q);
 *                 }
 * 
 *               local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
 *             }
 * 
 *         for (const auto &face : cell->face_iterators())
 *           if (face->at_boundary())
 *             {
 *               fe_face_values.reinit(cell, face);
 * 
 *               pressure_boundary_values.value_list(
 *                 fe_face_values.get_quadrature_points(), boundary_values);
 * 
 *               for (unsigned int q = 0; q < n_face_q_points; ++q)
 *                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                   local_rhs(i) += -(fe_face_values[velocities].value(i, q) * 
 *                                     fe_face_values.normal_vector(q) *        
 *                                     boundary_values[q] *                     
 *                                     fe_face_values.JxW(q));
 *             }
 * 
 * @endcode
 * 
 * The final step in the loop over all cells is to transfer local
 * contributions into the global matrix and right hand side
 * vector. Note that we use exactly the same interface as in previous
 * examples, although we now use block matrices and vectors instead of
 * the regular ones. In other words, to the outside world, block
 * objects have the same interface as matrices and vectors, but they
 * additionally allow to access individual blocks.
 * 
 * @code
 *         cell->get_dof_indices(local_dof_indices);
 *         for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *           for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *             system_matrix.add(local_dof_indices[i],
 *                               local_dof_indices[j],
 *                               local_matrix(i, j));
 *         for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *           system_rhs(local_dof_indices[i]) += local_rhs(i);
 *       }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Implementationoflinearsolversandpreconditioners"></a> 
 * <h3>Implementation of linear solvers and preconditioners</h3>
 * 
 
 * 
 * The linear solvers and preconditioners we use in this example have
 * been discussed in significant detail already in the introduction. We
 * will therefore not discuss the rationale for our approach here any
 * more, but rather only comment on some remaining implementational
 * aspects.
 * 
 
 * 
 * 
 * <a name="MixedLaplacesolve"></a> 
 * <h4>MixedLaplace::solve</h4>
 * 
 
 * 
 * As already outlined in the introduction, the solve function consists
 * essentially of two steps. First, we have to form the first equation
 * involving the Schur complement and solve for the pressure (component 1
 * of the solution). Then, we can reconstruct the velocities from the
 * second equation (component 0 of the solution).
 * 
 * @code
 *   template <int dim>
 *   void MixedLaplaceProblem<dim>::solve()
 *   {
 * @endcode
 * 
 * As a first step we declare references to all block components of the
 * matrix, the right hand side and the solution vector that we will
 * need.
 * 
 * @code
 *     const auto &M = system_matrix.block(0, 0);
 *     const auto &B = system_matrix.block(0, 1);
 * 
 *     const auto &F = system_rhs.block(0);
 *     const auto &G = system_rhs.block(1);
 * 
 *     auto &U = solution.block(0);
 *     auto &P = solution.block(1);
 * 
 * @endcode
 * 
 * Then, we will create corresponding LinearOperator objects and create
 * the <code>op_M_inv</code> operator:
 * 
 * @code
 *     const auto op_M = linear_operator(M);
 *     const auto op_B = linear_operator(B);
 * 
 *     ReductionControl         reduction_control_M(2000, 1.0e-18, 1.0e-10);
 *     SolverCG<Vector<double>> solver_M(reduction_control_M);
 *     PreconditionJacobi<SparseMatrix<double>> preconditioner_M;
 * 
 *     preconditioner_M.initialize(M);
 * 
 *     const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
 * 
 * @endcode
 * 
 * This allows us to declare the Schur complement <code>op_S</code> and
 * the approximate Schur complement <code>op_aS</code>:
 * 
 * @code
 *     const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
 *     const auto op_aS =
 *       transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;
 * 
 * @endcode
 * 
 * We now create a preconditioner out of <code>op_aS</code> that
 * applies a fixed number of 30 (inexpensive) CG iterations:
 * 
 * @code
 *     IterationNumberControl   iteration_number_control_aS(30, 1.e-18);
 *     SolverCG<Vector<double>> solver_aS(iteration_number_control_aS);
 * 
 *     const auto preconditioner_S =
 *       inverse_operator(op_aS, solver_aS, PreconditionIdentity());
 * 
 * @endcode
 * 
 * Now on to the first equation. The right hand side of it is
 * @f$B^TM^{-1}F-G@f$, which is what we compute in the first few lines. We
 * then solve the first equation with a CG solver and the
 * preconditioner we just declared.
 * 
 * @code
 *     const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
 * 
 *     SolverControl            solver_control_S(2000, 1.e-12);
 *     SolverCG<Vector<double>> solver_S(solver_control_S);
 * 
 *     const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);
 * 
 *     P = op_S_inv * schur_rhs;
 * 
 *     std::cout << solver_control_S.last_step()
 *               << " CG Schur complement iterations to obtain convergence."
 *               << std::endl;
 * 
 * @endcode
 * 
 * After we have the pressure, we can compute the velocity. The equation
 * reads @f$MU=-BP+F@f$, and we solve it by first computing the right hand
 * side, and then multiplying it with the object that represents the
 * inverse of the mass matrix:
 * 
 * @code
 *     U = op_M_inv * (F - op_B * P);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="MixedLaplaceProblemclassimplementationcontinued"></a> 
 * <h3>MixedLaplaceProblem class implementation (continued)</h3>
 * 
 
 * 
 * 
 * <a name="MixedLaplacecompute_errors"></a> 
 * <h4>MixedLaplace::compute_errors</h4>
 * 
 
 * 
 * After we have dealt with the linear solver and preconditioners, we
 * continue with the implementation of our main class. In particular, the
 * next task is to compute the errors in our numerical solution, in both the
 * pressures as well as velocities.
 *   
 
 * 
 * To compute errors in the solution, we have already introduced the
 * <code>VectorTools::integrate_difference</code> function in @ref step_7 "step-7" and
 * @ref step_11 "step-11". However, there we only dealt with scalar solutions, whereas here
 * we have a vector-valued solution with components that even denote
 * different quantities and may have different orders of convergence (this
 * isn't the case here, by choice of the used finite elements, but is
 * frequently the case in mixed finite element applications). What we
 * therefore have to do is to `mask' the components that we are interested
 * in. This is easily done: the
 * <code>VectorTools::integrate_difference</code> function takes as one of its
 * arguments a pointer to a weight function (the parameter defaults to the
 * null pointer, meaning unit weights). What we have to do is to pass
 * a function object that equals one in the components we are interested in,
 * and zero in the other ones. For example, to compute the pressure error,
 * we should pass a function that represents the constant vector with a unit
 * value in component <code>dim</code>, whereas for the velocity the
 * constant vector should be one in the first <code>dim</code> components,
 * and zero in the location of the pressure.
 *   
 
 * 
 * In deal.II, the <code>ComponentSelectFunction</code> does exactly this:
 * it wants to know how many vector components the function it is to
 * represent should have (in our case this would be <code>dim+1</code>, for
 * the joint velocity-pressure space) and which individual or range of
 * components should be equal to one. We therefore define two such masks at
 * the beginning of the function, following by an object representing the
 * exact solution and a vector in which we will store the cellwise errors as
 * computed by <code>integrate_difference</code>:
 * 
 * @code
 *   template <int dim>
 *   void MixedLaplaceProblem<dim>::compute_errors() const
 *   {
 *     const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
 *     const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
 *                                                      dim + 1);
 * 
 *     PrescribedSolution::ExactSolution<dim> exact_solution;
 *     Vector<double> cellwise_errors(triangulation.n_active_cells());
 * 
 * @endcode
 * 
 * As already discussed in @ref step_7 "step-7", we have to realize that it is
 * impossible to integrate the errors exactly. All we can do is
 * approximate this integral using quadrature. This actually presents a
 * slight twist here: if we naively chose an object of type
 * <code>QGauss@<dim@>(degree+1)</code> as one may be inclined to do (this
 * is what we used for integrating the linear system), one realizes that
 * the error is very small and does not follow the expected convergence
 * curves at all. What is happening is that for the mixed finite elements
 * used here, the Gauss points happen to be superconvergence points in
 * which the pointwise error is much smaller (and converges with higher
 * order) than anywhere else. These are therefore not particularly good
 * points for integration. To avoid this problem, we simply use a
 * trapezoidal rule and iterate it <code>degree+2</code> times in each
 * coordinate direction (again as explained in @ref step_7 "step-7"):
 * 
 * @code
 *     QTrapezoid<1>  q_trapez;
 *     QIterated<dim> quadrature(q_trapez, degree + 2);
 * 
 * @endcode
 * 
 * With this, we can then let the library compute the errors and output
 * them to the screen:
 * 
 * @code
 *     VectorTools::integrate_difference(dof_handler,
 *                                       solution,
 *                                       exact_solution,
 *                                       cellwise_errors,
 *                                       quadrature,
 *                                       VectorTools::L2_norm,
 *                                       &pressure_mask);
 *     const double p_l2_error =
 *       VectorTools::compute_global_error(triangulation,
 *                                         cellwise_errors,
 *                                         VectorTools::L2_norm);
 * 
 *     VectorTools::integrate_difference(dof_handler,
 *                                       solution,
 *                                       exact_solution,
 *                                       cellwise_errors,
 *                                       quadrature,
 *                                       VectorTools::L2_norm,
 *                                       &velocity_mask);
 *     const double u_l2_error =
 *       VectorTools::compute_global_error(triangulation,
 *                                         cellwise_errors,
 *                                         VectorTools::L2_norm);
 * 
 *     std::cout << "Errors: ||e_p||_L2 = " << p_l2_error
 *               << ",   ||e_u||_L2 = " << u_l2_error << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="MixedLaplaceoutput_results"></a> 
 * <h4>MixedLaplace::output_results</h4>
 * 
 
 * 
 * The last interesting function is the one in which we generate graphical
 * output. Note that all velocity components get the same solution name
 * "u". Together with using
 * DataComponentInterpretation::component_is_part_of_vector this will
 * cause DataOut<dim>::write_vtu() to generate a vector representation of
 * the individual velocity components, see @ref step_22 "step-22" or the
 * @ref VVOutput "Generating graphical output"
 * section of the
 * @ref vector_valued
 * module for more information. Finally, it seems inappropriate for higher
 * order elements to only show a single bilinear quadrilateral per cell in
 * the graphical output. We therefore generate patches of size
 * (degree+1)x(degree+1) to capture the full information content of the
 * solution. See the @ref step_7 "step-7" tutorial program for more information on this.
 * 
 * @code
 *   template <int dim>
 *   void MixedLaplaceProblem<dim>::output_results() const
 *   {
 *     std::vector<std::string> solution_names(dim, "u");
 *     solution_names.emplace_back("p");
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *       interpretation(dim,
 *                      DataComponentInterpretation::component_is_part_of_vector);
 *     interpretation.push_back(DataComponentInterpretation::component_is_scalar);
 * 
 *     DataOut<dim> data_out;
 *     data_out.add_data_vector(dof_handler,
 *                              solution,
 *                              solution_names,
 *                              interpretation);
 * 
 *     data_out.build_patches(degree + 1);
 * 
 *     std::ofstream output("solution.vtu");
 *     data_out.write_vtu(output);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="MixedLaplacerun"></a> 
 * <h4>MixedLaplace::run</h4>
 * 
 
 * 
 * This is the final function of our main class. It's only job is to call
 * the other functions in their natural order:
 * 
 * @code
 *   template <int dim>
 *   void MixedLaplaceProblem<dim>::run()
 *   {
 *     make_grid_and_dofs();
 *     assemble_system();
 *     solve();
 *     compute_errors();
 *     output_results();
 *   }
 * } // namespace Step20
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * The main function we stole from @ref step_6 "step-6" instead of @ref step_4 "step-4". It is almost
 * equal to the one in @ref step_6 "step-6" (apart from the changed class names, of course),
 * the only exception is that we pass the degree of the finite element space
 * to the constructor of the mixed Laplace problem (here, we use zero-th order
 * elements).
 * 
 * @code
 * int main()
 * {
 *   try
 *     {
 *       using namespace Step20;
 * 
 *       const unsigned int     fe_degree = 0;
 *       MixedLaplaceProblem<2> mixed_laplace_problem(fe_degree);
 *       mixed_laplace_problem.run();
 *     }
 *   catch (std::exception &exc)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Exception on processing: " << std::endl
 *                 << exc.what() << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 * 
 *       return 1;
 *     }
 *   catch (...)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Unknown exception!" << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 * 
 *   return 0;
 * }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
<a name="Outputoftheprogramandgraphicalvisualization"></a><h3>Output of the program and graphical visualization</h3>
 
 
 
If we run the program as is, we get this output for the @f$32\times 32@f$
mesh we use (for a total of 1024 cells with 1024 pressure degrees of
freedom since we use piecewise constants, and 2112 velocities because
the Raviart-Thomas element defines one degree per freedom per face and
there are @f$1024 + 32 = 1056@f$ faces parallel to the @f$x@f$-axis and the same
number parallel to the @f$y@f$-axis):
@verbatim
\$ make run
[ 66%] Built target step-20
Scanning dependencies of target run
[100%] Run step-20 with Release configuration
Number of active cells: 1024
Total number of cells: 1365
Number of degrees of freedom: 3136 (2112+1024)
24 CG Schur complement iterations to obtain convergence.
Errors: ||e_p||_L2 = 0.0445032,   ||e_u||_L2 = 0.010826
[100%] Built target run
@endverbatim
 
The fact that the number of iterations is so small, of course, is due to
the good (but expensive!) preconditioner we have developed. To get
confidence in the solution, let us take a look at it. The following three
images show (from left to right) the x-velocity, the y-velocity, and the
pressure:
 
<table style="width:60%" align="center">
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.u_new.jpg" width="400" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.v_new.jpg" width="400" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.p_new.jpg" width="400" alt=""></td>
  </tr>
</table>
 
 
 
Let us start with the pressure: it is highest at the left and lowest at the
right, so flow will be from left to right. In addition, though hardly visible
in the graph, we have chosen the pressure field such that the flow left-right
flow first channels towards the center and then outward again. Consequently,
the x-velocity has to increase to get the flow through the narrow part,
something that can easily be seen in the left image. The middle image
represents inward flow in y-direction at the left end of the domain, and
outward flow in y-direction at the right end of the domain.
 
 
 
As an additional remark, note how the x-velocity in the left image is only
continuous in x-direction, whereas the y-velocity is continuous in
y-direction. The flow fields are discontinuous in the other directions. This
very obviously reflects the continuity properties of the Raviart-Thomas
elements, which are, in fact, only in the space H(div) and not in the space
@f$H^1@f$. Finally, the pressure field is completely discontinuous, but
that should not surprise given that we have chosen <code>FE_DGQ(0)</code> as
the finite element for that solution component.
 
 
 
<a name="Convergence"></a><h3>Convergence</h3>
 
 
 
The program offers two obvious places where playing and observing convergence
is in order: the degree of the finite elements used (passed to the constructor
of the <code>MixedLaplaceProblem</code> class from <code>main()</code>), and
the refinement level (determined in
<code>MixedLaplaceProblem::make_grid_and_dofs</code>). What one can do is to
change these values and observe the errors computed later on in the course of
the program run.
 
 
 
If one does this, one finds the following pattern for the @f$L_2@f$ error
in the pressure variable:
<table align="center" class="doxtable">
  <tr>
    <th></th>
    <th colspan="3" align="center">Finite element order</th>
  </tr>
  <tr>
    <th>Refinement level</th>
    <th>0</th>
    <th>1</th>
    <th>2</th>
  </tr>
  <tr>
    <th>0</th>  <td>1.45344</td>  <td>0.0831743</td>  <td>0.0235186</td>
  </tr>
  <tr>
    <th>1</th>  <td>0.715099</td>  <td>0.0245341</td>  <td>0.00293983</td>
  </tr>
  <tr>
    <th>2</th>  <td>0.356383</td>  <td>0.0063458</td>  <td>0.000367478</td>
  </tr>
  <tr>
    <th>3</th>  <td>0.178055</td>  <td>0.00159944</td>  <td>4.59349e-05</td>
  </tr>
  <tr>
    <th>4</th>  <td>0.0890105</td>  <td>0.000400669</td>  <td>5.74184e-06</td>
  </tr>
  <tr>
    <th>5</th>  <td>0.0445032</td>  <td>0.000100218</td>  <td>7.17799e-07</td>
  </tr>
  <tr>
    <th>6</th>  <td>0.0222513</td>  <td>2.50576e-05</td>  <td>9.0164e-08</td>
  </tr>
  <tr>
    <th></th>  <th>@f$O(h)@f$</th>  <th>@f$O(h^2)@f$</th>  <th>@f$O(h^3)@f$</th>
  </tr>
</table>
 
The theoretically expected convergence orders are very nicely reflected by the
experimentally observed ones indicated in the last row of the table.
 
 
 
One can make the same experiment with the @f$L_2@f$ error
in the velocity variables:
<table align="center" class="doxtable">
  <tr>
    <th></th>
    <th colspan="3" align="center">Finite element order</th>
  </tr>
  <tr>
    <th>Refinement level</th>
    <th>0</th>
    <th>1</th>
    <th>2</th>
  </tr>
  <tr>
    <th>0</th> <td>0.367423</td> <td>0.127657</td> <td>5.10388e-14</td>
  </tr>
  <tr>
    <th>1</th> <td>0.175891</td> <td>0.0319142</td> <td>9.04414e-15</td>
  </tr>
  <tr>
    <th>2</th> <td>0.0869402</td> <td>0.00797856</td> <td>1.23723e-14</td>
  </tr>
  <tr>
    <th>3</th> <td>0.0433435</td> <td>0.00199464</td> <td>1.86345e-07</td>
  </tr>
  <tr>
    <th>4</th> <td>0.0216559</td> <td>0.00049866</td> <td>2.72566e-07</td>
  </tr>
  <tr>
    <th>5</th> <td>0.010826</td> <td>0.000124664</td> <td>3.57141e-07</td>
  </tr>
  <tr>
    <th>6</th> <td>0.00541274</td> <td>3.1166e-05</td> <td>4.46124e-07</td>
  </tr>
  <tr>
    <th></th>  <td>@f$O(h)@f$</td>  <td>@f$O(h^2)@f$</td>  <td>@f$O(h^3)@f$</td>
  </tr>
</table>
The result concerning the convergence order is the same here.
 
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
<a name="Morerealisticpermeabilityfields"></a><h4>More realistic permeability fields</h4>
 
 
Realistic flow computations for ground water or oil reservoir simulations will
not use a constant permeability. Here's a first, rather simple way to change
this situation: we use a permeability that decays very rapidly away from a
central flowline until it hits a background value of 0.001. This is to mimic
the behavior of fluids in sandstone: in most of the domain, the sandstone is
homogeneous and, while permeable to fluids, not overly so; on the other stone,
the stone has cracked, or faulted, along one line, and the fluids flow much
easier along this large crack. Here is how we could implement something like
this:
@code
template <int dim>
void
KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
                           std::vector<Tensor<2,dim> >    &values) const
{
  AssertDimension (points.size(), values.size());
 
  for (unsigned int p=0; p<points.size(); ++p)
    {
      values[p].clear ();
 
      const double distance_to_flowline
        = std::fabs(points[p][1]-0.2*std::sin(10*points[p][0]));
 
      const double permeability = std::max(std::exp(-(distance_to_flowline*
                                                      distance_to_flowline)
                                                    / (0.1 * 0.1)),
                                           0.001);
 
      for (unsigned int d=0; d<dim; ++d)
        values[p][d][d] = 1./permeability;
    }
}
@endcode
Remember that the function returns the inverse of the permeability tensor.
 
 
 
With a significantly higher mesh resolution, we can visualize this, here with
x- and y-velocity:
 
<table style="width:60%" align="center">
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-wiggle.png" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.v-wiggle.png" alt=""></td>
  </tr>
</table>
 
It is obvious how fluids flow essentially only along the middle line, and not
anywhere else.
 
 
 
Another possibility would be to use a random permeability field. A simple way
to achieve this would be to scatter a number of centers around the domain and
then use a permeability field that is the sum of (negative) exponentials for
each of these centers. Flow would then try to hop from one center of high
permeability to the next one. This is an entirely unscientific attempt at
describing a random medium, but one possibility to implement this behavior
would look like this:
@code
template <int dim>
class KInverse : public TensorFunction<2,dim>
{
  public:
    KInverse ();
 
    virtual void value_list (const std::vector<Point<dim> > &points,
                             std::vector<Tensor<2,dim> >    &values) const;
 
  private:
    std::vector<Point<dim> > centers;
};
 
 
template <int dim>
KInverse<dim>::KInverse ()
{
  const unsigned int N = 40;
  centers.resize (N);
  for (unsigned int i=0; i<N; ++i)
    for (unsigned int d=0; d<dim; ++d)
      centers[i][d] = 2.*rand()/RAND_MAX-1;
}
 
 
template <int dim>
void
KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
                           std::vector<Tensor<2,dim> >    &values) const
{
  AssertDimension (points.size(), values.size());
 
  for (unsigned int p=0; p<points.size(); ++p)
    {
      values[p].clear ();
 
      double permeability = 0;
      for (unsigned int i=0; i<centers.size(); ++i)
        permeability += std::exp(-(points[p] - centers[i]).norm_square() / (0.1 * 0.1));
 
      const double normalized_permeability
        = std::max(permeability, 0.005);
 
      for (unsigned int d=0; d<dim; ++d)
        values[p][d][d] = 1./normalized_permeability;
    }
}
@endcode
 
A piecewise constant interpolation of the diagonal elements of the
inverse of this tensor (i.e., of <code>normalized_permeability</code>)
looks as follows:
 
<img src="https://www.dealii.org/images/steps/developer/step-20.k-random.png" alt="">
 
 
With a permeability field like this, we would get x-velocities and pressures as
follows:
 
<table style="width:60%" align="center">
  <tr>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.u-random.png" alt=""></td>
    <td><img src="https://www.dealii.org/images/steps/developer/step-20.p-random.png" alt=""></td>
  </tr>
</table>
 
We will use these permeability fields again in @ref step_21 "step-21" and @ref step_43 "step-43".
 
 
<a name="Betterlinearsolvers"></a><h4>Better linear solvers</h4>
 
 
As mentioned in the introduction, the Schur complement solver used here is not
the best one conceivable (nor is it intended to be a particularly good
one). Better ones can be found in the literature and can be built using the
same block matrix techniques that were introduced here. We pick up on this
theme again in @ref step_22 "step-22", where we first build a Schur complement solver for the
Stokes equation as we did here, and then in the <a
href="step_22.html#improved-solver">Improved Solvers</a> section discuss better
ways based on solving the system as a whole but preconditioning based on
individual blocks. We will also come back to this in @ref step_43 "step-43".
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-20.cc"
*/