Improve ElementSubdomainModifier moving boundary and reinitialization

- moving_boundaries and moving_boundary_subdomain_pairs for sideset/nodeset modification
- active_subdomains and amorphous_activation for initialization
- Update tests for moving boundary
- Update gold for changed exodus output
- reinitialize_subdomains and previous_subdomain_reinitialized to control reinitialization of changed elements
- Remove confusion from previous use of 'active'

refs #25736

framework/doc/content/source/meshmodifiers/CoupledVarThresholdElementSubdomainModifier.md

@@ -13,12 +13,12 @@ The `CoupledVarThresholdElementSubdomainModifier` can model
 The `CoupledVarThresholdElementSubdomainModifier` changes the element subdomain based on the given criterion. It also handles the corresponding
 
 - Moving boundary/interface nodeset/sideset modification,
-- Solution initialization, and
-- Stateful material property initialization,
+- Solution reinitialization, and
+- Stateful material property reinitialization,
 
-all of which are demonstrated using the following example.
+all of which are demonstrated below.
 
-Consider a unit square domain, and an auxiliary variable defined by the function $\phi(x,y,t) = (x-t)^2+y^2-0.5^2$. The function represents a signed-distance function of a circle of radius $0.5$ whose center is moving along the x-axis towards right.
+Consider a unit square domain, and an auxiliary variable defined by the function $\phi(x,y,t) = (x-t)^2+y^2-0.5^2$. The function represents a signed-distance function of a circle of radius $0.5$ whose center is moving along the x-axis towards the right.
 
 Initially, the domain is decomposed by a vertical line $x=0.25$. The elements on the left side of the vertical line have subdomain ID of 1, and the elements on the right side have subdomain ID of 2. The `CoupledVarThresholdElementSubdomainModifier` is used to change the subdomain ID from 2 to 1 for elements within the circle.
 
@@ -28,42 +28,77 @@ If the `CoupledVarThresholdElementSubdomainModifier` is applied onto the entire
 
 !listing test/tests/meshmodifiers/element_subdomain_modifier/reversible.i start=[moving_circle] end=[] include-end=true
 
-!media media/meshmodifiers/esm_reversible.jpg style=float:center;width:100%; caption=The result of a reversible element subdomain modifier at three different time steps
+!media media/meshmodifiers/rev.png style=float:center;width:100%; caption=The result of a reversible element subdomain modifier at three different time steps
 
 However, in many applications, e.g. element death and activation, the equivalent movement of element subdomains is not reversible. In this case, omitting the parameter `complement_subdomain_id` will make the subdomain modification irreversible:
 
-!listing test/tests/meshmodifiers/element_subdomain_modifier/block_restricted.i start=[moving_circle] end=[] include-end=true
+!listing test/tests/meshmodifiers/element_subdomain_modifier/irreversible.i start=[moving_circle] end=[] include-end=true
 
-!media media/meshmodifiers/esm_irreversible.jpg style=float:center;width:100%; caption=The result of an irreversible element subdomain modifier at three different time steps
+!media media/meshmodifiers/irrev.png style=float:center;width:100%; caption=The result of an irreversible element subdomain modifier at three different time steps
 
 ## Moving boundary/interface nodeset/sideset modification
 
-The change of element subdomains will alter the definition of certain sidesets and nodesets. The `CoupledVarThresholdElementSubdomainModifier` optionally takes the parameter `moving_boundary_name` to help modify the corresponding sideset/nodeset. If the boundary provided through the `moving_boundary_name` parameter already exists, the modifier will attempt to modify the provided sideset/nodeset whenever an element changes subdomain. If the boundary does not exist, the modifier will create a sideset and a nodeset with the provided name.
+The change of element subdomains will alter the definitions of certain sidesets and nodesets. The parameters `moving_boundaries` and `moving_boundary_subdomain_pairs` can optionally be used to modify the corresponding sidesets/nodesets. The pair of subdomains that each boundary in `moving_boundaries` lies between must be specified in the corresponding `moving_boundary_subdomain_pairs`, with the element side from the first subdomain added to the boundary.
 
-!media media/meshmodifiers/esm_sideset.jpg style=float:center;width:100%; caption=The evolving sideset (green) at three different time steps
+If the boundaries provided through `moving_boundaries` already exist, the modifier will attempt to modify the provided sidesets/nodesets whenever an element changes subdomain ID. If the boundaries do not exist, the modifier will create sidesets and nodesets with the provided names.
 
