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().size() > 0 || _mesh->boundaryLowerDBlocks().size() > 0)
             mooseError("HFEM does not support mesh uniform refinement currently.");
 
           Adaptivity::uniformRefine(_mesh.get());

framework/src/bcs/ArrayHFEMDirichletBC.C

@@ -34,10 +34,11 @@ 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)
+        paramError("uhat",
+                   "Must be defined on lower-dimensional boundary subdomain that is added by "
+                   "Mesh/build_all_side_lowerd_mesh=true");
 
     if (_uhat_var->count() != _count)
       paramError("uhat",

framework/src/bcs/ArrayLowerDIntegratedBC.C

@@ -35,15 +35,15 @@ 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 : lower_domains)
+    if (_mesh.boundaryLowerDBlocks().count(id) == 0)
+      paramError("lowerd_variable",
+                 "Must be only defined on the boundary lower-dimensional subdomain that is "
+                 "added by Mesh/build_all_side_lowerd_mesh=true");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.boundaryLowerDBlocks().count(id) > 0)
+      paramError("variable", "Must not be defined on the boundary lower-dimensional subdomain");
 
   if (_lowerd_var.count() != _count)
     paramError("lowerd_variable",

framework/src/bcs/HFEMDirichletBC.C

@@ -29,10 +29,11 @@ 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)
+        paramError("uhat",
+                   "Must be defined on lower-dimensional boundary subdomain that is added by "
+                   "Mesh/build_all_side_lowerd_mesh=true");
 
     if (isParamValid("value"))
       paramError("uhat", "'uhat' and 'value' can not be both provided");

framework/src/bcs/LowerDIntegratedBC.C

@@ -34,15 +34,15 @@ 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 : lower_domains)
+    if (_mesh.boundaryLowerDBlocks().count(id) == 0)
+      paramError("lowerd_variable",
+                 "Must be only defined on the boundary lower-dimensional subdomain that is "
+                 "added by Mesh/build_all_side_lowerd_mesh=true");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.boundaryLowerDBlocks().count(id) > 0)
+      paramError("variable", "Must not be defined on the boundary lower-dimensional subdomain");
 
   // Note: the above two conditions also ensure that the variable and lower-d variable are
   // different.

framework/src/dgkernels/ArrayDGLowerDKernel.C

@@ -42,15 +42,15 @@ 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 : lower_domains)
+    if (_mesh.interiorLowerDBlocks().count(id) == 0)
+      paramError("lowerd_variable",
+                 "Must be only defined on the interior lower-dimensional subdomain that is "
+                 "added by Mesh/build_all_side_lowerd_mesh=true");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.interiorLowerDBlocks().count(id) > 0)
+      paramError("variable", "Must not be defined on the interior lower-dimensional subdomain");
 
   if (_lowerd_var.count() != _count)
     paramError("lowerd_variable",

framework/src/dgkernels/DGLowerDKernel.C

@@ -41,15 +41,15 @@ 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 : lower_domains)
+    if (_mesh.interiorLowerDBlocks().count(id) == 0)
+      paramError("lowerd_variable",
+                 "Must be only defined on the interior lower-dimensional subdomain that is "
+                 "added by Mesh/build_all_side_lowerd_mesh=true");
+
+  for (const auto & id : _var.activeSubdomains())
+    if (_mesh.interiorLowerDBlocks().count(id) > 0)
+      paramError("variable", "Must not be defined on the interior lower-dimensional subdomain");
 
   // Note: the above two conditions also ensure that the variable and lower-d variable are
   // different.

framework/src/interfaces/BlockRestrictable.C

@@ -338,8 +338,9 @@ 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)
-    return;
+  for (const auto & id : variable.activeSubdomains())
+    if (_blk_mesh->interiorLowerDBlocks().count(id) > 0)
+      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

@@ -297,6 +297,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 +634,65 @@ 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);
+  _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 +729,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 +761,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();

framework/src/problems/DisplacedProblem.C

@@ -875,7 +875,7 @@ DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem,
   reinitNeighbor(elem, side, tid);
 
   const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
