diff --git a/framework/doc/content/source/userobject/ActivateElementsByPath.md b/framework/doc/content/source/userobject/ActivateElementsByPath.md
new file mode 100644
index 000000000000..df8a6cb19ef1
--- /dev/null
+++ b/framework/doc/content/source/userobject/ActivateElementsByPath.md
@@ -0,0 +1,9 @@
+# ActivateElementsByPath
+
+This user object uses the +function path+ as the metric to activate (add) an element by moving the element from an "inactive" subdomain to the "active" subdomain. It uses the user provided points $(x(t), y(t), z(t))$ with components defined by the functions specified by the parameters `function_x`, `function_y`, and `function_z` in the input. An element is activated at time $t_0$ if this element is close (distance < `activate_distance`) to the point $(x(t_0), y(t_0), z(t_0))$.
+
+!syntax parameters /UserObjects/ActivateElementsByPath
+
+!syntax inputs /UserObjects/ActivateElementsByPath
+
+!syntax children /UserObjects/ActivateElementsByPath
diff --git a/framework/doc/content/source/userobject/ActivateElementsCoupled.md b/framework/doc/content/source/userobject/ActivateElementsCoupled.md
new file mode 100644
index 000000000000..2fcde7e96e8c
--- /dev/null
+++ b/framework/doc/content/source/userobject/ActivateElementsCoupled.md
@@ -0,0 +1,9 @@
+# ActivateElementsCoupled
+
+This user object uses the +coupled variable value+ as the metric to activate (add) an element by moving the element from an "inactive" subdomain to the "active" subdomain. It uses a coupled variable to decide whether to activate an element. The `coupled_var`, `activate_value` and the `activate_type` needs to be provided in the input to form the activation criterion. By default, the element is activated if the averaged value of the coupled variable in the element is `above` the `activate_value`. User can set `activate_type = 'below'` or `'equal'` to activate element when the averaged coupled variable value is below or equal to the `activate_value`.
+
+!syntax parameters /UserObjects/ActivateElementsCoupled
+
+!syntax inputs /UserObjects/ActivateElementsCoupled
+
+!syntax children /UserObjects/ActivateElementsCoupled
diff --git a/framework/doc/content/source/userobject/ActivateElementsUserObjectBase.md b/framework/doc/content/source/userobject/ActivateElementsUserObjectBase.md
new file mode 100644
index 000000000000..6161b1e1e93f
--- /dev/null
+++ b/framework/doc/content/source/userobject/ActivateElementsUserObjectBase.md
@@ -0,0 +1,5 @@
+# ActivateElementsUserObjectBase
+
+The `ActivateElementsUserObjectBase` class is for activating (adding) an element by moving the element from an "inactive" subdomain to the "active" subdomain. There are two metrics for activating an element, i.e., by +function path+ and by +coupled variable value+, which are implemented in the classes [ActivateElementsByPath](/ActivateElementsByPath.md) and [ActivateElementsCoupled](/ActivateElementsCoupled.md), respectively.
+
+This user object updates the boundary information due to the addition of elements across subdomains. User can save this boundary information by passing the changed boundary name to  `expand_boundary_name` in the input block.  Note that current implementation does not de-activate an element once the element is activated.
diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h
index 6f0e22ef57c2..6ef952da48e0 100644
--- a/framework/include/problems/FEProblemBase.h
+++ b/framework/include/problems/FEProblemBase.h
@@ -694,6 +694,14 @@ class FEProblemBase : public SubProblem, public Restartable
 
   void projectSolution();
 
