diff --git a/src/base/dof_map_constraints.C b/src/base/dof_map_constraints.C index c174e2cef0d..fdba3fb92de 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); @@ -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; } @@ -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; @@ -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 = @@ -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)) @@ -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; @@ -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) @@ -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; } 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; 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