Add support for declaring sidesets in input file that do not exist initially

[GiudGiud]

Reason

See discussion below.
The input errors because the sideset does not exist at t=0, but will exist later in the simulation

Design

Same as add-subdomain-id parameter in MooseMesh
A parameter, then simple logic adding the sideset id to the list of sideset ids and boundary ids

Impact

Facilitate simulations with changing meshes

Discussed in #27857

^{Originally posted by Traiwit June 10, 2024}

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hi guys,

I have 3 questions RE: ElementSubdomainModifier (more specifically CoupledVarThresholdElementSubdomainModifier)

First question, how do I make BC wait for the boundary to be generated? lets say i'm generating a boundary at every time step

 [new_side]
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'state_M0'
    criterion_type = EQUAL
    threshold = 1
    subdomain_id = 3
    complement_subdomain_id = 1
    moving_boundary_name = moving_boundary_tet
    execute_on = 'TIMESTEP_BEGIN'
[]

but during the initialization, it just gives this error

*** ERROR ***
(BCs/bc_drain0001_Pit/boundary):
    the following side sets (ids) do not exist on the mesh: moving_boundary_tet (21)
    
    MOOSE distinguishes between "node sets" and "side sets" depending on whether 
    you are using "Nodal" or "Integrated" BCs respectively. Node sets corresponding 
    to your side sets are constructed for you by default.
    
    Try setting "Mesh/construct_side_list_from_node_list=true" if you see this error.
    Note: If you are running with adaptivity you should prefer using side sets.

application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

I tried to ignore the checks, but still doesnt work

[Problem]
  boundary_restricted_elem_integrity_check = false  # default: true
  boundary_restricted_node_integrity_check = false   # default: true
  check_uo_aux_state = false                        # default: false
  error_on_jacobian_nonzero_reallocation = false    # default: false
  fv_bcs_integrity_check = false                     # default: true
  kernel_coverage_check = false                      # default: true
  material_coverage_check = false                    # default: true
  material_dependency_check = false                  # default: true
  skip_nl_system_check = false                      # default: false
[]

---

Question2: if I call CoupledVarThresholdElementSubdomainModifier twice in 1 input file, the surface/side generated is overwritten the previous one, for example if I do this

[tet]
  type = CoupledVarThresholdElementSubdomainModifier
  coupled_var = 'state_M0'
  criterion_type = EQUAL
  block = '1 3'
  threshold = 1
  subdomain_id = 3
  complement_subdomain_id = 1
  moving_boundary_name = moving_boundary_tet
  execute_on = 'TIMESTEP_BEGIN'
[]


[wedge]
  type = CoupledVarThresholdElementSubdomainModifier
  coupled_var = 'state_M0'
  block = '2 4'
  criterion_type = EQUAL
  threshold = 1
  subdomain_id = 4
  complement_subdomain_id = 2
  moving_boundary_name = moving_boundary_wedge
  execute_on = 'TIMESTEP_BEGIN'
[]

only moving_boundary_wedge will be produced in the mesh, how to overcome this problem? as later on, i might want to generate surface at the top of the mesh for surface recharge and stuff as well as the drain surface

---

Last question, this is probs the most headache one

I have modified ElementSubdomainModifier, so that it only produces surface(s) based on the variable, (i dont want anything to do with Subdomain whatsoever). I'm not so good with C++ so I mostly just commented out lines that related to subdomain

now I have a problem, the surface I generated is almost perfect, but it still generated artificial surfaces inside the main surface, I dig deep enough to realise that this 'inside surface' is caused by the parallelisation? the more mpi cores > the more inside surfaces

compare to the one I want (pre-made, not generated)

below is my modified ElementSubdomainModifier.C could you help me to improve this?

or should I make another function that remove the surfaces inside any surface (scan all directions, etc)?

#include "ElementSubdomainModifier.h"
#include "DisplacedProblem.h"

#include "libmesh/parallel_algebra.h"
#include "libmesh/parallel.h"
#include "libmesh/dof_map.h"
#include "libmesh/remote_elem.h"
#include "libmesh/parallel_ghost_sync.h"