-!media media/meshmodifiers/esm_nodeset.jpg style=float:center;width:100%; caption=The evolving nodeset (green) at three different time steps
+!listing test/tests/meshmodifiers/element_subdomain_modifier/moving_boundary.i start=[moving_circle] end=[] include-end=true
 
-Nodal and integrated BCs can be applied on the moving boundary.
+!media media/meshmodifiers/nodeset.png style=float:center;width:100%; caption=The evolving nodeset (green) between subdomains 1 and 2, as created by the modifier without an existing boundary.
+
+If only one boundary is provided but multiple pairs of subdomains are specified, then all the pairs are applied to the one boundary. Element sides on a subdomain's external boundary can be added by specifying only one subdomain.
+
+!listing test/tests/meshmodifiers/element_subdomain_modifier/external_moving_boundary.i start=[ext] end=[AuxVariables] include-end=false
 
-## Solution initialization
+!media media/meshmodifiers/ext.png style=float:center;width:100%; caption=The evolving sideset (green) around subdomain 1, including the external element sides, from an existing boundary.
+
+Nodal and integrated BCs can be applied on the moving boundary.
 
-Depending on the physics, one may or may not want to initialize the solution when an element and its related nodes change subdomain.
-The parameter `apply_initial_conditions` defaults to true and determines whether the initial conditions should be re-evaluated.
+## Solution reinitialization
 
-Suppose initially there is an auxiliary variable $u=1$ everywhere inside the domain, and the variable value in subdomain 1 (blue) doubles at each time step:
+By default, all elements that change subdomain ID are reinitialized to the new subdomain's initial condition. Suppose the auxiliary variable $u$ has an initial variable value of $1$ in subdomain 1 and $-0.5$ in subdomain 2, and the variable value doubles at each timestep in subdomain 1:
 
-!listing test/tests/meshmodifiers/element_subdomain_modifier/initial_condition.i start=[AuxVariables] end=[Executioner]
+!listing test/tests/meshmodifiers/element_subdomain_modifier/initial_condition.i start=[ICs] end=[Postprocessors]
 
-!media media/meshmodifiers/esm_ic.jpg style=float:center;width:100%; caption=The auxiliary variable $u$ at three different time steps
+!media media/meshmodifiers/init_cond.png style=float:center;width:100%; caption=The auxiliary variable $u$ is reinitialized to $1$, which doubles over the timestep to $2$, for all the elements that change subdomain ID to 1
 
-## Stateful material property initialization
+## Stateful material property reinitialization
 
-Similarly, all stateful material properties will be re-initialized when an element changes subdomain. Suppose initially the diffusivity is $0.5$ everywhere, and the diffusivity doubles at each time step:
+Similarly, all stateful material properties will be re-initialized when an element changes subdomain ID. Suppose initially the diffusivity is $0.5$ in subdomain 1 and $-1$ in subdomain 2, and the diffusivity doubles at each time step in subdomain 1:
 
 !listing test/tests/meshmodifiers/element_subdomain_modifier/stateful_property.i start=[Materials] end=[Executioner]
 
