step-72: use AffineConstraints only

Same as done in step-15 in #16981

diff --git a/examples/step-72/step-72.cc b/examples/step-72/step-72.cc
index b0479b77e5d0745c1ffe187a0dc227c7c5c23c92..b1fe79105cacb7a7d8c5ee1472c0c5659f14caeb 100644 (file)
--- a/examples/step-72/step-72.cc
+++ b/examples/step-72/step-72.cc
@@ -146,13 +146,12 @@ namespace Step72
     void run(const int formulation, const double tolerance);
 
   private:
-    void   setup_system(const bool initial_step);
+    void   setup_system();
     void   assemble_system_unassisted();
     void   assemble_system_with_residual_linearization();
     void   assemble_system_using_energy_functional();
     void   solve();
     void   refine_mesh();
-    void   set_boundary_values();
     double compute_residual(const double alpha) const;
     double determine_step_length() const;
     void   output_results(const unsigned int refinement_cycle) const;
@@ -163,7 +162,8 @@ namespace Step72
     const FE_Q<dim>   fe;
     const QGauss<dim> quadrature_formula;
 
-    AffineConstraints<double> hanging_node_constraints;
+    AffineConstraints<double> nonzero_constraints;
+    AffineConstraints<double> zero_constraints;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
@@ -212,26 +212,30 @@ namespace Step72
   // data structures, namely the DoFHandler, the hanging node constraints
   // applied to the problem, and the linear system.
   template <int dim>
-  void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
+  void MinimalSurfaceProblem<dim>::setup_system()
   {
-    if (initial_step)
-      {
-        dof_handler.distribute_dofs(fe);
-        current_solution.reinit(dof_handler.n_dofs());
-
-        hanging_node_constraints.clear();
-        DoFTools::make_hanging_node_constraints(dof_handler,
-                                                hanging_node_constraints);
-        hanging_node_constraints.close();
-      }
-
+    dof_handler.distribute_dofs(fe);
+    current_solution.reinit(dof_handler.n_dofs());
     newton_update.reinit(dof_handler.n_dofs());
     system_rhs.reinit(dof_handler.n_dofs());
 
-    DynamicSparsityPattern dsp(dof_handler.n_dofs());
-    DoFTools::make_sparsity_pattern(dof_handler, dsp);
+    zero_constraints.clear();
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             zero_constraints);
+    DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints);
+    zero_constraints.close();
 
-    hanging_node_constraints.condense(dsp);
+    nonzero_constraints.clear();
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             BoundaryValues<dim>(),
+                                             nonzero_constraints);
+    nonzero_constraints.close();
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, zero_constraints);
 
     sparsity_pattern.copy_from(dsp);
     system_matrix.reinit(sparsity_pattern);
@@ -385,15 +389,8 @@ namespace Step72
       const std::vector<types::global_dof_index> &local_dof_indices =
         copy_data.local_dof_indices[0];
 
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        {
-          for (unsigned int j = 0; j < dofs_per_cell; ++j)
-            system_matrix.add(local_dof_indices[i],
-                              local_dof_indices[j],
-                              cell_matrix(i, j));
-
-          system_rhs(local_dof_indices[i]) += cell_rhs(i);
-        }
+      zero_constraints.distribute_local_to_global(
+        cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
     };
 
     // We have all of the required functions definitions in place, so
@@ -410,22 +407,6 @@ namespace Step72
                           sample_scratch_data,
                           sample_copy_data,
                           MeshWorker::assemble_own_cells);
-
-    // And finally, as is done in step-15, we remove hanging nodes from the
-    // system and apply zero boundary values to the linear system that defines
-    // the Newton updates $\delta u^n$.
-    hanging_node_constraints.condense(system_matrix);
-    hanging_node_constraints.condense(system_rhs);
-
-    std::map<types::global_dof_index, double> boundary_values;
-    VectorTools::interpolate_boundary_values(dof_handler,
-                                             0,
-                                             Functions::ZeroFunction<dim>(),
-                                             boundary_values);
-    MatrixTools::apply_boundary_values(boundary_values,
-                                       system_matrix,
-                                       newton_update,
-                                       system_rhs);
   }
 
   // @sect5{Assembly via differentiation of the residual vector}
@@ -617,15 +598,8 @@ namespace Step72
       const std::vector<types::global_dof_index> &local_dof_indices =
         copy_data.local_dof_indices[0];
 
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        {
-          for (unsigned int j = 0; j < dofs_per_cell; ++j)
-            system_matrix.add(local_dof_indices[i],
-                              local_dof_indices[j],
-                              cell_matrix(i, j));
-
-          system_rhs(local_dof_indices[i]) += cell_rhs(i);
-        }
+      zero_constraints.distribute_local_to_global(
+        cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
     };
 
     MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
@@ -634,19 +608,6 @@ namespace Step72
                           sample_scratch_data,
                           sample_copy_data,
                           MeshWorker::assemble_own_cells);
