276 *
for (
unsigned int i = 0; i <
values.size(); ++i)
285 * SmoothSolution<dim>::gradient(
const Point<dim> &point,
286 *
const unsigned int )
const
294 *
return return_value;
301 * The corresponding right-hand side of the smooth function:
305 *
class SmoothRightHandSide :
public Function<dim>
308 * SmoothRightHandSide()
313 * std::vector<double> &values,
314 *
const unsigned int )
const override;
321 * SmoothRightHandSide<dim>::value_list(
const std::vector<
Point<dim>> &points,
322 * std::vector<double> &values,
323 *
const unsigned int )
const
326 *
for (
unsigned int i = 0; i <
values.size(); ++i)
327 * values[i] = 8. * PI * PI *
std::sin(2. * PI * points[i][0]) *
335 * The right-hand side that corresponds to the function
337 * assume that the diffusion coefficient @f$\nu = 1@f$:
341 *
class SingularRightHandSide :
public Function<dim>
344 * SingularRightHandSide()
349 * std::vector<double> &values,
350 *
const unsigned int )
const override;
360 * SingularRightHandSide<dim>::value_list(
const std::vector<
Point<dim>> &points,
361 * std::vector<double> &values,
362 *
const unsigned int )
const
364 *
for (
unsigned int i = 0; i <
values.size(); ++i)
373 * <a name=
"step_74-Auxiliaryfunctions"></a>
375 * This function computes the penalty @f$\sigma@f$.
378 *
double get_penalty_factor(
const unsigned int fe_degree,
379 *
const double cell_extent_left,
380 *
const double cell_extent_right)
382 *
const unsigned int degree =
std::max(1U, fe_degree);
383 *
return degree * (degree + 1.) * 0.5 *
384 * (1. / cell_extent_left + 1. / cell_extent_right);
391 * <a name=
"step_74-TheCopyData"></a>
392 * <h3>The CopyData</h3>
394 * which is essentially the same as @ref step_12
"step-12". Note that the
395 *
"Scratch" object is not defined here because we use
397 * objects is extensively explained in the
WorkStream namespace documentation.
400 *
struct CopyDataFace
403 * std::vector<types::global_dof_index> joint_dof_indices;
404 * std::array<double, 2>
values;
405 * std::array<unsigned int, 2> cell_indices;
414 * std::vector<types::global_dof_index> local_dof_indices;
415 * std::vector<CopyDataFace> face_data;
420 *
template <
class Iterator>
421 *
void reinit(
const Iterator &cell,
const unsigned int dofs_per_cell)
423 *
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
424 * cell_rhs.reinit(dofs_per_cell);
425 * local_dof_indices.resize(dofs_per_cell);
426 * cell->get_dof_indices(local_dof_indices);
435 * <a name=
"step_74-TheSIPGLaplaceclass"></a>
436 * <h3>The SIPGLaplace
class</h3>
437 * After these preparations, we proceed with the main
class of this program,
438 * called `SIPGLaplace`. The overall structure of the
class is as in many
439 * of the other tutorial programs. Major differences will only come up in the
448 * SIPGLaplace(
const TestCase &test_case);
452 *
void setup_system();
453 *
void assemble_system();
455 *
void refine_grid();
456 *
void output_results(
const unsigned int cycle)
const;
458 *
void compute_errors();
459 *
void compute_error_estimate();
460 *
double compute_energy_norm_error();
463 *
const unsigned int degree;
465 *
const QGauss<dim - 1> face_quadrature;
467 *
const QGauss<dim - 1> face_quadrature_overintegration;
482 * The remainder of the
class's members are used for the following:
483 * - Vectors to store error estimator square and energy norm square per
485 * - Print convergence rate and errors on the screen.
486 * - The fiffusion coefficient @f$\nu@f$ is set to 1.
487 * - Members that store information about the test case to be computed.
490 * Vector<double> estimated_error_square_per_cell;
491 * Vector<double> energy_norm_square_per_cell;
493 * ConvergenceTable convergence_table;
495 * const double diffusion_coefficient = 1.;
497 * const TestCase test_case;
498 * std::unique_ptr<const Function<dim>> exact_solution;
499 * std::unique_ptr<const Function<dim>> rhs_function;
504 * The constructor here takes the test case as input and then
505 * determines the correct solution and right-hand side classes. The
506 * remaining member variables are initialized in the obvious way.
510 * SIPGLaplace<dim>::SIPGLaplace(const TestCase &test_case)
512 * , quadrature(degree + 1)
513 * , face_quadrature(degree + 1)
514 * , quadrature_overintegration(degree + 2)
515 * , face_quadrature_overintegration(degree + 2)
518 * , dof_handler(triangulation)
519 * , test_case(test_case)
521 * if (test_case == TestCase::convergence_rate)
523 * exact_solution = std::make_unique<const SmoothSolution<dim>>();
524 * rhs_function = std::make_unique<const SmoothRightHandSide<dim>>();
527 * else if (test_case == TestCase::l_singularity)
530 * std::make_unique<const Functions::LSingularityFunction>();
531 * rhs_function = std::make_unique<const SingularRightHandSide<dim>>();
534 * AssertThrow(false, ExcNotImplemented());
540 * void SIPGLaplace<dim>::setup_system()
542 * dof_handler.distribute_dofs(fe);
543 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
544 * DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
545 * sparsity_pattern.copy_from(dsp);
547 * system_matrix.reinit(sparsity_pattern);
548 * solution.reinit(dof_handler.n_dofs());
549 * system_rhs.reinit(dof_handler.n_dofs());
557 * <a name="step_74-Theassemble_systemfunction"></a>
558 * <h3>The assemble_system function</h3>
559 * The assemble function here is similar to that in @ref step_12 "step-12" and @ref step_47 "step-47".
560 * Different from assembling by hand, we just need to focus
561 * on assembling on each cell, each boundary face, and each
562 * interior face. The loops over cells and faces are handled
563 * automatically by MeshWorker::mesh_loop().
567 * The function starts by defining a local (lambda) function that is
568 * used to integrate the cell terms:
572 * void SIPGLaplace<dim>::assemble_system()
574 * const auto cell_worker =
575 * [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
576 * ScratchData &scratch_data,
577 * CopyData ©_data) {
578 * const FEValues<dim> &fe_v = scratch_data.reinit(cell);
579 * const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
580 * copy_data.reinit(cell, dofs_per_cell);
582 * const std::vector<Point<dim>> &q_points =
583 * scratch_data.get_quadrature_points();
584 * const unsigned int n_q_points = q_points.size();
585 * const std::vector<double> &JxW = scratch_data.get_JxW_values();
587 * std::vector<double> rhs(n_q_points);
588 * rhs_function->value_list(q_points, rhs);
590 * for (unsigned int point = 0; point < n_q_points; ++point)
591 * for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
593 * for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
594 * copy_data.cell_matrix(i, j) +=
595 * diffusion_coefficient * // nu
596 * fe_v.shape_grad(i, point) * // grad v_h
597 * fe_v.shape_grad(j, point) * // grad u_h
600 * copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
608 * Next, we need a function that assembles face integrals on the boundary:
611 * const auto boundary_worker =
612 * [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
613 * const unsigned int &face_no,
614 * ScratchData &scratch_data,
615 * CopyData ©_data) {
616 * const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
618 * const std::vector<Point<dim>> &q_points =
619 * scratch_data.get_quadrature_points();
620 * const unsigned int n_q_points = q_points.size();
621 * const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
623 * const std::vector<double> &JxW = scratch_data.get_JxW_values();
624 * const std::vector<Tensor<1, dim>> &normals =
625 * scratch_data.get_normal_vectors();
627 * std::vector<double> g(n_q_points);
628 * exact_solution->value_list(q_points, g);
630 * const double extent1 = cell->measure() / cell->face(face_no)->measure();
631 * const double penalty = get_penalty_factor(degree, extent1, extent1);
633 * for (unsigned int point = 0; point < n_q_points; ++point)
635 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
636 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
637 * copy_data.cell_matrix(i, j) +=
638 * (-diffusion_coefficient * // - nu
639 * fe_fv.shape_value(i, point) * // v_h
640 * (fe_fv.shape_grad(j, point) * // (grad u_h .
641 * normals[point]) // n)
643 * - diffusion_coefficient * // - nu
644 * (fe_fv.shape_grad(i, point) * // (grad v_h .
645 * normals[point]) * // n)
646 * fe_fv.shape_value(j, point) // u_h
648 * + diffusion_coefficient * penalty * // + nu sigma
649 * fe_fv.shape_value(i, point) * // v_h
650 * fe_fv.shape_value(j, point) // u_h
655 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
656 * copy_data.cell_rhs(i) +=
657 * (-diffusion_coefficient * // - nu
658 * (fe_fv.shape_grad(i, point) * // (grad v_h .
659 * normals[point]) * // n)
663 * + diffusion_coefficient * penalty * // + nu sigma
664 * fe_fv.shape_value(i, point) * g[point] // v_h g
673 * Finally, a function that assembles face integrals on interior
674 * faces. To reinitialize FEInterfaceValues, we need to pass
675 * cells, face and subface indices (for adaptive refinement) to
676 * the reinit() function of FEInterfaceValues:
679 * const auto face_worker =
680 * [&](const typename DoFHandler<dim>::cell_iterator &cell,
681 * const unsigned int &f,
682 * const unsigned int &sf,
683 * const typename DoFHandler<dim>::cell_iterator &ncell,
684 * const unsigned int &nf,
685 * const unsigned int &nsf,
686 * ScratchData &scratch_data,
687 * CopyData ©_data) {
688 * const FEInterfaceValues<dim> &fe_iv =
689 * scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
691 * copy_data.face_data.emplace_back();
692 * CopyDataFace ©_data_face = copy_data.face_data.back();
693 * const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs();
694 * copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
695 * copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
697 * const std::vector<double> &JxW = fe_iv.get_JxW_values();
698 * const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
700 * const double extent1 = cell->measure() / cell->face(f)->measure();
701 * const double extent2 = ncell->measure() / ncell->face(nf)->measure();
702 * const double penalty = get_penalty_factor(degree, extent1, extent2);
704 * for (const unsigned int point : fe_iv.quadrature_point_indices())
706 * for (const unsigned int i : fe_iv.dof_indices())
707 * for (const unsigned int j : fe_iv.dof_indices())
708 * copy_data_face.cell_matrix(i, j) +=
709 * (-diffusion_coefficient * // - nu
710 * fe_iv.jump_in_shape_values(i, point) * // [v_h]
711 * (fe_iv.average_of_shape_gradients(j,
712 * point) * // ({grad u_h} .
713 * normals[point]) // n)
715 * - diffusion_coefficient * // - nu
716 * (fe_iv.average_of_shape_gradients(i,
717 * point) * // (grad v_h .
718 * normals[point]) * // n)
719 * fe_iv.jump_in_shape_values(j, point) // [u_h]
721 * + diffusion_coefficient * penalty * // + nu sigma
722 * fe_iv.jump_in_shape_values(i, point) * // [v_h]
723 * fe_iv.jump_in_shape_values(j, point) // [u_h]
732 * The following lambda function will then copy data into the
733 * global matrix and right-hand side. Though there are no hanging
734 * node constraints in DG discretization, we define an empty
735 * AffineConstraints object that allows us to use the
736 * AffineConstraints::distribute_local_to_global() functionality.
739 * AffineConstraints<double> constraints;
740 * constraints.close();
741 * const auto copier = [&](const CopyData &c) {
742 * constraints.distribute_local_to_global(c.cell_matrix,
744 * c.local_dof_indices,
750 * Copy data from interior face assembly to the global matrix.
753 * for (const CopyDataFace &cdf : c.face_data)
755 * constraints.distribute_local_to_global(cdf.cell_matrix,
756 * cdf.joint_dof_indices,
764 * With the assembly functions defined, we can now create
765 * ScratchData and CopyData objects, and pass them together with
766 * the lambda functions above to MeshWorker::mesh_loop(). In
767 * addition, we need to specify that we want to assemble on
768 * interior faces exactly once.
771 * const UpdateFlags cell_flags = update_values | update_gradients |
772 * update_quadrature_points | update_JxW_values;
773 * const UpdateFlags face_flags = update_values | update_gradients |
774 * update_quadrature_points |
775 * update_normal_vectors | update_JxW_values;
777 * ScratchData scratch_data(
778 * mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
779 * CopyData copy_data;
781 * MeshWorker::mesh_loop(dof_handler.begin_active(),
787 * MeshWorker::assemble_own_cells |
788 * MeshWorker::assemble_boundary_faces |
789 * MeshWorker::assemble_own_interior_faces_once,
799 * <a name="step_74-Thesolveandoutput_resultsfunction"></a>
800 * <h3>The solve() and output_results() function</h3>
801 * The following two functions are entirely standard and without difficulty.
805 * void SIPGLaplace<dim>::solve()
807 * SparseDirectUMFPACK A_direct;
808 * A_direct.initialize(system_matrix);
809 * A_direct.vmult(solution, system_rhs);
815 * void SIPGLaplace<dim>::output_results(const unsigned int cycle) const
817 * const std::string filename = "sol_Q" + Utilities::int_to_string(degree, 1) +
818 * "-" + Utilities::int_to_string(cycle, 2) +
820 * std::ofstream output(filename);
822 * DataOut<dim> data_out;
823 * data_out.attach_dof_handler(dof_handler);
824 * data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
825 * data_out.build_patches(mapping);
826 * data_out.write_vtu(output);
833 * <a name="step_74-Thecompute_error_estimatefunction"></a>
834 * <h3>The compute_error_estimate() function</h3>
835 * The assembly of the error estimator here is quite similar to
836 * that of the global matrix and right-had side and can be handled
837 * by the MeshWorker::mesh_loop() framework. To understand what
838 * each of the local (lambda) functions is doing, recall first that
839 * the local cell residual is defined as
840 * @f$h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2@f$:
844 * void SIPGLaplace<dim>::compute_error_estimate()
846 * const auto cell_worker =
847 * [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
848 * ScratchData &scratch_data,
849 * CopyData ©_data) {
850 * const FEValues<dim> &fe_v = scratch_data.reinit(cell);
852 * copy_data.cell_index = cell->active_cell_index();
854 * const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
855 * const unsigned int n_q_points = q_points.size();
856 * const std::vector<double> &JxW = fe_v.get_JxW_values();
858 * std::vector<Tensor<2, dim>> hessians(n_q_points);
859 * fe_v.get_function_hessians(solution, hessians);
861 * std::vector<double> rhs(n_q_points);
862 * rhs_function->value_list(q_points, rhs);
864 * const double hk = cell->diameter();
865 * double residual_norm_square = 0;
867 * for (unsigned int point = 0; point < n_q_points; ++point)
869 * const double residual =
870 * rhs[point] + diffusion_coefficient * trace(hessians[point]);
871 * residual_norm_square += residual * residual * JxW[point];
873 * copy_data.value = hk * hk * residual_norm_square;
878 * Next compute boundary terms @f$\sum_{f\in \partial K \cap \partial \Omega}
879 * \sigma \left\| [ u_h-g_D ] \right\|_f^2 @f$:
882 * const auto boundary_worker =
883 * [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
884 * const unsigned int &face_no,
885 * ScratchData &scratch_data,
886 * CopyData ©_data) {
887 * const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
889 * const std::vector<Point<dim>> &q_points = fe_fv.get_quadrature_points();
890 * const unsigned n_q_points = q_points.size();
892 * const std::vector<double> &JxW = fe_fv.get_JxW_values();
894 * std::vector<double> g(n_q_points);
895 * exact_solution->value_list(q_points, g);
897 * std::vector<double> sol_u(n_q_points);
898 * fe_fv.get_function_values(solution, sol_u);
900 * const double extent1 = cell->measure() / cell->face(face_no)->measure();
901 * const double penalty = get_penalty_factor(degree, extent1, extent1);
903 * double difference_norm_square = 0.;
904 * for (unsigned int point = 0; point < q_points.size(); ++point)
906 * const double diff = (g[point] - sol_u[point]);
907 * difference_norm_square += diff * diff * JxW[point];
909 * copy_data.value += penalty * difference_norm_square;
914 * And finally interior face terms @f$\sum_{f\in \partial K}\lbrace \sigma
915 * \left\| [u_h] \right\|_f^2 + h_f \left\| [\nu \nabla u_h \cdot
916 * \mathbf n ] \right\|_f^2 \rbrace@f$:
919 * const auto face_worker =
920 * [&](const typename DoFHandler<dim>::cell_iterator &cell,
921 * const unsigned int &f,
922 * const unsigned int &sf,
923 * const typename DoFHandler<dim>::cell_iterator &ncell,
924 * const unsigned int &nf,
925 * const unsigned int &nsf,
926 * ScratchData &scratch_data,
927 * CopyData ©_data) {
928 * const FEInterfaceValues<dim> &fe_iv =
929 * scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
931 * copy_data.face_data.emplace_back();
932 * CopyDataFace ©_data_face = copy_data.face_data.back();
934 * copy_data_face.cell_indices[0] = cell->active_cell_index();
935 * copy_data_face.cell_indices[1] = ncell->active_cell_index();
937 * const std::vector<double> &JxW = fe_iv.get_JxW_values();
938 * const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
940 * const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
941 * const unsigned int n_q_points = q_points.size();
943 * std::vector<double> jump(n_q_points);
944 * fe_iv.get_jump_in_function_values(solution, jump);
946 * std::vector<Tensor<1, dim>> grad_jump(n_q_points);
947 * fe_iv.get_jump_in_function_gradients(solution, grad_jump);
949 * const double h = cell->face(f)->diameter();
951 * const double extent1 = cell->measure() / cell->face(f)->measure();
952 * const double extent2 = ncell->measure() / ncell->face(nf)->measure();
953 * const double penalty = get_penalty_factor(degree, extent1, extent2);
955 * double flux_jump_square = 0;
956 * double u_jump_square = 0;
957 * for (unsigned int point = 0; point < n_q_points; ++point)
959 * u_jump_square += jump[point] * jump[point] * JxW[point];
960 * const double flux_jump = grad_jump[point] * normals[point];
961 * flux_jump_square +=
962 * diffusion_coefficient * flux_jump * flux_jump * JxW[point];
964 * copy_data_face.values[0] =
965 * 0.5 * h * (flux_jump_square + penalty * u_jump_square);
966 * copy_data_face.values[1] = copy_data_face.values[0];
971 * Having computed local contributions for each cell, we still
972 * need a way to copy these into the global vector that will hold
973 * the error estimators for all cells:
976 * const auto copier = [&](const CopyData ©_data) {
977 * if (copy_data.cell_index != numbers::invalid_unsigned_int)
978 * estimated_error_square_per_cell[copy_data.cell_index] +=
980 * for (const CopyDataFace &cdf : copy_data.face_data)
981 * for (unsigned int j = 0; j < 2; ++j)
982 * estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
987 * After all of this set-up, let's
do the actual work: We resize
988 * the vector into which the results will be written, and then
1001 * ScratchData scratch_data(
1002 * mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
1004 * CopyData copy_data;
1006 * dof_handler.end(),
1021 * <a name=
"step_74-Thecompute_energy_norm_errorfunction"></a>
1022 * <h3>The compute_energy_norm_error() function</h3>
1023 * Next, we evaluate the accuracy in terms of the energy norm.
1024 * This function is similar to the assembling of the error estimator above.
1025 * Here we compute the square of the energy norm defined by
1027 * \|u \|_{1,h}^2 = \sum_{K \in \Gamma_h} \nu\|\nabla u \|_K^2 +
1028 * \sum_{f \in F_i} \sigma \| [ u ] \|_f^2 +
1029 * \sum_{f \in F_b} \sigma \|u\|_f^2.
1031 * Therefore the corresponding error is
1033 * \|u -u_h \|_{1,h}^2 = \sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2
1034 * + \sum_{f \in F_i} \sigma \|[ u_h ] \|_f^2 + \sum_{f \in F_b}\sigma
1039 *
template <
int dim>
1040 *
double SIPGLaplace<dim>::compute_energy_norm_error()
1042 * energy_norm_square_per_cell.reinit(
triangulation.n_active_cells());
1046 * Assemble @f$\sum_{K \in \Gamma_h} \nu\|\nabla (u_h - u) \|_K^2 @f$.
1049 *
const auto cell_worker =
1051 * ScratchData &scratch_data,
1052 * CopyData ©_data) {
1055 * copy_data.cell_index = cell->active_cell_index();
1057 *
const std::vector<Point<dim>> &q_points = fe_v.get_quadrature_points();
1058 *
const unsigned int n_q_points = q_points.size();
1059 *
const std::vector<double> &JxW = fe_v.get_JxW_values();
1061 * std::vector<Tensor<1, dim>> grad_u(n_q_points);
1062 * fe_v.get_function_gradients(solution, grad_u);
1064 * std::vector<Tensor<1, dim>> grad_exact(n_q_points);
1065 * exact_solution->gradient_list(q_points, grad_exact);
1067 *
double norm_square = 0;
1068 *
for (
unsigned int point = 0;
point < n_q_points; ++
point)
1071 * (grad_u[point] - grad_exact[point]).norm_square() * JxW[
point];
1073 * copy_data.value = diffusion_coefficient * norm_square;
1078 * Assemble @f$\sum_{f \in F_b}\sigma \|u_h-g_D\|_f^2@f$.
1081 *
const auto boundary_worker =
1083 *
const unsigned int &face_no,
1084 * ScratchData &scratch_data,
1085 * CopyData ©_data) {
1089 *
const unsigned n_q_points = q_points.size();
1091 *
const std::vector<double> &JxW = fe_fv.get_JxW_values();
1093 * std::vector<double> g(n_q_points);
1094 * exact_solution->value_list(q_points, g);
1096 * std::vector<double> sol_u(n_q_points);
1097 * fe_fv.get_function_values(solution, sol_u);
1099 *
const double extent1 = cell->measure() / cell->face(face_no)->measure();
1100 *
const double penalty = get_penalty_factor(degree, extent1, extent1);
1102 *
double difference_norm_square = 0.;
1103 *
for (
unsigned int point = 0;
point < q_points.size(); ++
point)
1105 *
const double diff = (g[
point] - sol_u[
point]);
1106 * difference_norm_square += diff * diff * JxW[
point];
1108 * copy_data.value += penalty * difference_norm_square;
1113 * Assemble @f$\sum_{f \in F_i} \sigma \| [ u_h ] \|_f^2@f$.
1116 *
const auto face_worker =
1118 *
const unsigned int &f,
1119 *
const unsigned int &sf,
1121 *
const unsigned int &nf,
1122 *
const unsigned int &nsf,
1123 * ScratchData &scratch_data,
1124 * CopyData ©_data) {
1126 * scratch_data.
reinit(cell, f, sf, ncell, nf, nsf);
1128 * copy_data.face_data.emplace_back();
1129 * CopyDataFace ©_data_face = copy_data.face_data.back();
1131 * copy_data_face.cell_indices[0] = cell->active_cell_index();
1132 * copy_data_face.cell_indices[1] = ncell->active_cell_index();
1134 *
const std::vector<double> &JxW = fe_iv.get_JxW_values();
1136 *
const std::vector<Point<dim>> &q_points = fe_iv.get_quadrature_points();
1137 *
const unsigned int n_q_points = q_points.size();
1139 * std::vector<double> jump(n_q_points);
1140 * fe_iv.get_jump_in_function_values(solution, jump);
1142 *
const double extent1 = cell->measure() / cell->face(f)->measure();
1143 *
const double extent2 = ncell->measure() / ncell->face(nf)->measure();
1144 *
const double penalty = get_penalty_factor(degree, extent1, extent2);
1146 *
double u_jump_square = 0;
1147 *
for (
unsigned int point = 0;
point < n_q_points; ++
point)
1151 * copy_data_face.values[0] = 0.5 * penalty * u_jump_square;
1152 * copy_data_face.values[1] = copy_data_face.values[0];
1155 *
const auto copier = [&](
const CopyData ©_data) {
1157 * energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
1158 *
for (
const CopyDataFace &cdf : copy_data.face_data)
1159 * for (unsigned
int j = 0; j < 2; ++j)
1160 * energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
1168 *
const ScratchData scratch_data(mapping,
1170 * quadrature_overintegration,
1172 * face_quadrature_overintegration,
1175 * CopyData copy_data;
1177 * dof_handler.end(),
1187 *
const double energy_error =
1188 *
std::sqrt(energy_norm_square_per_cell.l1_norm());
1189 *
return energy_error;
1197 * <a name=
"step_74-Therefine_gridfunction"></a>
1198 * <h3>The refine_grid() function</h3>
1201 * template <
int dim>
1202 *
void SIPGLaplace<dim>::refine_grid()
1204 *
const double refinement_fraction = 0.1;
1207 *
triangulation, estimated_error_square_per_cell, refinement_fraction, 0.);
1217 * <a name=
"step_74-Thecompute_errorsfunction"></a>
1218 * <h3>The compute_errors() function</h3>
1219 * We compute three errors in the @f$L_2@f$ norm, @f$H_1@f$ seminorm, and
1220 * the energy norm, respectively. These are then printed to screen,
1221 * but also stored in a table that records how these errors decay
1222 * with mesh refinement and which can be output in one step at the
1223 * end of the program.
1226 * template <
int dim>
1227 *
void SIPGLaplace<dim>::compute_errors()
1229 *
double L2_error, H1_error, energy_error;
1236 * *(exact_solution.get()),
1237 * difference_per_cell,
1238 * quadrature_overintegration,
1242 * difference_per_cell,
1244 * convergence_table.add_value(
"L2", L2_error);
1252 * *(exact_solution.get()),
1253 * difference_per_cell,
1254 * quadrature_overintegration,
1258 * difference_per_cell,
1260 * convergence_table.add_value(
"H1", H1_error);
1264 * energy_error = compute_energy_norm_error();
1265 * convergence_table.add_value(
"Energy", energy_error);
1268 * std::cout <<
" Error in the L2 norm : " << L2_error << std::endl
1269 * <<
" Error in the H1 seminorm : " << H1_error << std::endl
1270 * <<
" Error in the energy norm : " << energy_error
1279 * <a name=
"step_74-Therunfunction"></a>
1280 * <h3>The
run() function</h3>
1283 * template <
int dim>
1284 *
void SIPGLaplace<dim>::run()
1286 *
const unsigned int max_cycle =
1287 * (test_case == TestCase::convergence_rate ? 6 : 20);
1288 *
for (
unsigned int cycle = 0; cycle < max_cycle; ++cycle)
1290 * std::cout <<
"Cycle " << cycle << std::endl;
1292 *
switch (test_case)
1294 *
case TestCase::convergence_rate:
1309 *
case TestCase::l_singularity:
1329 * std::cout <<
" Number of active cells : "
1333 * std::cout <<
" Number of degrees of freedom : " << dof_handler.n_dofs()
1336 * assemble_system();
1338 * output_results(cycle);
1340 * convergence_table.add_value(
"cycle", cycle);
1341 * convergence_table.add_value(
"cells",
triangulation.n_active_cells());
1342 * convergence_table.add_value(
"dofs", dof_handler.n_dofs());
1346 *
if (test_case == TestCase::l_singularity)
1348 * compute_error_estimate();
1349 * std::cout <<
" Estimated error : "
1350 * <<
std::sqrt(estimated_error_square_per_cell.l1_norm())
1353 * convergence_table.add_value(
1355 *
std::sqrt(estimated_error_square_per_cell.l1_norm()));
1357 * std::cout << std::endl;
1362 * Having
run all of our computations, let us tell the convergence
1363 * table how to format its data and output it to screen:
1366 * convergence_table.set_precision(
"L2", 3);
1367 * convergence_table.set_precision(
"H1", 3);
1368 * convergence_table.set_precision(
"Energy", 3);
1370 * convergence_table.set_scientific(
"L2",
true);
1371 * convergence_table.set_scientific(
"H1",
true);
1372 * convergence_table.set_scientific(
"Energy",
true);
1374 *
if (test_case == TestCase::convergence_rate)
1376 * convergence_table.evaluate_convergence_rates(
1378 * convergence_table.evaluate_convergence_rates(
1381 *
if (test_case == TestCase::l_singularity)
1383 * convergence_table.set_precision(
"Estimator", 3);
1384 * convergence_table.set_scientific(
"Estimator",
true);
1387 * std::cout <<
"degree = " << degree << std::endl;
1388 * convergence_table.write_text(
1398 * <a name=
"step_74-Themainfunction"></a>
1399 * <h3>The main() function</h3>
1400 * The following <code>main</code> function is similar to previous examples as
1401 * well, and need not be commented on.
1408 *
using namespace dealii;
1409 *
using namespace Step74;
1411 *
const TestCase test_case = TestCase::l_singularity;
1413 * SIPGLaplace<2> problem(test_case);
1416 *
catch (std::exception &exc)
1418 * std::cerr << std::endl
1420 * <<
"----------------------------------------------------"
1422 * std::cerr <<
"Exception on processing: " << std::endl
1423 * << exc.what() << std::endl
1424 * <<
"Aborting!" << std::endl
1425 * <<
"----------------------------------------------------"
1431 * std::cerr << std::endl
1433 * <<
"----------------------------------------------------"
1435 * std::cerr <<
"Unknown exception!" << std::endl
1436 * <<
"Aborting!" << std::endl
1437 * <<
"----------------------------------------------------"
1445<a name=
"step_74-Results"></a><h1>Results</h1>
1448The output of
this program consist of the console output and
1449solutions in
vtu format.
1451In the
first test
case, when you
run the program, the screen output should look like the following:
1454 Number of active cells : 16
1455 Number of degrees of freedom : 256
1456 Error in the
L2 norm : 0.00193285
1457 Error in the H1 seminorm : 0.106087
1458 Error in the energy
norm : 0.150625
1461 Number of active cells : 64
1462 Number of degrees of freedom : 1024
1463 Error in the
L2 norm : 9.60497e-05
1464 Error in the H1 seminorm : 0.0089954
1465 Error in the energy
norm : 0.0113265
1473When
using the smooth
case with polynomial degree 3, the convergence
1474table will look like
this:
1475<table align=
"center" class=
"doxtable">
1487 <td align=
"center">0</td>
1488 <td align=
"right">16</td>
1489 <td align=
"right">256</td>
1490 <td align=
"center">1.933e-03</td>
1492 <td align=
"center">1.061e-01</td>
1494 <td align=
"center">1.506e-01</td>
1497 <td align=
"center">1</td>
1498 <td align=
"right">64</td>
1499 <td align=
"right">1024</td>
1500 <td align=
"center">9.605e-05</td>
1501 <td align=
"center">4.33</td>
1502 <td align=
"center">8.995e-03</td>
1503 <td align=
"center">3.56</td>
1504 <td align=
"center">1.133e-02</td>
1507 <td align=
"center">2</td>
1508 <td align=
"right">256</td>
1509 <td align=
"right">4096</td>
1510 <td align=
"center">5.606e-06</td>
1511 <td align=
"center">4.10</td>
1512 <td align=
"center">9.018e-04</td>
1513 <td align=
"center">3.32</td>
1514 <td align=
"center">9.736e-04</td>
1517 <td align=
"center">3</td>
1518 <td align=
"right">1024</td>
1519 <td align=
"right">16384</td>
1520 <td align=
"center">3.484e-07</td>
1521 <td align=
"center">4.01</td>
1522 <td align=
"center">1.071e-04</td>
1523 <td align=
"center">3.07</td>
1524 <td align=
"center">1.088e-04</td>
1527 <td align=
"center">4</td>
1528 <td align=
"right">4096</td>
1529 <td align=
"right">65536</td>
1530 <td align=
"center">2.179e-08</td>
1531 <td align=
"center">4.00</td>
1532 <td align=
"center">1.327e-05</td>
1533 <td align=
"center">3.01</td>
1534 <td align=
"center">1.331e-05</td>
1537 <td align=
"center">5</td>
1538 <td align=
"right">16384</td>
1539 <td align=
"right">262144</td>
1540 <td align=
"center">1.363e-09</td>
1541 <td align=
"center">4.00</td>
1542 <td align=
"center">1.656e-06</td>
1543 <td align=
"center">3.00</td>
1544 <td align=
"center">1.657e-06</td>
1548Theoretically,
for polynomial degree @f$p@f$, the order of convergence in @f$L_2@f$
1549norm and @f$H^1@f$ seminorm should be @f$p+1@f$ and @f$p@f$, respectively. Our numerical
1550results are in good agreement with theory.
1552In the
second test
case, when you
run the program, the screen output should look like the following:
1555 Number of active cells : 192
1556 Number of degrees of freedom : 3072
1557 Error in the
L2 norm : 0.000323585
1558 Error in the H1 seminorm : 0.0296202
1559 Error in the energy
norm : 0.0420478
1560 Estimated error : 0.136067
1563 Number of active cells : 249
1564 Number of degrees of freedom : 3984
1565 Error in the
L2 norm : 0.000114739
1566 Error in the H1 seminorm : 0.0186571
1567 Error in the energy
norm : 0.0264879
1568 Estimated error : 0.0857186
1576The following figure provides a
log-
log plot of the errors versus
1577the number of degrees of freedom
for this test
case on the L-shaped
1578domain. In order to interpret it, let @f$n@f$ be the number of degrees of
1579freedom, then on uniformly refined meshes, @f$h@f$ is of order
1580@f$1/\
sqrt{n}@f$ in 2D. Combining the theoretical results in the previous
case,
1581we see that
if the solution is sufficiently smooth,
1582we can expect the error in the @f$L_2@f$
norm to be of order @f$O(n^{-\frac{p+1}{2}})@f$
1583and in @f$H^1@f$ seminorm to be @f$O(n^{-\frac{p}{2}})@f$. It is not a priori
1584clear that one would get the same kind of behavior as a function of
1585@f$n@f$ on adaptively refined meshes like the ones we use
for this second
1586test
case, but one can certainly hope. Indeed, from the figure, we see
1587that the SIPG with adaptive mesh refinement produces asymptotically
1588the kinds of hoped-
for results:
1590<img width=
"600px" src=
"https://www.dealii.org/images/steps/developer/step-74.log-log-plot.png" alt=
"">
1592In addition, we observe that the error estimator decreases
1593at almost the same rate as the errors in the energy
norm and @f$H^1@f$ seminorm,
1594and one order lower than the @f$L_2@f$ error. This suggests
1595its ability to predict regions with large errors.
1597While this tutorial is focused on the implementation, the @ref step_59
"step-59" tutorial program achieves an efficient
1598large-
scale solver in terms of computing time with
matrix-free solution techniques.
1599Note that the @ref step_59
"step-59" tutorial does not work with meshes containing hanging nodes at this moment,
1600because the multigrid interface matrices are not as easily determined,
1601but that is merely the lack of some interfaces in deal.II,
nothing fundamental.
1604<a name=
"step_74-PlainProg"></a>
1605<h1> The plain program</h1>
1606@include
"step-74.cc"
void reinit(const CellIteratorType &cell, const unsigned int face_no, const unsigned int sub_face_no, const CellNeighborIteratorType &cell_neighbor, const unsigned int face_no_neighbor, const unsigned int sub_face_no_neighbor, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int, const unsigned int fe_index_neighbor=numbers::invalid_unsigned_int)
const std::vector< Point< spacedim > > & get_quadrature_points() const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
virtual double laplacian(const Point< 2 > &p, const unsigned int component=0) const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
@ assemble_boundary_faces
@ assemble_own_interior_faces_once
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static constexpr double PI
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation