Move code to dof_map and compute on-the-fly

Author: lindsayad

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<subdomain_id_type, std::unordered_set<unsigned char>> &
-  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_type, std::unordered_set<unsigned char>> _subdomain_id_to_elem_dims;
-#endif
-
   /**
    * The partitioner class is a friend so that it can set
    * the number of partitions.

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<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)
@@ -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<unsigned char> 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.");

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 <algorithm> // 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<unsigned char>(elem->dim());
-    _elem_dims.insert(elem_dim);
+    _elem_dims.insert(cast_int<unsigned char>(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<Order>
         (std::min(static_cast<int>(_supported_nodal_order),
                   static_cast<int>(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<int>(_mesh_subdomains.size());
-    std::unordered_map<subdomain_id_type, subdomain_id_type> sub_to_dense_idx;
-    sub_to_dense_idx.reserve(n_sub);
-    std::vector<subdomain_id_type> 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<Mask> 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