-
-    hanging_node_constraints.condense(system_matrix);
-    hanging_node_constraints.condense(system_rhs);
-
-    std::map<types::global_dof_index, double> boundary_values;
-    VectorTools::interpolate_boundary_values(dof_handler,
-                                             0,
-                                             Functions::ZeroFunction<dim>(),
-                                             boundary_values);
-    MatrixTools::apply_boundary_values(boundary_values,
-                                       system_matrix,
-                                       newton_update,
-                                       system_rhs);
   }
 
   // @sect5{Assembly via differentiation of the energy functional}
@@ -785,15 +746,8 @@ namespace Step72
       const std::vector<types::global_dof_index> &local_dof_indices =
         copy_data.local_dof_indices[0];
 
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        {
-          for (unsigned int j = 0; j < dofs_per_cell; ++j)
-            system_matrix.add(local_dof_indices[i],
-                              local_dof_indices[j],
-                              cell_matrix(i, j));
-
-          system_rhs(local_dof_indices[i]) += cell_rhs(i);
-        }
+      zero_constraints.distribute_local_to_global(
+        cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
     };
 
     MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
@@ -802,19 +756,6 @@ namespace Step72
                           sample_scratch_data,
                           sample_copy_data,
                           MeshWorker::assemble_own_cells);
-
-    hanging_node_constraints.condense(system_matrix);
-    hanging_node_constraints.condense(system_rhs);
-
-    std::map<types::global_dof_index, double> boundary_values;
-    VectorTools::interpolate_boundary_values(dof_handler,
-                                             0,
-                                             Functions::ZeroFunction<dim>(),
-                                             boundary_values);
-    MatrixTools::apply_boundary_values(boundary_values,
-                                       system_matrix,
-                                       newton_update,
-                                       system_rhs);
   }
 
 
@@ -833,7 +774,7 @@ namespace Step72
 
     solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
 
-    hanging_node_constraints.distribute(newton_update);
+    zero_constraints.distribute(newton_update);
 
     const double alpha = determine_step_length();
     current_solution.add(alpha, newton_update);
@@ -862,46 +803,21 @@ namespace Step72
                                                     0.03);
 
     triangulation.prepare_coarsening_and_refinement();
-    SolutionTransfer<dim> solution_transfer(dof_handler);
-    solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
-    triangulation.execute_coarsening_and_refinement();
-
-    dof_handler.distribute_dofs(fe);
 
-    Vector<double> tmp(dof_handler.n_dofs());
-    solution_transfer.interpolate(current_solution, tmp);
-    current_solution = tmp;
+    SolutionTransfer<dim> solution_transfer(dof_handler);
+    const Vector<double>  coarse_solution = current_solution;
+    solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution);
 
-    hanging_node_constraints.clear();
-    DoFTools::make_hanging_node_constraints(dof_handler,
-                                            hanging_node_constraints);
-    hanging_node_constraints.close();
+    triangulation.execute_coarsening_and_refinement();
 
-    set_boundary_values();
+    setup_system();
 
-    setup_system(false);
+    solution_transfer.interpolate(coarse_solution, current_solution);
+    nonzero_constraints.distribute(current_solution);
   }
 
 
 
-  // @sect4{MinimalSurfaceProblem::set_boundary_values}
-
-  // The choice of boundary conditions remains identical to step-15...
-  template <int dim>
-  void MinimalSurfaceProblem<dim>::set_boundary_values()
-  {
-    std::map<types::global_dof_index, double> boundary_values;
-    VectorTools::interpolate_boundary_values(dof_handler,
-                                             0,
-                                             BoundaryValues<dim>(),
-                                             boundary_values);
-    for (auto &boundary_value : boundary_values)
-      current_solution(boundary_value.first) = boundary_value.second;
-
-    hanging_node_constraints.distribute(current_solution);
-  }
-
-
   // @sect4{MinimalSurfaceProblem::compute_residual}
 
   // ... as does the function used to compute the residual during the
@@ -951,16 +867,11 @@ namespace Step72
           }
 
         cell->get_dof_indices(local_dof_indices);
-        for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          residual(local_dof_indices[i]) += cell_residual(i);
+        zero_constraints.distribute_local_to_global(cell_residual,
+                                                    local_dof_indices,
+                                                    residual);
       }
 
-    hanging_node_constraints.condense(residual);
-
-    for (const types::global_dof_index i :
-         DoFTools::extract_boundary_dofs(dof_handler))
-      residual(i) = 0;
-
     return residual.l2_norm();
   }
 
@@ -1033,8 +944,8 @@ namespace Step72
     GridGenerator::hyper_ball(triangulation);
     triangulation.refine_global(2);
 
-    setup_system(/*first time=*/true);
-    set_boundary_values();
+    setup_system();
+    nonzero_constraints.distribute(current_solution);
 
     double       last_residual_norm = std::numeric_limits<double>::max();
     unsigned int refinement_cycle   = 0;