Merge pull request #32419 from lynnmunday/xfem_autoMeshNodes

automatically determine 3D mesh cutter crack front nodes

Author: bwspenc

modules/solid_mechanics/doc/content/source/actions/DomainIntegralAction.md

@@ -31,4 +31,8 @@ Details on the theory behind the computation of the various fracture integrals a
 | [DomainIntegralQFunction](DomainIntegralQFunction.md) | AuxKernel | Optionally populates AuxVariables with values of the geometry-based q function when `output_q=true` |
 | [DomainIntegralTopologicalQFunction](DomainIntegralTopologicalQFunction.md) | AuxKernel | Optionally populates AuxVariables with values of the topology-based q function when `output_q=true` |
 
+### Output with `crack_front_points_provider`
+
+The [!param](/DomainIntegral/DomainIntegralAction/crack_front_points_provider) parameter references a UserObject that defines a set of surface points at the edge of a cutting plane used by the XFEM module to define the crack. One of its primary purposes is to define the discrete locations where fracture integrals are computed when XFEM is used to define the crack. In XFEM, a crack can be described purely by geometric primitives or by a lower-dimensional mesh.  When [!param](/DomainIntegral/DomainIntegralAction/crack_front_points_provider) is specified, the `DomainIntegral` action can create output at the crack front points, which are not at node locations.  By default (`output_vpp=true`), the vectorpostprocessors will always be created and will be dynamically sized as the number of crack front points changes throughout the simulation. The Postprocessor and AuxKernel output will only be created if [!param](/DomainIntegral/DomainIntegralAction/number_points_from_provider) is specified.
+
 !syntax parameters /DomainIntegral/DomainIntegralAction

modules/solid_mechanics/include/userobjects/CrackFrontDefinition.h

@@ -23,9 +23,6 @@ namespace libMesh
 class QBase;
 }
 
-// fixme lynn this needs to be removed from the globalname space in a seperate commit after Grizzly
-// is fixed to use includeCrackFrontDefinitionParams
-void addCrackFrontDefinitionParams(InputParameters & params);
 /**
  * Class used in fracture integrals to define geometric characteristics of the crack front
  */

modules/solid_mechanics/include/userobjects/CrackFrontPointsProvider.h

@@ -21,6 +21,11 @@ class CrackFrontPointsProvider : public ElementUserObject
 
   CrackFrontPointsProvider(const InputParameters & parameters, const bool uses_mesh = false);
 
+  /** Get the current number of crack front points
+   * @return The current number of crack front points
+   */
+  virtual unsigned int getNumberOfCrackFrontPoints() const = 0;
+
   /** get a set of points along a crack front from a XFEM GeometricCutUserObject
    * @return A vector which contains all crack front points
    */

modules/solid_mechanics/src/actions/DomainIntegralAction.C

@@ -113,7 +113,7 @@ DomainIntegralAction::validParams()
       "equivalent_k",
       false,
       "Calculate an equivalent K from KI, KII and KIII, assuming self-similar crack growth.");
