From a90acfabcd585c9aefbebf5f85ae8b3debe870c2 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Fri, 14 Sep 2018 16:46:55 -0600
Subject: [PATCH] Deal better with boundary values in step-52.

---
 examples/step-52/step-52.cc | 41 +++++++++++++++++++++++++++++--------
 1 file changed, 32 insertions(+), 9 deletions(-)

diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc
index 77ef88d5e1..a9f3e8d1c0 100644
--- a/examples/step-52/step-52.cc
+++ b/examples/step-52/step-52.cc
@@ -466,12 +466,23 @@ namespace Step52
 
   // @sect4{<code>Diffusion::explicit_method</code>}
   //
-  // This function is the driver for all the explicit methods. It calls
-  // <code>evolve_one_time_step</code> which performs one time step.  For
-  // explicit methods, <code>evolve_one_time_step</code> needs to evaluate
-  // $M^{-1}(f(t,y))$, i.e, it needs <code>evaluate_diffusion</code>. Because
-  // <code>evaluate_diffusion</code> is a member function, it needs to be bound
-  // to <code>this</code>. Finally, the solution is output every 10 time steps.
+  // This function is the driver for all the explicit methods. At the
+  // top it initializes the time stepping and the solution (by setting
+  // it to zero and then ensuring that boundary value and hanging node
+  // constraints are respected; of course, with the mesh we use here,
+  // hanging node constraints are not in fact an issue). It then calls
+  // <code>evolve_one_time_step</code> which performs one time step.
+  //
+  // For explicit methods, <code>evolve_one_time_step</code> needs to
+  // evaluate $M^{-1}(f(t,y))$, i.e, it needs
+  // <code>evaluate_diffusion</code>. Because
+  // <code>evaluate_diffusion</code> is a member function, it needs to
+  // be bound to <code>this</code>. After each evolution step, we
+  // again apply the correct boundary values and hanging node
+  // constraints.
+  //
+  // Finally, the solution is output
+  // every 10 time steps.
   void Diffusion::explicit_method(const TimeStepping::runge_kutta_method method,
                                   const unsigned int n_time_steps,
                                   const double       initial_time,
@@ -480,7 +491,9 @@ namespace Step52
     const double time_step =
       (final_time - initial_time) / static_cast<double>(n_time_steps);
     double time = initial_time;
-    solution    = 0.;
+
+    solution = 0.;
+    constraint_matrix.distribute(solution);
 
     TimeStepping::ExplicitRungeKutta<Vector<double>> explicit_runge_kutta(
       method);
@@ -496,6 +509,8 @@ namespace Step52
           time_step,
           solution);
 
+        constraint_matrix.distribute(solution);
+
         if ((i + 1) % 10 == 0)
           output_results(i + 1, method);
       }
@@ -517,7 +532,9 @@ namespace Step52
     const double time_step =
       (final_time - initial_time) / static_cast<double>(n_time_steps);
     double time = initial_time;
-    solution    = 0.;
+
+    solution = 0.;
+    constraint_matrix.distribute(solution);
 
     TimeStepping::ImplicitRungeKutta<Vector<double>> implicit_runge_kutta(
       method);
@@ -538,6 +555,8 @@ namespace Step52
           time_step,
           solution);
 
+        constraint_matrix.distribute(solution);
+
         if ((i + 1) % 10 == 0)
           output_results(i + 1, method);
       }
@@ -576,7 +595,9 @@ namespace Step52
     const double max_delta     = 10 * time_step;
     const double refine_tol    = 1e-1;
     const double coarsen_tol   = 1e-5;
-    solution                   = 0.;
+
+    solution = 0.;
+    constraint_matrix.distribute(solution);
 
     TimeStepping::EmbeddedExplicitRungeKutta<Vector<double>>
       embedded_explicit_runge_kutta(method,
@@ -605,6 +626,8 @@ namespace Step52
           time_step,
           solution);
 
+        constraint_matrix.distribute(solution);
+
         if ((n_steps + 1) % 10 == 0)
           output_results(n_steps + 1, method);
 
-- 
2.39.5