InputParameters
ElementSubdomainModifier::validParams()
{
  InputParameters params = ElementUserObject::validParams();
  params.addClassDescription(
      "Modify element subdomain ID. This userobject only runs on the undisplaced mesh, and it will "
      "modify both the undisplaced and the displaced mesh.");
  params.addParam<bool>("apply_initial_conditions",
                        true,
                        "Whether to apply initial conditions on the moved nodes and elements");
  params.addParam<BoundaryName>(
      "moving_boundary_name",
      "Boundary to modify when an element is moved. A boundary with the provided name will be "
      "created if not already exists on the mesh.");
  params.addParam<BoundaryName>("complement_moving_boundary_name",
                                "Boundary that associated with the unmoved domain when neighbor "
                                "element(s) is moved. A boundary with the provided name will be "
                                "created if not already exists on the mesh.");
  params.set<bool>("use_displaced_mesh") = false;
  params.suppressParameter<bool>("use_displaced_mesh");
  return params;
}

ElementSubdomainModifier::ElementSubdomainModifier(const InputParameters & parameters)
  : ElementUserObject(parameters),
    _displaced_problem(_fe_problem.getDisplacedProblem().get()),
    _apply_ic(getParam<bool>("apply_initial_conditions")),
    _moving_boundary_specified(isParamValid("moving_boundary_name")),
    _complement_moving_boundary_specified(isParamValid("complement_moving_boundary_name")),
    _moving_boundary_id(-1)
{
}

void
ElementSubdomainModifier::initialSetup()
{
  if (_moving_boundary_specified)
  {
    _moving_boundary_name = getParam<BoundaryName>("moving_boundary_name");
    setMovingBoundaryName(_mesh);
    if (_displaced_problem)
      setMovingBoundaryName(_displaced_problem->mesh());
  }

  if (_complement_moving_boundary_specified)
  {
    _complement_moving_boundary_name = getParam<BoundaryName>("complement_moving_boundary_name");
    setComplementMovingBoundaryName(_mesh);
    if (_displaced_problem)
      setComplementMovingBoundaryName(_displaced_problem->mesh());
  }
}

void
ElementSubdomainModifier::setMovingBoundaryName(MooseMesh & mesh)
{
  // We only need one boundary to modify. Create a dummy vector just to use the API.
  const std::vector<BoundaryID> boundary_ids = mesh.getBoundaryIDs({{_moving_boundary_name}}, true);
  mooseAssert(boundary_ids.size() == 1, "Expect exactly one boundary ID.");
  _moving_boundary_id = boundary_ids[0];
  mesh.setBoundaryName(_moving_boundary_id, _moving_boundary_name);
  mesh.getMesh().get_boundary_info().sideset_name(_moving_boundary_id) = _moving_boundary_name;
  mesh.getMesh().get_boundary_info().nodeset_name(_moving_boundary_id) = _moving_boundary_name;
}

void
ElementSubdomainModifier::setComplementMovingBoundaryName(MooseMesh & mesh)
{
  // We only need one boundary to modify. Create a dummy vector just to use the API.
  const std::vector<BoundaryID> boundary_ids =
      mesh.getBoundaryIDs({{_complement_moving_boundary_name}}, true);
  mooseAssert(boundary_ids.size() == 1, "Expect exactly one boundary ID.");
  _complement_moving_boundary_id = boundary_ids[0];
  mesh.setBoundaryName(_complement_moving_boundary_id, _complement_moving_boundary_name);
  mesh.getMesh().get_boundary_info().sideset_name(_complement_moving_boundary_id) =
      _complement_moving_boundary_name;
  mesh.getMesh().get_boundary_info().nodeset_name(_complement_moving_boundary_id) =
      _complement_moving_boundary_name;
}

void
ElementSubdomainModifier::initialize()
{
  _moved_elems.clear();
  _moved_displaced_elems.clear();
  _moved_nodes.clear();
  _cached_subdomain_assignments.clear();
  _moving_boundary_subdomains.clear();
}

