Merge pull request #16145 from bangerth/add-constraints

Use AffineConstraints::add_constraint() throughout.

examples/step-16/step-16.cc

@@ -436,10 +436,13 @@ namespace Step16
         boundary_constraints[level].reinit(
           dof_handler.locally_owned_mg_dofs(level),
           DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_refinement_edge_indices(level));
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_boundary_indices(level));
+
+        for (const types::global_dof_index dof_index :
+             mg_constrained_dofs.get_refinement_edge_indices(level))
+          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+        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].close();
       }
 

examples/step-37/step-37.cc

@@ -852,8 +852,9 @@ namespace Step37
           AffineConstraints<double> level_constraints(
             dof_handler.locally_owned_mg_dofs(level),
             DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
-          level_constraints.add_lines(
-            mg_constrained_dofs.get_boundary_indices(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.close();
 
           typename MatrixFree<dim, float>::AdditionalData additional_data;

examples/step-42/step-42.cc

@@ -1275,9 +1275,9 @@ namespace Step42
                            0) &&
                           !constraints_hanging_nodes.is_constrained(index_z))
                         {
-                          all_constraints.add_line(index_z);
-                          all_constraints.set_inhomogeneity(index_z,
-                                                            undeformed_gap);
+                          all_constraints.add_constraint(index_z,
+                                                         {},
+                                                         undeformed_gap);
                           distributed_solution(index_z) = undeformed_gap;
 
                           active_set.add_index(index_z);

examples/step-50/step-50.cc

@@ -596,8 +596,9 @@ void LaplaceProblem<dim, degree>::setup_multigrid()
                 dof_handler.locally_owned_mg_dofs(level),
                 DoFTools::extract_locally_relevant_level_dofs(dof_handler,
                                                               level));
-              level_constraints.add_lines(
-                mg_constrained_dofs.get_boundary_indices(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.close();
 
               typename MatrixFree<dim, float>::AdditionalData additional_data;
@@ -826,11 +827,13 @@ void LaplaceProblem<dim, degree>::assemble_multigrid()
       boundary_constraints[level].reinit(
         dof_handler.locally_owned_mg_dofs(level),
         DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
-      boundary_constraints[level].add_lines(
-        mg_constrained_dofs.get_refinement_edge_indices(level));
-      boundary_constraints[level].add_lines(
-        mg_constrained_dofs.get_boundary_indices(level));
 
+      for (const types::global_dof_index dof_index :
+           mg_constrained_dofs.get_refinement_edge_indices(level))
+        boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+      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].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_line(n_u);
+        constraints.add_constraint(n_u, {}, 0.);
 
       constraints.close();
     }
@@ -732,16 +732,22 @@ namespace Step56
       triangulation.n_levels());
     for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
       {
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_refinement_edge_indices(level));
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_boundary_indices(level));
+        for (const types::global_dof_index dof_index :
+             mg_constrained_dofs.get_refinement_edge_indices(level))
+          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+        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].close();
 
-        IndexSet idx = mg_constrained_dofs.get_refinement_edge_indices(level) &
-                       mg_constrained_dofs.get_boundary_indices(level);
+        const IndexSet idx =
+          mg_constrained_dofs.get_refinement_edge_indices(level) &
+          mg_constrained_dofs.get_boundary_indices(level);
 
-        boundary_interface_constraints[level].add_lines(idx);
+        for (const types::global_dof_index dof_index : idx)
+          boundary_interface_constraints[level].add_constraint(dof_index,
+                                                               {},
+                                                               0.);
         boundary_interface_constraints[level].close();
       }
 

examples/step-63/step-63.cc

@@ -817,10 +817,13 @@ namespace Step63
         boundary_constraints[level].reinit(
           dof_handler.locally_owned_mg_dofs(level),
           DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_refinement_edge_indices(level));
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_boundary_indices(level));
+
+        for (const types::global_dof_index dof_index :
+             mg_constrained_dofs.get_refinement_edge_indices(level))
+          boundary_constraints[level].add_constraint(dof_index, {}, 0.);
+        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].close();
       }
 

examples/step-66/step-66.cc

