From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 3 May 2021 23:58:45 +0000 (-0600)
Subject: Fix a potential bug in step-15.
X-Git-Tag: v9.3.0-rc1~146^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2dfa0211a1e4407c0bd0803580f1a03fa72af147;p=dealii.git

Fix a potential bug in step-15.
---

diff --git a/examples/step-15/doc/results.dox b/examples/step-15/doc/results.dox
index bfd4d38c99..f939657c18 100644
--- a/examples/step-15/doc/results.dox
+++ b/examples/step-15/doc/results.dox
@@ -3,7 +3,7 @@
 
 The output of the program looks as follows:
 @code
-Mesh refinement step
+Mesh refinement step 0
   Initial residual: 1.53143
   Residual: 1.08746
   Residual: 0.966748
@@ -11,7 +11,7 @@ Mesh refinement step
   Residual: 0.766462
   Residual: 0.685475
 
-Mesh refinement step
+Mesh refinement step 1
   Initial residual: 0.868959
   Residual: 0.762125
   Residual: 0.677792
@@ -19,7 +19,7 @@ Mesh refinement step
   Residual: 0.542748
   Residual: 0.48704
 
-Mesh refinement step
+Mesh refinement step 2
   Initial residual: 0.426445
   Residual: 0.382731
   Residual: 0.343865
@@ -27,7 +27,7 @@ Mesh refinement step
   Residual: 0.278147
   Residual: 0.250327
 
-Mesh refinement step
+Mesh refinement step 3
   Initial residual: 0.282026
   Residual: 0.253146
   Residual: 0.227414
@@ -35,7 +35,7 @@ Mesh refinement step
   Residual: 0.183803
   Residual: 0.165319
 
-Mesh refinement step
+Mesh refinement step 4
   Initial residual: 0.154404
   Residual: 0.138723
   Residual: 0.124694
diff --git a/examples/step-15/step-15.cc b/examples/step-15/step-15.cc
index afb3f16ee3..766488b9bd 100644
--- a/examples/step-15/step-15.cc
+++ b/examples/step-15/step-15.cc
@@ -404,18 +404,11 @@ namespace Step15
     // values of the old solution, but rather indices, we need to preserve the
     // old solution vector until we have gotten the new interpolated
     // values. Thus, we have the new values written into a temporary vector,
-    // and only afterwards write them into the solution vector object. Once we
-    // have this solution we have to make sure that the $u^n$ we now have
-    // actually has the correct boundary values. As explained at the end of
-    // the introduction, this is not automatically the case even if the
-    // solution before refinement had the correct boundary values, and so we
-    // have to explicitly make sure that it now has:
+    // and only afterwards write them into the solution vector object:
     Vector<double> tmp(dof_handler.n_dofs());
     solution_transfer.interpolate(current_solution, tmp);
     current_solution = tmp;
 
-    set_boundary_values();
-
     // On the new mesh, there are different hanging nodes, which we have to
     // compute again. To ensure there are no hanging nodes of the old mesh in
     // the object, it's first cleared.  To be on the safe side, we then also
@@ -430,6 +423,14 @@ namespace Step15
 
     hanging_node_constraints.distribute(current_solution);
 
+    // Once we have the interpolated solution and all information about
+    // hanging nodes, we have to make sure that the $u^n$ we now have
+    // actually has the correct boundary values. As explained at the end of
+    // the introduction, this is not automatically the case even if the
+    // solution before refinement had the correct boundary values, and so we
+    // have to explicitly make sure that it now has:
+    set_boundary_values();
+
     // We end the function by updating all the remaining data structures,
     // indicating to <code>setup_dofs()</code> that this is not the first
     // go-around and that it needs to preserve the content of the solution
@@ -445,9 +446,15 @@ namespace Step15
   // boundary values for our problem.  Having refined the mesh (or just
   // started computations), there might be new nodal points on the
   // boundary. These have values that are simply interpolated from the
-  // previous mesh (or are just zero), instead of the correct boundary
-  // values. This is fixed up by setting all boundary nodes explicit to the
-  // right value:
+  // previous mesh in `refine_mesh()`, instead of the correct boundary
+  // values. This is fixed up by setting all boundary nodes of the current
+  // solution vector explicit to the right value.
+  //
+  // There is one issue we have to pay attention to, though: If we have
+  // a hanging node right next to a new boundary node, then its value
+  // must also be adjusted to make sure that the finite element field
+  // remains continuous. This is what the call in the last line of this
+  // function does.
   template <int dim>
   void MinimalSurfaceProblem<dim>::set_boundary_values()
   {
@@ -458,6 +465,8 @@ namespace Step15
                                              boundary_values);
     for (auto &boundary_value : boundary_values)
       current_solution(boundary_value.first) = boundary_value.second;
+
+    hanging_node_constraints.distribute(current_solution);
   }