step-39.h

Go to the documentation of this file.

 
 *         [&](const CellIterator     &cell,
 *             MatrixScratchData<dim> &scratch_data,
 *             MatrixCopyData<dim>    &copy) {
 *           copy.reinit(cell);
 *           scratch_data.fe_values.reinit(cell);
 *           MatrixIntegrator::cell<dim>(scratch_data.fe_values, copy.cell_matrix);
 *         },
 *         /* copier: */
 *         [&](const MatrixCopyData<dim> &copy) {
 *           matrix.add(copy.local_dof_indices, copy.cell_matrix);
 *           for (const auto &face : copy.face_data)
 *             {
 *               matrix.add(face.dof_indices_1, face.dof_indices_1, face.matrix_11);
 *               matrix.add(face.dof_indices_1, face.dof_indices_2, face.matrix_12);
 *               matrix.add(face.dof_indices_2, face.dof_indices_1, face.matrix_21);
 *               matrix.add(face.dof_indices_2, face.dof_indices_2, face.matrix_22);
 *             }
 *         },
 *         /* scratch and copy objects: */
 *         scratch,
 *         copy_data,
 *         /* where and how we want to integrate: */
 *         MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces |
 *           MeshWorker::assemble_own_interior_faces_once,
 *         /* boundary face worker: */
 *         [&](const CellIterator     &cell,
 *             const unsigned int      face_no,
 *             MatrixScratchData<dim> &scratch_data,
 *             MatrixCopyData<dim>    &copy) {
 *           scratch_data.boundary_fe_values.reinit(cell, face_no);
 *           MatrixIntegrator::boundary<dim>(cell,
 *                                           face_no,
 *                                           scratch_data.boundary_fe_values,
 *                                           copy.cell_matrix);
 *         },
 *         /* interior face worker: */
 *         [&](const CellIterator     &cell,
 *             const unsigned int      face_no,
 *             const unsigned int      subface_no,
 *             const CellIterator     &neighbor,
 *             const unsigned int      neighbor_face_no,
 *             const unsigned int      neighbor_subface_no,
 *             MatrixScratchData<dim> &scratch_data,
 *             MatrixCopyData<dim>    &copy) {
 *           FaceCopyData &face_copy = copy.emplace_face_data(cell, neighbor);
 *           if (subface_no == numbers::invalid_unsigned_int)
 *             {
 *               scratch_data.face_fe_values.reinit(cell, face_no);
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   MatrixIntegrator::face<dim>(cell,
 *                                               face_no,
 *                                               scratch_data.face_fe_values,
 *                                               neighbor,
 *                                               neighbor_face_no,
 *                                               scratch_data.neighbor_face_values,
 *                                               face_copy.matrix_11,
 *                                               face_copy.matrix_12,
 *                                               face_copy.matrix_21,
 *                                               face_copy.matrix_22);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   MatrixIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     scratch_data.face_fe_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_subface_values,
 *                     face_copy.matrix_11,
 *                     face_copy.matrix_12,
 *                     face_copy.matrix_21,
 *                     face_copy.matrix_22);
 *                 }
 *             }
 *           else
 *             {
 *               scratch_data.subface_values.reinit(cell, face_no, subface_no);
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   MatrixIntegrator::face<dim>(cell,
 *                                               face_no,
 *                                               scratch_data.subface_values,
 *                                               neighbor,
 *                                               neighbor_face_no,
 *                                               scratch_data.neighbor_face_values,
 *                                               face_copy.matrix_11,
 *                                               face_copy.matrix_12,
 *                                               face_copy.matrix_21,
 *                                               face_copy.matrix_22);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   MatrixIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     scratch_data.subface_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_subface_values,
 *                     face_copy.matrix_11,
 *                     face_copy.matrix_12,
 *                     face_copy.matrix_21,
 *                     face_copy.matrix_22);
 *                 }
 *             }
 *         });
 *     }
 *   
 *   
 * @endcode
 * 
 * Now, we do the same for the level matrices. Not too surprisingly, this
 * function looks like a twin of the previous one. The mesh loop again
 * traverses cells, boundary faces, and interior faces, and the worker
 * functions compute local matrix contributions without touching global data.
 * The scratch object provides the FEValues-like data needed for those local
 * computations, and the copy object collects the local cell matrix together
 * with the face blocks.
 *   
 
 * 
 * The main difference lies in the copier. Rather than assembling into a
 * single global matrix, it dispatches the local data into the appropriate
 * level matrix and, on refinement edges, into the up- and down-transfer
 * matrices that are required by the multigrid algorithm.
 * 
 * @code
 *     template <int dim>
 *     void InteriorPenaltyProblem<dim>::assemble_mg_matrix()
 *     {
 *       using CellIterator = typename DoFHandler<dim>::level_cell_iterator;
 *   
 *       const MatrixScratchData<dim> scratch(mapping,
 *                                            fe,
 *                                            QGauss<dim>(fe.degree + 1),
 *                                            QGauss<dim - 1>(fe.degree + 1),
 *                                            update_gradients | update_JxW_values,
 *                                            update_values | update_gradients |
 *                                              update_JxW_values |
 *                                              update_normal_vectors);
 *       const MatrixCopyData<dim>    copy_data(fe.n_dofs_per_cell());
 *   
 *       MeshWorker::mesh_loop(
 *         dof_handler.begin_mg(),
 *         dof_handler.end_mg(),
 *         /* cell worker: */
 *         [&](const CellIterator     &cell,
 *             MatrixScratchData<dim> &scratch_data,
 *             MatrixCopyData<dim>    &copy) {
 *           copy.reinit(cell, true);
 *           scratch_data.fe_values.reinit(cell);
 *           MatrixIntegrator::cell<dim>(scratch_data.fe_values, copy.cell_matrix);
 *         },
 *         /* copier: */
 *         [&](const MatrixCopyData<dim> &copy) {
 *           mg_matrix[copy.level].add(copy.local_dof_indices, copy.cell_matrix);
 *   
 *           for (const auto &face : copy.face_data)
 *             if (face.level_1 == face.level_2)
 *               {
 *                 mg_matrix[face.level_1].add(face.dof_indices_1,
 *                                             face.dof_indices_1,
 *                                             face.matrix_11);
 *                 mg_matrix[face.level_1].add(face.dof_indices_1,
 *                                             face.dof_indices_2,
 *                                             face.matrix_12);
 *                 mg_matrix[face.level_1].add(face.dof_indices_2,
 *                                             face.dof_indices_1,
 *                                             face.matrix_21);
 *                 mg_matrix[face.level_1].add(face.dof_indices_2,
 *                                             face.dof_indices_2,
 *                                             face.matrix_22);
 *               }
 *             else
 *               {
 *                 Assert(face.level_1 > face.level_2, ExcInternalError());
 *   
 *                 mg_matrix[face.level_1].add(face.dof_indices_1,
 *                                             face.dof_indices_1,
 *                                             face.matrix_11);
 *   
 *                 for (unsigned int j = 0; j < face.dof_indices_2.size(); ++j)
 *                   for (unsigned int k = 0; k < face.dof_indices_1.size(); ++k)
 *                     {
 *                       mg_matrix_dg_up[face.level_1].add(face.dof_indices_2[j],
 *                                                         face.dof_indices_1[k],
 *                                                         face.matrix_12(k, j));
 *                       mg_matrix_dg_down[face.level_1].add(face.dof_indices_2[j],
 *                                                           face.dof_indices_1[k],
 *                                                           face.matrix_21(j, k));
 *                     }
 *               }
 *         },
 *         /* scratch and copy objects: */
 *         scratch,
 *         copy_data,
 *         /* where and how we want to integrate: */
 *         MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces |
 *           MeshWorker::assemble_own_interior_faces_once,
 *         /* boundary face worker: */
 *         [&](const CellIterator     &cell,
 *             const unsigned int      face_no,
 *             MatrixScratchData<dim> &scratch_data,
 *             MatrixCopyData<dim>    &copy) {
 *           scratch_data.boundary_fe_values.reinit(cell, face_no);
 *           MatrixIntegrator::boundary<dim>(cell,
 *                                           face_no,
 *                                           scratch_data.boundary_fe_values,
 *                                           copy.cell_matrix);
 *         },
 *         /* interior face worker: */
 *         [&](const CellIterator     &cell,
 *             const unsigned int      face_no,
 *             const unsigned int      subface_no,
 *             const CellIterator     &neighbor,
 *             const unsigned int      neighbor_face_no,
 *             const unsigned int      neighbor_subface_no,
 *             MatrixScratchData<dim> &scratch_data,
 *             MatrixCopyData<dim>    &copy) {
 *           FaceCopyData &face_copy = copy.emplace_face_data(cell, neighbor, true);
 *           if (subface_no == numbers::invalid_unsigned_int)
 *             {
 *               scratch_data.face_fe_values.reinit(cell, face_no);
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   MatrixIntegrator::face<dim>(cell,
 *                                               face_no,
 *                                               scratch_data.face_fe_values,
 *                                               neighbor,
 *                                               neighbor_face_no,
 *                                               scratch_data.neighbor_face_values,
 *                                               face_copy.matrix_11,
 *                                               face_copy.matrix_12,
 *                                               face_copy.matrix_21,
 *                                               face_copy.matrix_22);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   MatrixIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     scratch_data.face_fe_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_subface_values,
 *                     face_copy.matrix_11,
 *                     face_copy.matrix_12,
 *                     face_copy.matrix_21,
 *                     face_copy.matrix_22);
 *                 }
 *             }
 *           else
 *             {
 *               scratch_data.subface_values.reinit(cell, face_no, subface_no);
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   MatrixIntegrator::face<dim>(cell,
 *                                               face_no,
 *                                               scratch_data.subface_values,
 *                                               neighbor,
 *                                               neighbor_face_no,
 *                                               scratch_data.neighbor_face_values,
 *                                               face_copy.matrix_11,
 *                                               face_copy.matrix_12,
 *                                               face_copy.matrix_21,
 *                                               face_copy.matrix_22);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   MatrixIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     scratch_data.subface_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_subface_values,
 *                     face_copy.matrix_11,
 *                     face_copy.matrix_12,
 *                     face_copy.matrix_21,
 *                     face_copy.matrix_22);
 *                 }
 *             }
 *         });
 *     }
 *   
 *   
 * @endcode
 * 
 * Here we have another clone of the assembly function. The difference to
 * assembling the system matrix consists in that we assemble a vector here.
 *   
 
 * 
 * The mesh loop still uses the same worker/copier split. The cell worker is
 * only responsible for initializing the copy object for the current cell,
 * whereas the actual local work happens on boundary faces: there, the worker
 * evaluates the inhomogeneous boundary values and accumulates the associated
 * Nitsche terms into the local right-hand-side vector. The copier then adds
 * that local vector to the global right-hand side.
 * 
 * @code
 *     template <int dim>
 *     void InteriorPenaltyProblem<dim>::assemble_right_hand_side()
 *     {
 *       using CellIterator = typename DoFHandler<dim>::active_cell_iterator;
 *   
 *       const RightHandSideScratchData<dim> scratch(mapping,
 *                                                   fe,
 *                                                   QGauss<dim - 1>(fe.degree + 1),
 *                                                   update_quadrature_points |
 *                                                     update_values |
 *                                                     update_gradients |
 *                                                     update_JxW_values |
 *                                                     update_normal_vectors);
 *       const RightHandSideCopyData         copy_data(fe.n_dofs_per_cell());
 *   
 *       MeshWorker::mesh_loop(
 *         dof_handler.begin_active(),
 *         dof_handler.end(),
 *         /* cell worker: */
 *         [&](const CellIterator &cell,
 *             RightHandSideScratchData<dim> &,
 *             RightHandSideCopyData &copy) { copy.reinit(cell); },
 *         /* copier: */
 *         [&](const RightHandSideCopyData &copy) {
 *           right_hand_side.add(copy.local_dof_indices, copy.cell_rhs);
 *         },
 *         /* scratch and copy objects: */
 *         scratch,
 *         copy_data,
 *         /* where and how we want to integrate: */
 *         MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces,
 *         /* boundary face worker: */
 *         [&](const CellIterator            &cell,
 *             const unsigned int             face_no,
 *             RightHandSideScratchData<dim> &scratch_data,
 *             RightHandSideCopyData         &copy) {
 *           scratch_data.boundary_fe_values.reinit(cell, face_no);
 *           exact_solution.value_list(
 *             scratch_data.boundary_fe_values.get_quadrature_points(),
 *             scratch_data.boundary_values);
 *           RHSIntegrator::boundary<dim>(cell,
 *                                        face_no,
 *                                        scratch_data.boundary_fe_values,
 *                                        scratch_data.boundary_values,
 *                                        copy.cell_rhs);
 *         });
 *   
 *       right_hand_side *= -1.;
 *     }
 *   
 *   
 * @endcode
 * 
 * Now that we have coded all functions building the discrete linear system,
 * it is about time that we actually solve it.
 * 
 * @code
 *     template <int dim>
 *     void InteriorPenaltyProblem<dim>::solve()
 *     {
 * @endcode
 * 
 * The solver of choice is conjugate gradient.
 * 
 * @code
 *       SolverControl            control(1000, 1.e-12);
 *       SolverCG<Vector<double>> solver(control);
 *   
 * @endcode
 * 
 * Now we are setting up the components of the multilevel
 * preconditioner. First, we need transfer between grid levels. The object
 * we are using here generates sparse matrices for these transfers.
 * 
 * @code
 *       MGTransferPrebuilt<Vector<double>> mg_transfer;
 *       mg_transfer.build(dof_handler);
 *   
 * @endcode
 * 
 * Then, we need an exact solver for the matrix on the coarsest level.
 * 
 * @code
 *       FullMatrix<double> coarse_matrix;
 *       coarse_matrix.copy_from(mg_matrix[0]);
 *       MGCoarseGridHouseholder<double, Vector<double>> mg_coarse;
 *       mg_coarse.initialize(coarse_matrix);
 *   
 * @endcode
 * 
 * While transfer and coarse grid solver are pretty much generic, more
 * flexibility is offered for the smoother. First, we choose Gauss-Seidel
 * as our smoothing method.
 * 
 * @code
 *       GrowingVectorMemory<Vector<double>> mem;
 *       using RELAXATION = PreconditionSOR<SparseMatrix<double>>;
 *       mg::SmootherRelaxation<RELAXATION, Vector<double>> mg_smoother;
 *       RELAXATION::AdditionalData                         smoother_data(1.);
 *       mg_smoother.initialize(mg_matrix, smoother_data);
 *   
 * @endcode
 * 
 * Do two smoothing steps on each level.
 * 
 * @code
 *       mg_smoother.set_steps(2);
 * @endcode
 * 
 * Since the SOR method is not symmetric, but we use conjugate gradient
 * iteration below, here is a trick to make the multilevel preconditioner
 * a symmetric operator even for nonsymmetric smoothers.
 * 
 * @code
 *       mg_smoother.set_symmetric(true);
 * @endcode
 * 
 * The smoother class optionally implements the variable V-cycle, which we
 * do not want here.
 * 
 * @code
 *       mg_smoother.set_variable(false);
 *   
 * @endcode
 * 
 * Finally, we must wrap our matrices in an object having the required
 * multiplication functions.
 * 
 * @code
 *       mg::Matrix<Vector<double>> mgmatrix(mg_matrix);
 *       mg::Matrix<Vector<double>> mgdown(mg_matrix_dg_down);
 *       mg::Matrix<Vector<double>> mgup(mg_matrix_dg_up);
 *   
 * @endcode
 * 
 * Now, we are ready to set up the V-cycle operator and the multilevel
 * preconditioner.
 * 
 * @code
 *       Multigrid<Vector<double>> mg(
 *         mgmatrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
 * @endcode
 * 
 * Let us not forget the edge matrices needed because of the adaptive
 * refinement.
 * 
 * @code
 *       mg.set_edge_flux_matrices(mgdown, mgup);
 *   
 * @endcode
 * 
 * After all preparations, wrap the Multigrid object into another object,
 * which can be used as a regular preconditioner,
 * 
 * @code
 *       PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<Vector<double>>>
 *         preconditioner(dof_handler, mg, mg_transfer);
 * @endcode
 * 
 * and use it to solve the system.
 * 
 * @code
 *       solver.solve(matrix, solution, right_hand_side, preconditioner);
 *   
 *       std::cout << "Converged in " << control.last_step() << " iterations"
 *                 << std::endl;
 *     }
 *   
 *   
 * @endcode
 * 
 * The next function estimates the error. The big difference to the previous
 * mesh loop functions is that we now also read from the discrete solution
 * vector. The results of the estimator are stored in a vector with one entry
 * per cell.
 *   
 
 * 
 * As before, MeshWorker::mesh_loop() separates the local work from the
 * accumulation into the global output vector. The workers evaluate the
 * current solution on cells and faces, compute the cell, boundary, and jump
 * contributions to the estimator, and store them in the copy object. The
 * copier then distributes these local indicators to the per-cell entries of
 * the global estimator vector, splitting face terms evenly between the two
 * cells that share the face.
 * 
 * @code
 *     template <int dim>
 *     double InteriorPenaltyProblem<dim>::estimate()
 *     {
 *       estimates.block(0).reinit(triangulation.n_active_cells());
 *       using CellIterator = typename DoFHandler<dim>::active_cell_iterator;
 *   
 *       const unsigned int n_gauss_points =
 *         dof_handler.get_fe().tensor_degree() + 1;
 *       const EstimatorScratchData<dim> scratch(mapping,
 *                                               fe,
 *                                               QGauss<dim>(n_gauss_points),
 *                                               QGauss<dim - 1>(n_gauss_points + 1),
 *                                               QGauss<dim - 1>(n_gauss_points),
 *                                               update_hessians | update_JxW_values,
 *                                               update_quadrature_points |
 *                                                 update_values | update_gradients |
 *                                                 update_JxW_values |
 *                                                 update_normal_vectors,
 *                                               update_values | update_gradients |
 *                                                 update_JxW_values |
 *                                                 update_normal_vectors);
 *       const ErrorCopyData<1>          copy_data;
 *   
 *       MeshWorker::mesh_loop(
 *         dof_handler.begin_active(),
 *         dof_handler.end(),
 *         /* cell worker: */
 *         [&](const CellIterator        &cell,
 *             EstimatorScratchData<dim> &scratch_data,
 *             ErrorCopyData<1>          &copy) {
 *           copy.reinit(cell);
 *           scratch_data.fe_values.reinit(cell);
 *           const FEValues<dim> &fe_values = scratch_data.fe_values;
 *           fe_values.get_function_hessians(solution, scratch_data.cell_hessians);
 *           copy.cell_values[0] =
 *             Estimator::cell<dim>(cell, fe_values, scratch_data.cell_hessians);
 *         },
 *         /* copier: */
 *         [&](const ErrorCopyData<1> &copy) {
 *           estimates.block(0)(copy.cell_index) += copy.cell_values[0];
 *           for (const auto &face : copy.face_data)
 *             {
 *               estimates.block(0)(face.cell_index_1) += 0.5 * face.values[0];
 *               estimates.block(0)(face.cell_index_2) += 0.5 * face.values[0];
 *             }
 *         },
 *         /* scratch and copy objects: */
 *         scratch,
 *         copy_data,
 *         /* where and how we want to integrate: */
 *         MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces |
 *           MeshWorker::assemble_own_interior_faces_once,
 *         /* boundary face worker: */
 *         [&](const CellIterator        &cell,
 *             const unsigned int         face_no,
 *             EstimatorScratchData<dim> &scratch_data,
 *             ErrorCopyData<1>          &copy) {
 *           scratch_data.boundary_fe_values.reinit(cell, face_no);
 *           const FEFaceValues<dim> &fe_face_values =
 *             scratch_data.boundary_fe_values;
 *           fe_face_values.get_function_values(
 *             solution, scratch_data.boundary_solution_values);
 *           exact_solution.value_list(fe_face_values.get_quadrature_points(),
 *                                     scratch_data.boundary_exact_values);
 *           copy.cell_values[0] +=
 *             Estimator::boundary<dim>(cell,
 *                                      face_no,
 *                                      fe_face_values,
 *                                      scratch_data.boundary_solution_values,
 *                                      scratch_data.boundary_exact_values);
 *         },
 *         /* interior face worker: */
 *         [&](const CellIterator        &cell,
 *             const unsigned int         face_no,
 *             const unsigned int         subface_no,
 *             const CellIterator        &neighbor,
 *             const unsigned int         neighbor_face_no,
 *             const unsigned int         neighbor_subface_no,
 *             EstimatorScratchData<dim> &scratch_data,
 *             ErrorCopyData<1>          &copy) {
 *           auto &face_data = copy.emplace_face_data(cell, neighbor);
 *   
 *           if (subface_no == numbers::invalid_unsigned_int)
 *             {
 *               scratch_data.face_fe_values.reinit(cell, face_no);
 *               const FEFaceValuesBase<dim> &fe_face_values =
 *                 scratch_data.face_fe_values;
 *   
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_face_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *                   fe_face_values.get_function_gradients(
 *                     solution, scratch_data.face_solution_gradients);
 *                   neighbor_fe_face_values.get_function_gradients(
 *                     solution, scratch_data.neighbor_face_solution_gradients);
 *   
 *                   face_data.values[0] = Estimator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     scratch_data.face_solution_gradients,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values,
 *                     scratch_data.neighbor_face_solution_gradients);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_subface_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *                   fe_face_values.get_function_gradients(
 *                     solution, scratch_data.face_solution_gradients);
 *                   neighbor_fe_face_values.get_function_gradients(
 *                     solution, scratch_data.neighbor_face_solution_gradients);
 *   
 *                   face_data.values[0] = Estimator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     scratch_data.face_solution_gradients,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values,
 *                     scratch_data.neighbor_face_solution_gradients);
 *                 }
 *             }
 *           else
 *             {
 *               scratch_data.subface_values.reinit(cell, face_no, subface_no);
 *               const FEFaceValuesBase<dim> &fe_face_values =
 *                 scratch_data.subface_values;
 *   
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_face_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *                   fe_face_values.get_function_gradients(
 *                     solution, scratch_data.face_solution_gradients);
 *                   neighbor_fe_face_values.get_function_gradients(
 *                     solution, scratch_data.neighbor_face_solution_gradients);
 *   
 *                   face_data.values[0] = Estimator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     scratch_data.face_solution_gradients,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values,
 *                     scratch_data.neighbor_face_solution_gradients);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_subface_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *                   fe_face_values.get_function_gradients(
 *                     solution, scratch_data.face_solution_gradients);
 *                   neighbor_fe_face_values.get_function_gradients(
 *                     solution, scratch_data.neighbor_face_solution_gradients);
 *   
 *                   face_data.values[0] = Estimator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     scratch_data.face_solution_gradients,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values,
 *                     scratch_data.neighbor_face_solution_gradients);
 *                 }
 *             }
 *         });
 *   
 *       return estimates.block(0).l2_norm();
 *     }
 *   
 * @endcode
 * 
 * Here we compare our finite element solution with the known exact solution
 * and compute the mean quadratic error of the gradient and the function
 * itself. This function is a close relative of the estimation function right
 * above: the mesh loop again visits cells, boundary faces, and interior
 * faces; the workers evaluate local quantities with the help of the scratch
 * object; and the copier writes the resulting indicators into global data
 * structures only after the local computation is finished.
 *   
 
 * 
 * Since we compute the error in the energy and the
 * <i>L<sup>2</sup></i>-norm, respectively, our block vector needs two
 * blocks here. Consequently, the copy object stores two local values per
 * cell and per face contribution, and the copier accumulates each of them
 * into the corresponding block.
 * 
 * @code
 *     template <int dim>
 *     void InteriorPenaltyProblem<dim>::error()
 *     {
 *       BlockVector<double> errors(2);
 *       errors.block(0).reinit(triangulation.n_active_cells());
 *       errors.block(1).reinit(triangulation.n_active_cells());
 *       using CellIterator = typename DoFHandler<dim>::active_cell_iterator;
 *   
 *       const unsigned int n_gauss_points =
 *         dof_handler.get_fe().tensor_degree() + 1;
 *       const ErrorScratchData<dim> scratch(mapping,
 *                                           fe,
 *                                           QGauss<dim>(n_gauss_points),
 *                                           QGauss<dim - 1>(n_gauss_points + 1),
 *                                           QGauss<dim - 1>(n_gauss_points),
 *                                           update_quadrature_points |
 *                                             update_values | update_gradients |
 *                                             update_JxW_values,
 *                                           update_quadrature_points |
 *                                             update_values | update_JxW_values,
 *                                           update_values | update_JxW_values);
 *       const ErrorCopyData<2>      copy_data;
 *   
 *       MeshWorker::mesh_loop(
 *         dof_handler.begin_active(),
 *         dof_handler.end(),
 *         /* cell worker: */
 *         [&](const CellIterator    &cell,
 *             ErrorScratchData<dim> &scratch_data,
 *             ErrorCopyData<2>      &copy) {
 *           copy.reinit(cell);
 *           scratch_data.fe_values.reinit(cell);
 *           const FEValues<dim> &fe_values = scratch_data.fe_values;
 *           fe_values.get_function_values(solution,
 *                                         scratch_data.cell_solution_values);
 *           fe_values.get_function_gradients(solution,
 *                                            scratch_data.cell_solution_gradients);
 *           exact_solution.value_list(fe_values.get_quadrature_points(),
 *                                     scratch_data.cell_exact_values);
 *           exact_solution.gradient_list(fe_values.get_quadrature_points(),
 *                                        scratch_data.cell_exact_gradients);
 *           copy.cell_values =
 *             ErrorIntegrator::cell<dim>(fe_values,
 *                                        scratch_data.cell_solution_values,
 *                                        scratch_data.cell_solution_gradients,
 *                                        scratch_data.cell_exact_values,
 *                                        scratch_data.cell_exact_gradients);
 *         },
 *         /* copier: */
 *         [&](const ErrorCopyData<2> &copy) {
 *           errors.block(0)(copy.cell_index) += copy.cell_values[0];
 *           errors.block(1)(copy.cell_index) += copy.cell_values[1];
 *           for (const auto &face : copy.face_data)
 *             {
 *               errors.block(0)(face.cell_index_1) += 0.5 * face.values[0];
 *               errors.block(0)(face.cell_index_2) += 0.5 * face.values[0];
 *             }
 *         },
 *         /* scratch and copy objects: */
 *         scratch,
 *         copy_data,
 *         /* where and how we want to integrate: */
 *         MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces |
 *           MeshWorker::assemble_own_interior_faces_once,
 *         /* boundary face worker: */
 *         [&](const CellIterator    &cell,
 *             const unsigned int     face_no,
 *             ErrorScratchData<dim> &scratch_data,
 *             ErrorCopyData<2>      &copy) {
 *           scratch_data.boundary_fe_values.reinit(cell, face_no);
 *           const FEFaceValues<dim> &fe_face_values =
 *             scratch_data.boundary_fe_values;
 *           fe_face_values.get_function_values(
 *             solution, scratch_data.boundary_solution_values);
 *           exact_solution.value_list(fe_face_values.get_quadrature_points(),
 *                                     scratch_data.boundary_exact_values);
 *           copy.cell_values[0] +=
 *             ErrorIntegrator::boundary<dim>(cell,
 *                                            face_no,
 *                                            fe_face_values,
 *                                            scratch_data.boundary_solution_values,
 *                                            scratch_data.boundary_exact_values);
 *         },
 *         /* interior face worker: */
 *         [&](const CellIterator    &cell,
 *             const unsigned int     face_no,
 *             const unsigned int     subface_no,
 *             const CellIterator    &neighbor,
 *             const unsigned int     neighbor_face_no,
 *             const unsigned int     neighbor_subface_no,
 *             ErrorScratchData<dim> &scratch_data,
 *             ErrorCopyData<2>      &copy) {
 *           auto &face_data = copy.emplace_face_data(cell, neighbor);
 *   
 *           if (subface_no == numbers::invalid_unsigned_int)
 *             {
 *               scratch_data.face_fe_values.reinit(cell, face_no);
 *               const FEFaceValuesBase<dim> &fe_face_values =
 *                 scratch_data.face_fe_values;
 *   
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_face_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *   
 *                   face_data.values[0] = ErrorIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_subface_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *   
 *                   face_data.values[0] = ErrorIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values);
 *                 }
 *             }
 *           else
 *             {
 *               scratch_data.subface_values.reinit(cell, face_no, subface_no);
 *               const FEFaceValuesBase<dim> &fe_face_values =
 *                 scratch_data.subface_values;
 *   
 *               if (neighbor_subface_no == numbers::invalid_unsigned_int)
 *                 {
 *                   scratch_data.neighbor_face_values.reinit(neighbor,
 *                                                            neighbor_face_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_face_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *   
 *                   face_data.values[0] = ErrorIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values);
 *                 }
 *               else
 *                 {
 *                   scratch_data.neighbor_subface_values.reinit(
 *                     neighbor, neighbor_face_no, neighbor_subface_no);
 *                   const FEFaceValuesBase<dim> &neighbor_fe_face_values =
 *                     scratch_data.neighbor_subface_values;
 *   
 *                   fe_face_values.get_function_values(
 *                     solution, scratch_data.face_solution_values);
 *                   neighbor_fe_face_values.get_function_values(
 *                     solution, scratch_data.neighbor_face_solution_values);
 *   
 *                   face_data.values[0] = ErrorIntegrator::face<dim>(
 *                     cell,
 *                     face_no,
 *                     fe_face_values,
 *                     scratch_data.face_solution_values,
 *                     neighbor,
 *                     neighbor_face_no,
 *                     scratch_data.neighbor_face_solution_values);
 *                 }
 *             }
 *         });
 *   
 *       std::cout << "energy-error: " << errors.block(0).l2_norm() << std::endl;
 *       std::cout << "L2-error:     " << errors.block(1).l2_norm() << std::endl;
 *     }
 *   
 *   
 * @endcode
 * 
 * Create graphical output. We produce the filename by collating the
 * name from its various components, including the refinement cycle
 * that we output with two digits.
 * 
 * @code
 *     template <int dim>
 *     void
 *     InteriorPenaltyProblem<dim>::output_results(const unsigned int cycle) const
 *     {
 *       const std::string filename =
 *         "sol-" + Utilities::int_to_string(cycle, 2) + ".gnuplot";
 *   
 *       std::cout << "Writing solution to <" << filename << ">..." << std::endl
 *                 << std::endl;
 *       std::ofstream gnuplot_output(filename);
 *   
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "u");
 *       data_out.add_data_vector(estimates.block(0), "est");
 *   
 *       data_out.build_patches();
 *   
 *       data_out.write_gnuplot(gnuplot_output);
 *     }
 *   
 * @endcode
 * 
 * And finally the adaptive loop, more or less like in previous examples.
 * 
 * @code
 *     template <int dim>
 *     void InteriorPenaltyProblem<dim>::run(unsigned int n_steps)
 *     {
 *       std::cout << "Element: " << fe.get_name() << std::endl;
 *       for (unsigned int s = 0; s < n_steps; ++s)
 *         {
 *           std::cout << "Step " << s << std::endl;
 *           if (estimates.block(0).empty())
 *             triangulation.refine_global(1);
 *           else
 *             {
 *               GridRefinement::refine_and_coarsen_fixed_fraction(
 *                 triangulation, estimates.block(0), 0.5, 0.0);
 *               triangulation.execute_coarsening_and_refinement();
 *             }
 *   
 *           std::cout << "Triangulation " << triangulation.n_active_cells()
 *                     << " cells, " << triangulation.n_levels() << " levels"
 *                     << std::endl;
 *   
 *           setup_system();
 *           std::cout << "DoFHandler " << dof_handler.n_dofs()
 *                     << " dofs, level dofs";
 *           for (unsigned int l = 0; l < triangulation.n_levels(); ++l)
 *             std::cout << ' ' << dof_handler.n_dofs(l);
 *           std::cout << std::endl;
 *   
 *           std::cout << "Assemble matrix" << std::endl;
 *           assemble_matrix();
 *           std::cout << "Assemble multilevel matrix" << std::endl;
 *           assemble_mg_matrix();
 *           std::cout << "Assemble right hand side" << std::endl;
 *           assemble_right_hand_side();
 *           std::cout << "Solve" << std::endl;
 *           solve();
 *           error();
 *           std::cout << "Estimate " << estimate() << std::endl;
 *           output_results(s);
 *         }
 *     }
 *   } // namespace Step39
 *   
 *   
 *   
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step39;
 *   
 *         InteriorPenaltyProblem<2> test1;
 *         test1.run(12);
 *       }
 *     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="step_39-Results"></a><h1>Results</h1>
 
 
