From fec2a8ab3b24c166a7f28924c6fff0f7f2be2e49 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 3 Jan 2021 12:45:32 +0100 Subject: [PATCH] Minor changes --- examples/step-74/doc/intro.dox | 4 +-- examples/step-74/step-74.cc | 48 ++++++++++++++++------------------ 2 files changed, 25 insertions(+), 27 deletions(-) diff --git a/examples/step-74/doc/intro.dox b/examples/step-74/doc/intro.dox index 585a87421b..20b2e24b38 100644 --- a/examples/step-74/doc/intro.dox +++ b/examples/step-74/doc/intro.dox @@ -73,8 +73,8 @@ The discretization using the SIPG is given by the following weak formula @f}

The penalty parameter

-The parameter is defined as $\sigma = \gamma/h_f$, where $h_f$ a local length scale associated -with the cell face; here we choose the length of the cell in the direction normal to the face, +The penalty parameter is defined as $\sigma = \gamma/h_f$, where $h_f$ a local length scale associated +with the cell face; here we choose the approximation of the length of the cell in the direction normal to the face, and $\gamma$ is the penalization constant. To ensure the discrete coercivity, the penalization constant has to be large enough @cite ainsworth2007posteriori. People do not really have consensus on which precise formula to choose, among what was proposed in the literature. diff --git a/examples/step-74/step-74.cc b/examples/step-74/step-74.cc index ccdd6e5f6f..d8beb3cab9 100644 --- a/examples/step-74/step-74.cc +++ b/examples/step-74/step-74.cc @@ -202,7 +202,7 @@ namespace Step74 const double cell_extend_left, const double cell_extend_right) { - const double degree = std::max(1., static_cast(fe_degree)); + const double degree = std::max(1, fe_degree); return degree * (degree + 1.) * 0.5 * (1. / cell_extend_left + 1. / cell_extend_right); } @@ -242,7 +242,7 @@ namespace Step74 // @sect3{The SIPGLaplace class} // After this preparations, we proceed with the main class of this program // called SIPGLaplace. Major differences will only come up in the - // implementation of the assemble functions, since use FEInterfaceValues to + // implementation of the assemble functions, since we use FEInterfaceValues to // assemble face terms. template class SIPGLaplace @@ -262,16 +262,18 @@ namespace Step74 void compute_error_estimate(); double compute_energy_norm(); - Triangulation triangulation; - const unsigned degree; - QGauss quadrature; - QGauss face_quadrature; - const MappingQ1 mapping; + Triangulation triangulation; + const unsigned degree; + const QGauss quadrature; + const QGauss face_quadrature; + const QGauss quadrature_overintegration; + const QGauss face_quadrature_overintegration; + const MappingQ1 mapping; using ScratchData = MeshWorker::ScratchData; - FE_DGQ fe; - DoFHandler dof_handler; + const FE_DGQ fe; + DoFHandler dof_handler; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -303,6 +305,8 @@ namespace Step74 : degree(3) , quadrature(degree + 1) , face_quadrature(degree + 1) + , quadrature_overintegration(degree + 2) + , face_quadrature_overintegration(degree + 2) , mapping() , fe(degree) , dof_handler(triangulation) @@ -616,8 +620,7 @@ namespace Step74 std::vector sol_u(n_q_points); fe_fv.get_function_values(solution, sol_u); - const double extent1 = cell->extent_in_direction( - GeometryInfo::unit_normal_direction[face_no]); + const double extent1 = cell->measure() / cell->face(face_no)->measure(); const double penalty = compute_penalty(degree, extent1, extent1); double difference_norm_square = 0.; @@ -663,10 +666,8 @@ namespace Step74 const double h = cell->face(f)->diameter(); - const double extent1 = - cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); - const double extent2 = ncell->extent_in_direction( - GeometryInfo::unit_normal_direction[nf]); + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); const double penalty = compute_penalty(degree, extent1, extent2); double flux_jump_square = 0; @@ -764,8 +765,7 @@ namespace Step74 std::vector sol_u(n_q_points); fe_fv.get_function_values(solution, sol_u); - const double extent1 = cell->extent_in_direction( - GeometryInfo::unit_normal_direction[face_no]); + const double extent1 = cell->measure() / cell->face(face_no)->measure(); const double penalty = compute_penalty(degree, extent1, extent1); double difference_norm_square = 0.; @@ -802,10 +802,8 @@ namespace Step74 std::vector jump(n_q_points); get_function_jump(fe_iv, solution, jump); - const double extent1 = - cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); - const double extent2 = ncell->extent_in_direction( - GeometryInfo::unit_normal_direction[nf]); + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); const double penalty = compute_penalty(degree, extent1, extent2); double u_jump_square = 0; @@ -832,9 +830,9 @@ namespace Step74 ScratchData scratch_data(mapping, fe, - QGauss(fe.degree + 2), + quadrature_overintegration, cell_flags, - QGauss(fe.degree + 2), + face_quadrature_overintegration, face_flags); CopyData cd; @@ -879,7 +877,7 @@ namespace Step74 solution, *(exact_solution.get()), difference_per_cell, - QGauss(fe.degree + 2), + quadrature_overintegration, VectorTools::L2_norm); L2_error = VectorTools::compute_global_error(triangulation, @@ -894,7 +892,7 @@ namespace Step74 solution, *(exact_solution.get()), difference_per_cell, - QGauss(fe.degree + 2), + quadrature_overintegration, VectorTools::H1_seminorm); H1_error = VectorTools::compute_global_error(triangulation, -- 2.39.5