void
ElementSubdomainModifier::execute()
{
  // First, compute the desired subdomain ID for the current element.
  SubdomainID subdomain_id = computeSubdomainID();

  // If the current element's subdomain ID isn't what we want
  if (subdomain_id != std::numeric_limits<SubdomainID>::max() &&
      _current_elem->subdomain_id() != subdomain_id)
  {
    _moving_boundary_subdomains.insert(subdomain_id);
    _moving_boundary_subdomains.insert(_current_elem->subdomain_id());
    // Current element ID, used to index both the element on the displaced and undisplaced meshes.
    _mesh.getMesh().allow_renumbering(false);
    dof_id_type elem_id = _current_elem->id();
    // Change the element's subdomain
    Elem * elem = _mesh.elemPtr(elem_id);
    // Elem * displaced_elem =
    //     _displaced_problem ? _displaced_problem->mesh().elemPtr(elem_id) : nullptr;
    _mesh.getMesh().allow_renumbering(false);

    // Save the affected nodes so that we can later update/initialize the solution
    for (unsigned int i = 0; i < elem->n_nodes(); ++i)
      _moved_nodes.insert(elem->node_id(i));
    // _cached_subdomain_assignments.emplace_back(elem, subdomain_id);
    _moved_elems.push_back(elem);

  }
}


void
ElementSubdomainModifier::threadJoin(const UserObject & in_uo)
{
  // Join the data from uo into _this_ object:
  const auto & uo = static_cast<const ElementSubdomainModifier &>(in_uo);

  _moved_elems.insert(_moved_elems.end(), uo._moved_elems.begin(), uo._moved_elems.end());

  _moved_displaced_elems.insert(_moved_displaced_elems.end(),
                                uo._moved_displaced_elems.begin(),
                                uo._moved_displaced_elems.end());

  _moved_nodes.insert(uo._moved_nodes.begin(), uo._moved_nodes.end());
  // _cached_subdomain_assignments.insert(_cached_subdomain_assignments.end(),
  //                                      uo._cached_subdomain_assignments.begin(),
  //                                      uo._cached_subdomain_assignments.end());
}

void
ElementSubdomainModifier::finalize()
{
  // Do not apply cached subdomain changes
  // for (auto & [elem, subdomain_id] : _cached_subdomain_assignments)
  //   elem->subdomain_id() = subdomain_id;

  _ghost_sides_to_add.clear();
  _complement_ghost_sides_to_add.clear();
  _complement_ghost_sides_to_remove.clear();
  _complement_ghost_nodes_to_remove.clear();

  /*
    Synchronize ghost element subdomain ID
    Note: this needs to be done before updating boundary info because
    updating boundary requires the updated element subdomain ids
  */
  SyncSubdomainIds sync(_mesh.getMesh());
  Parallel::sync_dofobject_data_by_id(_mesh.getMesh().comm(),
                                      _mesh.getMesh().elements_begin(),
                                      _mesh.getMesh().elements_end(),
                                      sync);


  if (_moving_boundary_specified)
    updateMovingBoundaryInfo(_mesh, _moved_elems);

  // Synchronize boundary information after mesh update
  synchronizeBoundaryInfo(_mesh);

  // // Reinit equation systems
  // _fe_problem.meshChanged();

  _mesh.getMesh().allow_renumbering(false);

  // // Apply initial condition for the newly moved elements and boundary nodes
  buildMovedElemsRange();
  buildMovedBndNodesRange();



}




