From: Wolfgang Bangerth Date: Thu, 22 Jun 2023 16:47:33 +0000 (-0600) Subject: Fix step-33. X-Git-Tag: v9.5.0-rc1~53^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15437%2Fhead;p=dealii.git Fix step-33. --- diff --git a/examples/step-33/doc/intro.dox b/examples/step-33/doc/intro.dox index 58ee54687b..4f5b09ccec 100644 --- a/examples/step-33/doc/intro.dox +++ b/examples/step-33/doc/intro.dox @@ -304,7 +304,7 @@ the continuous and DG methodologies. It also simplifies the generation of the Jacobian because we do not have to track constrained degrees of freedom through the automatic differentiation used to compute it. -@note Whereas this program was written in 2008, we were unaware of any +@note When this program was written in 2008, we were unaware of any publication that would actually have used this approach. However, a more recent paper by A. Dedner, R. Klöfkorn, and M. Kränkel ("Continuous Finite-Elements on Non-Conforming Grids Using @@ -312,6 +312,26 @@ Discontinuous Galerkin Stabilization", Proceedings of Finite Volumes for Complex Applications VII - Methods and Theoretical Aspects, Springer, 2014) comes close. +@note In hindsight, dealing with hanging nodes in this way is perhaps + not the best choice. deal.II contains many places that assume that + an element that has degrees of freedom on vertices, edges, and faces + (like the FE_Q element used here) represents finite element functions + on a triangulation that have certain continuity properties. Calling + DoFTools::make_hanging_node_constraints() and putting the corresponding + output into an AffineConstraints object allows one to enforce + these assumptions when one calls AffineConstraints::distribute() + at the end of the function that solves linear systems -- this call + can be found in all other tutorial programs that use meshes with + hanging nodes, and ensures that the resulting function's value at + hanging nodes is correct with regard to the values of the nodes on + the adjacent parent cell. On the other hand, this program does not do + this, and this runs afoul of an assertion that was later added to the + library and that tested that the function is indeed continuous. + To address the resulting failure due to the assertion not being + satisfied, we *do* after all enforce continuity at hanging nodes + just before we transfer the solution from one mesh to another, + in the `ConservationLaw::refine_grid()` function. + Further, we enforce a maximum number of refinement levels to keep refinement under check. It is the author's experience that for adaptivity for a time dependent problem, refinement can easily lead the simulation to a screeching halt, because of time step restrictions if the mesh diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 0e5b75d23f..6401e674eb 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -2258,6 +2258,29 @@ namespace Step33 cell->set_coarsen_flag(); } + // The next step addresses a problem mentioned in a remark in the + // introduction: The SolutionTransfer class we want to use later on tests + // the assumption that the solution function is continuous at hanging + // nodes. This is not actually the case in this program because we chose + // (perhaps unwisely) to enforce hanging node constraints in a weak way, + // as one would for example do with discontinuous elements. But the elements + // we use here are continuous (namely, multiple copies of FE_Q), and so + // the assertion would fail and the program abort. To avoid the issue + // (without having to rewrite the whole program), we simply ensure that the + // solution *does* satisfy the hanging node constraints, but creating + // an AffineConstraint object that contains the hanging node constraints and + // applying the constraints to the two solution vectors we want the + // SolutionTransfer class to transfer to the next mesh: + { + AffineConstraints hanging_node_constraints; + DoFTools::make_hanging_node_constraints(dof_handler, + hanging_node_constraints); + hanging_node_constraints.close(); + + hanging_node_constraints.distribute(old_solution); + hanging_node_constraints.distribute(predictor); + } + // Then we need to transfer the various solution vectors from the old to // the new grid while we do the refinement. The SolutionTransfer class is // our friend here; it has a fairly extensive documentation, including