From 49db35d86bb39249c59db7540669e01f84007178 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 13 May 2024 15:02:22 -0600 Subject: [PATCH] Fix a bug in the computation of the IP penalty. --- examples/step-39/step-39.cc | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/examples/step-39/step-39.cc b/examples/step-39/step-39.cc index 37f663bd1f..ddb40b553b 100644 --- a/examples/step-39/step-39.cc +++ b/examples/step-39/step-39.cc @@ -114,10 +114,12 @@ namespace Step39 // in two steps: First, we compute on each cell // $K_i$ the value $P_i = p_i(p_i+1)/h_i$, where // $p_i$ is the polynomial degree on cell $K_i$ and $h_i$ is the length of - // $K_i$ orthogonal to the current face. Second, if one of the two + // $K_i$ orthogonal to the current face. Second, if exactly one of the two // cells adjacent to the face has children, its penalty is multiplied // by two (to account for the fact that the mesh size $h_i$ there is - // only half that previously computed). Finally, we return the average + // only half that previously computed); it is possible that both adjacent + // cells are refined, in which case we are integrating over a non-active + // face and no adjustment is necessary. Finally, we return the average // of the two penalty values. template double ip_penalty_factor(const MeshWorker::DoFInfo &dinfo1, @@ -134,12 +136,11 @@ namespace Step39 double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1); double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2); - if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children()) - { - Assert(dinfo1.face == dinfo2.face, ExcInternalError()); - Assert(dinfo1.face->has_children(), ExcInternalError()); - penalty1 *= 2; - } + if (dinfo1.cell->has_children() && !dinfo2.cell->has_children()) + penalty1 *= 2; + else if (!dinfo1.cell->has_children() && dinfo2.cell->has_children()) + penalty2 *= 2; + const double penalty = 0.5 * (penalty1 + penalty2); return penalty; } -- 2.39.5