-  params.addParam<bool>("output_q", true, "Output q");
+  params.addParam<bool>("output_q", false, "Output q");
   params.addRequiredParam<bool>(
       "incremental", "Flag to indicate whether an incremental or total model is being used.");
   params.addParam<std::vector<MaterialPropertyName>>(
@@ -229,10 +229,6 @@ DomainIntegralAction::DomainIntegralAction(const InputParameters & params)
 
   if (isParamValid("crack_front_points_provider"))
   {
-    if (!isParamValid("number_points_from_provider"))
-      paramError("number_points_from_provider",
-                 "DomainIntegral error: when crack_front_points_provider is used, "
-                 "number_points_from_provider must be provided.");
     _use_crack_front_points_provider = true;
     _crack_front_points_provider = getParam<UserObjectName>("crack_front_points_provider");
   }
@@ -370,12 +366,10 @@ DomainIntegralAction::act()
     const std::string uo_type_name("CrackFrontDefinition");
 
     InputParameters params = _factory.getValidParams(uo_type_name);
+    // The CrackFrontDefinition updates the vpps and MUST execute before them
+    params.set<int>("execution_order_group") = -1;
     if (_use_crack_front_points_provider)
-    {
-      // The CrackFrontDefinition updates the vpps and MUST execute before them
-      params.set<int>("execution_order_group") = -1;
       params.set<ExecFlagEnum>("execute_on") = xfem_exec_flags;
-    }
     else
       params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
 
@@ -400,8 +394,12 @@ DomainIntegralAction::act()
     if (_crack_front_points.size() != 0)
       params.set<std::vector<Point>>("crack_front_points") = _crack_front_points;
     if (_use_crack_front_points_provider)
-      params.applyParameters(parameters(),
-                             {"crack_front_points_provider, number_points_from_provider"});
+    {
+      params.set<UserObjectName>("crack_front_points_provider") = _crack_front_points_provider;
+      if (isParamValid("number_points_from_provider"))
+        params.set<unsigned int>("number_points_from_provider") =
+            getParam<unsigned int>("number_points_from_provider");
+    }
     if (_closed_loop)
       params.set<bool>("closed_loop") = _closed_loop;
     params.set<bool>("use_displaced_mesh") = _use_displaced_mesh;
@@ -436,6 +434,14 @@ DomainIntegralAction::act()
   }
   else if (_current_task == "add_aux_variable" && _output_q)
   {
+    if (isParamValid("number_points_from_provider") && num_crack_front_points == 0)
+    {
+      paramError("output_q",
+                 "Requesting AuxVariable output of q functions but the number of crack fronts for "
+                 "output is zero.  AuxVariable output for XFEM cutter objects requires "
+                 "number_points_from_provider to be set.");
+    }
+
     for (unsigned int ring_index = 0; ring_index < _ring_vec.size(); ++ring_index)
     {
       std::string aux_var_type;
@@ -530,6 +536,26 @@ DomainIntegralAction::act()
 
   else if (_current_task == "add_postprocessor")
   {
+    // Check that the number specified by the XFEM cutter object matches the number of points
+    // specified in the DomainIntegralAction block.  This is being done in teh add_postprocessor
+    // block because it must be done after all userObjects have been created.
+    if (_use_crack_front_points_provider && isParamValid("number_points_from_provider"))
+    {
+      auto crack_front_points_provider = &_problem->getUserObject<CrackFrontPointsProvider>(
+          getParam<UserObjectName>("crack_front_points_provider"));
+      if (crack_front_points_provider->usesMesh())
+      {
+        auto xfem_cutter_points = crack_front_points_provider->getNumberOfCrackFrontPoints();
+        if (xfem_cutter_points != num_crack_front_points)
+          paramError("number_points_from_provider",
+                     "This must match the number of points provided by the XFEM mesh cutter "
+                     "object."
+                     "\n   number_points_from_provider:",
+                     num_crack_front_points,
+                     "\n   XFEM Crack Front Points: ",
+                     xfem_cutter_points);
+      }
+    }
     for (std::set<INTEGRAL>::iterator sit = _integrals.begin(); sit != _integrals.end(); ++sit)
     {
       std::string pp_base_name;
@@ -1007,7 +1033,17 @@ DomainIntegralAction::calcNumCrackFrontPoints()
   else if (_crack_front_points.size() != 0)
     num_points = _crack_front_points.size();
   else if (_use_crack_front_points_provider)
-    num_points = getParam<unsigned int>("number_points_from_provider");
+  {
+    if (isParamValid("number_points_from_provider"))
+      num_points = getParam<unsigned int>("number_points_from_provider");
+    else
+      // Actual count determined at runtime by CrackFrontDefinition::initialSetup()
+      // which calls provider->getNumberOfCrackFrontPoints(). Use 0 here so the
+      // action does not create per-point objects that would access crack front data
+      // before it exists. The VectorPostprocessors (JIntegral, InteractionIntegral,
+      // etc.) dynamically size based on the runtime count.
+      num_points = 0;
+  }
   else
     mooseError("Must define either 'boundary' or 'crack_front_points'");
   return num_points;

modules/solid_mechanics/src/userobjects/CrackFrontDefinition.C

@@ -38,12 +38,6 @@ CrackFrontDefinition::validParams()
   return params;
 }
 
-void
-addCrackFrontDefinitionParams(InputParameters & params)
-{
-  CrackFrontDefinition::includeCrackFrontDefinitionParams(params);
-}
-
 void
 CrackFrontDefinition::includeCrackFrontDefinitionParams(InputParameters & params)
 {
@@ -108,9 +102,9 @@ CrackFrontDefinition::includeCrackFrontDefinitionParams(InputParameters & params
       "crack_front_points_provider",
       "The UserObject provides the crack front points from XFEM GeometricCutObject");
 
-  params.addParam<unsigned int>(
-      "number_points_from_provider",
-      "The number of crack front points, only needed if crack_front_points_provider is used.");
+  params.addParam<unsigned int>("number_points_from_provider",
+                                "The number of crack front points, only needed for "
+                                "crack_front_points_provider that do not use a cut mesh.");
   params.addParam<Real>(
       "crack_front_node_tolerance",
       1e-10,
@@ -164,12 +158,7 @@ CrackFrontDefinition::CrackFrontDefinition(const InputParameters & parameters)
       paramError("crack_front_points_provider",
                  "CrackFrontDefinition error: since boundary is defined, "
                  "crack_front_points_provider should not be added.");
-    if (!isParamValid("number_points_from_provider"))
-      paramError("number_points_from_provider",
-                 "CrackFrontDefinition error: When crack_front_points_provider is used, the "
-                 "number_points_from_provider must be provided.");
 
-    _num_points_from_provider = getParam<unsigned int>("number_points_from_provider");
     _geom_definition_method = CRACK_GEOM_DEFINITION::CRACK_FRONT_POINTS;
   }
   else if (isParamValid("number_points_from_provider"))
@@ -302,6 +291,13 @@ CrackFrontDefinition::initialSetup()
     if (_crack_front_points_provider->usesMesh())
     {
       _use_mesh_cutter = true;
+
+      // Automatically get number of crack front points from mesh-based provider
+      _num_points_from_provider = _crack_front_points_provider->getNumberOfCrackFrontPoints();
+      mooseInfo("CrackFrontDefinition: Automatically detected ",
+                _num_points_from_provider,
+                " crack front points from mesh-based provider");
+
       if (_direction_method != DIRECTION_METHOD::CURVED_CRACK_FRONT)
         paramError("crack_direction_method",
                    "Using a `crack_front_points_provider` that uses an XFEM cutter mesh also "
@@ -311,6 +307,17 @@ CrackFrontDefinition::initialSetup()
                    "'crack_mouth_boundary' cannot be set when using a "
                    "'crack_front_points_provider' that uses an XFEM cutter mesh");
     }
+    else if (isParamValid("number_points_from_provider"))
+    {
+      _num_points_from_provider = getParam<unsigned int>("number_points_from_provider");
+    }
+    else if (_num_points_from_provider == 0)
+    {
+      paramError(
+          "number_points_from_provider",
+          "CrackFrontDefinition error: When using a non-mesh-based crack_front_points_provider, "
+          "the number_points_from_provider parameter must be provided.");
+    }
   }
   if (_crack_front_points_provider != nullptr)
   {
@@ -365,6 +372,9 @@ CrackFrontDefinition::initialize()
   // cutter
   if (_use_mesh_cutter && _is_cutter_modified)
   {
+    // Automatically update the number of crack front points as crack grows
+    _num_points_from_provider = _crack_front_points_provider->getNumberOfCrackFrontPoints();
+
     _crack_front_points =
         _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider);
     updateCrackFrontGeometry();

modules/solid_mechanics/test/tests/domain_integral_thermal/c_integral_2d.i

@@ -101,6 +101,7 @@
   crack_direction_vector = '1 0 0'
   2d = true
   axis_2d = 2
+  output_q = true
   radius_inner = '60.0 80.0 100.0 120.0'
   radius_outer = '80.0 100.0 120.0 140.0'
   incremental = true

modules/solid_mechanics/test/tests/domain_integral_thermal/interaction_integral_2d_c.i

@@ -136,6 +136,7 @@
   radius_outer = '80.0 100.0 120.0 140.0'
   symmetry_plane = 1
   incremental = true
+  output_q = true
 
   # interaction integral parameters
   block = 1

modules/solid_mechanics/test/tests/domain_integral_thermal/j_integral_2d.i

@@ -34,6 +34,7 @@
   crack_direction_vector = '1 0 0'
   2d = true
   axis_2d = 2
+  output_q = true
   radius_inner = '60.0 80.0 100.0 120.0'
   radius_outer = '80.0 100.0 120.0 140.0'
   temperature = temp

modules/solid_mechanics/test/tests/interaction_integral/interaction_integral_2d.i

@@ -35,7 +35,6 @@
   block = 1
   youngs_modulus = 207000
   poissons_ratio = 0.3
-  output_q = false
   incremental = true
   equivalent_k = true
 []

modules/solid_mechanics/test/tests/interaction_integral/interaction_integral_2d_rot.i

@@ -35,7 +35,6 @@
   block = 1
   youngs_modulus = 207000
   poissons_ratio = 0.3
-  output_q = false
   incremental = true
   equivalent_k = true
 []