-!media media/meshmodifiers/esm_material.jpg style=float:center;width:100%; caption=The diffusivity at three different time steps
+!media media/meshmodifiers/stateful_prop.png style=float:center;width:100%; caption=The diffusivity is reinitialized to $0.5$, which doubles over the timestep to $1$, for all the elements that change subdomain ID to 1.
+
+## Reinitialization restrictions
+
+Depending on the physics, one may or may not want to reinitialize the solution when an element and its related nodes change subdomain ID. For the below examples, consider a unit square domain decomposed by vertical lines $x=0.3$ and $x=0.6$. The elements on the left have subdomain ID of 1, the elements in the middle have subdomain ID of 2, and the elements on the right have subdomain ID of 3.
+
+Two auxiliary variables are defined by the functions $\phi(x,y,t) = (x-t)^2+y^2-0.3^2$ and $\phi(x,y,t) = (x-t)^2+(y-1)^2-0.3^2$, which represent the signed-distance functions of circles of radius $0.3$ whose centers are moving along the bottom and top of the square respectively. The `CoupledVarThresholdElementSubdomainModifier` is used to change the subdomain ID to 1 for elements within the bottom circle, and to subdomain ID 2 for elements within the top circle.
+
+An auxiliary variable $u$ is defined over the domain, with initial values of 1, 2, and 3 in subdomains 1, 2, and 3 respectively:
+
+!listing test/tests/meshmodifiers/element_subdomain_modifier/reinitialization.i start=[MeshModifiers] end=[AuxVariables]
+
+!media media/meshmodifiers/orig.png style=float:center;width:100%; caption=The default behaviour of the modifier is to reinitializate all elements that change subdomain ID.
+
+However, if a list of subdomains (IDs or names) is provided through the parameter `reinitialize_subdomains`, the reinitialization only occurs if an element changes subdomain ID, and the new subdomain ID is in the list:
+
+!listing test/tests/meshmodifiers/element_subdomain_modifier/reinitialization_into.i start=[MeshModifiers] end=[AuxVariables]
+
+!media media/meshmodifiers/into.png style=float:center;width:100%; caption=Reinitialization of only the elements that change subdomain ID to 1.
+
+If an empty list is given in `reinitialize_subdomains`, then there is no reinitialization of any elements that change subdomain ID.
+
+!listing test/tests/meshmodifiers/element_subdomain_modifier/no_reinitialization.i start=[MeshModifiers] end=[AuxVariables]
+
+!media media/meshmodifiers/none.png style=float:center;width:100%; caption=No reinitialization of any elements that change subdomain ID.
+
+Reinitialization can be further restricted by setting the parameter `previous_subdomain_reinitialized` to `false`. The modifier will then additionally check the changing element's previous subdomain ID. Reinitialization will only occur if the previous subdomain ID was not in the list provided in the parameter `reinitialize_subdomains`.
+
+!listing test/tests/meshmodifiers/element_subdomain_modifier/reinitialization_from_into.i start=[MeshModifiers] end=[AuxVariables]
+
+!media media/meshmodifiers/from_into.png style=float:center;width:100%; caption=Reinitialization of only the elements which change subdomain ID from 3, to subdomain IDs 1 or 2.
 
 !syntax parameters /MeshModifiers/CoupledVarThresholdElementSubdomainModifier
 

framework/include/meshmodifiers/ElementSubdomainModifier.h

@@ -9,149 +9,33 @@
 
 #pragma once
 
-#include "ElementUserObject.h"
-#include "NonlinearSystemBase.h"
-#include "AuxiliarySystem.h"
+#include "ElementSubdomainModifierBase.h"
 