void
ElementSubdomainModifier::updateMovingBoundaryInfo(MooseMesh & mesh, const std::vector<const Elem *> & moved_elems)
{


  if (_moving_boundary_subdomains.size())
    mooseAssert(_moving_boundary_subdomains.size() == 2,
                "The number of moving subdomains should be two");

  BoundaryInfo & bnd_info = mesh.getMesh().get_boundary_info();

  bnd_info.sideset_name(_moving_boundary_id) = _moving_boundary_name;
  bnd_info.nodeset_name(_moving_boundary_id) = _moving_boundary_name;

  auto & elem_side_bnd_ids = bnd_info.get_sideset_map();
  std::set<const Elem *> boundary_elem_candidates(moved_elems.begin(), moved_elems.end());
  std::vector<std::pair<const Elem *, unsigned int>> elem_sides_to_be_cleared;

  std::cout << "Number of moved elements: " << moved_elems.size() << std::endl;

  // Collect current boundary elements and sides to be cleared
  for (const auto & [elem, side_bnd] : elem_side_bnd_ids)
  {
    auto side = side_bnd.first;
    auto boundary_id = side_bnd.second;
    if (boundary_id == _moving_boundary_id)
    {
      elem_sides_to_be_cleared.emplace_back(elem, side);
      // std::cout << "Marking old boundary side for removal: Elem ID " << elem->id() << ", Side " << side << std::endl;
    }
  }

  // Clear old boundary sides
  for (auto & [elem, side] : elem_sides_to_be_cleared)
  {
    bnd_info.remove_side(elem, side, _moving_boundary_id);
    // std::cout << "Removed old boundary side: Elem ID " << elem->id() << ", Side " << side << std::endl;
  }

  std::cout << "Old boundary sides removed." << std::endl;

  // Add moved elements to the candidate list
  for (auto elem : moved_elems)
    boundary_elem_candidates.insert(elem);

  // Set to collect new boundary nodes
  std::set<const Node *> boundary_nodes;

  // Process the boundary element candidates
  for (auto elem : boundary_elem_candidates)
  {
    for (auto side : elem->side_index_range())
    {
      const Elem * neighbor = elem->neighbor_ptr(side);

      // Debug statement to print element and side info
      // std::cout << "Processing Element ID: " << elem->id() << ", Side: " << side << std::endl;

      if (neighbor && neighbor->active() && neighbor != libMesh::remote_elem)
      {
        // Only add side if the neighbor is not part of the moved elements
        if (std::find(moved_elems.begin(), moved_elems.end(), neighbor) == moved_elems.end())
        {
          bnd_info.add_side(elem, side, _moving_boundary_id);
          // std::cout << "Adding new boundary side: Elem ID " << elem->id() << ", Side " << side << std::endl;

          // Collect nodes for this side
          auto nodes = elem->nodes_on_side(side);
          for (auto node : nodes)
            boundary_nodes.insert(&(elem->node_ref(node)));
        }
      }
      else if (neighbor && !neighbor->active() && neighbor != libMesh::remote_elem)
      {
        // std::cout << "Neighbor is not active, checking family members." << std::endl;
        std::vector<const Elem *> active_family;
        auto top_parent = neighbor->top_parent();
        top_parent->active_family_tree_by_neighbor(active_family, elem);
        for (auto felem : active_family)
        {
          if (std::find(moved_elems.begin(), moved_elems.end(), felem) == moved_elems.end())
          {
            auto cside = felem->which_neighbor_am_i(elem);
            bnd_info.add_side(felem, cside, _moving_boundary_id);
            // std::cout << "Adding new boundary side from family: Elem ID " << felem->id() << ", Side " << cside << std::endl;

            if (felem->processor_id() != this->processor_id())
              _ghost_sides_to_add[felem->processor_id()].emplace_back(felem->id(), cside);

            // Collect nodes for this side
            auto nodes = felem->nodes_on_side(cside);
            for (auto node : nodes)
              boundary_nodes.insert(&(felem->node_ref(node)));
          }
        }
      }
    }
  }

  // Clear old boundary nodes
  auto & nodeset_map = bnd_info.get_nodeset_map();
  std::vector<const Node *> nodes_elem_sides_to_be_cleared;
  for (const auto & pair : bnd_info.get_nodeset_map())
    if (pair.second == _moving_boundary_id)
      nodes_elem_sides_to_be_cleared.push_back(pair.first);

  for (const auto node : nodes_elem_sides_to_be_cleared)
  {
    bnd_info.remove_node(node, _moving_boundary_id);
    // std::cout << "Removed old boundary node: Node ID " << node->id() << std::endl;
  }

    /* Reconstruct a new nodeset from the updated sideset */
  // std::set<const Node *> boundary_nodes;
  for (const auto & [elem, side_bnd] : elem_side_bnd_ids)
  {
    auto side = side_bnd.first;
    auto boundary_id = side_bnd.second;
    if (boundary_id == _moving_boundary_id)
    {
      auto nodes = elem->nodes_on_side(side);
      for (auto node : nodes)
        boundary_nodes.insert(&(elem->node_ref(node)));
    }
  }

  // Add new boundary nodes
  for (const auto node : boundary_nodes)
  {
    bnd_info.add_node(node, _moving_boundary_id);
    // std::cout << "Added new boundary node: Node ID " << node->id() << std::endl;
  }

  std::cout << "Boundary nodes updated." << std::endl;

  // Ensure synchronization
  // synchronizeBoundaryInfo(mesh);
}


