split lower-d mesh blocks based on side types so that exodus output c…

…an work properly #29142

framework/include/loops/ThreadedElementLoopBase.h

@@ -255,8 +255,8 @@ ThreadedElementLoopBase<RangeType>::operator()(const RangeType & range, bool byp
 
         onElement(elem);
 
-        if (elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID ||
-            elem->subdomain_id() == Moose::BOUNDARY_SIDE_LOWERD_ID)
+        if (_mesh.interiorLowerDBlocks().count(elem->subdomain_id()) > 0 ||
+            _mesh.boundaryLowerDBlocks().count(elem->subdomain_id()) > 0)
         {
           postElement(elem);
           continue;

framework/include/mesh/MooseMesh.h

@@ -1365,6 +1365,15 @@ class MooseMesh : public MooseObject, public Restartable, public PerfGraphInterf
    */
   bool hasLowerD() const { return _has_lower_d; }
 
+  /**
+   * @return The set of lower-dimensional blocks for interior sides
+   */
+  const std::set<SubdomainID> & interiorLowerDBlocks() const { return _lower_d_interior_blocks; }
+  /**
+   * @return The set of lower-dimensional blocks for boundary sides
+   */
+  const std::set<SubdomainID> & boundaryLowerDBlocks() const { return _lower_d_boundary_blocks; }
+
 protected:
   /// Deprecated (DO NOT USE)
   std::vector<std::unique_ptr<GhostingFunctor>> _ghosting_functors;
@@ -1752,6 +1761,10 @@ class MooseMesh : public MooseObject, public Restartable, public PerfGraphInterf
   /// Holds a map from neighbor subomdain ids to the boundary ids that are attached to it
   std::unordered_map<SubdomainID, std::set<BoundaryID>> _neighbor_subdomain_boundary_ids;
 
+  /// Mesh blocks for interior lower-d elements in different types
+  std::set<SubdomainID> _lower_d_interior_blocks;
+  /// Mesh blocks for boundary lower-d elements in different types
+  std::set<SubdomainID> _lower_d_boundary_blocks;
   /// Holds a map from a high-order element side to its corresponding lower-d element
   std::unordered_map<std::pair<const Elem *, unsigned short int>, const Elem *>
       _higher_d_elem_side_to_lower_d_elem;

framework/include/utils/MooseTypes.h

@@ -659,8 +659,6 @@ namespace Moose
 {
 extern const processor_id_type INVALID_PROCESSOR_ID;
 extern const SubdomainID ANY_BLOCK_ID;
-extern const SubdomainID INTERNAL_SIDE_LOWERD_ID;
-extern const SubdomainID BOUNDARY_SIDE_LOWERD_ID;
 extern const SubdomainID INVALID_BLOCK_ID;
 extern const BoundaryID ANY_BOUNDARY_ID;
 extern const BoundaryID INVALID_BOUNDARY_ID;

framework/src/actions/SetupMeshCompleteAction.C

@@ -71,8 +71,7 @@ SetupMeshCompleteAction::act()
 
         if (_mesh->uniformRefineLevel())
         {
-          if (_mesh->meshSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) ||
-              _mesh->meshSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
+          if (!_mesh->interiorLowerDBlocks().empty() || !_mesh->boundaryLowerDBlocks().empty())
             mooseError("HFEM does not support mesh uniform refinement currently.");
 
           Adaptivity::uniformRefine(_mesh.get());

framework/src/bcs/ArrayHFEMDirichletBC.C

@@ -34,10 +34,14 @@ ArrayHFEMDirichletBC::ArrayHFEMDirichletBC(const InputParameters & parameters)
 
   if (_uhat_var)
   {
-    if (!_uhat_var->activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
-      paramError("uhat",
-                 "Must be defined on BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is added by "
-                 "Mesh/build_all_side_lowerd_mesh=true");
+    for (const auto & id : _uhat_var->activeSubdomains())
+      if (_mesh.boundaryLowerDBlocks().count(id) == 0)
+        mooseDocumentedError("moose",
+                             29151,
+                             "'uhat' must be defined on lower-dimensional boundary subdomain '" +
+                                 _mesh.getSubdomainName(id) +
+                                 "' that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe "
+                                 "check could be overly restrictive.");
 
     if (_uhat_var->count() != _count)
       paramError("uhat",

framework/src/bcs/ArrayLowerDIntegratedBC.C

@@ -35,15 +35,21 @@ ArrayLowerDIntegratedBC::ArrayLowerDIntegratedBC(const InputParameters & paramet
     _work_vector(_count)
 {
   const auto & lower_domains = _lowerd_var.activeSubdomains();
-  if (!lower_domains.count(Moose::BOUNDARY_SIDE_LOWERD_ID) && lower_domains.size() != 1)
-    paramError(
-        "lowerd_variable",
-        "Must be only defined on the subdomain BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is "
-        "added by Mesh/build_all_side_lowerd_mesh=true");
-
-  if (_var.activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
-    paramError("variable",
-               "Must not be defined on the subdomain BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain");
+  for (const auto & id : _mesh.boundaryLowerDBlocks())
+    if (lower_domains.count(id) == 0)
+      mooseDocumentedError(
+          "moose",
+          29151,
+          "'lowerd_variable' must be defined on the boundary lower-dimensional subdomain '" +
+              _mesh.getSubdomainName(id) +
+              "' that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe check could be overly "
+              "restrictive.");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.boundaryLowerDBlocks().count(id) > 0)
+      paramError("variable",
+                 "Must not be defined on the boundary lower-dimensional subdomain '" +
+                     _mesh.getSubdomainName(id) + "'");
 
   if (_lowerd_var.count() != _count)
     paramError("lowerd_variable",

framework/src/bcs/HFEMDirichletBC.C

@@ -29,10 +29,14 @@ HFEMDirichletBC::HFEMDirichletBC(const InputParameters & parameters)
 {
   if (_uhat_var)
   {
-    if (!_uhat_var->activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
-      paramError("uhat",
-                 "Must be defined on BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is added by "
-                 "Mesh/build_all_side_lowerd_mesh=true");
+    for (const auto & id : _uhat_var->activeSubdomains())
+      if (_mesh.boundaryLowerDBlocks().count(id) == 0)
+        mooseDocumentedError("moose",
+                             29151,
+                             "'uhat' must be defined on lower-dimensional boundary subdomain '" +
+                                 _mesh.getSubdomainName(id) +
+                                 "' that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe "
+                                 "check could be overly restrictive.");
 
     if (isParamValid("value"))
       paramError("uhat", "'uhat' and 'value' can not be both provided");

framework/src/bcs/LowerDIntegratedBC.C

@@ -34,15 +34,21 @@ LowerDIntegratedBC::LowerDIntegratedBC(const InputParameters & parameters)
     _test_lambda(_lowerd_var.phiLower())
 {
   const auto & lower_domains = _lowerd_var.activeSubdomains();
-  if (!lower_domains.count(Moose::BOUNDARY_SIDE_LOWERD_ID) && lower_domains.size() != 1)
-    paramError(
-        "lowerd_variable",
-        "Must be only defined on the subdomain BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is "
-        "added by Mesh/build_all_side_lowerd_mesh=true");
-
-  if (_var.activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
-    paramError("variable",
-               "Must not be defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain");
+  for (const auto & id : _mesh.boundaryLowerDBlocks())
+    if (lower_domains.count(id) == 0)
+      mooseDocumentedError(
+          "moose",
+          29151,
+          "'lowerd_variable' must be defined on the boundary lower-dimensional subdomain '" +
+              _mesh.getSubdomainName(id) +
+              "' that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe check could be overly "
+              "restrictive.");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.boundaryLowerDBlocks().count(id) > 0)
+      paramError("variable",
+                 "Must not be defined on the boundary lower-dimensional subdomain '" +
+                     _mesh.getSubdomainName(id) + "'");
 
   // Note: the above two conditions also ensure that the variable and lower-d variable are
   // different.

framework/src/dgkernels/ArrayDGLowerDKernel.C

@@ -42,15 +42,21 @@ ArrayDGLowerDKernel::ArrayDGLowerDKernel(const InputParameters & parameters)
     _work_vector(_count)
 {
   const auto & lower_domains = _lowerd_var.activeSubdomains();
-  if (!lower_domains.count(Moose::INTERNAL_SIDE_LOWERD_ID) && lower_domains.size() != 1)
-    paramError(
-        "lowerd_variable",
-        "Must be only defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain that is "
-        "added by Mesh/build_all_side_lowerd_mesh=true");
-
-  if (_var.activeSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID))
-    paramError("variable",
-               "Must not be defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain");
+  for (const auto & id : _mesh.interiorLowerDBlocks())
+    if (lower_domains.count(id) == 0)
+      mooseDocumentedError(
+          "moose",
+          29151,
+          "'lowerd_variable' must be defined on the interior lower-dimensional subdomain '" +
+              _mesh.getSubdomainName(id) +
+              "'that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe check could be overly "
+              "restrictive.");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.interiorLowerDBlocks().count(id) > 0)
+      paramError("variable",
+                 "Must not be defined on the interior lower-dimensional subdomain '" +
+                     _mesh.getSubdomainName(id) + "'");
 
   if (_lowerd_var.count() != _count)
     paramError("lowerd_variable",

framework/src/dgkernels/DGLowerDKernel.C

@@ -41,15 +41,21 @@ DGLowerDKernel::DGLowerDKernel(const InputParameters & parameters)
     _test_lambda(_lowerd_var.phiLower())
 {
   const auto & lower_domains = _lowerd_var.activeSubdomains();
-  if (!lower_domains.count(Moose::INTERNAL_SIDE_LOWERD_ID) && lower_domains.size() != 1)
-    paramError(
-        "lowerd_variable",
-        "Must be only defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain that is "
-        "added by Mesh/build_all_side_lowerd_mesh=true");
-
-  if (_var.activeSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID))
-    paramError("variable",
-               "Must not be defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain");
+  for (const auto & id : _mesh.interiorLowerDBlocks())
+    if (lower_domains.count(id) == 0)
+      mooseDocumentedError(
+          "moose",
+          29151,
+          "'lowerd_variable' must be defined on the interior lower-dimensional subdomain '" +
+              _mesh.getSubdomainName(id) +
+              "'that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe check could be overly "
+              "restrictive.");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.interiorLowerDBlocks().count(id) > 0)
+      paramError("variable",
+                 "Must not be defined on the interior lower-dimensional subdomain'" +
+                     _mesh.getSubdomainName(id) + "'");
 
   // Note: the above two conditions also ensure that the variable and lower-d variable are
   // different.

framework/src/interfaces/BlockRestrictable.C

@@ -338,7 +338,7 @@ BlockRestrictable::checkVariable(const MooseVariableFieldBase & variable) const
 {
   // a variable defined on all internal sides does not need this check because
   // it can be coupled with other variables in DG kernels
-  if (variable.activeSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) > 0)
+  if (variable.activeOnSubdomains(_blk_mesh->interiorLowerDBlocks()))
     return;
 
   if (!isBlockSubset(variable.activeSubdomains()))

framework/src/loops/ComputeFullJacobianThread.C

@@ -300,7 +300,7 @@ ComputeFullJacobianThread::computeOnInternalFace(const Elem * neighbor)
         if (dg->variable().number() == ivar && dg->isImplicit() &&
             dg->hasBlocks(neighbor->subdomain_id()) &&
             (jvariable.activeOnSubdomain(_subdomain) ||
-             jvariable.activeOnSubdomain(Moose::INTERNAL_SIDE_LOWERD_ID)))
+             jvariable.activeOnSubdomains(_fe_problem.mesh().interiorLowerDBlocks())))
         {
           dg->prepareShapes(jvar);
           dg->prepareNeighborShapes(jvar);

framework/src/mesh/MooseMesh.C

@@ -60,6 +60,7 @@
 #include "libmesh/default_coupling.h"
 #include "libmesh/ghost_point_neighbors.h"
 #include "libmesh/fe_type.h"
+#include "libmesh/enum_to_string.h"
 
 static const int GRAIN_SIZE =
     1; // the grain_size does not have much influence on our execution speed
@@ -297,6 +298,8 @@ MooseMesh::MooseMesh(const MooseMesh & other_mesh)
     _patch_update_strategy(other_mesh._patch_update_strategy),
     _regular_orthogonal_mesh(false),
     _is_split(other_mesh._is_split),
+    _lower_d_interior_blocks(other_mesh._lower_d_interior_blocks),
+    _lower_d_boundary_blocks(other_mesh._lower_d_boundary_blocks),
     _has_lower_d(other_mesh._has_lower_d),
     _allow_recovery(other_mesh._allow_recovery),
     _construct_node_list_from_side_list(other_mesh._construct_node_list_from_side_list),
@@ -632,19 +635,69 @@ MooseMesh::buildLowerDMesh()
   // remove existing lower-d element first
   std::set<Elem *> deleteable_elems;
   for (auto & elem : mesh.element_ptr_range())
-    if (elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID ||
-        elem->subdomain_id() == Moose::BOUNDARY_SIDE_LOWERD_ID)
+    if (_lower_d_interior_blocks.count(elem->subdomain_id()) ||
+        _lower_d_boundary_blocks.count(elem->subdomain_id()))
       deleteable_elems.insert(elem);
     else if (elem->n_sides() > max_n_sides)
       max_n_sides = elem->n_sides();
 
   for (auto & elem : deleteable_elems)
     mesh.delete_elem(elem);
+  for (const auto & id : _lower_d_interior_blocks)
+    _mesh_subdomains.erase(id);
+  for (const auto & id : _lower_d_boundary_blocks)
+    _mesh_subdomains.erase(id);
+  _lower_d_interior_blocks.clear();
+  _lower_d_boundary_blocks.clear();
 
   mesh.comm().max(max_n_sides);
 
   deleteable_elems.clear();
 
+  // get all side types
+  std::set<int> interior_side_types;
+  std::set<int> boundary_side_types;
+  for (const auto & elem : mesh.active_element_ptr_range())
+    for (const auto side : elem->side_index_range())
+    {
+      Elem * neig = elem->neighbor_ptr(side);
+      std::unique_ptr<Elem> side_elem(elem->build_side_ptr(side));
+      if (neig)
+        interior_side_types.insert(side_elem->type());
+      else
+        boundary_side_types.insert(side_elem->type());
+    }
+  mesh.comm().set_union(interior_side_types);
+  mesh.comm().set_union(boundary_side_types);
+
+  // assign block ids for different side types
+  std::map<ElemType, SubdomainID> interior_block_ids;
+  std::map<ElemType, SubdomainID> boundary_block_ids;
+  // we assume this id is not used by the mesh
+  auto id = libMesh::Elem::invalid_subdomain_id - 2;
+  for (const auto & tpid : interior_side_types)
+  {
+    const auto type = ElemType(tpid);
+    mesh.subdomain_name(id) = "INTERNAL_SIDE_LOWERD_SUBDOMAIN_" + Utility::enum_to_string(type);
+    interior_block_ids[type] = id;
+    _lower_d_interior_blocks.insert(id);
+    if (_mesh_subdomains.count(id) > 0)
+      mooseError("Trying to add a mesh block with id ", id, " that has existed in the mesh");
+    _mesh_subdomains.insert(id);
+    --id;
+  }
+  for (const auto & tpid : boundary_side_types)
+  {
+    const auto type = ElemType(tpid);
+    mesh.subdomain_name(id) = "BOUNDARY_SIDE_LOWERD_SUBDOMAIN_" + Utility::enum_to_string(type);
+    boundary_block_ids[type] = id;
+    _lower_d_boundary_blocks.insert(id);
+    if (_mesh_subdomains.count(id) > 0)
+      mooseError("Trying to add a mesh block with id ", id, " that has existed in the mesh");
+    _mesh_subdomains.insert(id);
+    --id;
+  }
+
   dof_id_type max_elem_id = mesh.max_elem_id();
   unique_id_type max_unique_id = mesh.parallel_max_unique_id();
 
@@ -681,9 +734,9 @@ MooseMesh::buildLowerDMesh()
 
         // Add subdomain ID
         if (neig)
-          side_elem->subdomain_id() = Moose::INTERNAL_SIDE_LOWERD_ID;
+          side_elem->subdomain_id() = interior_block_ids.at(side_elem->type());
         else
-          side_elem->subdomain_id() = Moose::BOUNDARY_SIDE_LOWERD_ID;
+          side_elem->subdomain_id() = boundary_block_ids.at(side_elem->type());
 
         // set ids consistently across processors (these ids will be temporary)
         side_elem->set_id(max_elem_id + elem->id() * max_n_sides + side);
@@ -713,11 +766,6 @@ MooseMesh::buildLowerDMesh()
   for (auto & elem : side_elems)
     mesh.add_elem(elem);
 
-  _mesh_subdomains.insert(Moose::INTERNAL_SIDE_LOWERD_ID);
-  mesh.subdomain_name(Moose::INTERNAL_SIDE_LOWERD_ID) = "INTERNAL_SIDE_LOWERD_SUBDOMAIN";
-  _mesh_subdomains.insert(Moose::BOUNDARY_SIDE_LOWERD_ID);
-  mesh.subdomain_name(Moose::BOUNDARY_SIDE_LOWERD_ID) = "BOUNDARY_SIDE_LOWERD_SUBDOMAIN";
-
   // we do all the stuff in prepare_for_use such as renumber_nodes_and_elements(),
   // update_parallel_id_counts(), cache_elem_dims(), etc. except partitioning here.
   const bool skip_partitioning_old = mesh.skip_partitioning();
@@ -1326,9 +1374,13 @@ MooseMesh::cacheInfo()
   _block_node_list.clear();
   _higher_d_elem_side_to_lower_d_elem.clear();
   _lower_d_elem_to_higher_d_elem_side.clear();
+  _lower_d_interior_blocks.clear();
+  _lower_d_boundary_blocks.clear();
+
+  auto & mesh = getMesh();
 
   // TODO: Thread this!
-  for (const auto & elem : getMesh().element_ptr_range())
+  for (const auto & elem : mesh.element_ptr_range())
   {
     const Elem * ip_elem = elem->interior_parent();
 
@@ -1346,6 +1398,18 @@ MooseMesh::cacheInfo()
             std::pair<std::pair<const Elem *, unsigned short int>, const Elem *>(pair, elem));
         _lower_d_elem_to_higher_d_elem_side.insert(
             std::pair<const Elem *, unsigned short int>(elem, ip_side));
+
+        auto id = elem->subdomain_id();
+        if (ip_elem->neighbor_ptr(ip_side))
+        {
+          if (mesh.subdomain_name(id).find("INTERNAL_SIDE_LOWERD_SUBDOMAIN_") != std::string::npos)
+            _lower_d_interior_blocks.insert(id);
+        }
+        else
+        {
+          if (mesh.subdomain_name(id).find("BOUNDARY_SIDE_LOWERD_SUBDOMAIN_") != std::string::npos)
+            _lower_d_boundary_blocks.insert(id);
+        }
       }
     }
 
@@ -1356,7 +1420,7 @@ MooseMesh::cacheInfo()
     }
   }
 
-  for (const auto & elem : getMesh().active_local_element_ptr_range())
+  for (const auto & elem : mesh.active_local_element_ptr_range())
   {
     SubdomainID subdomain_id = elem->subdomain_id();
     auto & sub_data = _sub_to_data[subdomain_id];