<a name="step_39-Logfileoutput"></a><h3>Logfile output</h3>
 
First, the program produces the usual logfile here stored in <tt>deallog</tt>. It reads (with omission of intermediate steps)
 
@code
Element: FE_DGQ<2>(3)
Step 0
Triangulation 16 cells, 2 levels
DoFHandler 256 dofs, level dofs 64 256
Assemble matrix
Assemble multilevel matrix
Assemble right hand side
Solve
Converged in 13 iterations
energy-error: 0.297419
L2-error:     0.00452447
Estimate 0.990460
Writing solution to <sol-00.gnuplot>...
 
Step 1
Triangulation 25 cells, 3 levels
DoFHandler 400 dofs, level dofs 64 256 192
Assemble matrix
Assemble multilevel matrix
Assemble right hand side
Solve
Converged in 14 iterations
energy-error: 0.258559
L2-error:     0.00288510
Estimate 0.738624
Writing solution to <sol-01.gnuplot>...
 
...
 
Step 10
Triangulation 232 cells, 11 levels
DoFHandler 3712 dofs, level dofs 64 256 896 768 768 640 512 256 256 256 256
Assemble matrix
Assemble multilevel matrix
Assemble right hand side
Solve
Converged in 15 iterations
energy-error: 0.0132475
L2-error:     1.00423e-05
Estimate 0.0470724
Writing solution to <sol-10.gnuplot>...
 
