remove the property storage of XFEM deleted element and side

ref idaholab#17063

framework/include/materials/MaterialData.h

@@ -191,6 +191,13 @@ class MaterialData
    */
   bool isOnlyResizeIfSmaller() const { return _resize_only_if_smaller; };
 
+  /**
+   * Remove the property storage and element pointer from MaterialPropertyStorage data structures
+   * Use this when elements are deleted so we don't end up with invalid elem pointers (for e.g.
+   * stateful properties) hanging around in our data structures
+   */
+  void eraseProperty(const Elem * elem) { _storage.eraseProperty(elem); };
+
 protected:
   /// Reference to the MaterialStorage class
   MaterialPropertyStorage & _storage;

framework/include/materials/MaterialPropertyStorage.h

@@ -255,13 +255,17 @@ class MaterialPropertyStorage
     return _prop_names.count(retrievePropertyId(prop_name)) > 0;
   }
 
+  /**
+   * Remove the property storage and element pointer from internal data structures
+   * Use this when elements are deleted so we don't end up with invalid elem pointers (for e.g.
+   * stateful properties) hanging around in our data structures
+   */
+  void eraseProperty(const Elem * elem);
+
 protected:
   /// Release all internal data structures
   void releaseProperties();
 
-  /// Remove the property storage and element pointer from internal data structures
-  void eraseProperty(const Elem * elem);
-
   /// Internal property storage release helper
   void releasePropertyMap(HashMap<unsigned int, MaterialProperties> & inner_map);
 

modules/xfem/src/base/XFEM.C

@@ -1430,6 +1430,11 @@ XFEM::cutMeshWithEFA(NonlinearSystemBase & nl, AuxiliarySystem & aux)
     }
 
     elem_to_delete->nullify_neighbors();
+
+    // remove the property storage of deleted element/side
+    (*_material_data)[0]->eraseProperty(elem_to_delete);
+    (*_bnd_material_data)[0]->eraseProperty(elem_to_delete);
+
     _mesh->boundary_info->remove(elem_to_delete);
     unsigned int deleted_elem_id = elem_to_delete->id();
     _mesh->delete_elem(elem_to_delete);

modules/xfem/test/tests/checkpoint/checkpoint.i

