Merge pull request #16163 from bangerth/constrain-to-zero

Add AffineConstraints::constrain_dof_to_zero().

doc/news/changes/minor/20231002Bangerth_2

@@ -2,6 +2,8 @@ New: The new function AffineConstraints::add_constraint() adds an
 entire constraint all at once, rather than splitting the task across
 AffineConstraint::add_line(), multiple calls to
 AffineConstraint::add_entry(), and
-AffineConstraint::set_inhomogeneity().
+AffineConstraint::set_inhomogeneity(). The new function
+AffineConstraints::constrain_dof_to_zero() is a short-cut that adds a
+constraint that requires a specific degree of freedom to be zero.
 <br>
 (Wolfgang Bangerth, 2023/10/05)

examples/step-16/step-16.cc

@@ -439,10 +439,10 @@ namespace Step16
 
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_refinement_edge_indices(level))
-          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+          boundary_constraints[level].constrain_dof_to_zero(dof_index);
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_boundary_indices(level))
-          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+          boundary_constraints[level].constrain_dof_to_zero(dof_index);
         boundary_constraints[level].close();
       }
 

examples/step-37/step-37.cc

@@ -854,7 +854,7 @@ namespace Step37
             DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
           for (const types::global_dof_index dof_index :
                mg_constrained_dofs.get_boundary_indices(level))
-            level_constraints.add_constraint(dof_index, {}, 0.);
+            level_constraints.constrain_dof_to_zero(dof_index);
           level_constraints.close();
 
           typename MatrixFree<dim, float>::AdditionalData additional_data;

examples/step-46/doc/intro.dox

@@ -494,13 +494,12 @@ for (const auto &cell: dof_handler.active_cell_iterators())
              cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
              for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
              if (stokes_fe.face_system_to_component_index(i).first < dim)
-               constraints.add_constraint (local_face_dof_indices[i], {}, 0.);
+               constraints.constrain_dof_to_zero (local_face_dof_indices[i]);
            }
         }
 @endcode
 
-The last line calls
-AffineConstraints::add_constraint() and adds the constraint
+The last line adds the constraint
 <i>x<sub>local_face_dof_indices[i]</sub>=0</i>, which is exactly
 what we need in the current context. The call to
 FiniteElement::face_system_to_component_index() makes sure that we only set

examples/step-46/step-46.cc

@@ -385,9 +385,8 @@ namespace Step46
                          ++i)
                       if (stokes_fe.face_system_to_component_index(i).first <
                           dim)
-                        constraints.add_constraint(local_face_dof_indices[i],
-                                                   {},
-                                                   0.);
+                        constraints.constrain_dof_to_zero(
+                          local_face_dof_indices[i]);
                   }
               }
     }

examples/step-50/step-50.cc

@@ -598,7 +598,7 @@ void LaplaceProblem<dim, degree>::setup_multigrid()
                                                               level));
               for (const types::global_dof_index dof_index :
                    mg_constrained_dofs.get_boundary_indices(level))
-                level_constraints.add_constraint(dof_index, {}, 0.);
+                level_constraints.constrain_dof_to_zero(dof_index);
               level_constraints.close();
 
               typename MatrixFree<dim, float>::AdditionalData additional_data;
@@ -830,10 +830,10 @@ void LaplaceProblem<dim, degree>::assemble_multigrid()
 
       for (const types::global_dof_index dof_index :
            mg_constrained_dofs.get_refinement_edge_indices(level))
-        boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+        boundary_constraints[level].constrain_dof_to_zero(dof_index);
       for (const types::global_dof_index dof_index :
            mg_constrained_dofs.get_boundary_indices(level))
-        boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+        boundary_constraints[level].constrain_dof_to_zero(dof_index);
       boundary_constraints[level].close();
     }
 

examples/step-56/step-56.cc

@@ -572,7 +572,7 @@ namespace Step56
       // this here by marking the first pressure dof, which has index n_u as a
       // constrained dof.
       if (solver_type == SolverType::UMFPACK)
-        constraints.add_constraint(n_u, {}, 0.);
+        constraints.constrain_dof_to_zero(n_u);
 
       constraints.close();
     }
@@ -734,10 +734,10 @@ namespace Step56
       {
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_refinement_edge_indices(level))
-          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+          boundary_constraints[level].constrain_dof_to_zero(dof_index);
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_boundary_indices(level))
-          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+          boundary_constraints[level].constrain_dof_to_zero(dof_index);
         boundary_constraints[level].close();
 
         const IndexSet idx =

examples/step-63/step-63.cc

@@ -820,10 +820,10 @@ namespace Step63
 
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_refinement_edge_indices(level))
-          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+          boundary_constraints[level].constrain_dof_to_zero(dof_index);
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_boundary_indices(level))
-          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+          boundary_constraints[level].constrain_dof_to_zero(dof_index);
         boundary_constraints[level].close();
       }
 

examples/step-66/step-66.cc

