[WIP] Debugging check that variables of discontinuous family live on just one elem->dim()

src/base/dof_map.C

@@ -3139,6 +3139,47 @@ void DofMap::reinit_static_condensation()
     _sc->reinit();
 }
 
+#ifdef DEBUG
+namespace
+{
+unsigned char determine_elem_dims(const MeshBase & mesh,
+                                  const std::set<subdomain_id_type> & subdomains)
+{
+  unsigned char elem_dims = 0;
+  for (const auto * const elem : mesh.active_subdomain_set_element_ptr_range(subdomains))
+    elem_dims |= static_cast<unsigned char>(1u << cast_int<unsigned char>(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<subdomain_id_type> * const active_subdomains)
+{
+  const auto continuity = FEInterface::get_continuity(fe_type);
+  if (continuity != DISCONTINUOUS && continuity != SIDE_DISCONTINUOUS)
+    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 & var_subdomains = (active_subdomains && !active_subdomains->empty())
+                                    ? *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
+  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.");
+}
+}
+#endif
+
 unsigned int DofMap::add_variable(System & sys,
                                   std::string_view var,
                                   const FEType & type,
@@ -3153,6 +3194,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 +3307,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)

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<System> ("SimpleSystem");
-    sys.add_variable("u", CONSTANT, MONOMIAL);
+    const std::set<subdomain_id_type> u_subs = {0};
+    sys.add_variable("u", CONSTANT, MONOMIAL, &u_subs);
     es.init();
     es.reinit();
   }