Handle stateful properties in xfem moving interface

framework/include/materials/MaterialData.h

@@ -79,6 +79,22 @@ class MaterialData
   /// copy material properties from one element to another
   void copy(const Elem & elem_to, const Elem & elem_from, unsigned int side);
 
+  /**
+   * Copy material properties to elem_to.
+   *
+   * @param elem_to Element to copy data to
+   * @param props material property data to copy
+   * @param props_old old material property data to copy
+   * @param props_older older material property data to copy, can be empty if the target
+   * element has no older properties
+   * @param side Side number (elemental material properties have this equal to zero)
+   */
+  void copy(const Elem & elem_to,
+            const HashMap<unsigned int, MaterialProperties> & props,
+            const HashMap<unsigned int, MaterialProperties> & props_old,
+            const HashMap<unsigned int, MaterialProperties> & props_older,
+            unsigned int side);
+
   /// material properties for given element (and possible side)
   void swap(const Elem & elem, unsigned int side = 0);
 

framework/include/materials/MaterialProperty.h

@@ -254,7 +254,7 @@ rawValueEqualityHelper(std::vector<std::vector<T1>> & out, const std::vector<std
   for (MooseIndex(in) i = 0; i < in.size(); ++i)
   {
     out[i].resize(in[i].size());
-    for (MooseIndex(in[i].size()) j = 0; j < in[i].size(); ++j)
+    for (MooseIndex(in[i]) j = 0; j < in[i].size(); ++j)
       out[i][j] = MetaPhysicL::raw_value(in[i][j]);
   }
 }

framework/include/materials/MaterialPropertyStorage.h

@@ -138,13 +138,34 @@ class MaterialPropertyStorage
    * @param elem_to Element to copy data to
    * @param elem_from Element to copy data from
    * @param side Side number (elemental material properties have this equal to zero)
+   * @param n_qpoints number of quadrature points to work with
    */
   void copy(MaterialData & material_data,
             const Elem & elem_to,
             const Elem & elem_from,
             unsigned int side,
             unsigned int n_qpoints);
 
+  /**
+   * Copy material properties to elem_to. Thread safe.
+   *
+   * @param material_data MaterialData object to work with
+   * @param elem_to Element to copy data to
+   * @param props_from material property data to copy
+   * @param props_from_old old material property data to copy
+   * @param props_from_older older material property data to copy, can be empty if the target
+   * element has no older properties
+   * @param side Side number (elemental material properties have this equal to zero)
+   * @param n_qpoints number of quadrature points to work with
+   */
+  void copy(MaterialData & material_data,
+            const Elem & elem_to,
+            const HashMap<unsigned int, MaterialProperties> & props_from,
+            const HashMap<unsigned int, MaterialProperties> & props_from_old,
+            const HashMap<unsigned int, MaterialProperties> & props_from_older,
+            unsigned int side,
+            unsigned int n_qpoints);
+
   /**
    * Swap (shallow copy) material properties in MaterialData and MaterialPropertyStorage
    * Thread safe

framework/src/materials/MaterialData.C

@@ -56,6 +56,16 @@ MaterialData::copy(const Elem & elem_to, const Elem & elem_from, unsigned int si
   _storage.copy(*this, elem_to, elem_from, side, _n_qpoints);
 }
 
+void
+MaterialData::copy(const Elem & elem_to,
+                   const HashMap<unsigned int, MaterialProperties> & props,
+                   const HashMap<unsigned int, MaterialProperties> & props_old,
+                   const HashMap<unsigned int, MaterialProperties> & props_older,
+                   unsigned int side)
+{
+  _storage.copy(*this, elem_to, props, props_old, props_older, side, _n_qpoints);
+}
+
 void
 MaterialData::swap(const Elem & elem, unsigned int side /* = 0*/)
 {

framework/src/materials/MaterialPropertyStorage.C

@@ -355,6 +355,31 @@ MaterialPropertyStorage::copy(MaterialData & material_data,
   }
 }
 
