Add ActivateElementsUserObjectBase, ActivateElementsCoupled, and Acti…

…vateElementsByPath classes

Add methods in FEProblemBase in support of element activation

Add tests

Refs. issue #15795

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

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

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.

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.

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;
+};

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;
+};

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;
+};

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()
 {

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);
+}