Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

AffineConstraints::merge(): allow mixed numbers #15753

Merged
merged 1 commit into from
Jul 22, 2023
Merged

AffineConstraints::merge(): allow mixed numbers #15753

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
157 changes: 156 additions & 1 deletion include/deal.II/lac/affine_constraints.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -821,9 +821,10 @@ class AffineConstraints : public Subscriptor
* Merging a AffineConstraints that is initialized with an IndexSet
* and one that is not initialized with an IndexSet is not yet implemented.
*/
template <typename other_number>
void
merge(
const AffineConstraints & other_constraints,
const AffineConstraints<other_number> &other_constraints,
const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
const bool allow_different_local_lines = false);

Expand Down Expand Up @@ -2468,6 +2469,160 @@ AffineConstraints<number>::copy_from(
}


template <typename number>
template <typename other_number>
void
AffineConstraints<number>::merge(
const AffineConstraints<other_number> &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<size_type, number> &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<size_type, number> &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 <typename number>
template <typename MatrixType>
Expand Down
143 changes: 0 additions & 143 deletions include/deal.II/lac/affine_constraints.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1014,149 +1014,6 @@ AffineConstraints<number>::is_closed(const MPI_Comm comm) const
}


template <typename number>
void
AffineConstraints<number>::merge(
const AffineConstraints<number> &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<size_type, number> &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<size_type, number> &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 <typename number>
void
Expand Down