-class ElementSubdomainModifier : public ElementUserObject
+class ElementSubdomainModifier : public ElementSubdomainModifierBase
 {
 public:
   static InputParameters validParams();
 
   ElementSubdomainModifier(const InputParameters & parameters);
 
-  virtual void initialSetup() override;
   virtual void initialize() override;
   virtual void execute() override;
   virtual void threadJoin(const UserObject & /*uo*/) override;
   virtual void finalize() override;
 
 protected:
-  /// Compute the subdomain ID of the current element
+  /**
+   * Compute the subdomain ID of the current element
+   *
+   * If the computed subdomain ID is different from the current element's current subdomain ID, the
+   * current element WILL be "moved" from its current subdomain to the computed subdomain.
+   *
+   * If the computed subdomain ID is (1) the same as the current element's current subdomain ID or
+   * (2) Moose::INVALID_BLOCK_ID, the current element does not participate in subdomain modification
+   * activities.
+   */
   virtual SubdomainID computeSubdomainID() = 0;
 
-  /// The ID of the moving boundary that this object creates/modifies.
-  BoundaryID movingBoundaryID() const
-  {
-    if (!_moving_boundary_specified)
-      mooseError("no moving boundary specified");
-    return _moving_boundary_id;
-  }
-
-  /// The ID of the complement moving boundary that this object creates/modifies.
-  BoundaryID complementMovingBoundaryID() const
-  {
-    if (!_complement_moving_boundary_specified)
-      mooseError("no complement moving boundary specified");
-    return _complement_moving_boundary_id;
-  }
-
-  /// The name of the moving boundary that this object creates/modifies.
-  const BoundaryName movingBoundaryName() const
-  {
-    if (!_moving_boundary_specified)
-      mooseError("no moving boundary specified");
-    return _moving_boundary_name;
-  }
-
-  /// The name of the complement moving boundary that this object creates/modifies.
-  const BoundaryName complementMovingBoundaryName() const
-  {
-    if (!_complement_moving_boundary_specified)
-      mooseError("no complement moving boundary specified");
-    return _complement_moving_boundary_name;
-  }
-
-  /// Range of the elements who changed their subdomain ID
-  ConstElemRange & movedElemsRange() const { return *_moved_elems_range; }
-
-  /// Range of the boundary nodes on moved elements
-  ConstBndNodeRange & movedBndNodesRange() const { return *_moved_bnd_nodes_range; }
-
-  /// Pointer to the displaced problem
-  DisplacedProblem * _displaced_problem;
-
-private:
-  /// Set the name of the moving boundary. Create the nodeset/sideset if not exist.
-  void setMovingBoundaryName(MooseMesh & mesh);
-
-  /// Set the name of the complement moving boundary. Create the nodeset/sideset if not exist.
-  void setComplementMovingBoundaryName(MooseMesh & mesh);
-
-  /// Update the moving boundary (both the underlying sideset and nodeset)
-  void updateMovingBoundaryInfo(MooseMesh & mesh, const std::vector<const Elem *> & moved_elems);
-
-  /// Update the complement boundary (both the underlying sideset and nodeset)
-  void updateComplementBoundaryInfo(MooseMesh & mesh,
-                                    const std::vector<const Elem *> & moved_elems);
-
-  /// Remove ghosted element sides
-  void pushBoundarySideInfo(MooseMesh & mesh);
-
-  /// Remove ghosted boundary nodes
-  void pushBoundaryNodeInfo(MooseMesh & mesh);
-
-  /// synchronize boundary information across processors
-  void synchronizeBoundaryInfo(MooseMesh & mesh);
-
-  /// Change the subdomain ID of all ancestor elements
-  void setAncestorsSubdomainIDs(const SubdomainID & subdomain_id, const dof_id_type & elem_id);
-
-  /// Helper function to build the range of moved elements
-  void buildMovedElemsRange();
-
-  /// Helper function to build the range of boundary nodes on moved elements
-  void buildMovedBndNodesRange();
-
-  /// Set old and older solutions to be the same as the current solution
-  void setOldAndOlderSolutionsForMovedNodes(SystemBase & sys);
-
-  /// Elements on the undisplaced mesh whose subdomain IDs have changed
-  std::vector<const Elem *> _moved_elems;
-
-  /// Elements on the displaced mesh whose subdomain IDs have changed
-  std::vector<const Elem *> _moved_displaced_elems;
-
-  /// Nodes on the moved elements
-  std::set<dof_id_type> _moved_nodes;
-
-  /// Range of the moved elements
-  std::unique_ptr<ConstElemRange> _moved_elems_range;
-
-  /// Range of the boundary nodes on the moved elements
-  std::unique_ptr<ConstBndNodeRange> _moved_bnd_nodes_range;
-
-  /// Whether to re-apply ICs on moved elements and moved nodes
-  const bool _apply_ic;
-
-  /// Whether a moving boundary name is provided
-  const bool _moving_boundary_specified;
-
-  /// Whether a complement moving boundary name is provided
-  const bool _complement_moving_boundary_specified;
-
-  /// The name of the moving boundary
-  BoundaryName _moving_boundary_name;
-
-  /// The name of the complement moving boundary
-  BoundaryName _complement_moving_boundary_name;
-
-  /// The Id of the moving boundary
-  BoundaryID _moving_boundary_id;
-
-  /// The Id of the complement moving boundary
-  BoundaryID _complement_moving_boundary_id;
-
-  /// save the added/removed ghost sides to sync across processors
-  std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned short int>>>
-      _complement_ghost_sides_to_remove, _complement_ghost_sides_to_add, _ghost_sides_to_add;
-
-  /// save the added/removed ghost nodes to sync across processors
-  std::unordered_map<processor_id_type, std::vector<dof_id_type>> _complement_ghost_nodes_to_remove;
-
-  /// Any subdomain change is stored in this map and only applied in finalize to avoid messing up other UOs
-  std::vector<std::pair<Elem *, SubdomainID>> _cached_subdomain_assignments;
-
-  /// Subdomains between that the moving boundary is
-  std::set<SubdomainID> _moving_boundary_subdomains;
+  /// Element subdomain assignments
+  std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> _moved_elems;
 };