+  /**
+   * Project initial conditions for custom \p elem_range and \p bnd_node_range
+   * This is needed when elements/boundary nodes are added to a specific subdomain
+   * at an intermediate step
+   */
+  void projectInitialConditionOnCustomRange(ConstElemRange & elem_range,
+                                            ConstBndNodeRange & bnd_node_range);
+
   // Materials /////
   virtual void addMaterial(const std::string & kernel_name,
                            const std::string & name,
@@ -1450,6 +1458,13 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   void notifyWhenMeshChanges(MeshChangedInterface * mci);
 
+  /**
+   * Initialize stateful properties for elements in a specific \p elem_range
+   * This is needed when elements/boundary nodes are added to a specific subdomain
+   * at an intermediate step
+   */
+  void initElementStatefulProps(const ConstElemRange & elem_range);
+
   /**
    * Method called to perform a series of sanity checks before a simulation is run. This method
    * doesn't return when errors are found, instead it generally calls mooseError() directly.
diff --git a/framework/include/userobject/ActivateElementsByPath.h b/framework/include/userobject/ActivateElementsByPath.h
new file mode 100644
index 000000000000..701433cb4d3b
--- /dev/null
+++ b/framework/include/userobject/ActivateElementsByPath.h
@@ -0,0 +1,32 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "ActivateElementsUserObjectBase.h"
+#include "Function.h"
+
+class ActivateElementsByPath : public ActivateElementsUserObjectBase
+{
+public:
+  static InputParameters validParams();
+
+  ActivateElementsByPath(const InputParameters & parameters);
+
+  virtual bool isElementActivated() override;
+
+protected:
+  /// path of the heat source, x, y, z components
+  const Function * _function_x;
+  const Function * _function_y;
+  const Function * _function_z;
+  /// define the distance of the element to the point on the path,
+  /// below which the element will be activated
+  const Real _activate_distance;
+};
diff --git a/framework/include/userobject/ActivateElementsCoupled.h b/framework/include/userobject/ActivateElementsCoupled.h
new file mode 100644
index 000000000000..12391e69e606
--- /dev/null
+++ b/framework/include/userobject/ActivateElementsCoupled.h
@@ -0,0 +1,31 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "ActivateElementsUserObjectBase.h"
+#include "MooseEnum.h"
+
+class ActivateElementsCoupled : public ActivateElementsUserObjectBase
+{
+public:
+  static InputParameters validParams();
+
+  ActivateElementsCoupled(const InputParameters & parameters);
+
+  virtual bool isElementActivated() override;
+
+protected:
+  /// variable value to decide wether an element whould be activated
+  const VariableValue & _coupled_var;
+  /// variable value to decide wether an element whould be activated
+  const Real _activate_value;
+  /// type of activation - blow or above
+  const enum class ActivateType { BELOW, EQUAL, ABOVE } _activate_type;
+};
diff --git a/framework/include/userobject/ActivateElementsUserObjectBase.h b/framework/include/userobject/ActivateElementsUserObjectBase.h
new file mode 100644
index 000000000000..211df3436e39
--- /dev/null
+++ b/framework/include/userobject/ActivateElementsUserObjectBase.h
@@ -0,0 +1,96 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "ElementUserObject.h"
+#include "NonlinearSystemBase.h"
+#include "AuxiliarySystem.h"
+
+class ActivateElementsUserObjectBase : public ElementUserObject
+{
+public:
+  static InputParameters validParams();
+
+  ActivateElementsUserObjectBase(const InputParameters & parameters);
+
+  const std::set<dof_id_type> & getNewlyActivatedElements() const { return _newly_activated_elem; };
+
+  BoundaryID getExpandedBoundaryID()
+  {
+    mooseAssert(!_boundary_ids.empty(), "Boundary ID is empty");
+    return _boundary_ids[0];
+  }
+
+  virtual bool isElementActivated() = 0;
+
+  void initialize() override{};
+  void execute() override;
+  void threadJoin(const UserObject & /*uo*/) override{};
+  void finalize() override;
+
+protected:
+  void setNewBoundayName();
+
+  void updateBoundaryInfo(MooseMesh & mesh);
+
+  void push_boundary_side_info(
+      MooseMesh & mesh,
+      std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned int>>> &
+          elems_to_push);
+
+  void push_boundary_node_info(
+      MooseMesh & mesh,
+      std::unordered_map<processor_id_type, std::vector<dof_id_type>> & nodes_to_push);
+  /**
+   * Initialize solutions for the nodes
+   */
+  void initSolutions(ConstElemRange & elem_range, ConstBndNodeRange & bnd_node_range);
+  /**
+   * Returns true if all the connected elements are in the _newly_activated_elem
+   */
+  bool isNewlyActivated(const Node * node);
+
+  void getNodesToRemoveFromBnd(std::set<dof_id_type> & remove_set, std::set<dof_id_type> & add_set);
+
+  void insertNodeIdsOnSide(const Elem * ele,
+                           const unsigned short int side,
+                           std::set<dof_id_type> & node_ids);
+
+  /**
+   * Get ranges for use with threading.
+   */
+  ConstElemRange * getNewlyActivatedElementRange();
+  ConstBndNodeRange * getNewlyActivatedBndNodeRange();
+  ConstNodeRange * getNewlyActivatedNodeRange();
+
+  std::set<dof_id_type> _newly_activated_elem;
+  std::set<dof_id_type> _newly_activated_node;
+  /**
+   * Somes nodes are to be removed from the boundary
+   * when adding/removing sides
+   */
+  std::set<dof_id_type> _node_to_remove_from_bnd;
+
+  /**
+   * Ranges for use with threading.
+   */
+  std::unique_ptr<ConstElemRange> _activated_elem_range;
+  std::unique_ptr<ConstBndNodeRange> _activated_bnd_node_range;
+  std::unique_ptr<ConstNodeRange> _activated_node_range;
+
+  /// activate subdomain ID
+  const subdomain_id_type _active_subdomain_id;
+  /// inactivate subdomain ID (the subdomain that you want to keep the same)
+  const subdomain_id_type _inactive_subdomain_id;
+  /// expanded boundary name
+  const std::vector<BoundaryName> _expand_boundary_name;
+  /// expanded boundary IDs
+  std::vector<BoundaryID> _boundary_ids, _disp_boundary_ids;
+};
diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C
index b3214ea448ce..70d16f3da587 100644
--- a/framework/src/problems/FEProblemBase.C
+++ b/framework/src/problems/FEProblemBase.C
@@ -2885,6 +2885,54 @@ FEProblemBase::projectSolution()
   _aux->solution().localize(*_aux->sys().current_local_solution, _aux->dofMap().get_send_list());
 }
 
+void
+FEProblemBase::projectInitialConditionOnCustomRange(ConstElemRange & elem_range,
+                                                    ConstBndNodeRange & bnd_nodes)
+{
+  ComputeInitialConditionThread cic(*this);
+  Threads::parallel_reduce(elem_range, cic);
+
+  // Need to close the solution vector here so that boundary ICs take precendence
+  _nl->solution().close();
+  _aux->solution().close();
+
+  ComputeBoundaryInitialConditionThread cbic(*this);
+  Threads::parallel_reduce(bnd_nodes, cbic);
+
+  _nl->solution().close();
+  _aux->solution().close();
+
+  // Also, load values into the SCALAR dofs
+  // Note: We assume that all SCALAR dofs are on the
+  // processor with highest ID
+  if (processor_id() == (n_processors() - 1) && _scalar_ics.hasActiveObjects())
+  {
+    const auto & ics = _scalar_ics.getActiveObjects();
+    for (const auto & ic : ics)
+    {
+      MooseVariableScalar & var = ic->variable();
+      var.reinit();
+
+      DenseVector<Number> vals(var.order());
+      ic->compute(vals);
+
+      const unsigned int n_SCALAR_dofs = var.dofIndices().size();
+      for (unsigned int i = 0; i < n_SCALAR_dofs; i++)
+      {
+        const dof_id_type global_index = var.dofIndices()[i];
+        var.sys().solution().set(global_index, vals(i));
+        var.setValue(i, vals(i));
+      }
+    }
+  }
+
+  _nl->solution().close();
+  _nl->solution().localize(*_nl->system().current_local_solution, _nl->dofMap().get_send_list());
+
+  _aux->solution().close();
+  _aux->solution().localize(*_aux->sys().current_local_solution, _aux->dofMap().get_send_list());
+}
+
 std::shared_ptr<MaterialBase>
 FEProblemBase::getMaterial(std::string name,
                            Moose::MaterialDataType type,
@@ -6182,6 +6230,20 @@ FEProblemBase::notifyWhenMeshChanges(MeshChangedInterface * mci)
   _notify_when_mesh_changes.push_back(mci);
 }
 
