step-82.h

Go to the documentation of this file.

) const
 *     {
 *       double return_value = 0.0;
 *   
 *       if (dim == 2)
 *         {
 *           return_value = 24.0 * std::pow(p(1) * (1.0 - p(1)), 2) +
 *                          +24.0 * std::pow(p(0) * (1.0 - p(0)), 2) +
 *                          2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
 *                            (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1));
 *         }
 *       else if (dim == 3)
 *         {
 *           return_value =
 *             24.0 * std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2) +
 *             24.0 * std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2) +
 *             24.0 * std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2) +
 *             2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
 *               (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
 *               std::pow(p(2) * (1.0 - p(2)), 2) +
 *             2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
 *               (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
 *               std::pow(p(1) * (1.0 - p(1)), 2) +
 *             2.0 * (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
 *               (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
 *               std::pow(p(0) * (1.0 - p(0)), 2);
 *         }
 *       else
 *         Assert(false, ExcNotImplemented());
 *   
 *       return return_value;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * This class implement the manufactured (exact) solution @f$u@f$. To compute the
 * errors, we need the value of @f$u@f$ as well as its gradient and its Hessian.
 * 
 * @code
 *     template <int dim>
 *     class ExactSolution : public Function<dim>
 *     {
 *     public:
 *       ExactSolution()
 *         : Function<dim>()
 *       {}
 *   
 *       virtual double value(const Point<dim> & p,
 *                            const unsigned int component = 0) const override;
 *   
 *       virtual Tensor<1, dim>
 *       gradient(const Point<dim> & p,
 *                const unsigned int component = 0) const override;
 *   
 *       virtual SymmetricTensor<2, dim>
 *       hessian(const Point<dim> & p,
 *               const unsigned int component = 0) const override;
 *     };
 *   
 *   
 *   
 *     template <int dim>
 *     double ExactSolution<dim>::value(const Point<dim> &p,
 *                                      const unsigned int /*component*/) const
 *     {
 *       double return_value = 0.0;
 *   
 *       if (dim == 2)
 *         {
 *           return_value = std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2);
 *         }
 *       else if (dim == 3)
 *         {
 *           return_value = std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)) *
 *                                     p(2) * (1.0 - p(2)),
 *                                   2);
 *         }
 *       else
 *         Assert(false, ExcNotImplemented());
 *   
 *       return return_value;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     Tensor<1, dim>
 *     ExactSolution<dim>::gradient(const Point<dim> &p,
 *                                  const unsigned int /*component*/) const
 *     {
 *       Tensor<1, dim> return_gradient;
 *   
 *       if (dim == 2)
 *         {
 *           return_gradient[0] =
 *             (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
 *             std::pow(p(1) * (1.0 - p(1)), 2);
 *           return_gradient[1] =
 *             (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
 *             std::pow(p(0) * (1.0 - p(0)), 2);
 *         }
 *       else if (dim == 3)
 *         {
 *           return_gradient[0] =
 *             (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
 *             std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2);
 *           return_gradient[1] =
 *             (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
 *             std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2);
 *           return_gradient[2] =
 *             (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) *
 *             std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2);
 *         }
 *       else
 *         Assert(false, ExcNotImplemented());
 *   
 *       return return_gradient;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     SymmetricTensor<2, dim>
 *     ExactSolution<dim>::hessian(const Point<dim> &p,
 *                                 const unsigned int /*component*/) const
 *     {
 *       SymmetricTensor<2, dim> return_hessian;
 *   
 *       if (dim == 2)
 *         {
 *           return_hessian[0][0] = (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
 *                                  std::pow(p(1) * (1.0 - p(1)), 2);
 *           return_hessian[0][1] =
 *             (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
 *             (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3));
 *           return_hessian[1][1] = (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
 *                                  std::pow(p(0) * (1.0 - p(0)), 2);
 *         }
 *       else if (dim == 3)
 *         {
 *           return_hessian[0][0] =
 *             (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
 *             std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2);
 *           return_hessian[0][1] =
 *             (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
 *             (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
 *             std::pow(p(2) * (1.0 - p(2)), 2);
 *           return_hessian[0][2] =
 *             (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
 *             (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) *
 *             std::pow(p(1) * (1.0 - p(1)), 2);
 *           return_hessian[1][1] =
 *             (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
 *             std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2);
 *           return_hessian[1][2] =
 *             (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
 *             (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) *
 *             std::pow(p(0) * (1.0 - p(0)), 2);
 *           return_hessian[2][2] =
 *             (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
 *             std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2);
 *         }
 *       else
 *         Assert(false, ExcNotImplemented());
 *   
 *       return return_hessian;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ImplementationofthecodeBiLaplacianLDGLiftcodeclass"></a> 
 * <h3>Implementation of the <code>BiLaplacianLDGLift</code> class</h3>
 * 
 
 * 
 * 
 * <a name="BiLaplacianLDGLiftBiLaplacianLDGLift"></a> 
 * <h4>BiLaplacianLDGLift::BiLaplacianLDGLift</h4>
 * 
 
 * 
 * In the constructor, we set the polynomial degree of the two finite element
 * spaces, we associate the corresponding DoF handlers to the triangulation,
 * and we set the two penalty coefficients.
 * 
 * @code
 *     template <int dim>
 *     BiLaplacianLDGLift<dim>::BiLaplacianLDGLift(const unsigned int n_refinements,
 *                                                 const unsigned int fe_degree,
 *                                                 const double penalty_jump_grad,
 *                                                 const double penalty_jump_val)
 *       : n_refinements(n_refinements)
 *       , fe(fe_degree)
 *       , dof_handler(triangulation)
 *       , fe_lift(FE_DGQ<dim>(fe_degree), dim * dim)
 *       , penalty_jump_grad(penalty_jump_grad)
 *       , penalty_jump_val(penalty_jump_val)
 *     {}
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftmake_grid"></a> 
 * <h4>BiLaplacianLDGLift::make_grid</h4>
 * 
 
 * 
 * To build a mesh for @f$\Omega=(0,1)^d@f$, we simply call the function
 * <code>GridGenerator::hyper_cube</code> and then refine it using
 * <code>refine_global</code>. The number of refinements is hard-coded
 * here.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::make_grid()
 *     {
 *       std::cout << "Building the mesh............." << std::endl;
 *   
 *       GridGenerator::hyper_cube(triangulation, 0.0, 1.0);
 *   
 *       triangulation.refine_global(n_refinements);
 *   
 *       std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftsetup_system"></a> 
 * <h4>BiLaplacianLDGLift::setup_system</h4>
 * 
 
 * 
 * In the following function, we set up the degrees of freedom, the sparsity
 * pattern, the size of the matrix @f$A@f$, and the size of the solution and
 * right-hand side vectors
 * @f$\boldsymbol{U}@f$ and @f$\boldsymbol{F}@f$. For the sparsity pattern, we cannot
 * directly use the function <code>DoFTools::make_flux_sparsity_pattern</code>
 * (as we would do for instance for the SIPG method) because we need to take
 * into account the interactions of a neighboring cell with another
 * neighboring cell as described in the introduction. The extended sparsity
 * pattern is built by iterating over all the active cells. For the current
 * cell, we collect all its degrees of freedom as well as the degrees of
 * freedom of all its neighboring cells, and then couple everything with
 * everything.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::setup_system()
 *     {
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl;
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *   
 *       const auto dofs_per_cell = fe.dofs_per_cell;
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           std::vector<types::global_dof_index> dofs(dofs_per_cell);
 *           cell->get_dof_indices(dofs);
 *   
 *           for (unsigned int f = 0; f < cell->n_faces(); ++f)
 *             if (!cell->face(f)->at_boundary())
 *               {
 *                 const auto neighbor_cell = cell->neighbor(f);
 *   
 *                 std::vector<types::global_dof_index> tmp(dofs_per_cell);
 *                 neighbor_cell->get_dof_indices(tmp);
 *   
 *                 dofs.insert(std::end(dofs), std::begin(tmp), std::end(tmp));
 *               }
 *   
 *           for (const auto i : dofs)
 *             for (const auto j : dofs)
 *               {
 *                 dsp.add(i, j);
 *                 dsp.add(j, i);
 *               }
 *         }
 *   
 *       sparsity_pattern.copy_from(dsp);
 *   
 *   
 *       matrix.reinit(sparsity_pattern);
 *       rhs.reinit(dof_handler.n_dofs());
 *   
 *       solution.reinit(dof_handler.n_dofs());
 *   
 * @endcode
 * 
 * At the end of the function, we output this sparsity pattern as
 * a scalable vector graphic. You can visualize it by loading this
 * file in most web browsers:
 * 
 * @code
 *       std::ofstream out("sparsity-pattern.svg");
 *       sparsity_pattern.print_svg(out);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftassemble_system"></a> 
 * <h4>BiLaplacianLDGLift::assemble_system</h4>
 * 
 
 * 
 * This function simply calls the two functions responsible
 * for the assembly of the matrix and the right-hand side.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::assemble_system()
 *     {
 *       std::cout << "Assembling the system............." << std::endl;
 *   
 *       assemble_matrix();
 *       assemble_rhs();
 *   
 *       std::cout << "Done. " << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftassemble_matrix"></a> 
 * <h4>BiLaplacianLDGLift::assemble_matrix</h4>
 * 
 
 * 
 * This function assembles the matrix @f$A@f$ whose entries are defined
 * by @f$A_{ij}=A_h(\varphi_j,\varphi_i)@f$ which involves the product of
 * discrete Hessians and the penalty terms.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::assemble_matrix()
 *     {
 *       matrix = 0;
 *   
 *       QGauss<dim>     quad(fe.degree + 1);
 *       QGauss<dim - 1> quad_face(fe.degree + 1);
 *   
 *       const unsigned int n_q_points      = quad.size();
 *       const unsigned int n_q_points_face = quad_face.size();
 *   
 *       FEValues<dim> fe_values(fe, quad, update_hessians | update_JxW_values);
 *   
 *       FEFaceValues<dim> fe_face(
 *         fe, quad_face, update_values | update_gradients | update_normal_vectors);
 *   
 *       FEFaceValues<dim> fe_face_neighbor(
 *         fe, quad_face, update_values | update_gradients | update_normal_vectors);
 *   
 *       const unsigned int n_dofs = fe_values.dofs_per_cell;
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(n_dofs);
 *       std::vector<types::global_dof_index> local_dof_indices_neighbor(n_dofs);
 *       std::vector<types::global_dof_index> local_dof_indices_neighbor_2(n_dofs);
 *   
 * @endcode
 * 
 * As indicated in the introduction, the following matrices are used for
 * the contributions of the products of the discrete Hessians.
 * 
 * @code
 *       FullMatrix<double> stiffness_matrix_cc(n_dofs,
 *                                              n_dofs); // interactions cell / cell
 *       FullMatrix<double> stiffness_matrix_cn(
 *         n_dofs, n_dofs); // interactions cell / neighbor
 *       FullMatrix<double> stiffness_matrix_nc(
 *         n_dofs, n_dofs); // interactions neighbor / cell
 *       FullMatrix<double> stiffness_matrix_nn(
 *         n_dofs, n_dofs); // interactions neighbor / neighbor
 *       FullMatrix<double> stiffness_matrix_n1n2(
 *         n_dofs, n_dofs); // interactions neighbor1 / neighbor2
 *       FullMatrix<double> stiffness_matrix_n2n1(
 *         n_dofs, n_dofs); // interactions neighbor2 / neighbor1
 *   
 * @endcode
 * 
 * The following matrices are used for the contributions of the two
 * penalty terms:
 * 
 * @code
 *       FullMatrix<double> ip_matrix_cc(n_dofs, n_dofs); // interactions cell / cell
 *       FullMatrix<double> ip_matrix_cn(n_dofs,
 *                                       n_dofs); // interactions cell / neighbor
 *       FullMatrix<double> ip_matrix_nc(n_dofs,
 *                                       n_dofs); // interactions neighbor / cell
 *       FullMatrix<double> ip_matrix_nn(n_dofs,
 *                                       n_dofs); // interactions neighbor / neighbor
 *   
 *       std::vector<std::vector<Tensor<2, dim>>> discrete_hessians(
 *         n_dofs, std::vector<Tensor<2, dim>>(n_q_points));
 *       std::vector<std::vector<std::vector<Tensor<2, dim>>>>
 *         discrete_hessians_neigh(GeometryInfo<dim>::faces_per_cell,
 *                                 discrete_hessians);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *           cell->get_dof_indices(local_dof_indices);
 *   
 * @endcode
 * 
 * We now compute all the discrete Hessians that are not vanishing
 * on the current cell, i.e., the discrete Hessian of all the basis
 * functions with support on the current cell or on one of its
 * neighbors.
 * 
 * @code
 *           compute_discrete_hessians(cell,
 *                                     discrete_hessians,
 *                                     discrete_hessians_neigh);
 *   
 * @endcode
 * 
 * First, we compute and add the interactions of the degrees of freedom
 * of the current cell.
 * 
 * @code
 *           stiffness_matrix_cc = 0;
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               const double dx = fe_values.JxW(q);
 *   
 *               for (unsigned int i = 0; i < n_dofs; ++i)
 *                 for (unsigned int j = 0; j < n_dofs; ++j)
 *                   {
 *                     const Tensor<2, dim> &H_i = discrete_hessians[i][q];
 *                     const Tensor<2, dim> &H_j = discrete_hessians[j][q];
 *   
 *                     stiffness_matrix_cc(i, j) += scalar_product(H_j, H_i) * dx;
 *                   }
 *             }
 *   
 *           for (unsigned int i = 0; i < n_dofs; ++i)
 *             for (unsigned int j = 0; j < n_dofs; ++j)
 *               {
 *                 matrix(local_dof_indices[i], local_dof_indices[j]) +=
 *                   stiffness_matrix_cc(i, j);
 *               }
 *   
 * @endcode
 * 
 * Next, we compute and add the interactions of the degrees of freedom
 * of the current cell with those of its neighbors. Note that the
 * interactions of the degrees of freedom of a neighbor with those of
 * the same neighbor are included here.
 * 
 * @code
 *           for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
 *             {
 *               const typename DoFHandler<dim>::face_iterator face =
 *                 cell->face(face_no);
 *   
 *               const bool at_boundary = face->at_boundary();
 *               if (!at_boundary)
 *                 {
 * @endcode
 * 
 * There is nothing to be done if boundary face (the liftings of
 * the Dirichlet BCs are accounted for in the assembly of the
 * RHS; in fact, nothing to be done in this program since we
 * prescribe homogeneous BCs).
 * 
 
 * 
 * 
 * @code
 *                   const typename DoFHandler<dim>::active_cell_iterator
 *                     neighbor_cell = cell->neighbor(face_no);
 *                   neighbor_cell->get_dof_indices(local_dof_indices_neighbor);
 *   
 *                   stiffness_matrix_cn = 0;
 *                   stiffness_matrix_nc = 0;
 *                   stiffness_matrix_nn = 0;
 *                   for (unsigned int q = 0; q < n_q_points; ++q)
 *                     {
 *                       const double dx = fe_values.JxW(q);
 *   
 *                       for (unsigned int i = 0; i < n_dofs; ++i)
 *                         {
 *                           for (unsigned int j = 0; j < n_dofs; ++j)
 *                             {
 *                               const Tensor<2, dim> &H_i = discrete_hessians[i][q];
 *                               const Tensor<2, dim> &H_j = discrete_hessians[j][q];
 *   
 *                               const Tensor<2, dim> &H_i_neigh =
 *                                 discrete_hessians_neigh[face_no][i][q];
 *                               const Tensor<2, dim> &H_j_neigh =
 *                                 discrete_hessians_neigh[face_no][j][q];
 *   
 *                               stiffness_matrix_cn(i, j) +=
 *                                 scalar_product(H_j_neigh, H_i) * dx;
 *                               stiffness_matrix_nc(i, j) +=
 *                                 scalar_product(H_j, H_i_neigh) * dx;
 *                               stiffness_matrix_nn(i, j) +=
 *                                 scalar_product(H_j_neigh, H_i_neigh) * dx;
 *                             }
 *                         }
 *                     }
 *   
 *                   for (unsigned int i = 0; i < n_dofs; ++i)
 *                     {
 *                       for (unsigned int j = 0; j < n_dofs; ++j)
 *                         {
 *                           matrix(local_dof_indices[i],
 *                                  local_dof_indices_neighbor[j]) +=
 *                             stiffness_matrix_cn(i, j);
 *                           matrix(local_dof_indices_neighbor[i],
 *                                  local_dof_indices[j]) +=
 *                             stiffness_matrix_nc(i, j);
 *                           matrix(local_dof_indices_neighbor[i],
 *                                  local_dof_indices_neighbor[j]) +=
 *                             stiffness_matrix_nn(i, j);
 *                         }
 *                     }
 *   
 *                 } // boundary check
 *             }     // for face
 *   
 * @endcode
 * 
 * We now compute and add the interactions of the degrees of freedom of
 * a neighboring cells with those of another neighboring cell (this is
 * where we need the extended sparsity pattern).
 * 
 * @code
 *           for (unsigned int face_no = 0; face_no < cell->n_faces() - 1; ++face_no)
 *             {
 *               const typename DoFHandler<dim>::face_iterator face =
 *                 cell->face(face_no);
 *   
 *               const bool at_boundary = face->at_boundary();
 *               if (!at_boundary)
 *                 { // nothing to be done if boundary face (the liftings of the
 * @endcode
 * 
 * Dirichlet BCs are accounted for in the assembly of the RHS;
 * in fact, nothing to be done in this program since we
 * prescribe homogeneous BCs)
 * 
 
 * 
 * 
 
 * 
 * 
 * @code
 *                   for (unsigned int face_no_2 = face_no + 1;
 *                        face_no_2 < cell->n_faces();
 *                        ++face_no_2)
 *                     {
 *                       const typename DoFHandler<dim>::face_iterator face_2 =
 *                         cell->face(face_no_2);
 *   
 *                       const bool at_boundary_2 = face_2->at_boundary();
 *                       if (!at_boundary_2)
 *                         {
 *                           const typename DoFHandler<dim>::active_cell_iterator
 *                             neighbor_cell = cell->neighbor(face_no);
 *                           neighbor_cell->get_dof_indices(
 *                             local_dof_indices_neighbor);
 *                           const typename DoFHandler<dim>::active_cell_iterator
 *                             neighbor_cell_2 = cell->neighbor(face_no_2);
 *                           neighbor_cell_2->get_dof_indices(
 *                             local_dof_indices_neighbor_2);
 *   
 *                           stiffness_matrix_n1n2 = 0;
 *                           stiffness_matrix_n2n1 = 0;
 *   
 *                           for (unsigned int q = 0; q < n_q_points; ++q)
 *                             {
 *                               const double dx = fe_values.JxW(q);
 *   
 *                               for (unsigned int i = 0; i < n_dofs; ++i)
 *                                 for (unsigned int j = 0; j < n_dofs; ++j)
 *                                   {
 *                                     const Tensor<2, dim> &H_i_neigh =
 *                                       discrete_hessians_neigh[face_no][i][q];
 *                                     const Tensor<2, dim> &H_j_neigh =
 *                                       discrete_hessians_neigh[face_no][j][q];
 *   
 *                                     const Tensor<2, dim> &H_i_neigh2 =
 *                                       discrete_hessians_neigh[face_no_2][i][q];
 *                                     const Tensor<2, dim> &H_j_neigh2 =
 *                                       discrete_hessians_neigh[face_no_2][j][q];
 *   
 *                                     stiffness_matrix_n1n2(i, j) +=
 *                                       scalar_product(H_j_neigh2, H_i_neigh) * dx;
 *                                     stiffness_matrix_n2n1(i, j) +=
 *                                       scalar_product(H_j_neigh, H_i_neigh2) * dx;
 *                                   }
 *                             }
 *   
 *                           for (unsigned int i = 0; i < n_dofs; ++i)
 *                             for (unsigned int j = 0; j < n_dofs; ++j)
 *                               {
 *                                 matrix(local_dof_indices_neighbor[i],
 *                                        local_dof_indices_neighbor_2[j]) +=
 *                                   stiffness_matrix_n1n2(i, j);
 *                                 matrix(local_dof_indices_neighbor_2[i],
 *                                        local_dof_indices_neighbor[j]) +=
 *                                   stiffness_matrix_n2n1(i, j);
 *                               }
 *                         } // boundary check face_2
 *                     }     // for face_2
 *                 }         // boundary check face_1
 *             }             // for face_1
 *   
 *   
 * @endcode
 * 
 * Finally, we compute and add the two penalty terms.
 * 
 * @code
 *           for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
 *             {
 *               const typename DoFHandler<dim>::face_iterator face =
 *                 cell->face(face_no);
 *   
 *               const double mesh_inv = 1.0 / face->diameter(); // h_e^{-1}
 *               const double mesh3_inv =
 *                 1.0 / std::pow(face->diameter(), 3); // ĥ_e^{-3}
 *   
 *               fe_face.reinit(cell, face_no);
 *   
 *               ip_matrix_cc = 0; // filled in any case (boundary or interior face)
 *   
 *               const bool at_boundary = face->at_boundary();
 *               if (at_boundary)
 *                 {
 *                   for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                     {
 *                       const double dx = fe_face.JxW(q);
 *   
 *                       for (unsigned int i = 0; i < n_dofs; ++i)
 *                         for (unsigned int j = 0; j < n_dofs; ++j)
 *                           {
 *                             ip_matrix_cc(i, j) += penalty_jump_grad * mesh_inv *
 *                                                   fe_face.shape_grad(j, q) *
 *                                                   fe_face.shape_grad(i, q) * dx;
 *                             ip_matrix_cc(i, j) += penalty_jump_val * mesh3_inv *
 *                                                   fe_face.shape_value(j, q) *
 *                                                   fe_face.shape_value(i, q) * dx;
 *                           }
 *                     }
 *                 }
 *               else
 *                 { // interior face
 *   
 *                   const typename DoFHandler<dim>::active_cell_iterator
 *                                      neighbor_cell = cell->neighbor(face_no);
 *                   const unsigned int face_no_neighbor =
 *                     cell->neighbor_of_neighbor(face_no);
 *   
 * @endcode
 * 
 * In the next step, we need to have a global way to compare the
 * cells in order to not calculate the same jump term twice:
 * 
 * @code
 *                   if (neighbor_cell->id() < cell->id())
 *                     continue; // skip this face (already considered)
 *                   else
 *                     {
 *                       fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
 *                       neighbor_cell->get_dof_indices(local_dof_indices_neighbor);
 *   
 *                       ip_matrix_cn = 0;
 *                       ip_matrix_nc = 0;
 *                       ip_matrix_nn = 0;
 *   
 *                       for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                         {
 *                           const double dx = fe_face.JxW(q);
 *   
 *                           for (unsigned int i = 0; i < n_dofs; ++i)
 *                             {
 *                               for (unsigned int j = 0; j < n_dofs; ++j)
 *                                 {
 *                                   ip_matrix_cc(i, j) +=
 *                                     penalty_jump_grad * mesh_inv *
 *                                     fe_face.shape_grad(j, q) *
 *                                     fe_face.shape_grad(i, q) * dx;
 *                                   ip_matrix_cc(i, j) +=
 *                                     penalty_jump_val * mesh3_inv *
 *                                     fe_face.shape_value(j, q) *
 *                                     fe_face.shape_value(i, q) * dx;
 *   
 *                                   ip_matrix_cn(i, j) -=
 *                                     penalty_jump_grad * mesh_inv *
 *                                     fe_face_neighbor.shape_grad(j, q) *
 *                                     fe_face.shape_grad(i, q) * dx;
 *                                   ip_matrix_cn(i, j) -=
 *                                     penalty_jump_val * mesh3_inv *
 *                                     fe_face_neighbor.shape_value(j, q) *
 *                                     fe_face.shape_value(i, q) * dx;
 *   
 *                                   ip_matrix_nc(i, j) -=
 *                                     penalty_jump_grad * mesh_inv *
 *                                     fe_face.shape_grad(j, q) *
 *                                     fe_face_neighbor.shape_grad(i, q) * dx;
 *                                   ip_matrix_nc(i, j) -=
 *                                     penalty_jump_val * mesh3_inv *
 *                                     fe_face.shape_value(j, q) *
 *                                     fe_face_neighbor.shape_value(i, q) * dx;
 *   
 *                                   ip_matrix_nn(i, j) +=
 *                                     penalty_jump_grad * mesh_inv *
 *                                     fe_face_neighbor.shape_grad(j, q) *
 *                                     fe_face_neighbor.shape_grad(i, q) * dx;
 *                                   ip_matrix_nn(i, j) +=
 *                                     penalty_jump_val * mesh3_inv *
 *                                     fe_face_neighbor.shape_value(j, q) *
 *                                     fe_face_neighbor.shape_value(i, q) * dx;
 *                                 }
 *                             }
 *                         }
 *                     } // face not visited yet
 *   
 *                 } // boundary check
 *   
 *               for (unsigned int i = 0; i < n_dofs; ++i)
 *                 {
 *                   for (unsigned int j = 0; j < n_dofs; ++j)
 *                     {
 *                       matrix(local_dof_indices[i], local_dof_indices[j]) +=
 *                         ip_matrix_cc(i, j);
 *                     }
 *                 }
 *   
 *               if (!at_boundary)
 *                 {
 *                   for (unsigned int i = 0; i < n_dofs; ++i)
 *                     {
 *                       for (unsigned int j = 0; j < n_dofs; ++j)
 *                         {
 *                           matrix(local_dof_indices[i],
 *                                  local_dof_indices_neighbor[j]) +=
 *                             ip_matrix_cn(i, j);
 *                           matrix(local_dof_indices_neighbor[i],
 *                                  local_dof_indices[j]) += ip_matrix_nc(i, j);
 *                           matrix(local_dof_indices_neighbor[i],
 *                                  local_dof_indices_neighbor[j]) +=
 *                             ip_matrix_nn(i, j);
 *                         }
 *                     }
 *                 }
 *   
 *             } // for face
 *         }     // for cell
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftassemble_rhs"></a> 
 * <h4>BiLaplacianLDGLift::assemble_rhs</h4>
 * 
 
 * 
 * This function assemble the right-hand side of the system. Since we consider
 * homogeneous Dirichlet boundary conditions, imposed weakly in the bilinear
 * form using the Nitsche approach, it only involves the contribution of the
 * forcing term @f$\int_{\Omega}fv_h@f$.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::assemble_rhs()
 *     {
 *       rhs = 0;
 *   
 *       const QGauss<dim> quad(fe.degree + 1);
 *       FEValues<dim>     fe_values(
 *         fe, quad, update_values | update_quadrature_points | update_JxW_values);
 *   
 *       const unsigned int n_dofs     = fe_values.dofs_per_cell;
 *       const unsigned int n_quad_pts = quad.size();
 *   
 *       const RightHandSide<dim> right_hand_side;
 *   
 *       Vector<double>                       local_rhs(n_dofs);
 *       std::vector<types::global_dof_index> local_dof_indices(n_dofs);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           local_rhs = 0;
 *           for (unsigned int q = 0; q < n_quad_pts; ++q)
 *             {
 *               const double dx = fe_values.JxW(q);
 *   
 *               for (unsigned int i = 0; i < n_dofs; ++i)
 *                 {
 *                   local_rhs(i) +=
 *                     right_hand_side.value(fe_values.quadrature_point(q)) *
 *                     fe_values.shape_value(i, q) * dx;
 *                 }
 *             }
 *   
 *           for (unsigned int i = 0; i < n_dofs; ++i)
 *             rhs(local_dof_indices[i]) += local_rhs(i);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftsolve"></a> 
 * <h4>BiLaplacianLDGLift::solve</h4>
 * 
 
 * 
 * To solve the linear system @f$A\boldsymbol{U}=\boldsymbol{F}@f$,
 * we proceed as in @ref step_74 "step-74" and use a direct method. We could
 * as well use an iterative method, for instance the conjugate
 * gradient method as in @ref step_3 "step-3".
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::solve()
 *     {
 *       SparseDirectUMFPACK A_direct;
 *       A_direct.initialize(matrix);
 *       A_direct.vmult(solution, rhs);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftcompute_errors"></a> 
 * <h4>BiLaplacianLDGLift::compute_errors</h4>
 * 
 
 * 
 * This function computes the discrete @f$H^2@f$, @f$H^1@f$ and @f$L^2@f$ norms of
 * the error @f$u-u_h@f$, where @f$u@f$ is the exact solution and @f$u_h@f$ is
 * the approximate solution. See the introduction for the definition
 * of the norms.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::compute_errors()
 *     {
 *       double error_H2 = 0;
 *       double error_H1 = 0;
 *       double error_L2 = 0;
 *   
 *       QGauss<dim>     quad(fe.degree + 1);
 *       QGauss<dim - 1> quad_face(fe.degree + 1);
 *   
 *       FEValues<dim> fe_values(fe,
 *                               quad,
 *                               update_values | update_gradients | update_hessians |
 *                                 update_quadrature_points | update_JxW_values);
 *   
 *       FEFaceValues<dim> fe_face(fe,
 *                                 quad_face,
 *                                 update_values | update_gradients |
 *                                   update_quadrature_points | update_JxW_values);
 *   
 *       FEFaceValues<dim> fe_face_neighbor(fe,
 *                                          quad_face,
 *                                          update_values | update_gradients);
 *   
 *       const unsigned int n_q_points      = quad.size();
 *       const unsigned int n_q_points_face = quad_face.size();
 *   
 * @endcode
 * 
 * We introduce some variables for the exact solution...
 * 
 * @code
 *       const ExactSolution<dim> u_exact;
 *   
 * @endcode
 * 
 * ...and for the approximate solution:
 * 
 * @code
 *       std::vector<double>         solution_values_cell(n_q_points);
 *       std::vector<Tensor<1, dim>> solution_gradients_cell(n_q_points);
 *       std::vector<Tensor<2, dim>> solution_hessians_cell(n_q_points);
 *   
 *       std::vector<double>         solution_values(n_q_points_face);
 *       std::vector<double>         solution_values_neigh(n_q_points_face);
 *       std::vector<Tensor<1, dim>> solution_gradients(n_q_points_face);
 *       std::vector<Tensor<1, dim>> solution_gradients_neigh(n_q_points_face);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *   
 *           fe_values.get_function_values(solution, solution_values_cell);
 *           fe_values.get_function_gradients(solution, solution_gradients_cell);
 *           fe_values.get_function_hessians(solution, solution_hessians_cell);
 *   
 * @endcode
 * 
 * We first add the <i>bulk</i> terms.
 * 
 * @code
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               const double dx = fe_values.JxW(q);
 *   
 *               error_H2 += (u_exact.hessian(fe_values.quadrature_point(q)) -
 *                            solution_hessians_cell[q])
 *                             .norm_square() *
 *                           dx;
 *               error_H1 += (u_exact.gradient(fe_values.quadrature_point(q)) -
 *                            solution_gradients_cell[q])
 *                             .norm_square() *
 *                           dx;
 *               error_L2 += std::pow(u_exact.value(fe_values.quadrature_point(q)) -
 *                                      solution_values_cell[q],
 *                                    2) *
 *                           dx;
 *             } // for quadrature points
 *   
 * @endcode
 * 
 * We then add the face contributions.
 * 
 * @code
 *           for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
 *             {
 *               const typename DoFHandler<dim>::face_iterator face =
 *                 cell->face(face_no);
 *   
 *               const double mesh_inv = 1.0 / face->diameter(); // h^{-1}
 *               const double mesh3_inv =
 *                 1.0 / std::pow(face->diameter(), 3); // h^{-3}
 *   
 *               fe_face.reinit(cell, face_no);
 *   
 *               fe_face.get_function_values(solution, solution_values);
 *               fe_face.get_function_gradients(solution, solution_gradients);
 *   
 *               const bool at_boundary = face->at_boundary();
 *               if (at_boundary)
 *                 {
 *                   for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                     {
 *                       const double dx = fe_face.JxW(q);
 *                       const double u_exact_q =
 *                         u_exact.value(fe_face.quadrature_point(q));
 *                       const Tensor<1, dim> u_exact_grad_q =
 *                         u_exact.gradient(fe_face.quadrature_point(q));
 *   
 *                       error_H2 +=
 *                         mesh_inv *
 *                         (u_exact_grad_q - solution_gradients[q]).norm_square() *
 *                         dx;
 *                       error_H2 += mesh3_inv *
 *                                   std::pow(u_exact_q - solution_values[q], 2) *
 *                                   dx;
 *                       error_H1 += mesh_inv *
 *                                   std::pow(u_exact_q - solution_values[q], 2) *
 *                                   dx;
 *                     }
 *                 }
 *               else
 *                 { // interior face
 *   
 *                   const typename DoFHandler<dim>::active_cell_iterator
 *                                      neighbor_cell = cell->neighbor(face_no);
 *                   const unsigned int face_no_neighbor =
 *                     cell->neighbor_of_neighbor(face_no);
 *   
 * @endcode
 * 
 * In the next step, we need to have a global way to compare the
 * cells in order to not calculate the same jump term twice:
 * 
 * @code
 *                   if (neighbor_cell->id() < cell->id())
 *                     continue; // skip this face (already considered)
 *                   else
 *                     {
 *                       fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
 *   
 *                       fe_face.get_function_values(solution, solution_values);
 *                       fe_face_neighbor.get_function_values(solution,
 *                                                            solution_values_neigh);
 *                       fe_face.get_function_gradients(solution,
 *                                                      solution_gradients);
 *                       fe_face_neighbor.get_function_gradients(
 *                         solution, solution_gradients_neigh);
 *   
 *                       for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                         {
 *                           const double dx = fe_face.JxW(q);
 *   
 * @endcode
 * 
 * To compute the jump term, we use the fact that
 * @f$\jump{u}=0@f$ and
 * @f$\jump{\nabla u}=\mathbf{0}@f$ since @f$u\in
 * H^2(\Omega)@f$.
 * 
 * @code
 *                           error_H2 +=
 *                             mesh_inv *
 *                             (solution_gradients_neigh[q] - solution_gradients[q])
 *                               .norm_square() *
 *                             dx;
 *                           error_H2 += mesh3_inv *
 *                                       std::pow(solution_values_neigh[q] -
 *                                                  solution_values[q],
 *                                                2) *
 *                                       dx;
 *                           error_H1 += mesh_inv *
 *                                       std::pow(solution_values_neigh[q] -
 *                                                  solution_values[q],
 *                                                2) *
 *                                       dx;
 *                         }
 *                     } // face not visited yet
 *   
 *                 } // boundary check
 *   
 *             } // for face
 *   
 *         } // for cell
 *   
 *       error_H2 = std::sqrt(error_H2);
 *       error_H1 = std::sqrt(error_H1);
 *       error_L2 = std::sqrt(error_L2);
 *   
 *       std::cout << "DG H2 norm of the error: " << error_H2 << std::endl;
 *       std::cout << "DG H1 norm of the error: " << error_H1 << std::endl;
 *       std::cout << "   L2 norm of the error: " << error_L2 << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftoutput_results"></a> 
 * <h4>BiLaplacianLDGLift::output_results</h4>
 * 
 
 * 
 * This function, which writes the solution to a vtk file,
 * is copied from @ref step_3 "step-3".
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::output_results() const
 *     {
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "solution");
 *       data_out.build_patches();
 *   
 *       std::ofstream output("solution.vtk");
 *       data_out.write_vtk(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftassemble_local_matrix"></a> 
 * <h4>BiLaplacianLDGLift::assemble_local_matrix</h4>
 * 
 
 * 
 * As already mentioned above, this function is used to assemble
 * the (local) mass matrices needed for the computations of the
 * lifting terms. We reiterate that only the basis functions with
 * support on the current cell are considered.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::assemble_local_matrix(
 *       const FEValues<dim> &fe_values_lift,
 *       const unsigned int   n_q_points,
 *       FullMatrix<double> & local_matrix)
 *     {
 *       const FEValuesExtractors::Tensor<2> tau_ext(0);
 *   
 *       const unsigned int n_dofs = fe_values_lift.dofs_per_cell;
 *   
 *       local_matrix = 0;
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           const double dx = fe_values_lift.JxW(q);
 *   
 *           for (unsigned int m = 0; m < n_dofs; ++m)
 *             for (unsigned int n = 0; n < n_dofs; ++n)
 *               {
 *                 local_matrix(m, n) +=
 *                   scalar_product(fe_values_lift[tau_ext].value(n, q),
 *                                  fe_values_lift[tau_ext].value(m, q)) *
 *                   dx;
 *               }
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftcompute_discrete_hessians"></a> 
 * <h4>BiLaplacianLDGLift::compute_discrete_hessians</h4>
 * 
 
 * 
 * This function is the main novelty of this program. It computes
 * the discrete Hessian @f$H_h(\varphi)@f$ for all the basis functions
 * @f$\varphi@f$ of @f$\mathbb{V}_h@f$ supported on the current cell and
 * those supported on a neighboring cell. The first argument
 * indicates the current cell (referring to the global DoFHandler
 * object), while the other two arguments are output variables that
 * are filled by this function.
 *   
 
 * 
 * In the following, we need to evaluate finite element shape
 * functions for the `fe_lift` finite element on the current
 * cell. Like for example in @ref step_61 "step-61", this "lift" space is defined
 * on every cell individually; as a consequence, there is no global
 * DoFHandler associated with this because we simply have no need
 * for such a DoFHandler. That leaves the question of what we should
 * initialize the FEValues and FEFaceValues objects with when we ask
 * them to evaluate shape functions of `fe_lift` on a concrete
 * cell. If we simply provide the first argument to this function,
 * `cell`, to FEValues::reinit(), we will receive an error message
 * that the given `cell` belongs to a DoFHandler that has a
 * different finite element associated with it than the `fe_lift`
 * object we want to evaluate. Fortunately, there is a relatively
 * easy solution: We can call FEValues::reinit() with a cell that
 * points into a triangulation -- the same cell, but not associated
 * with a DoFHandler, and consequently no finite element space. In
 * that case, FEValues::reinit() will skip the check that would
 * otherwise lead to an error message. All we have to do is to convert
 * the DoFHandler cell iterator into a Triangulation cell iterator;
 * see the first couple of lines of the function below to see how
 * this can be done.
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::compute_discrete_hessians(
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       std::vector<std::vector<Tensor<2, dim>>> &            discrete_hessians,
 *       std::vector<std::vector<std::vector<Tensor<2, dim>>>>
 *         &discrete_hessians_neigh)
 *     {
 *       const typename Triangulation<dim>::cell_iterator cell_lift =
 *         static_cast<typename Triangulation<dim>::cell_iterator>(cell);
 *   
 *       QGauss<dim>     quad(fe.degree + 1);
 *       QGauss<dim - 1> quad_face(fe.degree + 1);
 *   
 *       const unsigned int n_q_points      = quad.size();
 *       const unsigned int n_q_points_face = quad_face.size();
 *   
 * @endcode
 * 
 * The information we need from the basis functions of
 * @f$\mathbb{V}_h@f$: <code>fe_values</code> is needed to add
 * the broken Hessian part of the discrete Hessian, while
 * <code>fe_face</code> and <code>fe_face_neighbor</code>
 * are used to compute the right-hand sides for the local
 * problems.
 * 
 * @code
 *       FEValues<dim> fe_values(fe, quad, update_hessians | update_JxW_values);
 *   
 *       FEFaceValues<dim> fe_face(
 *         fe, quad_face, update_values | update_gradients | update_normal_vectors);
 *   
 *       FEFaceValues<dim> fe_face_neighbor(
 *         fe, quad_face, update_values | update_gradients | update_normal_vectors);
 *   
 *       const unsigned int n_dofs = fe_values.dofs_per_cell;
 *   
 * @endcode
 * 
 * The information needed from the basis functions
 * of the finite element space for the lifting terms:
 * <code>fe_values_lift</code> is used for the (local)
 * mass matrix (see @f$\boldsymbol{M}_c@f$ in the introduction),
 * while <code>fe_face_lift</code> is used to compute the
 * right-hand sides (see @f$\boldsymbol{G}_c@f$ for @f$b_e@f$).
 * 
 * @code
 *       FEValues<dim> fe_values_lift(fe_lift,
 *                                    quad,
 *                                    update_values | update_JxW_values);
 *   
 *       FEFaceValues<dim> fe_face_lift(
 *         fe_lift, quad_face, update_values | update_gradients | update_JxW_values);
 *   
 *       const FEValuesExtractors::Tensor<2> tau_ext(0);
 *   
 *       const unsigned int n_dofs_lift = fe_values_lift.dofs_per_cell;
 *       FullMatrix<double> local_matrix_lift(n_dofs_lift, n_dofs_lift);
 *   
 *       Vector<double> local_rhs_re(n_dofs_lift), local_rhs_be(n_dofs_lift),
 *         coeffs_re(n_dofs_lift), coeffs_be(n_dofs_lift), coeffs_tmp(n_dofs_lift);
 *   
 *       SolverControl            solver_control(1000, 1e-12);
 *       SolverCG<Vector<double>> solver(solver_control);
 *   
 *       double factor_avg; // 0.5 for interior faces, 1.0 for boundary faces
 *   
 *       fe_values.reinit(cell);
 *       fe_values_lift.reinit(cell_lift);
 *   
 * @endcode
 * 
 * We start by assembling the (local) mass matrix used for the computation
 * of the lifting terms @f$r_e@f$ and @f$b_e@f$.
 * 
 * @code
 *       assemble_local_matrix(fe_values_lift, n_q_points, local_matrix_lift);
 *   
 *       for (unsigned int i = 0; i < n_dofs; ++i)
 *         for (unsigned int q = 0; q < n_q_points; ++q)
 *           {
 *             discrete_hessians[i][q] = 0;
 *   
 *             for (unsigned int face_no = 0;
 *                  face_no < discrete_hessians_neigh.size();
 *                  ++face_no)
 *               {
 *                 discrete_hessians_neigh[face_no][i][q] = 0;
 *               }
 *           }
 *   
 * @endcode
 * 
 * In this loop, we compute the discrete Hessian at each quadrature point
 * @f$x_q@f$ of <code>cell</code> for each basis function supported on
 * <code>cell</code>, namely we fill-in the variable
 * <code>discrete_hessians[i][q]</code>. For the lifting terms, we need to
 * add the contribution of all the faces of <code>cell</code>.
 * 
 * @code
 *       for (unsigned int i = 0; i < n_dofs; ++i)
 *         {
 *           coeffs_re = 0;
 *           coeffs_be = 0;
 *   
 *           for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
 *             {
 *               const typename DoFHandler<dim>::face_iterator face =
 *                 cell->face(face_no);
 *   
 *               const bool at_boundary = face->at_boundary();
 *   
 * @endcode
 * 
 * Recall that by convention, the average of a function across a
 * boundary face @f$e@f$ reduces to the trace of the function on the
 * only element adjacent to @f$e@f$, namely there is no factor
 * @f$\frac{1}{2}@f$. We distinguish between the two cases (the current
 * face lies in the interior or on the boundary of the domain) using
 * the variable <code>factor_avg</code>.
 * 
 * @code
 *               factor_avg = 0.5;
 *               if (at_boundary)
 *                 {
 *                   factor_avg = 1.0;
 *                 }
 *   
 *               fe_face.reinit(cell, face_no);
 *               fe_face_lift.reinit(cell_lift, face_no);
 *   
 *               local_rhs_re = 0;
 *               for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                 {
 *                   const double         dx     = fe_face_lift.JxW(q);
 *                   const Tensor<1, dim> normal = fe_face.normal_vector(
 *                     q); // same as fe_face_lift.normal_vector(q)
 *   
 *                   for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                     {
 *                       local_rhs_re(m) +=
 *                         factor_avg *
 *                         (fe_face_lift[tau_ext].value(m, q) * normal) *
 *                         fe_face.shape_grad(i, q) * dx;
 *                     }
 *                 }
 *   
 * @endcode
 * 
 * Here, <code>local_rhs_be(m)</code> corresponds to @f$G_m@f$
 * introduced in the comments about the implementation of the
 * lifting @f$b_e@f$ in the case
 * @f$\varphi=\varphi^c@f$.
 * 
 * @code
 *               local_rhs_be = 0;
 *               for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                 {
 *                   const double         dx     = fe_face_lift.JxW(q);
 *                   const Tensor<1, dim> normal = fe_face.normal_vector(
 *                     q); // same as fe_face_lift.normal_vector(q)
 *   
 *                   for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                     {
 *                       local_rhs_be(m) += factor_avg *
 *                                          fe_face_lift[tau_ext].divergence(m, q) *
 *                                          normal * fe_face.shape_value(i, q) * dx;
 *                     }
 *                 }
 *   
 *               coeffs_tmp = 0;
 *               solver.solve(local_matrix_lift,
 *                            coeffs_tmp,
 *                            local_rhs_re,
 *                            PreconditionIdentity());
 *               coeffs_re += coeffs_tmp;
 *   
 *               coeffs_tmp = 0;
 *               solver.solve(local_matrix_lift,
 *                            coeffs_tmp,
 *                            local_rhs_be,
 *                            PreconditionIdentity());
 *               coeffs_be += coeffs_tmp;
 *   
 *             } // for face
 *   
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               discrete_hessians[i][q] += fe_values.shape_hessian(i, q);
 *   
 *               for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                 {
 *                   discrete_hessians[i][q] -=
 *                     coeffs_re[m] * fe_values_lift[tau_ext].value(m, q);
 *                 }
 *   
 *               for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                 {
 *                   discrete_hessians[i][q] +=
 *                     coeffs_be[m] * fe_values_lift[tau_ext].value(m, q);
 *                 }
 *             }
 *         } // for dof i
 *   
 *   
 *   
 * @endcode
 * 
 * In this loop, we compute the discrete Hessian at each quadrature point
 * @f$x_q@f$ of <code>cell</code> for each basis function supported on a
 * neighboring <code>neighbor_cell</code> of <code>cell</code>, namely we
 * fill-in the variable <code>discrete_hessians_neigh[face_no][i][q]</code>.
 * For the lifting terms, we only need to add the contribution of the
 * face adjacent to <code>cell</code> and <code>neighbor_cell</code>.
 * 
 * @code
 *       for (unsigned int face_no = 0; face_no < cell->n_faces(); ++face_no)
 *         {
 *           const typename DoFHandler<dim>::face_iterator face =
 *             cell->face(face_no);
 *   
 *           const bool at_boundary = face->at_boundary();
 *   
 *           if (!at_boundary)
 *             {
 * @endcode
 * 
 * For non-homogeneous Dirichlet BCs, we would need to
 * compute the lifting of the prescribed BC (see the
 * "Possible Extensions" section for more details).
 * 
 
 * 
 * 
 * @code
 *               const typename DoFHandler<2, dim>::active_cell_iterator
 *                                  neighbor_cell = cell->neighbor(face_no);
 *               const unsigned int face_no_neighbor =
 *                 cell->neighbor_of_neighbor(face_no);
 *               fe_face_neighbor.reinit(neighbor_cell, face_no_neighbor);
 *   
 *               for (unsigned int i = 0; i < n_dofs; ++i)
 *                 {
 *                   coeffs_re = 0;
 *                   coeffs_be = 0;
 *   
 *                   fe_face_lift.reinit(cell_lift, face_no);
 *   
 *                   local_rhs_re = 0;
 *                   for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                     {
 *                       const double         dx = fe_face_lift.JxW(q);
 *                       const Tensor<1, dim> normal =
 *                         fe_face_neighbor.normal_vector(q);
 *   
 *                       for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                         {
 *                           local_rhs_re(m) +=
 *                             0.5 * (fe_face_lift[tau_ext].value(m, q) * normal) *
 *                             fe_face_neighbor.shape_grad(i, q) * dx;
 *                         }
 *                     }
 *   
 * @endcode
 * 
 * Here, <code>local_rhs_be(m)</code> corresponds to @f$G_m@f$
 * introduced in the comments about the implementation of the
 * lifting @f$b_e@f$ in the case
 * @f$\varphi=\varphi^n@f$.
 * 
 * @code
 *                   local_rhs_be = 0;
 *                   for (unsigned int q = 0; q < n_q_points_face; ++q)
 *                     {
 *                       const double         dx = fe_face_lift.JxW(q);
 *                       const Tensor<1, dim> normal =
 *                         fe_face_neighbor.normal_vector(q);
 *   
 *                       for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                         {
 *                           local_rhs_be(m) +=
 *                             0.5 * fe_face_lift[tau_ext].divergence(m, q) *
 *                             normal * fe_face_neighbor.shape_value(i, q) * dx;
 *                         }
 *                     }
 *   
 *                   solver.solve(local_matrix_lift,
 *                                coeffs_re,
 *                                local_rhs_re,
 *                                PreconditionIdentity());
 *                   solver.solve(local_matrix_lift,
 *                                coeffs_be,
 *                                local_rhs_be,
 *                                PreconditionIdentity());
 *   
 *                   for (unsigned int q = 0; q < n_q_points; ++q)
 *                     {
 *                       for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                         {
 *                           discrete_hessians_neigh[face_no][i][q] -=
 *                             coeffs_re[m] * fe_values_lift[tau_ext].value(m, q);
 *                         }
 *   
 *                       for (unsigned int m = 0; m < n_dofs_lift; ++m)
 *                         {
 *                           discrete_hessians_neigh[face_no][i][q] +=
 *                             coeffs_be[m] * fe_values_lift[tau_ext].value(m, q);
 *                         }
 *                     }
 *   
 *                 } // for dof i
 *             }     // boundary check
 *         }         // for face
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="BiLaplacianLDGLiftrun"></a> 
 * <h4>BiLaplacianLDGLift::run</h4>
 * 
 * @code
 *     template <int dim>
 *     void BiLaplacianLDGLift<dim>::run()
 *     {
 *       make_grid();
 *   
 *       setup_system();
 *       assemble_system();
 *   
 *       solve();
 *   
 *       compute_errors();
 *       output_results();
 *     }
 *   
 *   } // namespace Step82
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * This is the <code>main</code> function. We define here the number of mesh
 * refinements, the polynomial degree for the two finite element spaces
 * (for the solution and the two liftings) and the two penalty coefficients.
 * We can also change the dimension to run the code in 3d.
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         const unsigned int n_ref = 3; // number of mesh refinements
 *   
 *         const unsigned int degree =
 *           2; // FE degree for u_h and the two lifting terms
 *   
 *         const double penalty_grad =
 *           1.0; // penalty coefficient for the jump of the gradients
 *         const double penalty_val =
 *           1.0; // penalty coefficient for the jump of the values
 *   
 *         Step82::BiLaplacianLDGLift<2> problem(n_ref,
 *                                               degree,
 *                                               penalty_grad,
 *                                               penalty_val);
 *   
 *         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>
 
 
 
When running the program, the sparsity pattern is written to an svg file, the solution is written to a vtk file, and some results are printed to the console. With the current setup, the output should read
 
@code
 
Number of active cells: 64
Number of degrees of freedom: 576
Assembling the system.............
Done.
DG H2 norm of the error: 0.0151063
DG H1 norm of the error: 0.000399747
   L2 norm of the error: 5.33856e-05
 
@endcode
 
This corresponds to the bi-Laplacian problem with the manufactured solution mentioned above for @f$d=2@f$, 3 refinements of the mesh, degree @f$k=2@f$, and @f$\gamma_0=\gamma_1=1@f$ for the penalty coefficients. By changing the number of refinements, we get the following results:
 
<table align="center" class="doxtable">
  <tr>
    <th>n_ref</th>
    <th>n_cells</th>
    <th>n_dofs</th>
    <th>error H2 </th>
    <th>rate</th>
    <th>error H1</th>
    <th>rate</th>
    <th>error L2</th>
    <th>rate</th>
  </tr>
  <tr>
    <td align="center">1</td>
    <td align="right">4</td>
    <td align="right">36</td>
    <td align="center">5.651e-02</td>
    <td align="center">--</td>
    <td align="center">3.366e-03</td>
    <td align="center">--</td>
    <td align="center">3.473e-04</td>
    <td align="center">--</td>
  </tr>
  <tr>
    <td align="center">2</td>
    <td align="right">16</td>
    <td align="right">144</td>
    <td align="center">3.095e-02</td>
    <td align="center">0.87</td>
    <td align="center">1.284e-03</td>
    <td align="center">1.39</td>
    <td align="center">1.369e-04</td>
    <td align="center">1.34</td>
  </tr>
  <tr>
    <td align="center">3</td>
    <td align="right">64</td>
    <td align="right">576</td>
    <td align="center">1.511e-02</td>
    <td align="center">1.03</td>
    <td align="center">3.997e-04</td>
    <td align="center">1.68</td>
    <td align="center">5.339e-05</td>
    <td align="center">1.36</td>
  </tr>
  <tr>
    <td align="center">4</td>
    <td align="right">256</td>
    <td align="right">2304</td>
    <td align="center">7.353e-03</td>
    <td align="center">1.04</td>
    <td align="center">1.129e-04</td>
    <td align="center">1.82</td>
    <td align="center">1.691e-05</td>
    <td align="center">1.66</td>
  </tr>
  <tr>
    <td align="center">5</td>
    <td align="right">1024</td>
    <td align="right">9216</td>
    <td align="center">3.609e-03</td>
    <td align="center">1.03</td>
    <td align="center">3.024e-05</td>
    <td align="center">1.90</td>
    <td align="center">4.789e-06</td>
    <td align="center">1.82</td>
  </tr>
  <tr>
    <td align="center">6</td>
    <td align="right">4096</td>
    <td align="right">36864</td>
    <td align="center">1.785e-03</td>
    <td align="center">1.02</td>
    <td align="center">7.850e-06</td>
    <td align="center">1.95</td>
    <td align="center">1.277e-06</td>
    <td align="center">1.91</td>
  </tr>
</table>
 
This matches the expected optimal convergence rates for the @f$H^2@f$ and
@f$H^1@f$ norms, but is sub-optimal for the @f$L_2@f$ norm. Incidentally, this
also matches the results seen in @ref step_47 "step-47" when using polynomial degree
@f$k=2@f$.
 
Indeed, just like in @ref step_47 "step-47", we can regain the optimal convergence
order if we set the polynomial degree of the finite elements to @f$k=3@f$
or higher. Here are the numbers for @f$k=3@f$:
 
<table align="center" class="doxtable">
  <tr> <th> n_ref </th> <th> n_cells </th>      <th> n_dofs </th>       <th> error H2 </th>     <th> rate </th> <th> error H1 </th>     <th> rate </th> <th> error L2 </th>     <th> rate</th> </tr>
  <tr> <td> 1 </td>     <td> 4 </td>    <td> 36 </td>   <td> 1.451e-02 </td>    <td> -- </td>   <td> 5.494e-04 </td>    <td> -- </td>   <td> 3.035e-05 </td>    <td> --</td> </tr>
  <tr> <td> 2 </td>     <td> 16 </td>   <td> 144 </td>  <td> 3.565e-03 </td>    <td> 2.02 </td> <td> 6.870e-05 </td>    <td> 3.00 </td> <td> 2.091e-06 </td>    <td> 3.86</td> </tr>
  <tr> <td> 3 </td>     <td> 64 </td>   <td> 576 </td>  <td> 8.891e-04 </td>    <td> 2.00 </td> <td> 8.584e-06 </td>    <td> 3.00 </td> <td> 1.352e-07 </td>    <td> 3.95</td> </tr>
  <tr> <td> 4 </td>     <td> 256 </td>  <td> 2304 </td> <td> 2.223e-04 </td>    <td> 2.00 </td> <td> 1.073e-06 </td>    <td> 3.00 </td> <td> 8.594e-09 </td>    <td> 3.98</td> </tr>
  <tr> <td> 5 </td>     <td> 1024 </td> <td> 9216 </td> <td> 5.560e-05 </td>    <td> 2.00 </td> <td> 1.341e-07 </td>    <td> 3.00 </td> <td> 5.418e-10 </td>    <td> 3.99</td> </tr>
  <tr> <td> 6 </td>     <td> 4096 </td> <td> 36864 </td>        <td> 1.390e-05 </td>    <td> 2.00 </td> <td> 1.676e-08 </td>    <td> 3.00 </td> <td> 3.245e-11 </td>    <td> 4.06</td> </tr>
</table>
 
 
<a name="Possibleextensions"></a><h3>Possible extensions</h3>
 
 
The code can be easily adapted to deal with the following cases:
 
<ol>
  <li>Non-homogeneous Dirichlet boundary conditions on (part of) the boundary @f$\partial \Omega@f$ of @f$\Omega@f$.</li>
  <li>Hanging-nodes (proceed as in @ref step_14 "step-14" to not visit a sub-face twice when computing the lifting terms in <code>compute_discrete_hessians</code> and the penalty terms in <code>assemble_matrix</code>).</li>
  <li>LDG method for the Poisson problem (use the discrete gradient consisting of the broken gradient and the lifting of the jump of @f$u_h@f$).</li>
</ol>
 
We give below additional details for the first of these points.
 
 
<a name="NonhomogeneousDirichletboundaryconditions"></a><h4>Non-homogeneous Dirichlet boundary conditions</h4>
 
If we prescribe non-homogeneous Dirichlet conditions, say
@f[
\nabla u=\mathbf{g} \quad \mbox{and} \quad u=g \qquad \mbox{on } \partial \Omega,
@f]
then the right-hand side @f$\boldsymbol{F}@f$ of the linear system needs to be modified as follows
@f[
F_i:=\int_{\Omega}f\varphi_i-\sum_{e\in\mathcal{E}_h^b}\int_{\Omega}r_e(\mathbf{g}):H_h(\varphi_i)+\sum_{e\in\mathcal{E}_h^b}\int_{\Omega}b_e(g):H_h(\varphi_i)+\gamma_1\sum_{e\in\mathcal{E}_h^b}h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0\sum_{e\in\mathcal{E}_h^b}h_e^{-3}\int_e g\varphi_i, \qquad 1\leq i \leq N_h.
@f]
Note that for any given index @f$i@f$, many of the terms are zero. Indeed, for @f$e\in \mathcal{E}_h^b@f$ we have @f${\rm supp}\,(r_e(\mathbf{g}))={\rm supp}\,(b_e(g))=K@f$, where @f$K@f$ is the element for which @f$e\subset\partial K@f$. Therefore, the liftings @f$r_e(\mathbf{g})@f$ and @f$b_e(g)@f$ contribute to @f$F_i@f$ only if @f$\varphi_i@f$ has support on @f$K@f$ or a neighbor of @f$K@f$. In other words, when integrating on a cell @f$K@f$, we need to add
@f[
\int_{K}f\varphi_i+\sum_{e\in\mathcal{E}_h^b, e\subset\partial K}\left[-\int_{K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{K}b_e(g):H_h(\varphi_i)+\gamma_1h_e^{-1}\int_e\mathbf{g}\cdot\nabla\varphi_i+\gamma_0h_e^{-3}\int_e g\varphi_i\right]
@f]
to @f$F_i@f$ for the indices @f$i@f$ such that @f$\varphi_i@f$ has support on @f$K@f$ and
@f[
\sum_{e\in\mathcal{E}_h^b, e\subset\partial K}\left[-\int_{K}r_e(\mathbf{g}):H_h(\varphi_i)+\int_{K}b_e(g):H_h(\varphi_i)\right]
@f]
to @f$F_i@f$ for the indices @f$i@f$ such that @f$\varphi_i@f$ has support on a neighbor of @f$K@f$.
 
@note
Note that we can easily consider the case where Dirichlet boundary conditions are imposed only on a subset @f$\emptyset\neq\Gamma_D\subset\partial \Omega@f$. In this case, we simply need to replace @f$\mathcal{E}_h^b@f$ by @f$\mathcal{E}_h^D\subset\mathcal{E}_h^b@f$ consisting of the faces belonging to @f$\Gamma_D@f$. This also affects the matrix @f$A@f$ (simply replace @f$\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^b@f$ by @f$\mathcal{E}_h=\mathcal{E}_h^0\cup\mathcal{E}_h^D@f$).
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-82.cc"
*/