void
ElementSubdomainModifier::updateComplementBoundaryInfo(
    MooseMesh & mesh, const std::vector<const Elem *> & moved_elems)
{

}

void
ElementSubdomainModifier::synchronizeBoundaryInfo(MooseMesh & mesh)
{
  // std::cout << "Pushing boundary side info..." << std::endl;
  pushBoundarySideInfo(mesh);
  // std::cout << "Boundary side info pushed." << std::endl;

  // std::cout << "Pushing boundary node info..." << std::endl;
  pushBoundaryNodeInfo(mesh);
  // std::cout << "Boundary node info pushed." << std::endl;

  // // std::cout << "Synchronizing side ids..." << std::endl;
  mesh.getMesh().get_boundary_info().parallel_sync_side_ids();
  // // std::cout << "Side ids synchronized." << std::endl;

  // // std::cout << "Synchronizing node ids..." << std::endl;
  mesh.getMesh().get_boundary_info().parallel_sync_node_ids();

  std::cout << "SynchronizeBoundaryInfo is done." << std::endl;

  // std::cout << "Updating mesh..." << std::endl;
  // mesh.update();
  // std::cout << "Mesh updated." << std::endl;
}


void
ElementSubdomainModifier::pushBoundarySideInfo(MooseMesh & moose_mesh)
{
  auto & mesh = moose_mesh.getMesh();
  auto elem_remove_functor =
      [&mesh, this](processor_id_type,
                    const std::vector<std::pair<dof_id_type, unsigned short int>> & received_elem)
  {
    // remove the side
    for (const auto & pr : received_elem)
      mesh.get_boundary_info().remove_side(
          mesh.elem_ptr(pr.first), pr.second, _complement_moving_boundary_id);
  };

  // We create a tempalte functor to add with custom boundary IDs
  auto elem_add_functor_with_boundary_id =
      [&mesh](const std::vector<std::pair<dof_id_type, unsigned short int>> & received_elem,
              const BoundaryID boundary_id)
  {
    // add the side
    for (const auto & pr : received_elem)
      mesh.get_boundary_info().add_side(mesh.elem_ptr(pr.first), pr.second, boundary_id);
  };

  // Then we use it with the regular and the complement boundary id
  auto complement_elem_add_functor =
      [this, elem_add_functor_with_boundary_id](
          processor_id_type,
          const std::vector<std::pair<dof_id_type, unsigned short int>> & received_elem)
  { elem_add_functor_with_boundary_id(received_elem, _complement_moving_boundary_id); };

  auto elem_add_functor =
      [this, elem_add_functor_with_boundary_id](
          processor_id_type,
          const std::vector<std::pair<dof_id_type, unsigned short int>> & received_elem)
  { elem_add_functor_with_boundary_id(received_elem, _moving_boundary_id); };

  // Push/pull the ghost cell sides for the regular and complement boundaries
  Parallel::push_parallel_vector_data(
      mesh.get_boundary_info().comm(), _ghost_sides_to_add, elem_add_functor);
  Parallel::push_parallel_vector_data(
      mesh.get_boundary_info().comm(), _complement_ghost_sides_to_remove, elem_remove_functor);
  if (_complement_moving_boundary_specified)
    Parallel::push_parallel_vector_data(mesh.get_boundary_info().comm(),
                                        _complement_ghost_sides_to_add,
                                        complement_elem_add_functor);
}

