From 0551f9593e3ea9d6e35a3638e8d6866ade9613c7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 10 Oct 2023 10:28:41 -0600 Subject: [PATCH] Inline the action of the three individual functions into add_constraint(). --- include/deal.II/lac/affine_constraints.h | 57 ++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 3 deletions(-) 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::add_constraint( const ArrayView> &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 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(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; } -- 2.39.5