-  if (lower_d_elem && lower_d_elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID)
+  if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0)
     reinitLowerDElem(lower_d_elem, tid);
   else
   {
@@ -884,7 +884,7 @@ DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem,
     auto & neighbor_side = _assembly[tid][currentNlSysNum()]->neighborSide();
     const Elem * lower_d_elem_neighbor = _mesh.getLowerDElem(neighbor, neighbor_side);
     if (lower_d_elem_neighbor &&
-        lower_d_elem_neighbor->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID)
+        _mesh.interiorLowerDBlocks().count(lower_d_elem_neighbor->subdomain_id()) > 0)
     {
       auto qps = _assembly[tid][currentNlSysNum()]->qPointsFaceNeighbor().stdVector();
       std::vector<Point> reference_points;

framework/src/problems/FEProblemBase.C

@@ -2253,7 +2253,7 @@ FEProblemBase::reinitElemNeighborAndLowerD(const Elem * elem,
   reinitNeighbor(elem, side, tid);
 
   const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
-  if (lower_d_elem && lower_d_elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID)
+  if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0)
     reinitLowerDElem(lower_d_elem, tid);
   else
   {
@@ -2262,7 +2262,7 @@ FEProblemBase::reinitElemNeighborAndLowerD(const Elem * elem,
     auto & neighbor_side = _assembly[tid][0]->neighborSide();
     const Elem * lower_d_elem_neighbor = _mesh.getLowerDElem(neighbor, neighbor_side);
     if (lower_d_elem_neighbor &&
-        lower_d_elem_neighbor->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID)
+        _mesh.interiorLowerDBlocks().count(lower_d_elem_neighbor->subdomain_id()) > 0)
     {
       auto qps = _assembly[tid][0]->qPointsFaceNeighbor().stdVector();
       std::vector<Point> reference_points;
@@ -7711,8 +7711,7 @@ FEProblemBase::initialAdaptMesh()
   _cycles_completed = 0;
   if (n)
   {
-    if (_mesh.meshSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) ||
-        _mesh.meshSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
+    if (_mesh.interiorLowerDBlocks().size() > 0 || _mesh.boundaryLowerDBlocks().size() > 0)
       mooseError("HFEM does not support mesh adaptivity currently.");
 
     TIME_SECTION("initialAdaptMesh", 2, "Performing Initial Adaptivity");
@@ -7757,8 +7756,7 @@ FEProblemBase::adaptMesh()
 
   for (unsigned int i = 0; i < cycles_per_step; ++i)
   {
-    if (_mesh.meshSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) ||
-        _mesh.meshSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID))
+    if (_mesh.interiorLowerDBlocks().size() > 0 || _mesh.boundaryLowerDBlocks().size() > 0)
       mooseError("HFEM does not support mesh adaptivity currently.");
 
     // Markers were already computed once by Executioner

framework/src/systems/NonlinearSystemBase.C

@@ -3640,8 +3640,10 @@ NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subd
                         std::inserter(difference, difference.end()));
 
     // there supposed to be no kernels on this lower-dimensional subdomain
-    difference.erase(Moose::INTERNAL_SIDE_LOWERD_ID);
-    difference.erase(Moose::BOUNDARY_SIDE_LOWERD_ID);
+    for (const auto & id : _mesh.interiorLowerDBlocks())
+      difference.erase(id);
+    for (const auto & id : _mesh.boundaryLowerDBlocks())
+      difference.erase(id);
 
     if (!difference.empty())
     {
@@ -3681,8 +3683,9 @@ NonlinearSystemBase::checkKernelCoverage(const std::set<SubdomainID> & mesh_subd
   for (auto & var_name : vars)
   {
     auto blks = getSubdomainsForVar(var_name);
-    if (blks.count(Moose::INTERNAL_SIDE_LOWERD_ID) || blks.count(Moose::BOUNDARY_SIDE_LOWERD_ID))
-      difference.erase(var_name);
+    for (const auto & id : blks)
+      if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0)
+        difference.erase(var_name);
   }
 
   if (!difference.empty())