+void
+MaterialPropertyStorage::copy(MaterialData & material_data,
+                              const Elem & elem_to,
+                              const HashMap<unsigned int, MaterialProperties> & props_from,
+                              const HashMap<unsigned int, MaterialProperties> & props_from_old,
+                              const HashMap<unsigned int, MaterialProperties> & props_from_older,
+                              unsigned int side,
+                              unsigned int n_qpoints)
+{
+  initProps(material_data, elem_to, side, n_qpoints);
+  for (unsigned int i = 0; i < _stateful_prop_id_to_prop_id.size(); ++i)
+  {
+    for (unsigned int qp = 0; qp < n_qpoints; ++qp)
+    {
+      props(&elem_to, side)[i]->qpCopy(qp, props_from.at(side)[i], qp);
+      propsOld(&elem_to, side)[i]->qpCopy(qp, props_from_old.at(side)[i], qp);
+      if (hasOlderProperties())
+      {
+        mooseAssert(!props_from_older.empty(), "trying to access empty older properties");
+        propsOlder(&elem_to, side)[i]->qpCopy(qp, props_from_older.at(side)[i], qp);
+      }
+    }
+  }
+}
+
 void
 MaterialPropertyStorage::swap(MaterialData & material_data, const Elem & elem, unsigned int side)
 {

modules/tensor_mechanics/src/materials/ADComputeFiniteStrainElasticStress.C

@@ -26,8 +26,8 @@ ADComputeFiniteStrainElasticStress::ADComputeFiniteStrainElasticStress(
     _elasticity_tensor_name(_base_name + "elasticity_tensor"),
     _elasticity_tensor(getADMaterialProperty<RankFourTensor>(_elasticity_tensor_name)),
     _strain_increment(getADMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
-    _rotation_total(declareADProperty<RankTwoTensor>("rotation_total")),
-    _rotation_total_old(getMaterialPropertyOldByName<RankTwoTensor>("rotation_total")),
+    _rotation_total(declareADProperty<RankTwoTensor>(_base_name + "rotation_total")),
+    _rotation_total_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "rotation_total")),
     _rotation_increment(
         getADMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
     _stress_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "stress")),

modules/tensor_mechanics/src/materials/ComputeFiniteStrainElasticStress.C

@@ -25,8 +25,8 @@ ComputeFiniteStrainElasticStress::ComputeFiniteStrainElasticStress(
     GuaranteeConsumer(this),
     _elasticity_tensor_name(_base_name + "elasticity_tensor"),
     _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
-    _rotation_total(declareProperty<RankTwoTensor>("rotation_total")),
-    _rotation_total_old(getMaterialPropertyOldByName<RankTwoTensor>("rotation_total")),
+    _rotation_total(declareProperty<RankTwoTensor>(_base_name + "rotation_total")),
+    _rotation_total_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "rotation_total")),
     _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
     _rotation_increment(
         getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),

modules/xfem/include/base/XFEM.h

@@ -319,6 +319,27 @@ class XFEM : public XFEMInterface
    */
   std::map<unique_id_type, std::vector<Real>> _cached_aux_solution;
 
+  /**
+   * Data structures to store material properties of the children elements prior to heal. These
+   * material properties are copied back to the children elements if the healed element is re-cut.
+   */
+  std::map<const Elem *,
+           std::pair<HashMap<unsigned int, MaterialProperties>,
+                     HashMap<unsigned int, MaterialProperties>>>
+      _healed_material_properties;
+  std::map<const Elem *,
+           std::pair<HashMap<unsigned int, MaterialProperties>,
+                     HashMap<unsigned int, MaterialProperties>>>
+      _healed_material_properties_old;
+  std::map<const Elem *,
+           std::pair<HashMap<unsigned int, MaterialProperties>,
+                     HashMap<unsigned int, MaterialProperties>>>
+      _healed_material_properties_older;
+  std::map<const Elem *, std::pair<bool, bool>> _healed_material_properties_used;
+
+  /// healed geometric cuts
+  std::map<const Elem *, const GeometricCutUserObject *> _healed_cuts;
+
   /**
    * Store the solution in stored_solution for a given node
    * @param node_to_store_to   Node for which the solution will be stored
@@ -398,4 +419,46 @@ class XFEM : public XFEMInterface
    * @param sys  System for which the dof indices are found
    */
   std::vector<dof_id_type> getNodeSolutionDofs(const Node * node, SystemBase & sys) const;