@@ -593,8 +593,10 @@ namespace Step66
         AffineConstraints<double> level_constraints(
           dof_handler.locally_owned_mg_dofs(level),
           DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
-        level_constraints.add_lines(
-          mg_constrained_dofs.get_boundary_indices(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.close();
 
         typename MatrixFree<dim, float>::AdditionalData additional_data;

examples/step-79/step-79.cc

@@ -536,13 +536,16 @@ namespace SAND
     types::global_dof_index last_density_dof =
       density_dofs.nth_index_in_set(density_dofs.n_elements() - 1);
     constraints.clear();
-    constraints.add_line(last_density_dof);
-    for (unsigned int i = 0; i < density_dofs.n_elements() - 1; ++i)
-      constraints.add_entry(last_density_dof,
-                            density_dofs.nth_index_in_set(i),
-                            -1);
-    constraints.set_inhomogeneity(last_density_dof, 0);
-
+    {
+      std::vector<std::pair<types::global_dof_index, double>>
+        constraint_entries;
+      constraint_entries.reserve(density_dofs.n_elements() - 1);
+      for (const types::global_dof_index dof_index : density_dofs)
+        if (dof_index != last_density_dof)
+          constraint_entries.emplace_back(dof_index, -1.);
+
+      constraints.add_constraint(last_density_dof, constraint_entries, 0.);
+    }
     constraints.close();
 
     // We can now finally create the sparsity pattern for the

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

@@ -568,18 +568,7 @@ AffineConstraints<number>::make_consistent_in_parallel(
 
   // 3) refill this constraint matrix
   for (const auto &line : temporal_constraint_matrix)
-    {
-      // ... line
-      this->add_line(line.index);
-
-      // ... inhomogeneity
-      if (line.inhomogeneity != number())
-        this->set_inhomogeneity(line.index, line.inhomogeneity);
-
-      // ... entries
-      if (!line.entries.empty())
-        this->add_entries(line.index, line.entries);
-    }
+    this->add_constraint(line.index, line.entries, line.inhomogeneity);
 
 #ifdef DEBUG
   Assert(this->is_consistent_in_parallel(
@@ -1175,8 +1164,7 @@ AffineConstraints<number>::get_view(const IndexSet &mask) const
   for (const ConstraintLine &line : lines)
     if (mask.is_element(line.index))
       {
-        const size_type row = mask.index_within_set(line.index);
-        view.add_line(row);
+#ifdef DEBUG
         for (const std::pair<size_type, number> &entry : line.entries)
           {
             Assert(
@@ -1185,18 +1173,24 @@ AffineConstraints<number>::get_view(const IndexSet &mask) const
                 "In creating a view of an AffineConstraints "
                 "object, the constraint on degree of freedom " +
                 std::to_string(line.index) + " (which corresponds to the " +
-                std::to_string(row) +
+                std::to_string(mask.index_within_set(line.index)) +
                 "th degree of freedom selected in the mask) "
                 "is constrained against degree of freedom " +
                 std::to_string(entry.first) +
                 ", but this degree of freedom is not listed in the mask and "
                 "consequently cannot be transcribed into the index space "
                 "of the output object."));
-            view.add_entry(row,
-                           mask.index_within_set(entry.first),
-                           entry.second);
           }
-        view.set_inhomogeneity(row, line.inhomogeneity);
+#endif
+
+        std::vector<std::pair<size_type, number>> translated_entries =
+          line.entries;
+        for (auto &entry : translated_entries)
+          entry.first = mask.index_within_set(entry.first);
+
+        view.add_constraint(mask.index_within_set(line.index),
+                            translated_entries,
+                            line.inhomogeneity);
       }
 
   view.close();

include/deal.II/multigrid/mg_constrained_dofs.h

@@ -356,8 +356,9 @@ MGConstrainedDoFs::initialize(
                       !level_constraints[l].is_constrained(dofs_2[i]) &&
                       !level_constraints[l].is_constrained(dofs_1[i]))
                     {
-                      level_constraints[l].add_line(dofs_2[i]);
-                      level_constraints[l].add_entry(dofs_2[i], dofs_1[i], 1.);
+                      level_constraints[l].add_constraint(dofs_2[i],
+                                                          {{dofs_1[i], 1.}},
+                                                          0.);
                     }
               }
       level_constraints[l].close();
@@ -624,10 +625,12 @@ MGConstrainedDoFs::merge_constraints(AffineConstraints<Number> &constraints,
 
   // merge constraints
   if (add_boundary_indices && this->have_boundary_indices())
-    constraints.add_lines(this->get_boundary_indices(level));
+    for (const auto i : this->get_boundary_indices(level))
+      constraints.add_constraint(i, {}, 0.);
 
   if (add_refinement_edge_indices)
-    constraints.add_lines(this->get_refinement_edge_indices(level));
+    for (const auto i : this->get_refinement_edge_indices(level))
+      constraints.add_constraint(i, {}, 0.);
 
   if (add_level_constraints)
     constraints.merge(this->get_level_constraints(level),