From 1e949f539f25af79bcaef284afe229eceda56b30 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 5 Dec 2025 15:56:23 -0700 Subject: [PATCH 1/6] Debugging check that variables of discontinuous family live on just one elem->dim() Refs idaholab/moose#32003 --- include/mesh/mesh_base.h | 15 ++++++++++ src/base/dof_map.C | 39 ++++++++++++++++++++++++ src/mesh/mesh_base.C | 65 ++++++++++++++++++++++++++++++++++++++-- 3 files changed, 117 insertions(+), 2 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index b1c40a44ee8..dcd9a1652fa 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -1844,6 +1844,14 @@ class MeshBase : public ParallelObject const boundary_id_type b2); #endif +#ifdef DEBUG + /** + * Get the map from subdomain ID to the set of element dimensions contained within + * that subdomain + */ + const std::unordered_map> & + get_subdomain_to_elem_dims_map() const { return _subdomain_id_to_elem_dims; } +#endif protected: @@ -2154,6 +2162,13 @@ class MeshBase : public ParallelObject */ Real _point_locator_close_to_point_tol; +#ifdef DEBUG + /** + * Map from subdomain ID to the set of element dimensions present in that subdomain + */ + std::unordered_map> _subdomain_id_to_elem_dims; +#endif + /** * The partitioner class is a friend so that it can set * the number of partitions. diff --git a/src/base/dof_map.C b/src/base/dof_map.C index a988d15bd08..0f1955c5866 100644 --- a/src/base/dof_map.C +++ b/src/base/dof_map.C @@ -3139,6 +3139,37 @@ void DofMap::reinit_static_condensation() _sc->reinit(); } +#ifdef DEBUG +namespace +{ +void discontinuity_sanity_check(const System & sys, + const FEType & fe_type, + const std::set * const active_subdomains) +{ + const auto continuity = FEInterface::get_continuity(fe_type); + if (continuity != DISCONTINUOUS && continuity != SIDE_DISCONTINUOUS) + return; + + const auto & mesh = sys.get_mesh(); + const auto & sub_to_elem_dims = mesh.get_subdomain_to_elem_dims_map(); + const auto & var_subdomains = (active_subdomains && !active_subdomains->empty()) + ? *active_subdomains + : mesh.get_mesh_subdomains(); + std::unordered_set var_elem_dims; + for (const auto sub_id : var_subdomains) + { + const auto & sub_elem_dims = libmesh_map_find(sub_to_elem_dims, sub_id); + for (const auto dim : sub_elem_dims) + var_elem_dims.insert(dim); + } + libmesh_assert_msg(var_elem_dims.size() <= 1, + "Discontinuous finite element families cannot live on elements with " + "different dimensions because this violates the idea that the family is " + "discontinuous (undefined) at the interface between elements"); +} +} +#endif + unsigned int DofMap::add_variable(System & sys, std::string_view var, const FEType & type, @@ -3153,6 +3184,10 @@ unsigned int DofMap::add_variable(System & sys, if (active_subdomains) libmesh_assert(this->comm().verify(active_subdomains->size())); +#ifdef DEBUG + discontinuity_sanity_check(sys, type, active_subdomains); +#endif + // Make sure the variable isn't there already // or if it is, that it's the type we want for (auto v : make_range(this->n_vars())) @@ -3262,6 +3297,10 @@ unsigned int DofMap::add_variables(System & sys, if (active_subdomains) libmesh_assert(this->comm().verify(active_subdomains->size())); +#ifdef DEBUG + discontinuity_sanity_check(sys, type, active_subdomains); +#endif + // Make sure the variable isn't there already // or if it is, that it's the type we want for (auto ovar : vars) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 535117af6b5..52f4d5744de 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -41,6 +41,7 @@ #include "libmesh/enum_to_string.h" #include "libmesh/point_locator_nanoflann.h" #include "libmesh/elem_side_builder.h" +#include "libmesh/utility.h" // C++ includes #include // for std::min @@ -1824,13 +1825,18 @@ void MeshBase::cache_elem_data() for (const auto & elem : this->active_element_ptr_range()) { - _elem_dims.insert(cast_int(elem->dim())); + const auto elem_dim = cast_int(elem->dim()); + _elem_dims.insert(elem_dim); _elem_default_orders.insert(elem->default_order()); - _mesh_subdomains.insert(elem->subdomain_id()); + const auto sub_id = elem->subdomain_id(); + _mesh_subdomains.insert(sub_id); _supported_nodal_order = static_cast (std::min(static_cast(_supported_nodal_order), static_cast(elem->supported_nodal_order()))); +#ifdef DEBUG + _subdomain_id_to_elem_dims[sub_id].insert(elem_dim); +#endif } if (!this->is_serial()) @@ -1841,6 +1847,61 @@ void MeshBase::cache_elem_data() this->comm().set_union(_elem_default_orders); this->comm().min(_supported_nodal_order); this->comm().set_union(_mesh_subdomains); + +#ifdef DEBUG +#ifdef LIBMESH_HAVE_MPI + + using Mask = unsigned char; + + // Prepare sparse subdomain id to dense subdomain "id" + const auto n_sub = cast_int(_mesh_subdomains.size()); + std::unordered_map sub_to_dense_idx; + sub_to_dense_idx.reserve(n_sub); + std::vector dense_idx_to_sub; + dense_idx_to_sub.reserve(n_sub); + { + subdomain_id_type dense_idx = 0; + for (const auto sub_id : _mesh_subdomains) + { + sub_to_dense_idx.emplace(sub_id, dense_idx++); + dense_idx_to_sub.push_back(sub_id); + } + } + + // Pack local dimension sets into Mask + std::vector masks(n_sub, Mask(0)); + for (const auto & [sub_id, elem_dims] : _subdomain_id_to_elem_dims) + { + const auto dense_idx = libmesh_map_find(sub_to_dense_idx, sub_id); + Mask m = 0; + for (const auto dim : elem_dims) + m |= Mask(1u << dim); + + masks[dense_idx] = m; + } + _subdomain_id_to_elem_dims.clear(); + + // Communicate the mask data + MPI_Allreduce( + MPI_IN_PLACE, masks.data(), n_sub, MPI_UNSIGNED_CHAR, MPI_BOR, _communicator.get()); + + constexpr unsigned char max_elem_dim = LIBMESH_DIM; + + for (const auto dense_idx : make_range(n_sub)) + { + const auto mask = masks[dense_idx]; + if (mask == 0) + // No elements for this subdomain + continue; + + const auto sub_id = dense_idx_to_sub[dense_idx]; + auto & dim_set = _subdomain_id_to_elem_dims[sub_id]; + for (const auto dim : make_range(max_elem_dim)) + if (mask & Mask(1u << dim)) + dim_set.insert(dim); + } +#endif // LIBMESH_HAVE_MPI +#endif // DEBUG } // If the largest element dimension found is larger than the current From cdcae75ace013ef6e4387511111569cf1d673fb8 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 5 Dec 2025 17:27:42 -0700 Subject: [PATCH 2/6] Early return in sanity check if mesh isn't prepared --- src/base/dof_map.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/base/dof_map.C b/src/base/dof_map.C index 0f1955c5866..70fbb361284 100644 --- a/src/base/dof_map.C +++ b/src/base/dof_map.C @@ -3151,6 +3151,10 @@ void discontinuity_sanity_check(const System & sys, return; const auto & mesh = sys.get_mesh(); + // This check won't be meaninful if the mesh isn't prepared + if (!mesh.is_prepared()) + return; + const auto & sub_to_elem_dims = mesh.get_subdomain_to_elem_dims_map(); const auto & var_subdomains = (active_subdomains && !active_subdomains->empty()) ? *active_subdomains From ddbbb0cb8b5aa5aa3f82f672f344da086c88852d Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 5 Dec 2025 17:28:04 -0700 Subject: [PATCH 3/6] Associate var in NodeElem test with 1D elements --- tests/systems/equation_systems_test.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/systems/equation_systems_test.C b/tests/systems/equation_systems_test.C index c4ac2d4a191..79894ef075a 100644 --- a/tests/systems/equation_systems_test.C +++ b/tests/systems/equation_systems_test.C @@ -182,11 +182,13 @@ public: MeshTools::Generation::build_line (mesh, 10, 0., 1., EDGE2); auto node_elem = mesh.add_elem(Elem::build(NODEELEM)); node_elem->set_node(0, mesh.node_ptr(0)); + node_elem->subdomain_id() = 1; mesh.prepare_for_use(); EquationSystems es(mesh); System &sys = es.add_system ("SimpleSystem"); - sys.add_variable("u", CONSTANT, MONOMIAL); + const std::set u_subs = {0}; + sys.add_variable("u", CONSTANT, MONOMIAL, &u_subs); es.init(); es.reinit(); } From 51dee792289d1ef72b556d29e281c52896a68d54 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 5 Dec 2025 17:30:32 -0700 Subject: [PATCH 4/6] Update sanity check message --- src/base/dof_map.C | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/base/dof_map.C b/src/base/dof_map.C index 70fbb361284..1ee5623ac53 100644 --- a/src/base/dof_map.C +++ b/src/base/dof_map.C @@ -3166,10 +3166,11 @@ void discontinuity_sanity_check(const System & sys, for (const auto dim : sub_elem_dims) var_elem_dims.insert(dim); } - libmesh_assert_msg(var_elem_dims.size() <= 1, - "Discontinuous finite element families cannot live on elements with " - "different dimensions because this violates the idea that the family is " - "discontinuous (undefined) at the interface between elements"); + libmesh_assert_msg( + var_elem_dims.size() <= 1, + "Discontinuous finite element families are typically associated with discontinuous Galerkin " + "methods. If degrees of freedom are associated with different values of element dimension, " + "then generally this will result in a singular system after application of the DG method."); } } #endif From 315f6eb2db40103fb713f316d1357e9c2cb649cf Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 5 Dec 2025 21:46:53 -0700 Subject: [PATCH 5/6] Move code to dof_map and compute on-the-fly --- include/mesh/mesh_base.h | 15 ---------- src/base/dof_map.C | 23 ++++++++------ src/mesh/mesh_base.C | 65 ++-------------------------------------- 3 files changed, 16 insertions(+), 87 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index dcd9a1652fa..b1c40a44ee8 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -1844,14 +1844,6 @@ class MeshBase : public ParallelObject const boundary_id_type b2); #endif -#ifdef DEBUG - /** - * Get the map from subdomain ID to the set of element dimensions contained within - * that subdomain - */ - const std::unordered_map> & - get_subdomain_to_elem_dims_map() const { return _subdomain_id_to_elem_dims; } -#endif protected: @@ -2162,13 +2154,6 @@ class MeshBase : public ParallelObject */ Real _point_locator_close_to_point_tol; -#ifdef DEBUG - /** - * Map from subdomain ID to the set of element dimensions present in that subdomain - */ - std::unordered_map> _subdomain_id_to_elem_dims; -#endif - /** * The partitioner class is a friend so that it can set * the number of partitions. diff --git a/src/base/dof_map.C b/src/base/dof_map.C index 1ee5623ac53..4d7590aa20f 100644 --- a/src/base/dof_map.C +++ b/src/base/dof_map.C @@ -3142,6 +3142,16 @@ void DofMap::reinit_static_condensation() #ifdef DEBUG namespace { +unsigned char determine_elem_dims(const MeshBase & mesh, + const std::set & subdomains) +{ + unsigned char elem_dims = 0; + for (const auto * const elem : mesh.active_subdomain_set_element_ptr_range(subdomains)) + elem_dims |= static_cast(1u << cast_int(elem->dim())); + mesh.comm().bitwise_or(elem_dims); + return elem_dims; +} + void discontinuity_sanity_check(const System & sys, const FEType & fe_type, const std::set * const active_subdomains) @@ -3155,19 +3165,14 @@ void discontinuity_sanity_check(const System & sys, if (!mesh.is_prepared()) return; - const auto & sub_to_elem_dims = mesh.get_subdomain_to_elem_dims_map(); const auto & var_subdomains = (active_subdomains && !active_subdomains->empty()) ? *active_subdomains : mesh.get_mesh_subdomains(); - std::unordered_set var_elem_dims; - for (const auto sub_id : var_subdomains) - { - const auto & sub_elem_dims = libmesh_map_find(sub_to_elem_dims, sub_id); - for (const auto dim : sub_elem_dims) - var_elem_dims.insert(dim); - } + const auto var_elem_dims = determine_elem_dims(mesh, var_subdomains); + // Power of two trick. Note that unsigned char automatically promoted to int for arithmetic and bitwise operations + const bool more_than_one_bit_set = (var_elem_dims & (var_elem_dims - 1)) != 0; libmesh_assert_msg( - var_elem_dims.size() <= 1, + !more_than_one_bit_set, "Discontinuous finite element families are typically associated with discontinuous Galerkin " "methods. If degrees of freedom are associated with different values of element dimension, " "then generally this will result in a singular system after application of the DG method."); diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 52f4d5744de..535117af6b5 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -41,7 +41,6 @@ #include "libmesh/enum_to_string.h" #include "libmesh/point_locator_nanoflann.h" #include "libmesh/elem_side_builder.h" -#include "libmesh/utility.h" // C++ includes #include // for std::min @@ -1825,18 +1824,13 @@ void MeshBase::cache_elem_data() for (const auto & elem : this->active_element_ptr_range()) { - const auto elem_dim = cast_int(elem->dim()); - _elem_dims.insert(elem_dim); + _elem_dims.insert(cast_int(elem->dim())); _elem_default_orders.insert(elem->default_order()); - const auto sub_id = elem->subdomain_id(); - _mesh_subdomains.insert(sub_id); + _mesh_subdomains.insert(elem->subdomain_id()); _supported_nodal_order = static_cast (std::min(static_cast(_supported_nodal_order), static_cast(elem->supported_nodal_order()))); -#ifdef DEBUG - _subdomain_id_to_elem_dims[sub_id].insert(elem_dim); -#endif } if (!this->is_serial()) @@ -1847,61 +1841,6 @@ void MeshBase::cache_elem_data() this->comm().set_union(_elem_default_orders); this->comm().min(_supported_nodal_order); this->comm().set_union(_mesh_subdomains); - -#ifdef DEBUG -#ifdef LIBMESH_HAVE_MPI - - using Mask = unsigned char; - - // Prepare sparse subdomain id to dense subdomain "id" - const auto n_sub = cast_int(_mesh_subdomains.size()); - std::unordered_map sub_to_dense_idx; - sub_to_dense_idx.reserve(n_sub); - std::vector dense_idx_to_sub; - dense_idx_to_sub.reserve(n_sub); - { - subdomain_id_type dense_idx = 0; - for (const auto sub_id : _mesh_subdomains) - { - sub_to_dense_idx.emplace(sub_id, dense_idx++); - dense_idx_to_sub.push_back(sub_id); - } - } - - // Pack local dimension sets into Mask - std::vector masks(n_sub, Mask(0)); - for (const auto & [sub_id, elem_dims] : _subdomain_id_to_elem_dims) - { - const auto dense_idx = libmesh_map_find(sub_to_dense_idx, sub_id); - Mask m = 0; - for (const auto dim : elem_dims) - m |= Mask(1u << dim); - - masks[dense_idx] = m; - } - _subdomain_id_to_elem_dims.clear(); - - // Communicate the mask data - MPI_Allreduce( - MPI_IN_PLACE, masks.data(), n_sub, MPI_UNSIGNED_CHAR, MPI_BOR, _communicator.get()); - - constexpr unsigned char max_elem_dim = LIBMESH_DIM; - - for (const auto dense_idx : make_range(n_sub)) - { - const auto mask = masks[dense_idx]; - if (mask == 0) - // No elements for this subdomain - continue; - - const auto sub_id = dense_idx_to_sub[dense_idx]; - auto & dim_set = _subdomain_id_to_elem_dims[sub_id]; - for (const auto dim : make_range(max_elem_dim)) - if (mask & Mask(1u << dim)) - dim_set.insert(dim); - } -#endif // LIBMESH_HAVE_MPI -#endif // DEBUG } // If the largest element dimension found is larger than the current From 138e739fc3dacc67d9c33a136f31dd4fbb0e4d63 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 9 Dec 2025 10:50:04 -0700 Subject: [PATCH 6/6] Add variable name(s) in discontinuity error message --- src/base/dof_map.C | 45 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/src/base/dof_map.C b/src/base/dof_map.C index 4d7590aa20f..f213f952c20 100644 --- a/src/base/dof_map.C +++ b/src/base/dof_map.C @@ -3152,12 +3152,36 @@ unsigned char determine_elem_dims(const MeshBase & mesh, return elem_dims; } -void discontinuity_sanity_check(const System & sys, +void discontinuity_sanity_check(const std::vector & vars, + const System & sys, const FEType & fe_type, const std::set * const active_subdomains) { + auto make_error_message = [&vars]() { + std::ostringstream oss; + for (const auto i : index_range(vars)) + { + oss << "'" << vars[i] << "'"; + if (i + 1 < vars.size()) + oss << ", "; + } + const std::string joined = oss.str(); + + const bool plural = vars.size() > 1; + const std::string var_word = plural ? "variables" : "variable"; + const std::string is_are = plural ? "are" : "is"; + + // Build the final message + return "Discontinuous finite element families are typically associated " + "with discontinuous Galerkin methods. If degrees of freedom are " + "associated with different values of element dimension, then " + "generally this will result in a singular system after application " + "of the DG method. The discontinuous " + + var_word + " " + joined + " " + is_are + " defined on multiple element dimensions."; + }; + const auto continuity = FEInterface::get_continuity(fe_type); - if (continuity != DISCONTINUOUS && continuity != SIDE_DISCONTINUOUS) + if (fe_type.family == SCALAR || (continuity != DISCONTINUOUS && continuity != SIDE_DISCONTINUOUS)) return; const auto & mesh = sys.get_mesh(); @@ -3169,13 +3193,12 @@ void discontinuity_sanity_check(const System & sys, ? *active_subdomains : mesh.get_mesh_subdomains(); const auto var_elem_dims = determine_elem_dims(mesh, var_subdomains); - // Power of two trick. Note that unsigned char automatically promoted to int for arithmetic and bitwise operations + // Power of two trick. Note that unsigned char automatically promoted to int for arithmetic and + // bitwise operations const bool more_than_one_bit_set = (var_elem_dims & (var_elem_dims - 1)) != 0; - libmesh_assert_msg( - !more_than_one_bit_set, - "Discontinuous finite element families are typically associated with discontinuous Galerkin " - "methods. If degrees of freedom are associated with different values of element dimension, " - "then generally this will result in a singular system after application of the DG method."); + + if (more_than_one_bit_set) + libmesh_error_msg(make_error_message()); } } #endif @@ -3195,7 +3218,8 @@ unsigned int DofMap::add_variable(System & sys, libmesh_assert(this->comm().verify(active_subdomains->size())); #ifdef DEBUG - discontinuity_sanity_check(sys, type, active_subdomains); + discontinuity_sanity_check( + std::vector{std::string(var)}, sys, type, active_subdomains); #endif // Make sure the variable isn't there already @@ -3300,6 +3324,7 @@ unsigned int DofMap::add_variables(System & sys, libmesh_assert(!sys.is_initialized()); + libmesh_assert(vars.size()); libmesh_assert(this->comm().verify(vars.size())); libmesh_assert(this->comm().verify(type)); libmesh_assert(this->comm().verify((active_subdomains == nullptr))); @@ -3308,7 +3333,7 @@ unsigned int DofMap::add_variables(System & sys, libmesh_assert(this->comm().verify(active_subdomains->size())); #ifdef DEBUG - discontinuity_sanity_check(sys, type, active_subdomains); + discontinuity_sanity_check(vars, sys, type, active_subdomains); #endif // Make sure the variable isn't there already