Step 11
Triangulation 322 cells, 12 levels
DoFHandler 5152 dofs, level dofs 64 256 1024 1024 896 768 768 640 448 320 320 320
Assemble matrix
Assemble multilevel matrix
Assemble right hand side
Solve
Converged in 15 iterations
energy-error: 0.00934891
L2-error:     5.41095e-06
Estimate 0.0329102
Writing solution to <sol-11.gnuplot>...
@endcode
 
This log for instance shows that the number of conjugate gradient
iteration steps is constant at approximately 15.
This is the key qualitative result of the example: despite adaptive mesh
refinement and the discontinuous Galerkin discretization, the multigrid
preconditioner keeps the iteration count essentially mesh-independent.
 
 
<a name="step_39-Postprocessingofthelogfile"></a><h3>Postprocessing of the logfile</h3>
 
 
<img src="https://dealii.org/images/steps/developer/step-39-convergence.svg" alt="">
Using the perl script <tt>postprocess.pl</tt>, we extract relevant
data into <tt>output.dat</tt>, which can be used to plot graphs with
<tt>gnuplot</tt>. The graph above for instance was produced using the gnuplot
script <tt>plot_errors.gpl</tt> via
 
@code
./step-39 | perl postprocess.pl >output.dat
gnuplot plot_errors.gpl
@endcode
 
Reference data can be found in <tt>output.reference.dat</tt>.
 *
 *
<a name="step_39-PlainProg"></a>
<h1> The plain program</h1>
@include "step-39.cc"
*/