From c5edf23d7752fe90a8c576e072a646325850ab75 Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Thu, 23 Aug 2018 14:42:43 -0500 Subject: [PATCH 1/3] Add assertion in add_constraint_row Finding cycles of constraints is easier when there's only one in the cycle. This turned out to be a red herring for the bug I'm investigating but it might be a good idea later. --- src/base/dof_map_constraints.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/base/dof_map_constraints.C b/src/base/dof_map_constraints.C index c174e2cef0d..b3b2cb8defa 100644 --- a/src/base/dof_map_constraints.C +++ b/src/base/dof_map_constraints.C @@ -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 @@ -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 it = _dof_constraints.emplace(dof_number, constraint_row); From 94ccfe1960dcabda54bdea6da7a8588846a02d2f Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Mon, 1 Jun 2020 16:26:31 -0500 Subject: [PATCH 2/3] Add more const, assertions to constraints code This helped pin down a bug, and has been cleaned up enough we ought to be able to leave it in for good. --- src/base/dof_map_constraints.C | 22 +++++++++++++++++----- src/fe/fe_lagrange.C | 4 ++++ 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/base/dof_map_constraints.C b/src/base/dof_map_constraints.C index b3b2cb8defa..fdba3fb92de 100644 --- a/src/base/dof_map_constraints.C +++ b/src/base/dof_map_constraints.C @@ -3320,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; } @@ -3993,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; @@ -4044,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 = @@ -4109,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)) @@ -4289,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; @@ -4497,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) @@ -4605,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; } diff --git a/src/fe/fe_lagrange.C b/src/fe/fe_lagrange.C index 88879ce0b35..4e9239e4736 100644 --- a/src/fe/fe_lagrange.C +++ b/src/fe/fe_lagrange.C @@ -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 From b0a610b0320eb04476800299320021abbd127fe5 Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Mon, 1 Jun 2020 15:58:10 -0500 Subject: [PATCH 3/3] Actually fix the constraints bug! It turns out that, without ghosting of all point neighbors, we can't trust find_point_neighbors to return a complete set for partition boundary points on a distributed mesh, which means we can't necessarily look at a complete set of element ids to disambiguate which side should constrain which. Looking at the specific points themselves to disambiguate constraints fixes the problem for me, at least on the case where it manifested, systems_of_equations_ex9.C with 21 processors. --- src/fe/fe_base.C | 34 ++++++++++++---------------------- 1 file changed, 12 insertions(+), 22 deletions(-) diff --git a/src/fe/fe_base.C b/src/fe/fe_base.C index c505edcbbe0..ffaf530d87f 100644 --- a/src/fe/fe_base.C +++ b/src/fe/fe_base.C @@ -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; @@ -2331,6 +2326,14 @@ 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 @@ -2338,7 +2341,7 @@ compute_periodic_constraints (DofConstraints & constraints, // 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); @@ -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;