From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 10 Oct 2023 16:28:41 +0000 (-0600)
Subject: Inline the action of the three individual functions into add_constraint().
X-Git-Tag: relicensing~406^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0551f9593e3ea9d6e35a3638e8d6866ade9613c7;p=dealii.git

Inline the action of the three individual functions into add_constraint().
---

diff --git a/include/deal.II/lac/affine_constraints.h b/include/deal.II/lac/affine_constraints.h
index 470a405e91..4eb0e7fe68 100644
--- a/include/deal.II/lac/affine_constraints.h
+++ b/include/deal.II/lac/affine_constraints.h
@@ -2301,13 +2301,64 @@ AffineConstraints<number>::add_constraint(
   const ArrayView<const std::pair<size_type, number>> &dependencies,
   const number                                         inhomogeneity)
 {
+  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."));
 
-  add_line(constrained_dof);
-  add_entries(constrained_dof, dependencies);
-  set_inhomogeneity(constrained_dof, inhomogeneity);
+  // In debug mode, make sure that columns don't show up more than
+  // once. Do this by sorting an array of the column indices. Try to
+  // avoid memory allocation by using a vector type that allocates up
+  // to 25 entries on the stack -- this is enough for the constraints
+  // of Q4 elements in 3d, and so should cover the vast majority of
+  // cases. If we have a constraint with more dependencies, then
+  // that's just going to require a heap allocation.
+#ifdef DEBUG
+  {
+    boost::container::small_vector<size_type, 25> column_indices;
+    column_indices.reserve(dependencies.size());
+    for (const auto &d : dependencies)
+      column_indices.emplace_back(d.first);
+    std::sort(column_indices.begin(), column_indices.end());
+    Assert(std::adjacent_find(column_indices.begin(), column_indices.end()) ==
+             column_indices.end(),
+           ExcMessage(
+             "You are trying to insert a constraint that lists the same "
+             "degree of freedom more than once on the right hand side. This is "
+             "not allowed."));
+  }
+#endif
+
+
+  // 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.entries.reserve(dependencies.size());
+  for (const auto &[column, weight] : dependencies)
+    {
+      Assert(column != constrained_dof,
+             ExcMessage("You cannot constrain a degree of freedom against "
+                        "itself."));
+      constraint.entries.emplace_back(column, weight);
+    }
+  constraint.inhomogeneity = inhomogeneity;
+
+  // Record the new constraint in the cache:
+  lines_cache[line_index] = lines.size() - 1;
 }