]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use AffineConstraints::add_constraint() throughout.
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 12 Oct 2023 22:37:35 +0000 (16:37 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sun, 15 Oct 2023 17:24:49 +0000 (11:24 -0600)
13 files changed:
examples/step-16/step-16.cc
examples/step-37/step-37.cc
examples/step-42/step-42.cc
examples/step-50/step-50.cc
examples/step-56/step-56.cc
examples/step-63/step-63.cc
examples/step-66/step-66.cc
examples/step-79/step-79.cc
include/deal.II/lac/affine_constraints.templates.h
include/deal.II/multigrid/mg_constrained_dofs.h
include/deal.II/numerics/vector_tools_constraints.templates.h
source/dofs/dof_tools_constraints.cc
source/grid/grid_tools.cc

index 09c4302d7de23816bef67673c67dd9731f73b07f..a79e2f1fde5c0772b8f9a5c0b7f8fc18fd19920b 100644 (file)
@@ -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();
       }
 
index f428beb6c666a6d5dbe1eb9c9ef74154427134c0..a4e98ab431c0eeaaabacc07e0340d8fd73f8a292 100644 (file)
@@ -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;
index 3fcd08064db4b9277761b0e9ae48928877e815a6..9d1374c245be7c36c9ec042bd35b0bd16be42ac7 100644 (file)
@@ -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);
index 6492d805cffc44ada43e5e95c7f2f16e4f871341..30d0327651e92555c85a73f68e0cb3113097e9b5 100644 (file)
@@ -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();
     }
 
index 8231689a7bb40eb4219f2549cae62847552869f2..1dca4bc59e69f64c4321dcf62fa648e118db01a5 100644 (file)
@@ -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();
       }
 
index bfead4ce5f8e6ad73283856fb8d389a3884f5e4d..81c59e9f64d1bebf7012fff0d61c643deef6b903 100644 (file)
@@ -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();
       }
 