+
+  /**
+   * Get the GeometricCutUserObject associated with an element
+   * @param elem The element
+   * @return     A constant pointer to the GeometricCutUserObject, nullptr if nothing found
+   */
+  const GeometricCutUserObject * getGeometricCutForElem(const Elem * elem) const;
+
+  /**
+   * Helper function to store the material properties of a healed element
+   * @param parent_elem The parent element
+   * @param elem1       The first child element
+   * @param elem2       The second child element
+   */
+  void storeMaterialPropertiesForElements(const Elem * parent_elem,
+                                          const Elem * elem1,
+                                          const Elem * elem2);
+
+  /**
+   * Helper function to store the material properties of a healed element
+   * @param parent_elem The parent element
+   * @param cut_elem    The element being cut
+   * @param props       The material properties to be set
+   * @param props_old   The old material properties to be set
+   * @param props_older The older material properties to be set
+   */
+  void setMaterialPropertiesForElement(
+      const Elem * parent_elem,
+      const Elem * cut_elem,
+      const HashMap<unsigned int, MaterialProperties> & props,
+      const HashMap<unsigned int, MaterialProperties> & props_old,
+      const HashMap<unsigned int, MaterialProperties> & props_older) const;
+
+  /**
+   * Determine which side of the element belongs to relative to the cut
+   * @param parent_elem The parent element
+   * @param cut_elem    The element being cut
+   * @param gcuo        The GeometricCutUserObject for the cut
+   */
+  bool getElementSideRelativeToInterface(const Elem * parent_elem,
+                                         const Elem * cut_elem,
+                                         const GeometricCutUserObject * gcuo) const;
 };

modules/xfem/include/materials/LevelSetBiMaterialBase.h

@@ -27,6 +27,7 @@ class LevelSetBiMaterialBaseTempl : public Material
   LevelSetBiMaterialBaseTempl(const InputParameters & parameters);
 
 protected:
+  virtual void initQpStatefulProperties() override;
   virtual void computeProperties() override;
   virtual void computeQpProperties() override;
 

modules/xfem/include/userobjects/GeometricCutUserObject.h

@@ -182,6 +182,20 @@ class GeometricCutUserObject : public CrackFrontPointsProvider
    */
   bool shouldHealMesh() const { return _heal_always; };
 
+  /**
+   * Get an unsigned integer telling which side the node belongs to relative to the cut.
+   * The returned ID contains no physical meaning, but should be consistent through out the
+   * simulation.
+   * @param node Pointer to the node
+   * @return an unsigned int indicating the side
+   */
+  virtual unsigned int getCutSideID(const Node * /*node*/) const
+  {
+    mooseError("Objects that inherit from GeometricCutUserObject should override the "
+               "getCutSideID method");
+    return 0;
+  }
+
 protected:
   /// Pointer to the XFEM controller object
   std::shared_ptr<XFEM> _xfem;

modules/xfem/include/userobjects/LevelSetCutUserObject.h

@@ -38,6 +38,13 @@ class LevelSetCutUserObject : public GeometricCutUserObject
   virtual const std::vector<Point>
   getCrackFrontPoints(unsigned int num_crack_front_points) const override;
 
+  /**
+   * If the levelset value is positive, return 1, otherwise return 0.
+   * @param node Pointer to the node
+   * @return an unsigned int indicating the side
+   */
+  virtual unsigned int getCutSideID(const Node * node) const override;
+
 protected:
   /// The variable number of the level set variable we using to define the cuts
   const unsigned int _level_set_var_number;