From 0ab1b65e9e1ec3bcb330f18f28866f33bb53360b Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 15 Jul 2023 20:00:06 +0200 Subject: [PATCH] AffineConstraints::merge(): allow mixed numbers --- include/deal.II/lac/affine_constraints.h | 157 +++++++++++++++++- .../lac/affine_constraints.templates.h | 143 ---------------- 2 files changed, 156 insertions(+), 144 deletions(-) diff --git a/include/deal.II/lac/affine_constraints.h b/include/deal.II/lac/affine_constraints.h index 49c1b64259..b786924f4c 100644 --- a/include/deal.II/lac/affine_constraints.h +++ b/include/deal.II/lac/affine_constraints.h @@ -821,9 +821,10 @@ public: * Merging a AffineConstraints that is initialized with an IndexSet * and one that is not initialized with an IndexSet is not yet implemented. */ + template void merge( - const AffineConstraints & other_constraints, + const AffineConstraints &other_constraints, const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed, const bool allow_different_local_lines = false); @@ -2468,6 +2469,160 @@ AffineConstraints::copy_from( } +template +template +void +AffineConstraints::merge( + const AffineConstraints &other_constraints, + const MergeConflictBehavior merge_conflict_behavior, + const bool allow_different_local_lines) +{ + (void)allow_different_local_lines; + Assert(allow_different_local_lines || + local_lines == other_constraints.local_lines, + ExcMessage( + "local_lines for this and the other objects are not the same " + "although allow_different_local_lines is false.")); + + // store the previous state with respect to sorting + const bool object_was_sorted = sorted; + sorted = false; + + // first action is to fold into the present object possible constraints + // in the second object. we don't strictly need to do this any more since + // the AffineConstraints container has learned to deal with chains of + // constraints in the close() function, but we have traditionally done + // this and it's not overly hard to do. + // + // for this, loop over all constraints and replace the constraint lines + // with a new one where constraints are replaced if necessary. + typename ConstraintLine::Entries tmp; + for (ConstraintLine &line : lines) + { + tmp.clear(); + for (const std::pair &entry : line.entries) + { + // if the present dof is not stored, or not constrained, or if we + // won't take the constraint from the other object, then simply copy + // it over + if ((other_constraints.local_lines.size() != 0. && + other_constraints.local_lines.is_element(entry.first) == + false) || + other_constraints.is_constrained(entry.first) == false || + ((merge_conflict_behavior != right_object_wins) && + other_constraints.is_constrained(entry.first) && + this->is_constrained(entry.first))) + tmp.push_back(entry); + else + // otherwise resolve further constraints by replacing the old + // entry by a sequence of new entries taken from the other + // object, but with multiplied weights + { + const auto *other_entries = + other_constraints.get_constraint_entries(entry.first); + Assert(other_entries != nullptr, ExcInternalError()); + + const number weight = entry.second; + + for (const std::pair &other_entry : + *other_entries) + tmp.emplace_back(other_entry.first, + other_entry.second * weight); + + line.inhomogeneity += + other_constraints.get_inhomogeneity(entry.first) * weight; + } + } + // finally exchange old and newly resolved line + line.entries.swap(tmp); + } + + if (local_lines.size() != 0) + local_lines.add_indices(other_constraints.local_lines); + + { + // do not bother to resize the lines cache exactly since it is pretty + // cheap to adjust it along the way. + std::fill(lines_cache.begin(), + lines_cache.end(), + numbers::invalid_size_type); + + // reset lines_cache for our own constraints + size_type index = 0; + for (const ConstraintLine &line : lines) + { + const size_type local_line_no = calculate_line_index(line.index); + if (local_line_no >= lines_cache.size()) + lines_cache.resize(local_line_no + 1, numbers::invalid_size_type); + lines_cache[local_line_no] = index++; + } + + // Add other_constraints to lines cache and our list of constraints + for (const auto &line : other_constraints.lines) + { + const size_type local_line_no = calculate_line_index(line.index); + if (local_line_no >= lines_cache.size()) + { + lines_cache.resize(local_line_no + 1, numbers::invalid_size_type); + lines.emplace_back(line.index, + typename ConstraintLine::Entries( + line.entries.begin(), line.entries.end()), + line.inhomogeneity); + lines_cache[local_line_no] = index++; + } + else if (lines_cache[local_line_no] == numbers::invalid_size_type) + { + // there are no constraints for that line yet + lines.emplace_back(line.index, + typename ConstraintLine::Entries( + line.entries.begin(), line.entries.end()), + line.inhomogeneity); + AssertIndexRange(local_line_no, lines_cache.size()); + lines_cache[local_line_no] = index++; + } + else + { + // we already store that line + switch (merge_conflict_behavior) + { + case no_conflicts_allowed: + AssertThrow(false, + ExcDoFIsConstrainedFromBothObjects(line.index)); + break; + + case left_object_wins: + // ignore this constraint + break; + + case right_object_wins: + AssertIndexRange(local_line_no, lines_cache.size()); + lines[lines_cache[local_line_no]] = { + line.index, + typename ConstraintLine::Entries(line.entries.begin(), + line.entries.end()), + line.inhomogeneity}; + break; + + default: + Assert(false, ExcNotImplemented()); + } + } + } + + // check that we set the pointers correctly + for (size_type i = 0; i < lines_cache.size(); ++i) + if (lines_cache[i] != numbers::invalid_size_type) + Assert(i == calculate_line_index(lines[lines_cache[i]].index), + ExcInternalError()); + } + + // if the object was sorted before, then make sure it is so afterward as + // well. otherwise leave everything in the unsorted state + if (object_was_sorted == true) + close(); +} + + template template diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index d465c9f710..7f37bbdac1 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -1014,149 +1014,6 @@ AffineConstraints::is_closed(const MPI_Comm comm) const } -template -void -AffineConstraints::merge( - const AffineConstraints &other_constraints, - const MergeConflictBehavior merge_conflict_behavior, - const bool allow_different_local_lines) -{ - (void)allow_different_local_lines; - Assert(allow_different_local_lines || - local_lines == other_constraints.local_lines, - ExcMessage( - "local_lines for this and the other objects are not the same " - "although allow_different_local_lines is false.")); - - // store the previous state with respect to sorting - const bool object_was_sorted = sorted; - sorted = false; - - // first action is to fold into the present object possible constraints - // in the second object. we don't strictly need to do this any more since - // the AffineConstraints container has learned to deal with chains of - // constraints in the close() function, but we have traditionally done - // this and it's not overly hard to do. - // - // for this, loop over all constraints and replace the constraint lines - // with a new one where constraints are replaced if necessary. - typename ConstraintLine::Entries tmp; - for (ConstraintLine &line : lines) - { - tmp.clear(); - for (const std::pair &entry : line.entries) - { - // if the present dof is not stored, or not constrained, or if we - // won't take the constraint from the other object, then simply copy - // it over - if ((other_constraints.local_lines.size() != 0. && - other_constraints.local_lines.is_element(entry.first) == - false) || - other_constraints.is_constrained(entry.first) == false || - ((merge_conflict_behavior != right_object_wins) && - other_constraints.is_constrained(entry.first) && - this->is_constrained(entry.first))) - tmp.push_back(entry); - else - // otherwise resolve further constraints by replacing the old - // entry by a sequence of new entries taken from the other - // object, but with multiplied weights - { - const typename ConstraintLine::Entries *other_entries = - other_constraints.get_constraint_entries(entry.first); - Assert(other_entries != nullptr, ExcInternalError()); - - const number weight = entry.second; - - for (const std::pair &other_entry : - *other_entries) - tmp.emplace_back(other_entry.first, - other_entry.second * weight); - - line.inhomogeneity += - other_constraints.get_inhomogeneity(entry.first) * weight; - } - } - // finally exchange old and newly resolved line - line.entries.swap(tmp); - } - - if (local_lines.size() != 0) - local_lines.add_indices(other_constraints.local_lines); - - { - // do not bother to resize the lines cache exactly since it is pretty - // cheap to adjust it along the way. - std::fill(lines_cache.begin(), - lines_cache.end(), - numbers::invalid_size_type); - - // reset lines_cache for our own constraints - size_type index = 0; - for (const ConstraintLine &line : lines) - { - const size_type local_line_no = calculate_line_index(line.index); - if (local_line_no >= lines_cache.size()) - lines_cache.resize(local_line_no + 1, numbers::invalid_size_type); - lines_cache[local_line_no] = index++; - } - - // Add other_constraints to lines cache and our list of constraints - for (const ConstraintLine &line : other_constraints.lines) - { - const size_type local_line_no = calculate_line_index(line.index); - if (local_line_no >= lines_cache.size()) - { - lines_cache.resize(local_line_no + 1, numbers::invalid_size_type); - lines.push_back(line); - lines_cache[local_line_no] = index++; - } - else if (lines_cache[local_line_no] == numbers::invalid_size_type) - { - // there are no constraints for that line yet - lines.push_back(line); - AssertIndexRange(local_line_no, lines_cache.size()); - lines_cache[local_line_no] = index++; - } - else - { - // we already store that line - switch (merge_conflict_behavior) - { - case no_conflicts_allowed: - AssertThrow(false, - ExcDoFIsConstrainedFromBothObjects(line.index)); - break; - - case left_object_wins: - // ignore this constraint - break; - - case right_object_wins: - AssertIndexRange(local_line_no, lines_cache.size()); - lines[lines_cache[local_line_no]] = line; - break; - - default: - Assert(false, ExcNotImplemented()); - } - } - } - - // check that we set the pointers correctly - for (size_type i = 0; i < lines_cache.size(); ++i) - if (lines_cache[i] != numbers::invalid_size_type) - Assert(i == calculate_line_index(lines[lines_cache[i]].index), - ExcInternalError()); - } - - // if the object was sorted before, then make sure it is so afterward as - // well. otherwise leave everything in the unsorted state - if (object_was_sorted == true) - close(); -} - - template void -- 2.39.5