void
ElementSubdomainModifier::pushBoundaryNodeInfo(MooseMesh & moose_mesh)
{
  auto & mesh = moose_mesh.getMesh();
  auto node_remove_functor =
      [&mesh, this](processor_id_type, const std::vector<dof_id_type> & received_nodes)
  {
    for (const auto & pr : received_nodes)
    {
      if (_moving_boundary_specified)
        mesh.get_boundary_info().remove_node(mesh.node_ptr(pr), _moving_boundary_id);
      if (_complement_moving_boundary_specified)
        mesh.get_boundary_info().remove_node(mesh.node_ptr(pr), _complement_moving_boundary_id);
    }
  };

  Parallel::push_parallel_vector_data(
      mesh.get_boundary_info().comm(), _complement_ghost_nodes_to_remove, node_remove_functor);
}



void
ElementSubdomainModifier::buildMovedElemsRange()
{
  // Clear the object first
  _moved_elems_range.reset();

  // Make some fake element iterators defining this vector of elements
  Elem * const * elem_itr_begin = const_cast<Elem * const *>(_moved_elems.data());
  Elem * const * elem_itr_end = elem_itr_begin + _moved_elems.size();

  const auto elems_begin = MeshBase::const_element_iterator(
      elem_itr_begin, elem_itr_end, Predicates::NotNull<Elem * const *>());
  const auto elems_end = MeshBase::const_element_iterator(
      elem_itr_end, elem_itr_end, Predicates::NotNull<Elem * const *>());

  _moved_elems_range = std::make_unique<ConstElemRange>(elems_begin, elems_end);
}

void
ElementSubdomainModifier::buildMovedBndNodesRange()
{
  // This is more involved than building the element range, because not all moved nodes are
  // necessarily associated with a boundary initial condition. We need to first build a set of
  // boundary nodes. Clear the object first:
  _moved_bnd_nodes_range.reset();

  // create a vector of the newly activated nodes
  std::set<const BndNode *> moved_bnd_nodes_set;
  for (auto & bnd_node : *_mesh.getBoundaryNodeRange())
  {
    dof_id_type bnd_node_id = bnd_node->_node->id();
    if (_moved_nodes.find(bnd_node_id) != _moved_nodes.end())
      moved_bnd_nodes_set.insert(bnd_node);
  }

  // Dump all the boundary nodes into a vector so that we can build a range out of it
  std::vector<const BndNode *> moved_bnd_nodes;
  moved_bnd_nodes.assign(moved_bnd_nodes_set.begin(), moved_bnd_nodes_set.end());

  // Make some fake node iterators defining this vector of nodes
  BndNode * const * bnd_node_itr_begin = const_cast<BndNode * const *>(moved_bnd_nodes.data());
  BndNode * const * bnd_node_itr_end = bnd_node_itr_begin + moved_bnd_nodes.size();

  const auto bnd_nodes_begin = MooseMesh::const_bnd_node_iterator(
      bnd_node_itr_begin, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
  const auto bnd_nodes_end = MooseMesh::const_bnd_node_iterator(
      bnd_node_itr_end, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());

  _moved_bnd_nodes_range = std::make_unique<ConstBndNodeRange>(bnd_nodes_begin, bnd_nodes_end);
}

void
ElementSubdomainModifier::setOldAndOlderSolutionsForMovedNodes(SystemBase & sys)
{

}

void
ElementSubdomainModifier::setAncestorsSubdomainIDs(const SubdomainID & subdomain_id,
                                                   const dof_id_type & elem_id)
{

}

if I'm not mistaken @hugary1995 is the one who created this class, right?

I'm also tagging: @GiudGiud since you know the general framework pretty well.

Thank you very much!

Kind regards,
Trai