@@ -596,7 +596,7 @@ namespace Step66
 
         for (const types::global_dof_index dof_index :
              mg_constrained_dofs.get_boundary_indices(level))
-          level_constraints.add_constraint(dof_index, {}, 0.);
+          level_constraints.constrain_dof_to_zero(dof_index);
         level_constraints.close();
 
         typename MatrixFree<dim, float>::AdditionalData additional_data;

include/deal.II/lac/affine_constraints.h

@@ -793,13 +793,42 @@ class AffineConstraints : public Subscriptor
    * @code
    *   constraints.add_constraint (42, {}, 27.0);
    * @endcode
+   * If you want to constrain a degree of freedom to zero, i.e.,
+   * require that
+   * @f[
+   *   x_{42} = 0
+   * @f]
+   * you would call this function as follows:
+   * @code
+   *   constraints.add_constraint (42, {}, 0.0);
+   * @endcode
+   * That said, this special case can be achieved in a more obvious way by
+   * calling
+   * @code
+   *   constraints.constrain_dof_to_zero (42);
+   * @endcode
+   * instead.
    */
   void
   add_constraint(
     const size_type                                      constrained_dof,
     const ArrayView<const std::pair<size_type, number>> &dependencies,
     const number                                         inhomogeneity = 0);
 
+  /**
+   * Constrain the given degree of freedom to be zero, i.e.,
+   * require a constraint like
+   * @f[
+   *   x_{42} = 0.
+   * @f]
+   * Calling this function is equivalent to, but more readable than, saying
+   * @code
+   *   constraints.add_constraint (42, {}, 0.0);
+   * @endcode
+   */
+  void
+  constrain_dof_to_zero(const size_type constrained_dof);
+
   /**
    * Add a new line to the matrix. If the line already exists, then the
    * function simply returns without doing anything.

include/deal.II/lac/affine_constraints.templates.h

@@ -143,6 +143,41 @@ AffineConstraints<number>::add_constraint(
 
 
 
+template <typename number>
+void
+AffineConstraints<number>::constrain_dof_to_zero(
+  const size_type constrained_dof)
+{
+  Assert(sorted == false, ExcMatrixIsClosed());
+  Assert(is_constrained(constrained_dof) == false,
+         ExcMessage("You cannot add a constraint for a degree of freedom "
+                    "that is already constrained."));
+
+  // The following can happen when we compute with distributed meshes and dof
+  // handlers and we constrain a degree of freedom whose number we don't have
+  // locally. if we don't abort here the program will try to allocate several
+  // terabytes of memory to resize the various arrays below :-)
+  Assert(constrained_dof != numbers::invalid_size_type, ExcInternalError());
+
+  // if necessary enlarge vector of existing entries for cache
+  const size_type line_index = calculate_line_index(constrained_dof);
+  if (line_index >= lines_cache.size())
+    lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
+                                line_index + 1),
+                       numbers::invalid_size_type);
+
+  // Push a new line to the end of the list and fill it with the
+  // provided information:
+  ConstraintLine &constraint = lines.emplace_back();
+  constraint.index           = constrained_dof;
+  constraint.inhomogeneity   = 0.;
+
+  // Record the new constraint in the cache:
+  lines_cache[line_index] = lines.size() - 1;
+}
+
+
+
 template <typename number>
 typename AffineConstraints<number>::LineRange
 AffineConstraints<number>::get_lines() const

include/deal.II/multigrid/mg_constrained_dofs.h

@@ -613,11 +613,11 @@ MGConstrainedDoFs::merge_constraints(AffineConstraints<Number> &constraints,
   // merge constraints
   if (add_boundary_indices && this->have_boundary_indices())
     for (const auto i : this->get_boundary_indices(level))
-      constraints.add_constraint(i, {}, 0.);
+      constraints.constrain_dof_to_zero(i);
 
   if (add_refinement_edge_indices)
     for (const auto i : this->get_refinement_edge_indices(level))
-      constraints.add_constraint(i, {}, 0.);
+      constraints.constrain_dof_to_zero(i);
 
   if (add_level_constraints)
     constraints.merge(this->get_level_constraints(level),

source/dofs/dof_tools_constraints.cc

@@ -2190,7 +2190,7 @@ namespace DoFTools
           if (constraints_are_cyclic)
             {
               if (std::abs(cycle_constraint_factor - number(1.)) > eps)
-                affine_constraints.add_constraint(dof_left, {}, 0.);
+                affine_constraints.constrain_dof_to_zero(dof_left);
             }
           else
             {
@@ -3587,9 +3587,8 @@ namespace DoFTools
                             // inhomogeneity is zero:
                             if (zero_boundary_constraints.is_constrained(
                                   face_dof) == false)
-                              zero_boundary_constraints.add_constraint(face_dof,
-                                                                       {},
-                                                                       0.);
+                              zero_boundary_constraints.constrain_dof_to_zero(
+                                face_dof);
                             else
                               Assert(zero_boundary_constraints
                                          .is_inhomogeneously_constrained(