@@ -0,0 +1,314 @@
+# This test is for two layer materials with different youngs modulus
+# The global stress is determined by switching the stress based on level set values
+# The material interface is marked by a level set function
+# The two layer materials are glued together
+
+[GlobalParams]
+  order = FIRST
+  family = LAGRANGE
+  displacements = 'disp_x disp_y'
+[]
+
+[XFEM]
+  qrule = volfrac
+  output_cut_plane = true
+[]
+
+[UserObjects]
+  [./level_set_cut_uo]
+    type = LevelSetCutUserObject
+    level_set_var = ls
+  [../]
+[]
+
+[Mesh]
+  displacements = 'disp_x disp_y'
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 5
+    ny = 5
+    xmin = 0.0
+    xmax = 5.
+    ymin = 0.0
+    ymax = 5.
+    elem_type = QUAD4
+  []
+  [./left_bottom]
+    type = ExtraNodesetGenerator
+    new_boundary = 'left_bottom'
+    coord = '0.0 0.0'
+    input = gen
+  [../]
+  [./left_top]
+    type = ExtraNodesetGenerator
+    new_boundary = 'left_top'
+    coord = '0.0 5.'
+    input = left_bottom
+  [../]
+[]
+
+[AuxVariables]
+  [./ls]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[AuxKernels]
+  [./ls_function]
+    type = FunctionAux
+    variable = ls
+    function = ls_func
+  [../]
+[]
+
+[Variables]
+  [./disp_x]
+  [../]
+  [./disp_y]
+  [../]
+[]
+
+[Functions]
+  [./ls_func]
+    type = ParsedFunction
+    value = 'y-2.5'
+  [../]
+[]
+
+[AuxVariables]
+  [./stress_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./a_strain_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./a_strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./a_strain_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./b_strain_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./b_strain_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./b_strain_xy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+[]
+
+[Kernels]
+  [./TensorMechanics]
+  [../]
+[]
+
+[AuxKernels]
+  [./stress_xx]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    index_i = 0
+    index_j = 0
+    variable = stress_xx
+  [../]
+  [./stress_yy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    index_i = 1
+    index_j = 1
+    variable = stress_yy
+  [../]
+  [./stress_xy]
+    type = RankTwoAux
+    rank_two_tensor = stress
+    index_i = 0
+    index_j = 1
+    variable = stress_xy
+  [../]
+  [./a_strain_xx]
+    type = RankTwoAux
+    rank_two_tensor = A_total_strain
+    index_i = 0
+    index_j = 0
+    variable = a_strain_xx
+  [../]
+  [./a_strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = A_total_strain
+    index_i = 1
+    index_j = 1
+    variable = a_strain_yy
+  [../]
+  [./a_strain_xy]
+    type = RankTwoAux
+    rank_two_tensor = A_total_strain
+    index_i = 0
+    index_j = 1
+    variable = a_strain_xy
+  [../]
+  [./b_strain_xx]
+    type = RankTwoAux
+    rank_two_tensor = B_total_strain
+    index_i = 0
+    index_j = 0
+    variable = b_strain_xx
+  [../]
+  [./b_strain_yy]
+    type = RankTwoAux
+    rank_two_tensor = B_total_strain
+    index_i = 1
+    index_j = 1
+    variable = b_strain_yy
+  [../]
+  [./b_strain_xy]
+    type = RankTwoAux
+    rank_two_tensor = B_total_strain
+    index_i = 0
+    index_j = 1
+    variable = b_strain_xy
+  [../]
+[]
+
+[Constraints]
+  [./dispx_constraint]
+    type = XFEMSingleVariableConstraint
+    use_displaced_mesh = false
+    variable = disp_x
+    alpha = 1e8
+    geometric_cut_userobject = 'level_set_cut_uo'
+  [../]
+  [./dispy_constraint]
+    type = XFEMSingleVariableConstraint
+    use_displaced_mesh = false
+    variable = disp_y
+    alpha = 1e8
+    geometric_cut_userobject = 'level_set_cut_uo'
+  [../]
+[]
+
+[BCs]
+  [./bottomx]
+    type = DirichletBC
+    boundary = bottom
+    variable = disp_x
+    value = 0.0
+  [../]
+  [./bottomy]
+    type = DirichletBC
+    boundary = bottom
+    variable = disp_y
+    value = 0.0
+  [../]
+  [./topx]
+    type = FunctionDirichletBC
+    boundary = top
+    variable = disp_x
+    function = 0.03*t
+  [../]
+  [./topy]
+    type = FunctionDirichletBC
+    boundary = top
+    variable = disp_y
+    function = '0.03*t'
+  [../]
+[]
+
+[Materials]
+  [./elasticity_tensor_A]
+    type = ComputeIsotropicElasticityTensor
+    base_name = A
+    youngs_modulus = 1e9
+    poissons_ratio = 0.3
+  [../]
+  [./strain_A]
+    type = ComputeSmallStrain
+    base_name = A
+  [../]
+  [./stress_A]
+    type = ComputeLinearElasticStress
+    base_name = A
+  [../]
+  [./elasticity_tensor_B]
+    type = ComputeIsotropicElasticityTensor
+    base_name = B
+    youngs_modulus = 1e5
+    poissons_ratio = 0.3
+  [../]
+  [./strain_B]
+    type = ComputeSmallStrain
+    base_name = B
+  [../]
+  [./stress_B]
+    type = ComputeLinearElasticStress
+    base_name = B
+  [../]
+  [./combined_stress]
+    type = LevelSetBiMaterialRankTwo
+    levelset_positive_base = 'A'
+    levelset_negative_base = 'B'
+    level_set_var = ls
+    prop_name = stress
+  [../]
+  [./combined_dstressdstrain]
+    type = LevelSetBiMaterialRankFour
+    levelset_positive_base = 'A'
+    levelset_negative_base = 'B'
+    level_set_var = ls
+    prop_name = Jacobian_mult
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+
+  solve_type = 'PJFNK'
+  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
+  petsc_options_value = '201                hypre    boomeramg      8'
+
+  line_search = 'bt'
+
+# controls for linear iterations
+  l_max_its = 20
+  l_tol = 1e-3
+
+# controls for nonlinear iterations
+  nl_max_its = 15
+  nl_rel_tol = 1e-14
+  nl_abs_tol = 1e-7
+
+# time control
+  start_time = 0.0
+  dt = 0.1
+  num_steps = 2
+
+  max_xfem_update = 1
+[]
+
+[Outputs]
+  checkpoint = true
+  exodus = true
+  execute_on = timestep_end
+  csv = true
+  [./console]
+    type = Console
+    output_linear = true
+  [../]
+[]

modules/xfem/test/tests/checkpoint/tests

@@ -0,0 +1,19 @@
+[Tests]
+  [checkpoint]
+    type = 'CheckFiles'
+    input = checkpoint.i
+    max_parallel = 1
+    max_threads = 1
+    check_files ='checkpoint_out_cp/0001.xdr
+                  checkpoint_out_cp/0001.xdr.0000
+                  checkpoint_out_cp/0001.rd-0
+                  checkpoint_out_cp/0001_mesh.cpr/1/header.cpr
+                  checkpoint_out_cp/0002.xdr
+                  checkpoint_out_cp/0002.xdr.0000
+                  checkpoint_out_cp/0002.rd-0
+                  checkpoint_out_cp/0002_mesh.cpr/1/header.cpr'
+   requirement = 'The XFEM module shall support outputting of checkpoint files.'
+   design = 'outputs/Checkpoint.md'
+   issues = '#17063'
+  []
+[]