Skip to content

Fix distributed periodic constraints bug #2590

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

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Fix distributed periodic constraints bug #2590

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c5edf23
Add assertion in add_constraint_row
roystgnr Aug 23, 2018
94ccfe1
Add more const, assertions to constraints code
roystgnr Jun 1, 2020
b0a610b
Actually fix the constraints bug!
roystgnr Jun 1, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 21 additions & 5 deletions src/base/dof_map_constraints.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1703,6 +1703,7 @@ void DofMap::add_constraint_row (const dof_id_type dof_number,
if (this->is_constrained_dof(dof_number))
libmesh_error_msg("ERROR: DOF " << dof_number << " was already constrained!");

// We don't allow nonsensical constraints
libmesh_assert_less(dof_number, this->n_dofs());

// There is an implied "1" on the diagonal of the constraint row, and the user
Expand All @@ -1716,6 +1717,9 @@ void DofMap::add_constraint_row (const dof_id_type dof_number,
libmesh_assert_less(pr.first, this->n_dofs());
#endif

// We don't allow constraints-against-ourselves
libmesh_assert(!constraint_row.count(dof_number));

// We don't get insert_or_assign until C++17 so we make do.
std::pair<DofConstraints::iterator, bool> it =
_dof_constraints.emplace(dof_number, constraint_row);
Expand Down Expand Up @@ -3316,6 +3320,10 @@ void DofMap::allgather_recursive_constraints(MeshBase & mesh)
for (auto & kv : pushed_keys_vals_to_me[i])
{
libmesh_assert_less(kv.first, this->n_dofs());

// Don't create an obviously-cyclic constraint
libmesh_assert_not_equal_to(constrained, kv.first);

row[kv.first] = kv.second;
}

Expand Down Expand Up @@ -3989,8 +3997,8 @@ void DofMap::scatter_constraints(MeshBase & mesh)
if (constrained_proc_id != this->processor_id())
continue;

DofConstraintRow & row = i.second;
for (auto & j : row)
const DofConstraintRow & row = i.second;
for (const auto & j : row)
{
const dof_id_type constraining = j.first;

Expand Down Expand Up @@ -4040,7 +4048,7 @@ void DofMap::scatter_constraints(MeshBase & mesh)
it != pid_ids.end(); ++push_i, ++it)
{
const dof_id_type constrained = *it;
DofConstraintRow & row = _dof_constraints[constrained];
const DofConstraintRow & row = _dof_constraints[constrained];
keys_vals[push_i].assign(row.begin(), row.end());

DofConstraintValueMap::const_iterator rhsit =
Expand Down Expand Up @@ -4105,6 +4113,10 @@ void DofMap::scatter_constraints(MeshBase & mesh)
for (auto & key_val : keys_vals[i])
{
libmesh_assert_less(key_val.first, this->n_dofs());

// Don't create an obviously-cyclic constraint
libmesh_assert_not_equal_to(constrained, key_val.first);

row[key_val.first] = key_val.second;
}
if (ids_rhss[i].second != Number(0))
Expand Down Expand Up @@ -4285,7 +4297,7 @@ void DofMap::scatter_constraints(MeshBase & mesh)
for (auto & i : _dof_constraints)
{
const dof_id_type constrained = i.first;
DofConstraintRow & row = i.second;
const DofConstraintRow & row = i.second;
for (const auto & j : row)
{
const dof_id_type constraining = j.first;
Expand Down Expand Up @@ -4493,7 +4505,7 @@ void DofMap::gather_constraints (MeshBase & /*mesh*/,
dof_id_type constrained = ids[i];
if (_dof_constraints.count(constrained))
{
DofConstraintRow & row = _dof_constraints[constrained];
const DofConstraintRow & row = _dof_constraints[constrained];
std::size_t row_size = row.size();
data[i].reserve(row_size);
for (const auto & j : row)
Expand Down Expand Up @@ -4601,6 +4613,10 @@ void DofMap::gather_constraints (MeshBase & /*mesh*/,
for (auto & pair : data[i])
{
libmesh_assert_less(pair.first, this->n_dofs());

// Don't create an obviously-cyclic constraint
libmesh_assert_not_equal_to(constrained, pair.first);

row[pair.first] = pair.second;
}

Expand Down
34 changes: 12 additions & 22 deletions src/fe/fe_base.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2194,15 +2194,10 @@ compute_periodic_constraints (DofConstraints & constraints,
// other way around
if ((primary_neigh->level() > primary_elem->level()) ||

// For equal-level elements, the one with
// higher id gets constrained in terms of
// the one with lower id
(primary_neigh->level() == primary_elem->level() &&
primary_neigh->id() > primary_elem->id()) ||

// On a one-element-thick mesh, we compare
// For equal-level elements (or on a
// one-element-thick mesh), we compare
// points to see what side gets constrained
(primary_neigh == primary_elem &&
(primary_neigh->level() == primary_elem->level() &&
(neigh_pt > primary_pt)))
continue;

Expand Down Expand Up @@ -2331,14 +2326,22 @@ compute_periodic_constraints (DofConstraints & constraints,
main_pt2 = neigh_pt2;
}

// Otherwise:
// Finer elements will get constrained in
// terms of coarser ones, not the other way
// around
if (primary_neigh->level() > primary_elem->level())
continue;

// For equal-level elements, or if
// If we have a one-element thick mesh,
// we'll need to sort our points to get a
// consistent ordering rule
//
// Use >= in this test to make sure that,
// for angular constraints, no node gets
// constrained to itself.
if (primary_neigh == primary_elem)
if (primary_neigh->level() == primary_elem->level())
{
if (primary_pt1 > primary_pt2)
std::swap(primary_pt1, primary_pt2);
Expand All @@ -2349,19 +2352,6 @@ compute_periodic_constraints (DofConstraints & constraints,
continue;
}

// Otherwise:
// Finer elements will get constrained in
// terms of coarser ones, not the other way
// around
if ((primary_neigh->level() > primary_elem->level()) ||

// For equal-level elements, the one with
// higher id gets constrained in terms of
// the one with lower id
(primary_neigh->level() == primary_elem->level() &&
primary_neigh->id() > primary_elem->id()))
continue;

primary_elem = primary_neigh;
primary_pt1 = neigh_pt1;
primary_pt2 = neigh_pt2;
Expand Down
4 changes: 4 additions & 0 deletions src/fe/fe_lagrange.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -795,6 +795,10 @@ void lagrange_compute_constraints (DofConstraints & constraints,
if ((std::abs(their_dof_value) > 1.e-5) &&
(std::abs(their_dof_value) < .999))
{
// That should be enough to make sure it isn't a
// self-constraint
libmesh_assert_not_equal_to(my_dof_g, their_dof_g);

constraint_row->emplace(their_dof_g, their_dof_value);
}
#ifdef DEBUG
Expand Down