index d5c258506329c470daaec0f239442e3a217156e9..504042a676753a039ca0e0980cc5129beb4fda56 100644 (file)
@@ -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;
index faf392acd59ad03419c78831a8b3dda64e1b8ced..7d4a69c333161f9d8e9d5d2eb334a327fb76fc8c 100644 (file)
@@ -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
index 83d9847583d8627326f1c86383d3256192f57818..a6e7abaf1d2140956aaebc3b8011ec0ee9f57e6a 100644 (file)
@@ -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(
@@ -1193,8 +1182,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(
@@ -1203,18 +1191,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();
index dd3e9a2b394903a392f4a093c2319958450e5591..e176b208e107b833eb0d9181e648770b97a13e9e 100644 (file)
@@ -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),
index 7a28a5efb83c638d7e276ac9e7cb39e14ea22844..b5a2f681869c48a82914c9cfc7707fd4fb096c6e 100644 (file)
@@ -170,21 +170,24 @@ namespace VectorTools
                   if (!constraints.is_constrained(dof_indices.dof_indices[0]) &&
                       constraints.can_store_line(dof_indices.dof_indices[0]))
                     {
-                      constraints.add_line(dof_indices.dof_indices[0]);
+                      const double normalized_inhomogeneity =
+                        (std::fabs(inhomogeneity / constraining_vector[0]) >
+                             std::numeric_limits<double>::epsilon() ?
+                           inhomogeneity / constraining_vector[0] :
+                           0);
 
                       if (std::fabs(constraining_vector[1] /
                                     constraining_vector[0]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[0],
-                                              dof_indices.dof_indices[1],
-                                              -constraining_vector[1] /
-                                                constraining_vector[0]);
-
-                      if (std::fabs(inhomogeneity / constraining_vector[0]) >
-                          std::numeric_limits<double>::epsilon())
-                        constraints.set_inhomogeneity(
-                          dof_indices.dof_indices[0],
-                          inhomogeneity / constraining_vector[0]);
+                        constraints.add_constraint(dof_indices.dof_indices[0],
+                                                   {{dof_indices.dof_indices[1],
+                                                     -constraining_vector[1] /
+                                                       constraining_vector[0]}},
+                                                   normalized_inhomogeneity);
+                      else
+                        constraints.add_constraint(dof_indices.dof_indices[0],
+                                                   {},
+                                                   normalized_inhomogeneity);
                     }
                 }
               else
@@ -192,21 +195,24 @@ namespace VectorTools
                   if (!constraints.is_constrained(dof_indices.dof_indices[1]) &&
                       constraints.can_store_line(dof_indices.dof_indices[1]))
                     {
-                      constraints.add_line(dof_indices.dof_indices[1]);
+                      const double normalized_inhomogeneity =
+                        (std::fabs(inhomogeneity / constraining_vector[1]) >
+                             std::numeric_limits<double>::epsilon() ?
+                           inhomogeneity / constraining_vector[1] :
+                           0);
 
                       if (std::fabs(constraining_vector[0] /
                                     constraining_vector[1]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[1],
-                                              dof_indices.dof_indices[0],
-                                              -constraining_vector[0] /
-                                                constraining_vector[1]);
-
-                      if (std::fabs(inhomogeneity / constraining_vector[1]) >
-                          std::numeric_limits<double>::epsilon())
-                        constraints.set_inhomogeneity(
-                          dof_indices.dof_indices[1],
-                          inhomogeneity / constraining_vector[1]);
+                        constraints.add_constraint(dof_indices.dof_indices[1],
+                                                   {{dof_indices.dof_indices[0],
+                                                     -constraining_vector[0] /
+                                                       constraining_vector[1]}},
+                                                   normalized_inhomogeneity);
+                      else
+                        constraints.add_constraint(dof_indices.dof_indices[1],
+                                                   {},
+                                                   normalized_inhomogeneity);
                     }
                 }
               break;
@@ -214,6 +220,16 @@ namespace VectorTools
 
           case 3:
             {
+              // Store constraints. There are at most 2 entries per
+              // constraint, so a boost::container::small_vector with
+              // static size 2 will do just fine.
+              //
+              // TODO: This could use std::inplace_vector once available (in
+              // C++26?)
+              boost::container::
+                small_vector<std::pair<types::global_dof_index, double>, 2>
+                  constraint_entries;
+
               if ((std::fabs(constraining_vector[0]) >=
                    std::fabs(constraining_vector[1]) + 1e-10) &&
                   (std::fabs(constraining_vector[0]) >=
@@ -222,29 +238,29 @@ namespace VectorTools
                   if (!constraints.is_constrained(dof_indices.dof_indices[0]) &&
                       constraints.can_store_line(dof_indices.dof_indices[0]))
                     {
-                      constraints.add_line(dof_indices.dof_indices[0]);
-
                       if (std::fabs(constraining_vector[1] /
                                     constraining_vector[0]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[0],
-                                              dof_indices.dof_indices[1],
-                                              -constraining_vector[1] /
-                                                constraining_vector[0]);
+                        constraint_entries.emplace_back(
+                          dof_indices.dof_indices[1],
+                          -constraining_vector[1] / constraining_vector[0]);
 
                       if (std::fabs(constraining_vector[2] /
                                     constraining_vector[0]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[0],
-                                              dof_indices.dof_indices[2],
-                                              -constraining_vector[2] /
-                                                constraining_vector[0]);
+                        constraint_entries.emplace_back(
+                          dof_indices.dof_indices[2],
+                          -constraining_vector[2] / constraining_vector[0]);
 
-                      if (std::fabs(inhomogeneity / constraining_vector[0]) >
-                          std::numeric_limits<double>::epsilon())
-                        constraints.set_inhomogeneity(
-                          dof_indices.dof_indices[0],
-                          inhomogeneity / constraining_vector[0]);
+                      const double normalized_inhomogeneity =
+                        (std::fabs(inhomogeneity / constraining_vector[0]) >
+                             std::numeric_limits<double>::epsilon() ?
+                           inhomogeneity / constraining_vector[0] :
+                           0);
+
+                      constraints.add_constraint(dof_indices.dof_indices[0],
+                                                 constraint_entries,
+                                                 normalized_inhomogeneity);
                     }
                 }
               else if ((std::fabs(constraining_vector[1]) + 1e-10 >=
@@ -255,29 +271,29 @@ namespace VectorTools
                   if (!constraints.is_constrained(dof_indices.dof_indices[1]) &&
                       constraints.can_store_line(dof_indices.dof_indices[1]))
                     {
-                      constraints.add_line(dof_indices.dof_indices[1]);
-
                       if (std::fabs(constraining_vector[0] /
                                     constraining_vector[1]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[1],
-                                              dof_indices.dof_indices[0],
-                                              -constraining_vector[0] /
-                                                constraining_vector[1]);
+                        constraint_entries.emplace_back(
+                          dof_indices.dof_indices[0],
+                          -constraining_vector[0] / constraining_vector[1]);
 
                       if (std::fabs(constraining_vector[2] /
                                     constraining_vector[1]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[1],
-                                              dof_indices.dof_indices[2],
-                                              -constraining_vector[2] /
-                                                constraining_vector[1]);
+                        constraint_entries.emplace_back(
+                          dof_indices.dof_indices[2],
+                          -constraining_vector[2] / constraining_vector[1]);
 
-                      if (std::fabs(inhomogeneity / constraining_vector[1]) >
-                          std::numeric_limits<double>::epsilon())
-                        constraints.set_inhomogeneity(
-                          dof_indices.dof_indices[1],
-                          inhomogeneity / constraining_vector[1]);
+                      const double normalized_inhomogeneity =
+                        (std::fabs(inhomogeneity / constraining_vector[1]) >
+                             std::numeric_limits<double>::epsilon() ?
+                           inhomogeneity / constraining_vector[1] :
+                           0);
+
+                      constraints.add_constraint(dof_indices.dof_indices[1],
+                                                 constraint_entries,
+                                                 normalized_inhomogeneity);
                     }
                 }
               else
@@ -285,29 +301,29 @@ namespace VectorTools
                   if (!constraints.is_constrained(dof_indices.dof_indices[2]) &&
                       constraints.can_store_line(dof_indices.dof_indices[2]))
                     {
-                      constraints.add_line(dof_indices.dof_indices[2]);
-
                       if (std::fabs(constraining_vector[0] /
                                     constraining_vector[2]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[2],
-                                              dof_indices.dof_indices[0],
-                                              -constraining_vector[0] /
-                                                constraining_vector[2]);
+                        constraint_entries.emplace_back(
+                          dof_indices.dof_indices[0],
+                          -constraining_vector[0] / constraining_vector[2]);
 
                       if (std::fabs(constraining_vector[1] /
                                     constraining_vector[2]) >
                           std::numeric_limits<double>::epsilon())
-                        constraints.add_entry(dof_indices.dof_indices[2],
-                                              dof_indices.dof_indices[1],
-                                              -constraining_vector[1] /
-                                                constraining_vector[2]);
+                        constraint_entries.emplace_back(
+                          dof_indices.dof_indices[1],
+                          -constraining_vector[1] / constraining_vector[2]);
 
-                      if (std::fabs(inhomogeneity / constraining_vector[2]) >
-                          std::numeric_limits<double>::epsilon())
-                        constraints.set_inhomogeneity(
-                          dof_indices.dof_indices[2],
-                          inhomogeneity / constraining_vector[2]);
+                      const double normalized_inhomogeneity =
+                        (std::fabs(inhomogeneity / constraining_vector[2]) >
+                             std::numeric_limits<double>::epsilon() ?
+                           inhomogeneity / constraining_vector[2] :
+                           0);
+
+                      constraints.add_constraint(dof_indices.dof_indices[2],
+                                                 constraint_entries,
+                                                 normalized_inhomogeneity);
                     }
                 }
 
@@ -363,25 +379,29 @@ namespace VectorTools
           if (!constraints.is_constrained(dof_indices.dof_indices[d]) &&
               constraints.can_store_line(dof_indices.dof_indices[d]))
             {
-              constraints.add_line(dof_indices.dof_indices[d]);
-
-              if (std::fabs(tangent_vector[d] /
-                            tangent_vector[largest_component]) >
-                  std::numeric_limits<double>::epsilon())
-                constraints.add_entry(
-                  dof_indices.dof_indices[d],
-                  dof_indices.dof_indices[largest_component],
-                  tangent_vector[d] / tangent_vector[largest_component]);
-
               const double inhomogeneity =
                 (b_values(d) * tangent_vector[largest_component] -
                  b_values(largest_component) * tangent_vector[d]) /
                 tangent_vector[largest_component];
 
-              if (std::fabs(inhomogeneity) >
+              const double normalized_inhomogeneity =
+                (std::fabs(inhomogeneity) >
+                     std::numeric_limits<double>::epsilon() ?
+                   inhomogeneity :
+                   0);
+
+              if (std::fabs(tangent_vector[d] /
+                            tangent_vector[largest_component]) >
                   std::numeric_limits<double>::epsilon())
-                constraints.set_inhomogeneity(dof_indices.dof_indices[d],
-                                              inhomogeneity);
+                constraints.add_constraint(
+                  dof_indices.dof_indices[d],
+                  {{dof_indices.dof_indices[largest_component],
+                    tangent_vector[d] / tangent_vector[largest_component]}},
+                  normalized_inhomogeneity);
+              else
+                constraints.add_constraint(dof_indices.dof_indices[d],
+                                           {},
+                                           normalized_inhomogeneity);
             }
     }
 
index 716de8eca869550566432db35c34f1733cfb0254..6ec936f7b2c2bb1c26113ccd4a1f252b056c04c9 100644 (file)
@@ -509,8 +509,14 @@ namespace DoFTools
           Assert(primary_dofs[col] != numbers::invalid_dof_index,
                  ExcInternalError());
 
-        std::vector<
-          std::pair<typename AffineConstraints<number2>::size_type, number2>>
+        // Build constraints in a vector of pairs that can be
+        // arbitrarily large, but that holds up to 25 elements without
+        // external memory allocation. This is good enough for hanging
+        // node constraints of Q4 elements in 3d, so covers most
+        // common cases.
+        boost::container::small_vector<
+          std::pair<typename AffineConstraints<number2>::size_type, number2>,
+          25>
           entries;
         entries.reserve(n_primary_dofs);
         for (unsigned int row = 0; row != n_dependent_dofs; ++row)
@@ -650,6 +656,16 @@ namespace DoFTools
       std::vector<types::global_dof_index> dofs_on_mother;
       std::vector<types::global_dof_index> dofs_on_children;
 
+      // Build constraints in a vector of pairs that can be
+      // arbitrarily large, but that holds up to 25 elements without
+      // external memory allocation. This is good enough for hanging
+      // node constraints of Q4 elements in 3d, so covers most
+      // common cases.
+      boost::container::small_vector<
+        std::pair<typename AffineConstraints<number>::size_type, number>,
+        25>
+        constraint_entries;
+
       // loop over all lines; only on lines there can be constraints. We do so
       // by looping over all active cells and checking whether any of the faces
       // are refined which can only be from the neighboring cell because this
@@ -658,10 +674,7 @@ namespace DoFTools
       // note that even though we may visit a face twice if the neighboring
       // cells are equally refined, we can only visit each face with hanging
       // nodes once
-      typename DoFHandler<dim_, spacedim>::active_cell_iterator
-        cell = dof_handler.begin_active(),
-        endc = dof_handler.end();
-      for (; cell != endc; ++cell)
+      for (const auto &cell : dof_handler.active_cell_iterators())
         {
           // artificial cells can at best neighbor ghost cells, but we're not
           // interested in these interfaces
@@ -760,13 +773,15 @@ namespace DoFTools
                 for (unsigned int row = 0; row != dofs_on_children.size();
                      ++row)
                   {
-                    constraints.add_line(dofs_on_children[row]);
+                    constraint_entries.clear();
+                    constraint_entries.reserve(dofs_on_mother.size());
                     for (unsigned int i = 0; i != dofs_on_mother.size(); ++i)
-                      constraints.add_entry(dofs_on_children[row],
-                                            dofs_on_mother[i],
-                                            fe.constraints()(row, i));
+                      constraint_entries.emplace_back(dofs_on_mother[i],
+                                                      fe.constraints()(row, i));
 
-                    constraints.set_inhomogeneity(dofs_on_children[row], 0.);
+                    constraints.add_constraint(dofs_on_children[row],
+                                               constraint_entries,
+                                               0.);
                   }
               }
             else
@@ -802,6 +817,16 @@ namespace DoFTools
       std::vector<types::global_dof_index> dofs_on_mother;
       std::vector<types::global_dof_index> dofs_on_children;
 
+      // Build constraints in a vector of pairs that can be
+      // arbitrarily large, but that holds up to 25 elements without
+      // external memory allocation. This is good enough for hanging
+      // node constraints of Q4 elements in 3d, so covers most
+      // common cases.
+      boost::container::small_vector<
+        std::pair<typename AffineConstraints<number>::size_type, number>,
+        25>
+        constraint_entries;
+
       // loop over all quads; only on quads there can be constraints. We do so
       // by looping over all active cells and checking whether any of the faces
       // are refined which can only be from the neighboring cell because this
@@ -810,10 +835,7 @@ namespace DoFTools
       // note that even though we may visit a face twice if the neighboring
       // cells are equally refined, we can only visit each face with hanging
       // nodes once
-      typename DoFHandler<dim_, spacedim>::active_cell_iterator
-        cell = dof_handler.begin_active(),
-        endc = dof_handler.end();
-      for (; cell != endc; ++cell)
+      for (const auto &cell : dof_handler.active_cell_iterators())
         {
           // artificial cells can at best neighbor ghost cells, but we're not
           // interested in these interfaces
@@ -1004,13 +1026,15 @@ namespace DoFTools
                 for (unsigned int row = 0; row != dofs_on_children.size();
                      ++row)
                   {
-                    constraints.add_line(dofs_on_children[row]);
+                    constraint_entries.clear();
+                    constraint_entries.reserve(dofs_on_mother.size());
                     for (unsigned int i = 0; i != dofs_on_mother.size(); ++i)
-                      constraints.add_entry(dofs_on_children[row],
-                                            dofs_on_mother[i],
-                                            fe.constraints()(row, i));
+                      constraint_entries.emplace_back(dofs_on_mother[i],
+                                                      fe.constraints()(row, i));
 
-                    constraints.set_inhomogeneity(dofs_on_children[row], 0.);
+                    constraints.add_constraint(dofs_on_children[row],
+                                               constraint_entries,
+                                               0.);
                   }
               }
             else
@@ -1993,11 +2017,20 @@ namespace DoFTools
           cell_to_rotated_face_index[cell_index] = i;
         }
 
+      // Build constraints in a vector of pairs that can be
+      // arbitrarily large, but that holds up to 25 elements without
+      // external memory allocation. This is good enough for hanging
+      // node constraints of Q4 elements in 3d, so covers most
+      // common cases.
+      boost::container::small_vector<
+        std::pair<typename AffineConstraints<number>::size_type, number>,
+        25>
+        constraint_entries;
+
       //
       // Loop over all dofs on face 2 and constrain them against all
       // matching dofs on face 1:
       //
-
       for (unsigned int i = 0; i < dofs_per_face; ++i)
         {
           // Obey the component mask
@@ -2055,8 +2088,8 @@ namespace DoFTools
               if (affine_constraints.is_constrained(dofs_2[i]))
                 continue;
 
-              // Enter the constraint piece by piece:
-              affine_constraints.add_line(dofs_2[i]);
+              constraint_entries.clear();
+              constraint_entries.reserve(dofs_per_face);
 
               for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
                 {
@@ -2067,11 +2100,16 @@ namespace DoFTools
                       jj, 0, face_orientation, face_flip, face_rotation)];
 
                   if (std::abs(transformation(i, jj)) > eps)
-                    affine_constraints.add_entry(dofs_2[i],
-                                                 dofs_1[j],
-                                                 transformation(i, jj));
+                    constraint_entries.emplace_back(dofs_1[j],
+                                                    transformation(i, jj));
                 }
 
+              // Enter the constraint::
+              affine_constraints.add_constraint(dofs_2[i],
+                                                constraint_entries,
+                                                0.);
+
+
               // Continue with next dof.
               continue;
             }
@@ -2149,14 +2187,12 @@ namespace DoFTools
           if (constraints_are_cyclic)
             {
               if (std::abs(cycle_constraint_factor - number(1.)) > eps)
-                affine_constraints.add_line(dof_left);
+                affine_constraints.add_constraint(dof_left, {}, 0.);
             }
           else
             {
-              affine_constraints.add_line(dof_left);
-              affine_constraints.add_entry(dof_left,
-                                           dof_right,
-                                           constraint_factor);
+              affine_constraints.add_constraint(
+                dof_left, {{dof_right, constraint_factor}}, 0.);
               // The number 1e10 in the assert below is arbitrary. If the
               // absolute value of constraint_factor is too large, then probably
               // the absolute value of periodicity_factor is too large or too
@@ -3339,8 +3375,6 @@ namespace DoFTools
 
 
           // otherwise enter all constraints
-          constraints.add_line(global_dof);
-
           constraint_line.clear();
           for (types::global_dof_index row = first_used_row;
                row < n_coarse_dofs;
@@ -3352,7 +3386,7 @@ namespace DoFTools
                 constraint_line.emplace_back(representants[row], j->second);
             }
 
-          constraints.add_entries(global_dof, constraint_line);
+          constraints.add_constraint(global_dof, constraint_line, 0.);
         }
   }
 
@@ -3483,10 +3517,14 @@ namespace DoFTools
     std::vector<types::global_dof_index> cell_dofs;
     cell_dofs.reserve(dof.get_fe_collection().max_dofs_per_cell());
 
-    typename DoFHandler<dim, spacedim>::active_cell_iterator
-      cell = dof.begin_active(),
-      endc = dof.end();
-    for (; cell != endc; ++cell)
+    // In looping over faces, we will encounter some DoFs multiple
+    // times (namely, the ones on vertices and (in 3d) edges shared
+    // between multiple boundary faces. Keep track of which DoFs we
+    // have already encountered, so that we do not have to consider
+    // them a second time.
+    std::set<types::global_dof_index> dofs_already_treated;
+
+    for (const auto &cell : dof.active_cell_iterators())
       if (!cell->is_artificial() && cell->at_boundary())
         {
           const FiniteElement<dim, spacedim> &fe = cell->get_fe();
@@ -3513,30 +3551,38 @@ namespace DoFTools
                   // enter those dofs into the list that match the component
                   // signature.
                   for (const types::global_dof_index face_dof : face_dofs)
-                    {
-                      // Find out if a dof has a contribution in this component,
-                      // and if so, add it to the list
-                      const std::vector<types::global_dof_index>::iterator
-                        it_index_on_cell = std::find(cell_dofs.begin(),
-                                                     cell_dofs.end(),
-                                                     face_dof);
-                      Assert(it_index_on_cell != cell_dofs.end(),
-                             ExcInvalidIterator());
-                      const unsigned int index_on_cell =
-                        std::distance(cell_dofs.begin(), it_index_on_cell);
-                      const ComponentMask &nonzero_component_array =
-                        cell->get_fe().get_nonzero_components(index_on_cell);
-                      bool nonzero = false;
-                      for (unsigned int c = 0; c < n_components; ++c)
-                        if (nonzero_component_array[c] && component_mask[c])
-                          {
-                            nonzero = true;
-                            break;
-                          }
-
-                      if (nonzero)
-                        zero_boundary_constraints.add_line(face_dof);
-                    }
+                    if (dofs_already_treated.find(face_dof) ==
+                        dofs_already_treated.end())
+                      {
+                        // Find out if a dof has a contribution in this
+                        // component, and if so, add it to the list
+                        const std::vector<types::global_dof_index>::iterator
+                          it_index_on_cell = std::find(cell_dofs.begin(),
+                                                       cell_dofs.end(),
+                                                       face_dof);
+                        Assert(it_index_on_cell != cell_dofs.end(),
+                               ExcInvalidIterator());
+                        const unsigned int index_on_cell =
+                          std::distance(cell_dofs.begin(), it_index_on_cell);
+                        const ComponentMask &nonzero_component_array =
+                          cell->get_fe().get_nonzero_components(index_on_cell);
+
+                        bool nonzero = false;
+                        for (unsigned int c = 0; c < n_components; ++c)
+                          if (nonzero_component_array[c] && component_mask[c])
+                            {
+                              nonzero = true;
+                              break;
+                            }
+
+                        if (nonzero)
+                          zero_boundary_constraints.add_constraint(face_dof,
+                                                                   {},
+                                                                   0.);
+                        // We already dealt with this DoF. Make sure we
+                        // don't touch it again.
+                        dofs_already_treated.insert(face_dof);
+                      }
                 }
             }
         }
index 6abda211d979adf9d27935bd01065536bd0d5d5b..b62a661b9091ba3fd99d0073f55aa4e15796a5ef 100644 (file)
@@ -2297,9 +2297,9 @@ namespace GridTools
             if (map_iter != map_end)
               for (unsigned int i = 0; i < dim; ++i)
                 {
-                  constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
-                  constraints[i].set_inhomogeneity(
+                  constraints[i].add_constraint(
                     cell->vertex_dof_index(vertex_no, 0),
+                    {},
                     (solve_for_absolute_positions ?
                        map_iter->second(i) :
                        map_iter->second(i) - vertex_point[i]));

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.