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

src/base/dof_map.C

@@ -3139,6 +3139,70 @@ 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 std::vector<std::string> & vars,
+                                const System & sys,
+                                const FEType & fe_type,
+                                const std::set<subdomain_id_type> * 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 (fe_type.family == SCALAR || (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;
+
+  if (more_than_one_bit_set)
+    libmesh_error_msg(make_error_message());
+}
+}
+#endif
+
 unsigned int DofMap::add_variable(System & sys,
                                   std::string_view var,
                                   const FEType & type,
@@ -3153,6 +3217,11 @@ unsigned int DofMap::add_variable(System & sys,
   if (active_subdomains)
     libmesh_assert(this->comm().verify(active_subdomains->size()));
 
+#ifdef DEBUG
+  discontinuity_sanity_check(
+      std::vector<std::string>{std::string(var)}, 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()))
@@ -3255,13 +3324,18 @@ 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)));
 
   if (active_subdomains)
     libmesh_assert(this->comm().verify(active_subdomains->size()));
 
+#ifdef DEBUG
+  discontinuity_sanity_check(vars, 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();
   }