Skip to content

Commit

Permalink
Merge pull request #15753 from peterrum/af_merge_number
Browse files Browse the repository at this point in the history 
AffineConstraints::merge(): allow mixed numbers
  • Loading branch information
@tamiko
tamiko committed Jul 22, 2023
2 parents 48a8713 + 0ab1b65 commit 74b8db5
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 144 deletions.
157 changes: 156 additions & 1 deletion include/deal.II/lac/affine_constraints.h
View file
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
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

0 comments on commit 74b8db5

Please sign in to comment.