+void
+FEProblemBase::initElementStatefulProps(const ConstElemRange & elem_range)
+{
+  ComputeMaterialsObjectThread cmt(*this,
+                                   _material_data,
+                                   _bnd_material_data,
+                                   _neighbor_material_data,
+                                   _material_props,
+                                   _bnd_material_props,
+                                   _neighbor_material_props,
+                                   _assembly);
+  cmt(elem_range, true);
+}
+
 void
 FEProblemBase::checkProblemIntegrity()
 {
diff --git a/framework/src/userobject/ActivateElementsByPath.C b/framework/src/userobject/ActivateElementsByPath.C
new file mode 100644
index 000000000000..7132f3234f9e
--- /dev/null
+++ b/framework/src/userobject/ActivateElementsByPath.C
@@ -0,0 +1,53 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ActivateElementsByPath.h"
+#include "libmesh/point.h"
+
+registerMooseObject("MooseApp", ActivateElementsByPath);
+
+InputParameters
+ActivateElementsByPath::validParams()
+{
+  InputParameters params = ActivateElementsUserObjectBase::validParams();
+
+  params.addParam<FunctionName>(
+      "function_x", "0", "The x component of the heating spot travel path");
+  params.addParam<FunctionName>(
+      "function_y", "0", "The y component of the heating spot travel path");
+  params.addParam<FunctionName>(
+      "function_z", "0", "The z component of the heating spot travel path");
+  params.addParam<Real>("activate_distance",
+                        1e-4,
+                        "The maximum distance of the activated element to the point on the path.");
+  return params;
+}
+
+ActivateElementsByPath::ActivateElementsByPath(const InputParameters & parameters)
+  : ActivateElementsUserObjectBase(parameters),
+    _function_x(isParamValid("function_x") ? &getFunction("function_x") : nullptr),
+    _function_y(isParamValid("function_y") ? &getFunction("function_y") : nullptr),
+    _function_z(isParamValid("function_z") ? &getFunction("function_z") : nullptr),
+    _activate_distance(
+        declareRestartableData<Real>("activate_distance", getParam<Real>("activate_distance")))
+{
+}
+
+bool
+ActivateElementsByPath::isElementActivated()
+{
+  // activate center (assume position of the activate center is only time dependent)
+  const static Point dummy;
+  Real x_t = _function_x ? _function_x->value(_t, dummy) : 0.0;
+  Real y_t = _function_y ? _function_y->value(_t, dummy) : 0.0;
+  Real z_t = _function_z ? _function_z->value(_t, dummy) : 0.0;
+
+  // activate element when element is close to the point
+  return _current_elem->close_to_point(Point(x_t, y_t, z_t), _activate_distance);
+}
diff --git a/framework/src/userobject/ActivateElementsCoupled.C b/framework/src/userobject/ActivateElementsCoupled.C
new file mode 100644
index 000000000000..95fad4687a2c
--- /dev/null
+++ b/framework/src/userobject/ActivateElementsCoupled.C
@@ -0,0 +1,64 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ActivateElementsCoupled.h"
+
+registerMooseObject("MooseApp", ActivateElementsCoupled);
+
+InputParameters
+ActivateElementsCoupled::validParams()
+{
+  InputParameters params = ActivateElementsUserObjectBase::validParams();
+
+  params.addRequiredParam<Real>("activate_value", "The value above which to activate the element");
+  params.addRequiredCoupledVar(
+      "coupled_var",
+      "The variable value will be used to decide wether an element whould be activated.");
+  params.addParam<MooseEnum>("activate_type",
+                             MooseEnum("below equal above", "above"),
+                             "Activate element when below or above the activate_value");
+  return params;
+}
+
+ActivateElementsCoupled::ActivateElementsCoupled(const InputParameters & parameters)
+  : ActivateElementsUserObjectBase(parameters),
+    _coupled_var(coupledValue("coupled_var")),
+    _activate_value(
+        declareRestartableData<Real>("activate_value", getParam<Real>("activate_value"))),
+    _activate_type(getParam<MooseEnum>("activate_type").getEnum<ActivateType>())
+{
+}
+
+bool
+ActivateElementsCoupled::isElementActivated()
+{
+  bool is_activated = false;
+  Real avg_val = 0.0;
+
+  for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
+    avg_val += _coupled_var[qp];
+  avg_val /= _qrule->n_points();
+
+  switch (_activate_type)
+  {
+    case ActivateType::BELOW:
+      is_activated = (avg_val < _activate_value);
+      break;
+
+    case ActivateType::EQUAL:
+      is_activated = MooseUtils::absoluteFuzzyEqual(avg_val - _activate_value, 0.0);
+      break;
+
+    case ActivateType::ABOVE:
+      is_activated = (avg_val > _activate_value);
+      break;
+  }
+
+  return is_activated;
+}
diff --git a/framework/src/userobject/ActivateElementsUserObjectBase.C b/framework/src/userobject/ActivateElementsUserObjectBase.C
new file mode 100644
index 000000000000..eaa65b8f7d1c
--- /dev/null
+++ b/framework/src/userobject/ActivateElementsUserObjectBase.C
@@ -0,0 +1,463 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ActivateElementsUserObjectBase.h"
+#include "DisplacedProblem.h"
+
+#include "libmesh/quadrature.h"
+#include "libmesh/parallel_algebra.h"
+#include "libmesh/parallel.h"
+#include "libmesh/point.h"
+#include "libmesh/dof_map.h"
+
+#include "libmesh/parallel_ghost_sync.h"
+#include "libmesh/mesh_communication.h"
+
+InputParameters
+ActivateElementsUserObjectBase::validParams()
+{
+  InputParameters params = ElementUserObject::validParams();
+  params.addClassDescription("Determine activated elements.");
+  params.addRequiredParam<subdomain_id_type>("active_subdomain_id", "The active subdomain ID.");
+  params.addParam<subdomain_id_type>(
+      "inactive_subdomain_id",
+      libMesh::invalid_uint,
+      "The inactivate subdomain ID, i.e., the subdomain that you want to keep the same.");
+  params.addRequiredParam<std::vector<BoundaryName>>("expand_boundary_name",
+                                                     "The expanded boundary name.");
+  return params;
+}
+
+ActivateElementsUserObjectBase::ActivateElementsUserObjectBase(const InputParameters & parameters)
+  : ElementUserObject(parameters),
+    _active_subdomain_id(declareRestartableData<subdomain_id_type>(
+        "active_subdomain_id", getParam<subdomain_id_type>("active_subdomain_id"))),
+    _inactive_subdomain_id(declareRestartableData<subdomain_id_type>(
+        "inactive_subdomain_id", getParam<subdomain_id_type>("inactive_subdomain_id"))),
+    _expand_boundary_name(getParam<std::vector<BoundaryName>>("expand_boundary_name"))
+{
+  setNewBoundayName();
+}
+
+void
+ActivateElementsUserObjectBase::setNewBoundayName()
+{
+  mooseAssert(!_expand_boundary_name.empty(), "Expanded boundary name is empty");
+  // add the new boundary and get its boundary id
+  _boundary_ids = _mesh.getBoundaryIDs(_expand_boundary_name, true);
+  mooseAssert(!_boundary_ids.empty(), "Boundary ID is empty.");
+  _mesh.setBoundaryName(_boundary_ids[0], _expand_boundary_name[0]);
+  _mesh.getMesh().get_boundary_info().sideset_name(_boundary_ids[0]) = _expand_boundary_name[0];
+  _mesh.getMesh().get_boundary_info().nodeset_name(_boundary_ids[0]) = _expand_boundary_name[0];
+
+  auto displaced_problem = _fe_problem.getDisplacedProblem();
+  if (displaced_problem)
+  {
+    _disp_boundary_ids = displaced_problem->mesh().getBoundaryIDs(_expand_boundary_name, true);
+    mooseAssert(!_disp_boundary_ids.empty(), "Boundary ID in the displaced mesh is empty");
+
+    displaced_problem->mesh().setBoundaryName(_disp_boundary_ids[0], _expand_boundary_name[0]);
+    displaced_problem->mesh().getMesh().get_boundary_info().sideset_name(_disp_boundary_ids[0]) =
+        _expand_boundary_name[0];
+    displaced_problem->mesh().getMesh().get_boundary_info().nodeset_name(_disp_boundary_ids[0]) =
+        _expand_boundary_name[0];
+  }
+}
+
+void
+ActivateElementsUserObjectBase::execute()
+{
+  if (isElementActivated() && _current_elem->subdomain_id() != _active_subdomain_id &&
+      _current_elem->subdomain_id() != _inactive_subdomain_id)
+  {
+    /*
+      _current_elem subdomain id is not assignable
+      create a copy of this element from MooseMesh
+    */
+    dof_id_type ele_id = _current_elem->id();
+    Elem * ele = _mesh.elemPtr(ele_id);
+
+    // Add element to the activate subdomain
+    ele->subdomain_id() = _active_subdomain_id;
+
+    //  Reassign element in the reference mesh while using a displaced mesh
+    auto displaced_problem = _fe_problem.getDisplacedProblem();
+    if (displaced_problem)
+    {
+      Elem * disp_ele = displaced_problem->mesh().elemPtr(ele_id);
+      disp_ele->subdomain_id() = _active_subdomain_id;
+    }
+
+    // Save the newly activated element id and node for updating boundary info later
+    _newly_activated_elem.insert(ele_id);
+    for (unsigned int i = 0; i < ele->n_nodes(); ++i)
+      _newly_activated_node.insert(ele->node_id(i));
+  }
+}
+
+void
+ActivateElementsUserObjectBase::finalize()
+{
+  /*
+    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);
+  // Update boundary info
+  updateBoundaryInfo(_mesh);
+
+  // Similarly for the displaced mesh
+  auto displaced_problem = _fe_problem.getDisplacedProblem();
+  if (displaced_problem)
+  {
+    SyncSubdomainIds sync_mesh(displaced_problem->mesh().getMesh());
+    Parallel::sync_dofobject_data_by_id(displaced_problem->mesh().getMesh().comm(),
+                                        displaced_problem->mesh().getMesh().elements_begin(),
+                                        displaced_problem->mesh().getMesh().elements_end(),
+                                        sync_mesh);
+    updateBoundaryInfo(displaced_problem->mesh());
+  }
+
+  // Reinit equation systems
+  _fe_problem.meshChanged();
+
+  // Get storage ranges for the newly activated elements and boundary nodes
+  ConstElemRange & elem_range = *this->getNewlyActivatedElementRange();
+  ConstBndNodeRange & bnd_node_range = *this->getNewlyActivatedBndNodeRange();
+
+  // Apply initial condition for the newly activated elements
+  initSolutions(elem_range, bnd_node_range);
+
+  //  Initialize stateful material properties for the newly activated elements
+  _fe_problem.initElementStatefulProps(elem_range);
+
+  //  Clear the list
+  _newly_activated_elem.clear();
+  _newly_activated_node.clear();
+
+  _node_to_remove_from_bnd.clear();
+}
+
+void
+ActivateElementsUserObjectBase::getNodesToRemoveFromBnd(std::set<dof_id_type> & remove_set,
+                                                        std::set<dof_id_type> & add_set)
+{
+  // get the difference between the remove_set and the add_set,
+  // save the difference in _node_to_remove_from_bnd
+  int sz = remove_set.size() + add_set.size();
+  std::vector<dof_id_type> v(sz);
+  std::vector<dof_id_type>::iterator it = std::set_difference(
+      remove_set.begin(), remove_set.end(), add_set.begin(), add_set.end(), v.begin());
+  v.resize(it - v.begin());
+  _node_to_remove_from_bnd.clear();
+  for (auto id : v)
+    _node_to_remove_from_bnd.insert(id);
+}
+
+void
+ActivateElementsUserObjectBase::insertNodeIdsOnSide(const Elem * ele,
+                                                    const unsigned short int side,
+                                                    std::set<dof_id_type> & node_ids)
+{
+  for (unsigned int i = 0; i < ele->side_ptr(side)->n_nodes(); ++i)
+    node_ids.insert(ele->side_ptr(side)->node_id(i));
+}
+
+void
+ActivateElementsUserObjectBase::updateBoundaryInfo(MooseMesh & mesh)
+{
+  // save the removed ghost sides and associated nodes to sync across processors
+  std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned int>>>
+      ghost_sides_to_remove;
+  std::unordered_map<processor_id_type, std::vector<dof_id_type>> ghost_nodes_to_remove;
+
+  // save nodes are added and removed
+  std::set<dof_id_type> add_nodes, remove_nodes;
+
+  for (auto ele_id : _newly_activated_elem)
+  {
+    Elem * ele = mesh.elemPtr(ele_id);
+    for (auto s : ele->side_index_range())
+    {
+      Elem * neighbor_ele = ele->neighbor_ptr(s);
+      if (neighbor_ele == nullptr)
+      {
+        // add this side to boundary
+        mesh.getMesh().get_boundary_info().add_side(ele, s, _boundary_ids[0]);
+        insertNodeIdsOnSide(ele, s, add_nodes);
+      }
+      else
+      {
+        if (neighbor_ele->subdomain_id() != _active_subdomain_id &&
+            neighbor_ele->subdomain_id() != _inactive_subdomain_id)
+        {
+          // add this side to boundary
+          mesh.getMesh().get_boundary_info().add_side(ele, s, _boundary_ids[0]);
+          insertNodeIdsOnSide(ele, s, add_nodes);
+        }
+        else
+        {
+          // remove this side from the boundary
+          mesh.getMesh().get_boundary_info().remove_side(ele, s);
+          insertNodeIdsOnSide(ele, s, remove_nodes);
+
+          // remove the neighbor side from the boundary
+          unsigned int neighbor_s = neighbor_ele->which_neighbor_am_i(ele);
+          mesh.getMesh().get_boundary_info().remove_side(neighbor_ele, neighbor_s);
+          insertNodeIdsOnSide(neighbor_ele, neighbor_s, remove_nodes);
+
+          if (neighbor_ele->processor_id() != this->processor_id())
+            ghost_sides_to_remove[neighbor_ele->processor_id()].emplace_back(neighbor_ele->id(),
+                                                                             neighbor_s);
+        }
+      }
+    }
+  }
+  // make sure to remove nodes that are not in the add list
+  getNodesToRemoveFromBnd(remove_nodes, add_nodes);
+  for (auto node_id : _node_to_remove_from_bnd)
+    mesh.getMesh().get_boundary_info().remove_node(mesh.nodePtr(node_id), _boundary_ids[0]);
+
+  // synchronize boundary information across processors
+  push_boundary_side_info(mesh, ghost_sides_to_remove);
+  push_boundary_node_info(mesh, ghost_nodes_to_remove);
+  mesh.getMesh().get_boundary_info().parallel_sync_side_ids();
+  mesh.getMesh().get_boundary_info().parallel_sync_node_ids();
+  mesh.update();
+}
+
+void
+ActivateElementsUserObjectBase::push_boundary_side_info(
+    MooseMesh & mesh,
+    std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned int>>> &
+        elems_to_push)
+{
+  auto elem_action_functor =
+      [&mesh, this](processor_id_type,
+                    const std::vector<std::pair<dof_id_type, unsigned int>> & received_elem) {
+        // remove the side
+        for (const auto & pr : received_elem)
+          mesh.getMesh().get_boundary_info().remove_side(
+              mesh.getMesh().elem_ptr(pr.first), pr.second, this->getExpandedBoundaryID());
+      };
+
+  Parallel::push_parallel_vector_data(
+      mesh.getMesh().get_boundary_info().comm(), elems_to_push, elem_action_functor);
+}
+
+void
+ActivateElementsUserObjectBase::push_boundary_node_info(
+    MooseMesh & mesh,
+    std::unordered_map<processor_id_type, std::vector<dof_id_type>> & nodes_to_push)
+{
+  auto node_action_functor = [&mesh, this](processor_id_type,
+                                           const std::vector<dof_id_type> & received_nodes) {
+    for (const auto & pr : received_nodes)
+    {
+      // remove the node
+      mesh.getMesh().get_boundary_info().remove_node(mesh.getMesh().node_ptr(pr),
+                                                     this->getExpandedBoundaryID());
+    }
+  };
+
+  Parallel::push_parallel_vector_data(
+      mesh.getMesh().get_boundary_info().comm(), nodes_to_push, node_action_functor);
+}
+
+ConstElemRange *
+ActivateElementsUserObjectBase::getNewlyActivatedElementRange()
+{
+  // deletes the object first
+  _activated_elem_range.reset();
+
+  // create a vector of the newly activated elements
+  std::vector<Elem *> elems;
+  for (auto elem_id : _newly_activated_elem)
+    elems.push_back(_mesh.elemPtr(elem_id));
+
+  // Make some fake element iterators defining this vector of
+  // elements
+  Elem * const * elempp = const_cast<Elem * const *>(elems.data());
+  Elem * const * elemend = elempp + elems.size();
+
+  const auto elems_begin =
+      MeshBase::const_element_iterator(elempp, elemend, Predicates::NotNull<Elem * const *>());
+
+  const auto elems_end =
+      MeshBase::const_element_iterator(elemend, elemend, Predicates::NotNull<Elem * const *>());
+  if (!_activated_elem_range)
+    _activated_elem_range = libmesh_make_unique<ConstElemRange>(elems_begin, elems_end);
+
+  return _activated_elem_range.get();
+}
+
+ConstBndNodeRange *
+ActivateElementsUserObjectBase::getNewlyActivatedBndNodeRange()
+{
+  // deletes the object first
+  _activated_bnd_node_range.reset();
+
+  // create a vector of the newly activated nodes
+  std::vector<const BndNode *> nodes;
+  std::set<const BndNode *> set_nodes;
+  ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
+  for (auto & bnode : bnd_nodes)
+  {
+    dof_id_type bnode_id = bnode->_node->id();
+    auto it = _newly_activated_node.find(bnode_id);
+    if (it != _newly_activated_node.end())
+      set_nodes.insert(bnode);
+  }
+
+  nodes.assign(set_nodes.begin(), set_nodes.end());
+
+  // Make some fake element iterators defining this vector of
+  // nodes
+  BndNode * const * nodepp = const_cast<BndNode * const *>(nodes.data());
+  BndNode * const * nodeend = nodepp + nodes.size();
+
+  const auto nodes_begin =
+      MooseMesh::const_bnd_node_iterator(nodepp, nodeend, Predicates::NotNull<BndNode * const *>());
+
+  const auto nodes_end = MooseMesh::const_bnd_node_iterator(
+      nodeend, nodeend, Predicates::NotNull<BndNode * const *>());
+
+  if (!_activated_bnd_node_range)
+    _activated_bnd_node_range = libmesh_make_unique<ConstBndNodeRange>(nodes_begin, nodes_end);
+
+  return _activated_bnd_node_range.get();
+}
+
+ConstNodeRange *
+ActivateElementsUserObjectBase::getNewlyActivatedNodeRange()
+{
+  // deletes the object first
+  _activated_node_range.reset();
+
+  // create a vector of the newly activated nodes
+  std::vector<const Node *> nodes;
+  for (auto elem_id : _newly_activated_elem)
+  {
+    const Node * const * elem_nodes = _mesh.elemPtr(elem_id)->get_nodes();
+    unsigned int n_nodes = _mesh.elemPtr(elem_id)->n_nodes();
+    for (unsigned int n = 0; n < n_nodes; ++n)
+    {
+      // check if all the elements connected to this node are newly activated
+      const Node * nd = elem_nodes[n];
+      if (isNewlyActivated(nd))
+        nodes.push_back(nd);
+    }
+  }
+
+  // Make some fake node iterators defining this vector of
+  // nodes
+  Node * const * nodepp = const_cast<Node * const *>(nodes.data());
+  Node * const * nodeend = nodepp + nodes.size();
+
+  const auto nodes_begin =
+      MeshBase::const_node_iterator(nodepp, nodeend, Predicates::NotNull<Node * const *>());
+
+  const auto nodes_end =
+      MeshBase::const_node_iterator(nodeend, nodeend, Predicates::NotNull<Node * const *>());
+
+  if (!_activated_node_range)
+    _activated_node_range = libmesh_make_unique<ConstNodeRange>(nodes_begin, nodes_end);
+
+  return _activated_node_range.get();
+}
+
+bool
+ActivateElementsUserObjectBase::isNewlyActivated(const Node * nd)
+{
+  const auto & node_to_elem_map = _mesh.nodeToElemMap();
+  auto node_to_elem_pair = node_to_elem_map.find(nd->id());
+  if (node_to_elem_pair != node_to_elem_map.end())
+  {
+    const std::vector<dof_id_type> & connected_ele_ids = node_to_elem_pair->second;
+    for (auto connected_ele_id : connected_ele_ids)
+    {
+      // check the connected elements
+      if (_mesh.elemPtr(connected_ele_id)->subdomain_id() == _inactive_subdomain_id)
+        return false;
+      if (_mesh.elemPtr(connected_ele_id)->subdomain_id() == _active_subdomain_id &&
+          std::find(_newly_activated_elem.begin(), _newly_activated_elem.end(), connected_ele_id) ==
+              _newly_activated_elem.end())
+        return false;
+    }
+  }
+  return true;
+}
+
+void
+ActivateElementsUserObjectBase::initSolutions(ConstElemRange & elem_range,
+                                              ConstBndNodeRange & bnd_node_range)
+{
+  // project initial condition to the current solution
+  _fe_problem.projectInitialConditionOnCustomRange(elem_range, bnd_node_range);
+
+  NumericVector<Number> & current_solution =
+      *_fe_problem.getNonlinearSystemBase().system().current_local_solution;
+  NumericVector<Number> & old_solution = _fe_problem.getNonlinearSystemBase().solutionOld();
+  NumericVector<Number> & older_solution = _fe_problem.getNonlinearSystemBase().solutionOlder();
+
+  NumericVector<Number> & current_aux_solution =
+      *_fe_problem.getAuxiliarySystem().system().current_local_solution;
+  NumericVector<Number> & old_aux_solution = _fe_problem.getAuxiliarySystem().solutionOld();
+  NumericVector<Number> & older_aux_solution = _fe_problem.getAuxiliarySystem().solutionOlder();
+
+  DofMap & dof_map = _fe_problem.getNonlinearSystemBase().dofMap();
+  DofMap & dof_map_aux = _fe_problem.getAuxiliarySystem().dofMap();
+
+  std::set<dof_id_type> dofs, dofs_aux;
+  // get dofs for the newly added elements
+  for (auto & elem : elem_range)
+  {
+    std::vector<dof_id_type> di, di_aux;
+    dof_map.dof_indices(elem, di);
+    dof_map_aux.dof_indices(elem, di_aux);
+    for (unsigned int i = 0; i < di.size(); ++i)
+      dofs.insert(di[i]);
+    for (unsigned int i = 0; i < di_aux.size(); ++i)
+      dofs_aux.insert(di_aux[i]);
+
+    di.clear();
+    di_aux.clear();
+  }
+
+  // update solutions
+  for (auto dof : dofs)
+  {
+    old_solution.set(dof, current_solution(dof));
+    older_solution.set(dof, current_solution(dof));
+  }
+  // update aux solutions
+  for (auto dof_aux : dofs_aux)
+  {
+    old_aux_solution.set(dof_aux, current_aux_solution(dof_aux));
+    older_aux_solution.set(dof_aux, current_aux_solution(dof_aux));
+  }
+
+  dofs.clear();
+  dofs_aux.clear();
+
+  current_solution.close();
+  old_solution.close();
+  older_solution.close();
+
+  current_aux_solution.close();
+  old_aux_solution.close();
+  older_aux_solution.close();
+
+  _fe_problem.restoreSolutions();
+}
diff --git a/modules/combined/test/tests/additive_manufacturing/check_element_addition.i b/modules/combined/test/tests/additive_manufacturing/check_element_addition.i
new file mode 100644
index 000000000000..36320dd23086
--- /dev/null
+++ b/modules/combined/test/tests/additive_manufacturing/check_element_addition.i
@@ -0,0 +1,107 @@
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Mesh]
+  [./gen]
+    type = GeneratedMeshGenerator
+    dim = 3
+    xmin =0
+    xmax =10
+    ymin =0
+    ymax =10
+    zmin =0
+    zmax =0.5
+    nx=20
+    ny=20
+    nz=1
+  [../]
+  [./left_domain]
+    input = gen
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '0 0 0'
+    top_right = '5 10 0.5'
+    block_id = 1
+  [../]
+  [./right_domain]
+    input = left_domain
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '5 0 0'
+    top_right = '10 10 0.5'
+    block_id = 2
+  [../]
+  [./sidesets]
+    input = right_domain
+    type = SideSetsAroundSubdomainGenerator
+    normal = '1 0 0'
+    block = 1
+    new_boundary = 'moving_interface'
+  []
+[]
+
+[Variables]
+  [./temp]
+    block = '1'
+  [../]
+[]
+
+[Functions]
+  [./fx]
+    type = ParsedFunction
+    value= '5.25'
+  [../]
+  [./fy]
+    type = ParsedFunction
+    value= '2.5*t'
+  [../]
+  [./fz]
+    type = ParsedFunction
+    value= '0.25'
+  [../]
+[]
+
+[Preconditioning]
+  [./smp]
+    type = SMP
+    full = true
+  [../]
+[]
+
+
+[Executioner]
+  type = Transient
+
+  automatic_scaling = true
+
+  solve_type = 'NEWTON'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  line_search = 'none'
+
+  l_max_its = 10
+  nl_max_its = 20
+  nl_rel_tol = 1e-4
+
+  start_time = 0.0
+  end_time = 1.0
+  dt = 1e-1
+  dtmin = 1e-4
+[]
+
+[UserObjects]
+  [./activated_elem_uo]
+    type = ActivateElementsByPath
+    execute_on = timestep_begin
+    function_x= fx
+    function_y= fy
+    function_z= fz
+    active_subdomain_id = 1
+    expand_boundary_name = 'moving_interface'
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/modules/combined/test/tests/additive_manufacturing/check_element_addition_by_variable.i b/modules/combined/test/tests/additive_manufacturing/check_element_addition_by_variable.i
new file mode 100644
index 000000000000..f438f05c6ca1
--- /dev/null
+++ b/modules/combined/test/tests/additive_manufacturing/check_element_addition_by_variable.i
@@ -0,0 +1,149 @@
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Mesh]
+  [./gen]
+    type = GeneratedMeshGenerator
+    dim = 3
+    xmin =0
+    xmax =2.0
+    ymin =0
+    ymax =2.0
+    zmin =0
+    zmax =2.0
+    nx=10
+    ny=10
+    nz=10
+  [../]
+  [./left_domain]
+    input = gen
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '0 0 0'
+    top_right = '2 2 1'
+    block_id = 1
+  [../]
+  [./right_domain]
+    input = left_domain
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '0 0 1'
+    top_right = '2 2 2'
+    block_id = 2
+  [../]
+  [./sidesets]
+    input = right_domain
+    type = SideSetsAroundSubdomainGenerator
+    normal = '0 0 1'
+    block = 1
+    new_boundary = 'moving_interface'
+  []
+[]
+
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+  block = '1 2'
+[]
+
+[Modules/TensorMechanics/Master]
+  [./all]
+    # strain = FINITE
+    add_variables = true
+    generate_output = 'stress_zz strain_zz'
+    block = '1 2'
+    use_automatic_differentiation = true
+  [../]
+[]
+
+[Materials]
+  [./elasticity]
+    type = ADComputeVariableIsotropicElasticityTensor
+    poissons_ratio = 0.3
+    youngs_modulus = 1e3
+    block = '1 2'
+  [../]
+  [./stress]
+    type = ADComputeLinearElasticStress
+    block = '1 2'
+  [../]
+[]
+
+[Functions]
+   [./front_pull]
+     type = PiecewiseLinear
+     x = '0 1'
+     y = '0 1'
+     scale_factor = 0.5
+   [../]
+ []
+
+[BCs]
+  [./disp_front_pull]
+    type  = ADFunctionDirichletBC
+    variable = disp_z
+    boundary = front
+    function = front_pull
+  [../]
+  [./uz_back_fix]
+     type = ADDirichletBC
+     variable = disp_z
+     boundary = back
+     value = 0.0
+   [../]
+   [./u_yz_fix]
+     type = ADDirichletBC
+     variable = disp_x
+     boundary = left
+     value = 0.0
+   [../]
+   [./u_xz_fix]
+     type = ADDirichletBC
+     variable = disp_y
+     boundary = bottom
+     value = 0.0
+   [../]
+[]
+
+[Preconditioning]
+  [./smp]
+    type = SMP
+    full = true
+  [../]
+[]
+
+
+[Executioner]
+  type = Transient
+
+  automatic_scaling = true
+
+  solve_type = 'NEWTON'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  line_search = 'none'
+
+  l_max_its = 10
+  nl_max_its = 20
+  nl_rel_tol = 1e-4
+
+  start_time = 0.0
+  end_time = 1.0
+  dt = 1e-1
+  dtmin = 1e-4
+[]
+
+[UserObjects]
+  [./activated_elem_uo]
+    type = ActivateElementsCoupled
+    execute_on = timestep_begin
+    coupled_var = strain_zz
+    activate_value = 0.05
+    active_subdomain_id = 1
+    expand_boundary_name = 'moving_interface'
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/modules/combined/test/tests/additive_manufacturing/check_initial_condition.i b/modules/combined/test/tests/additive_manufacturing/check_initial_condition.i
new file mode 100644
index 000000000000..a8d4b0f98f04
--- /dev/null
+++ b/modules/combined/test/tests/additive_manufacturing/check_initial_condition.i
@@ -0,0 +1,129 @@
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Mesh]
+  [./gen]
+    type = GeneratedMeshGenerator
+    dim = 3
+    xmin =0
+    xmax =10
+    ymin =0
+    ymax =10
+    zmin =0
+    zmax =0.5
+    nx=20
+    ny=20
+    nz=1
+  [../]
+  [./left_domain]
+    input = gen
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '0 0 0'
+    top_right = '2.5 10 0.5'
+    block_id = 1
+  [../]
+  [./middle_domain]
+    input = left_domain
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '2.5 0 0'
+    top_right = '5 10 0.5'
+    block_id = 2
+  [../]
+  [./right_domain]
+    input = middle_domain
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '5 0 0'
+    top_right = '10 10 0.5'
+    block_id = 3
+  [../]
+  [./sidesets]
+    input = right_domain
+    type = SideSetsAroundSubdomainGenerator
+    normal = '1 0 0'
+    block = 2
+    new_boundary = 'moving_interface'
+  []
+[]
+
+[Variables]
+  [./temp]
+    block = '1 2'
+  [../]
+[]
+
+[ICs]
+  [./temp_block1]
+    type =  ConstantIC
+    variable = temp
+    value = 300
+    block = 1
+  [../]
+  [./temp_block2]
+    type =  ConstantIC
+    variable = temp
+    value = 1000
+    block = 2
+  [../]
+[]
+
+[Functions]
+  [./fx]
+    type = ParsedFunction
+    value= '5.25'
+  [../]
+  [./fy]
+    type = ParsedFunction
+    value= '2.5*t'
+  [../]
+  [./fz]
+    type = ParsedFunction
+    value= '0.25'
+  [../]
+[]
+
+[Preconditioning]
+  [./smp]
+    type = SMP
+    full = true
+  [../]
+[]
+
+
+[Executioner]
+  type = Transient
+
+  automatic_scaling = true
+
+  solve_type = 'NEWTON'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  line_search = 'none'
+
+  l_max_its = 10
+  nl_max_its = 20
+  nl_rel_tol = 1e-4
+
+  start_time = 0.0
+  end_time = 1.0
+  dt = 1e-1
+  dtmin = 1e-4
+[]
+
+[UserObjects]
+  [./activated_elem_uo]
+    type = ActivateElementsByPath
+    execute_on = timestep_begin
+    function_x= fx
+    function_y= fy
+    function_z= fz
+    active_subdomain_id = 2
+    expand_boundary_name = 'moving_interface'
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/modules/combined/test/tests/additive_manufacturing/check_stateful_properties.i b/modules/combined/test/tests/additive_manufacturing/check_stateful_properties.i
new file mode 100644
index 000000000000..84273dc0d342
--- /dev/null
+++ b/modules/combined/test/tests/additive_manufacturing/check_stateful_properties.i
@@ -0,0 +1,178 @@
+[Problem]
+  kernel_coverage_check = false
+  material_coverage_check = false
+[]
+
+[Mesh]
+  [./gen]
+    type = GeneratedMeshGenerator
+    dim = 3
+    xmin =0
+    xmax =10
+    ymin =0
+    ymax =10
+    zmin =0
+    zmax =0.5
+    nx=20
+    ny=20
+    nz=1
+  [../]
+  [./left_domain]
+    input = gen
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '0 0 0'
+    top_right = '5 10 0.5'
+    block_id = 1
+  [../]
+  [./right_domain]
+    input = left_domain
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '5 0 0'
+    top_right = '10 10 0.5'
+    block_id = 2
+  [../]
+  [./sidesets]
+    input = right_domain
+    type = SideSetsAroundSubdomainGenerator
+    normal = '1 0 0'
+    block = 1
+    new_boundary = 'moving_interface'
+  []
+[]
+
+[Variables]
+  [./temp]
+    initial_condition = 300
+    block = '1'
+  [../]
+[]
+
+# Output aux variables to check if stateful properties
+# are initialized properly for newly added elements
+[AuxVariables]
+  [./density_aux]
+    order = CONSTANT
+    family = MONOMIAL
+    block = '1'
+  [../]
+  [./specific_heat_aux]
+    order = CONSTANT
+    family = MONOMIAL
+    block = '1'
+  [../]
+  [./thermal_conductivity_aux]
+    order = CONSTANT
+    family = MONOMIAL
+    block = '1'
+  [../]
+[]
+
+[Kernels]
+  [./null]
+    type = NullKernel
+    variable = temp
+    jacobian_fill = 1e-5
+  [../]
+[]
+
+[AuxKernels]
+  [./density]
+    type = ADMaterialRealAux
+    property = density
+    variable = density_aux
+    block = 1
+  []
+  [./specific_heat]
+    type = ADMaterialRealAux
+    property = specific_heat
+    variable = specific_heat_aux
+    block = 1
+  []
+  [./thermal_conductivity]
+    type = ADMaterialRealAux
+    property = thermal_conductivity
+    variable = thermal_conductivity_aux
+    block = 1
+  []
+[]
+
+[Functions]
+  [./fx]
+    type = ParsedFunction
+    value= '5.25'
+  [../]
+  [./fy]
+    type = ParsedFunction
+    value= '2.5*t'
+  [../]
+  [./fz]
+    type = ParsedFunction
+    value= '0.25'
+  [../]
+[]
+
+[Materials]
+  [./density]
+    type = ADDensity
+    density = 4.43e-6
+    block = '1'
+  [../]
+  [./heat]
+    type = ADHeatConductionMaterial
+    specific_heat = 600
+    thermal_conductivity = 10e-3
+    block = '1'
+  [../]
+  [./volumetric_heat]
+    type = ADGenericConstantMaterial
+    prop_names = 'volumetric_heat'
+    prop_values = 100
+    block = '1'
+  [../]
+[]
+
+[Preconditioning]
+  [./smp]
+    type = SMP
+    full = true
+  [../]
+[]
+
+
+[Executioner]
+  type = Transient
+
+  automatic_scaling = true
+
+  solve_type = 'NEWTON'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
+
+  line_search = 'none'
+
+  l_max_its = 10
+  nl_max_its = 20
+  nl_rel_tol = 1e-4
+
+  start_time = 0.0
+  end_time = 1.0
+  dt = 1e-1
+  dtmin = 1e-4
+[]
+
+[UserObjects]
+  [./activated_elem_uo]
+    type = ActivateElementsByPath
+    execute_on = timestep_begin
+    function_x= fx
+    function_y= fy
+    function_z= fz
+    active_subdomain_id = 1
+    expand_boundary_name = 'moving_interface'
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/modules/combined/test/tests/additive_manufacturing/gold/check_element_addition_by_variable_out.e-s002 b/modules/combined/test/tests/additive_manufacturing/gold/check_element_addition_by_variable_out.e-s002
new file mode 100644
index 000000000000..9d128d1d52ac
Binary files /dev/null and b/modules/combined/test/tests/additive_manufacturing/gold/check_element_addition_by_variable_out.e-s002 differ
diff --git a/modules/combined/test/tests/additive_manufacturing/gold/check_element_addition_out.e-s002 b/modules/combined/test/tests/additive_manufacturing/gold/check_element_addition_out.e-s002
new file mode 100644
index 000000000000..a803842ffc6e
Binary files /dev/null and b/modules/combined/test/tests/additive_manufacturing/gold/check_element_addition_out.e-s002 differ
diff --git a/modules/combined/test/tests/additive_manufacturing/gold/check_initial_condition_out.e-s002 b/modules/combined/test/tests/additive_manufacturing/gold/check_initial_condition_out.e-s002
new file mode 100644
index 000000000000..63d0b77b418e
Binary files /dev/null and b/modules/combined/test/tests/additive_manufacturing/gold/check_initial_condition_out.e-s002 differ
diff --git a/modules/combined/test/tests/additive_manufacturing/gold/check_stateful_properties_out.e-s002 b/modules/combined/test/tests/additive_manufacturing/gold/check_stateful_properties_out.e-s002
new file mode 100644
index 000000000000..559243f456ce
Binary files /dev/null and b/modules/combined/test/tests/additive_manufacturing/gold/check_stateful_properties_out.e-s002 differ
diff --git a/modules/combined/test/tests/additive_manufacturing/tests b/modules/combined/test/tests/additive_manufacturing/tests
new file mode 100644
index 000000000000..831f393fc470
--- /dev/null
+++ b/modules/combined/test/tests/additive_manufacturing/tests
@@ -0,0 +1,32 @@
+[Tests]
+  issues = '#15795'
+  design = 'ActivateElementsUserObjectBase.md ActivateElementsByPath.md ActivateElementsCoupled.md'
+  [./check_element_addition]
+    type = 'Exodiff'
+    input = 'check_element_addition.i'
+    exodiff = 'check_element_addition_out.e-s002'
+    cli_args = "Executioner/end_time=0.5 Outputs/interval=5"
+    requirement = 'The subdomain 1 should be expanded while subdomain 2 should be contracted via ActivateElementsUserObject. The addition of elements is controlled by a function path.'
+  [../]
+  [./check_element_addition_by_variable]
+    type = 'Exodiff'
+    input = 'check_element_addition_by_variable.i'
+    exodiff = 'check_element_addition_by_variable_out.e-s002'
+    cli_args = "Executioner/end_time=0.5 Outputs/interval=5"
+    requirement = 'The subdomain 1 should be expanded via ActivateElementsUserObject. The addition of elements is controlled by a coupled variable strain_zz.'
+  [../]
+  [./check_stateful_properties]
+    type = 'Exodiff'
+    input = 'check_stateful_properties.i'
+    exodiff = 'check_stateful_properties_out.e-s002'
+    cli_args = "Executioner/end_time=0.5 Outputs/interval=5"
+    requirement = 'The newly activated elements should have the stateful material properties for heat conduction as specified in the input.'
+  [../]
+  [./check_initial_condition]
+    type = 'Exodiff'
+    input = 'check_initial_condition.i'
+    exodiff = 'check_initial_condition_out.e-s002'
+    cli_args = "Executioner/end_time=0.5 Outputs/interval=5"
+    requirement = 'The newly activated elements should have the IC of variables as specified in the input.'
+  [../]
+[]