Merge pull request #27055 from lynnmunday/xfem_surfaceCut

Xfem crack nucleation in bulk

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

@@ -8,7 +8,7 @@ The `DomainIntegral` action is used to set up all of the objects used in computi
 
 Meshed crack: The crack can be explicitly included by creating a mesh with a topology that conforms to the crack. The location of the crack tip is provided to the code by defining a nodeset that includes all nodes in the finite element mesh that are located along the crack front. For 2D analyses, this nodeset would only contain a single node at the crack tip.  For 3D analyses, the mesh connectivity is used to construct a set of line segments that connect these nodes, and this is used to order the crack nodes.
 
-XFEM: Rather than defining the topology of the crack through the mesh, XFEM can be used to cut the mesh. In this case, a set of points, which does not need to conform to points in the mesh, must be provided by the user, and is used to define the location of the crack for computing the fracture integrals. Fracture integrals are computed at the locations of these points, in the order provided by the user.
+XFEM: Rather than defining the topology of the crack through the mesh, XFEM can be used to cut the mesh. In this case, a set of points, which does not need to conform to points in the mesh, must be provided by the user, and is used to define the location of the crack for computing the fracture integrals. Fracture integrals are computed at the locations of these points, in the order provided by the user.  If the `DomainIntegral` is setting up an [InteractionIntegral.md] being used as a growth criterion for the XFEM crack, i.e. [!param](/DomainIntegral/DomainIntegralAction/used_by_xfem_to_grow_crack) = true, the [CrackFrontDefinition.md] must be executed before the `DomainIntegral` vectorpostprocessors.  This is done by default in the `DomainIntegral` action which sets the the `CrackFrontDefinition` [!param](/UserObjects/CrackFrontDefinition/execution_order_group) to -1.
 
 ## Theory
 

modules/solid_mechanics/src/actions/DomainIntegralAction.C

@@ -374,7 +374,11 @@ DomainIntegralAction::act()
 
     InputParameters params = _factory.getValidParams(uo_type_name);
     if (_use_crack_front_points_provider && _used_by_xfem_to_grow_crack)
+    {
       params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END, EXEC_NONLINEAR};
+      // The CrackFrontDefinition updates the vpps and MUST execute before them
+      params.set<int>("execution_order_group") = -1;
+    }
     else
       params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
 
@@ -735,7 +739,12 @@ DomainIntegralAction::act()
         params.set<std::vector<OutputName>>("outputs") = {"none"};
 
       if (_use_crack_front_points_provider && _used_by_xfem_to_grow_crack)
+      {
+        // The CrackFrontDefinition updates this vpp and MUST execute before it
+        // This is enforced by setting the execution_order_group = -1 for the CrackFrontDefinition
+        // The CrackFrontDefinition must execute on nonlinear to update with xfem updates
         params.set<ExecFlagEnum>("execute_on") = {EXEC_TIMESTEP_END, EXEC_NONLINEAR};
+      }
       else
         params.set<ExecFlagEnum>("execute_on") = {EXEC_TIMESTEP_END};
 

modules/xfem/doc/content/source/userobjects/MeshCut2DRankTwoTensorNucleation.md

@@ -11,14 +11,14 @@ such as a principal stress or a component of stress.
 
 If the scalar exceeds a threshold specified by [!param](/UserObjects/MeshCut2DRankTwoTensorNucleation/nucleation_threshold), a line segment with a length specified by [!param](/UserObjects/MeshCut2DRankTwoTensorNucleation/nucleation_length) will be inserted
 into the [MeshCut2DFractureUserObject.md] cutter mesh.
-Cracks are only allowed to initiate from elements on boundaries specified by [!param](/UserObjects/MeshCut2DRankTwoTensorNucleation/initiate_on_boundary).
+Cracks will nucleate on the domain interior or boundary.  Crack nucleation can be limited to elements on boundaries specified by [!param](/UserObjects/MeshCut2DRankTwoTensorNucleation/initiate_on_boundary).
 Once the nucleation criterion is reached, a line segment of the specified length is inserted into the cutter mesh,
 centered on the element centroid it nucleates from.
 The direction of the nucleated crack is normal to the direction returned by the `RankTwoTensor` scalar.
 For example, `MaxInPlanePrincipal` returns the direction of the maximum in-plane principal component and the crack direction is normal to this.
-The nucleation length should be at least the length of the element it nucleates in so that the nucleated crack will completely cut the element.
+It is an error for [!param](/UserObjects/MeshCut2DRankTwoTensorNucleation/nucleation_length) to be smaller than the length of the element the crack nucleates in, unless the user sets `nucleate_across_full_element=true`.
 A crack will only be nucleated if it is at least a distance specified by [!param](/UserObjects/MeshCut2DRankTwoTensorNucleation/nucleation_radius) away from existing or nucleated cracks.
-If multiple cracks nucleate in the same xfem update and are within the specified `nucleation_radius`, the crack nucleated from the element with the lowest id will be retained and no other cracks within the `nucleation_radius` will be nucleated.
+If multiple cracks nucleate in the same xfem update and are within the specified `nucleation_radius`, the crack nucleated from the element with the lowest id will be retained and no other cracks within the `nucleation_radius` will be nucleated.  The `nucleation_radius` is a good way to limit the number of cracks that nucleate within a region, which can help limit problems related to excessive numbers of XFEM cuts in a single element.
 
 `MeshCut2DRankTwoTensorNucleation` copies several features available in the [XFEMRankTwoTensorMarkerUserObject.md].
 These include the nucleation threshold being provided as a coupled variable and the computation of the maximum value of the scalar

modules/xfem/include/userobjects/MeshCut2DRankTwoTensorNucleation.h

@@ -21,6 +21,18 @@ class MeshCut2DRankTwoTensorNucleation : public MeshCut2DNucleationBase
   virtual ~MeshCut2DRankTwoTensorNucleation() {}
 
 protected:
+  /// scaling factor to extend the nucleated crack beyond the element edges.
+  const Real _edge_extension_factor;
+
+  /// should element be cut if the nucleation_length is smaller than the element length.
+  const bool _nucleate_across_full_element;
+
+  /// is the nucleation length provided in the input file.
+  const bool _is_nucleation_length_provided;
+
+  /// Length of the nucleated cracks
+  const Real _nucleation_length;
+
   /// The tensor from which the scalar quantity used as a nucleating criterion is extracted
   const MaterialProperty<RankTwoTensor> & _tensor;
 
@@ -34,13 +46,20 @@ class MeshCut2DRankTwoTensorNucleation : public MeshCut2DNucleationBase
   const Point _point1;
   const Point _point2;
 
-  /// Length of crack to be nucleated
-  const Real _nucleation_length;
-
   /// Transformed Jacobian weights
   const MooseArray<Real> & _JxW;
   const MooseArray<Real> & _coord;
 
   virtual bool
   doesElementCrack(std::pair<RealVectorValue, RealVectorValue> & cutterElemNodes) override;
+
+  // FIXME Lynn Copy from TraceRayTools.C in rayTracing module.  Remove once this function migrates
+  // to libmesh.
+  bool lineLineIntersect2D(const Point & start,
+                           const Point & direction,
+                           const Real length,
+                           const Point & v0,
+                           const Point & v1,
+                           Point & intersection_point,
+                           Real & intersection_distance);
 };

modules/xfem/src/userobjects/MeshCut2DNucleationBase.C

@@ -20,7 +20,7 @@ InputParameters
 MeshCut2DNucleationBase::validParams()
 {
   InputParameters params = ElementUserObject::validParams();
-  params.addRequiredParam<std::vector<BoundaryName>>(
+  params.addParam<std::vector<BoundaryName>>(
       "initiate_on_boundary",
       "Cracks can only initiate on elements adjacent to specified boundaries.");
   params.addParam<Real>("nucleation_radius",
@@ -56,7 +56,7 @@ MeshCut2DNucleationBase::MeshCut2DNucleationBase(const InputParameters & paramet
   {
     std::vector<BoundaryName> initiation_boundary_names =
         getParam<std::vector<BoundaryName>>("initiate_on_boundary");
-    _initiation_boundary_ids = _mesh.getBoundaryIDs(initiation_boundary_names, true);
+    _initiation_boundary_ids = _mesh.getBoundaryIDs(initiation_boundary_names);
   }
 }
 
@@ -74,21 +74,36 @@ MeshCut2DNucleationBase::execute()
   if (_current_elem->processor_id() != processor_id())
     return;
 
-  bool isOnBoundary = false;
   unsigned int current_eid = _current_elem->id();
   std::map<unsigned int, std::pair<RealVectorValue, RealVectorValue>>::iterator mit;
   mit = _nucleated_elems.find(current_eid);
 
-  for (unsigned int i = 0; i < _initiation_boundary_ids.size(); ++i)
-    if (_mesh.isBoundaryElem(current_eid, _initiation_boundary_ids[i]))
-      isOnBoundary = true;
-  // This does not currently allow for nucleation in an element that is already cut
-  if (!is_cut && isOnBoundary && doesElementCrack(cutterElemNodes))
+  bool cut_element = false;
+  if (_initiation_boundary_ids.empty())
+  {
+    // nucleating crack in bulk
+    //  This does not currently allow for nucleation in an element that is already cut
+    cut_element = (!is_cut && doesElementCrack(cutterElemNodes));
+  }
+  else
+  {
+    // nucleating crack on specified boundary
+    bool isOnBoundary = false;
+    for (unsigned int i = 0; i < _initiation_boundary_ids.size(); ++i)
+      if (_mesh.isBoundaryElem(current_eid, _initiation_boundary_ids[i]))
+      {
+        isOnBoundary = true;
+        break;
+      }
+    // This does not currently allow for nucleation in an element that is already cut
+    cut_element = (!is_cut && isOnBoundary && doesElementCrack(cutterElemNodes));
+  }
+
+  if (cut_element)
   {
     if (mit != _nucleated_elems.end())
-    {
       mooseError("ERROR: element ", current_eid, " already marked for crack nucleation.");
-    }
+
     _nucleated_elems[current_eid] = cutterElemNodes;
   }
 }
@@ -97,7 +112,6 @@ void
 MeshCut2DNucleationBase::threadJoin(const UserObject & y)
 {
   const auto & xmuo = static_cast<const MeshCut2DNucleationBase &>(y);
-
   for (std::map<unsigned int, std::pair<RealVectorValue, RealVectorValue>>::const_iterator mit =
            xmuo._nucleated_elems.begin();
        mit != xmuo._nucleated_elems.end();

modules/xfem/src/userobjects/MeshCut2DRankTwoTensorNucleation.C

@@ -9,11 +9,12 @@
 
 #include "MeshCut2DRankTwoTensorNucleation.h"
 
+#include "MooseError.h"
+#include "MooseTypes.h"
 #include "libmesh/quadrature.h"
 #include "RankTwoTensor.h"
 #include "RankTwoScalarTools.h"
 #include "Assembly.h"
-#include <limits>
 
 registerMooseObject("XFEMApp", MeshCut2DRankTwoTensorNucleation);
 
@@ -23,6 +24,18 @@ MeshCut2DRankTwoTensorNucleation::validParams()
   InputParameters params = MeshCut2DNucleationBase::validParams();
   params.addClassDescription(
       "Nucleate a crack in MeshCut2D UO based on a scalar extracted from a RankTwoTensor");
+  params.addRangeCheckedParam<Real>(
+      "edge_extension_factor",
+      1e-5,
+      "edge_extension_factor >= 0",
+      "Factor by which the length of cracks that have been adjusted to coincide with element edges "
+      "are increased to avoid missing an intersection with the edge due to numerical tolerance "
+      "issues.");
+  params.addRangeCheckedParam<Real>(
+      "nucleation_length",
+      "nucleation_length >= 0",
+      "Size of crack to nucleate.  Must be larger than the length of the element in which the "
+      "crack is nucleated, unless 'nucleate_across_full_element' is set to 'true'.");
   params.addParam<MooseEnum>(
       "scalar_type",
       RankTwoScalarTools::scalarOptions(),
@@ -31,7 +44,6 @@ MeshCut2DRankTwoTensorNucleation::validParams()
   params.addRequiredCoupledVar(
       "nucleation_threshold",
       "Threshold for the scalar quantity of the RankTwoTensor to nucleate new cracks");
-  params.addRequiredParam<Real>("nucleation_length", "Nucleated crack length.");
   params.addParam<Point>(
       "point1",
       Point(0, 0, 0),
@@ -40,18 +52,29 @@ MeshCut2DRankTwoTensorNucleation::validParams()
       "point2",
       Point(0, 1, 0),
       "End point for axis used to calculate some cylindrical material tensor quantities");
+  params.addParam<bool>(
+      "nucleate_across_full_element",
+      false,
+      "Controls the behavior of nucleating cracks if 'nucleation_length' is smaller than the "
+      "length required to travserse an element with a nucleating crack. If this is set to 'false', "
+      "that condition will result in an error, but if set to 'true', a crack with the length "
+      "needed "
+      "to traverse the element will be inserted.");
   return params;
 }
 
 MeshCut2DRankTwoTensorNucleation::MeshCut2DRankTwoTensorNucleation(
     const InputParameters & parameters)
   : MeshCut2DNucleationBase(parameters),
+    _edge_extension_factor(getParam<Real>("edge_extension_factor")),
+    _nucleate_across_full_element(getParam<bool>("nucleate_across_full_element")),
+    _is_nucleation_length_provided(isParamValid("nucleation_length")),
+    _nucleation_length(_is_nucleation_length_provided ? getParam<Real>("nucleation_length") : 0),
     _tensor(getMaterialProperty<RankTwoTensor>(getParam<std::string>("tensor"))),
     _nucleation_threshold(coupledValue("nucleation_threshold")),
     _scalar_type(getParam<MooseEnum>("scalar_type")),
     _point1(parameters.get<Point>("point1")),
     _point2(parameters.get<Point>("point2")),
-    _nucleation_length(getParam<Real>("nucleation_length")),
     _JxW(_assembly.JxW()),
     _coord(_assembly.coordTransformation())
 {
@@ -84,13 +107,125 @@ MeshCut2DRankTwoTensorNucleation::doesElementCrack(
 
   if (tensor_quantity > average_threshold)
   {
-    const Point elem_center(_current_elem->vertex_average());
     does_it_crack = true;
+
+    const Point elem_center(_current_elem->vertex_average());
     Point crack_dir = point_dir.cross(Point(0, 0, 1));
-    RealVectorValue point_0 = elem_center - _nucleation_length / 2.0 * crack_dir.unit();
-    RealVectorValue point_1 = elem_center + _nucleation_length / 2.0 * crack_dir.unit();
+
+    // get max element length to use for ray length
+    std::unique_ptr<const Elem> edge;
+    Real circumference = 0;
+    for (const auto e : _current_elem->edge_index_range())
+    {
+      _current_elem->build_edge_ptr(edge, e);
+      circumference += edge->length(0, 1);
+    }
+
+    // Temproraries for doing things below
+    Real intersection_distance;
+
+    Point point_0;
+    Point point_1;
+    bool is_point_0_on_external_boundary = false;
+    bool is_point_1_on_external_boundary = false;
+    for (const auto e : _current_elem->edge_index_range())
+    {
+      _current_elem->build_edge_ptr(edge, e);
+      const auto & edge0 = edge->point(0);
+      const auto & edge1 = edge->point(1);
+      if (lineLineIntersect2D(
+              elem_center, -crack_dir, circumference, edge0, edge1, point_0, intersection_distance))
+        is_point_0_on_external_boundary = (_current_elem->neighbor_ptr(e) == nullptr);
+      else if (lineLineIntersect2D(elem_center,
+                                   crack_dir,
+                                   circumference,
+                                   edge0,
+                                   edge1,
+                                   point_1,
+                                   intersection_distance))
+        is_point_1_on_external_boundary = (_current_elem->neighbor_ptr(e) == nullptr);
+    }
+
+    // traverse_length is the length of a cut that goes across a single elment
+    // extend_length is the amount that needs to be added to each side of the traverse_length cut to
+    // make its length equal to the nucleation length
+    // We also add a small amount to the crack length to make sure it fully cuts the element edge.
+    // This factor is based on the element length.
+    Real traverse_length = (point_0 - point_1).norm();
+    Real extend_length = (_nucleation_length - traverse_length) / 2.0;
+    if (_nucleation_length < traverse_length)
+    {
+      if (_is_nucleation_length_provided && !_nucleate_across_full_element)
+        mooseError(
+            "Trying to nucleate crack smaller than element length, increase nucleation_length.\n  "
+            "Location of crack being nucleated: ",
+            point_0,
+            "\n  nucleation_length: ",
+            _nucleation_length,
+            "\n  Length to traverse element: ",
+            traverse_length,
+            "\n This error can be suppressed by setting nucleate_across_full_element=True which "
+            "will nucleate a crack the size of the element.");
+      else
+        //_edge_extension_factor is used to make sure cut will cut both edges of the element.
+        extend_length = (traverse_length * _edge_extension_factor) / 2.0;
+    }
+
+    // First create a cut assuming bulk nucleation
+    point_0 = point_0 - extend_length * crack_dir.unit();
+    point_1 = point_1 + extend_length * crack_dir.unit();
+
+    // modify edges of cut should go so that they are on the edge of the domain and have a length
+    // equal to nucleation_length
+    if (is_point_0_on_external_boundary)
+    {
+      point_0 = point_0 - (traverse_length * _edge_extension_factor) / 2.0 * crack_dir.unit();
+    }
+    else if (is_point_1_on_external_boundary)
+    {
+      point_1 = point_1 + (traverse_length * _edge_extension_factor) / 2.0 * crack_dir.unit();
+    }
     cutterElemNodes = {point_0, point_1};
   }
 
   return does_it_crack;
 }
+
+bool
+MeshCut2DRankTwoTensorNucleation::lineLineIntersect2D(const Point & start,
+                                                      const Point & direction,
+                                                      const Real length,
+                                                      const Point & v0,
+                                                      const Point & v1,
+                                                      Point & intersection_point,
+                                                      Real & intersection_distance)
+{
+  /// The standard "looser" tolerance to use in ray tracing when having difficulty finding intersection
+  const Real TRACE_TOLERANCE = 1e-5;
+
+  const auto r = direction * length;
+  const auto s = v1 - v0;
+  const auto rxs = r(0) * s(1) - r(1) * s(0);
+
+  // Lines are parallel or colinear
+  if (std::abs(rxs) < TRACE_TOLERANCE)
+    return false;
+
+  const auto v0mu0 = v0 - start;
+  const auto t = (v0mu0(0) * s(1) - v0mu0(1) * s(0)) / rxs;
+  if (0 >= t + TRACE_TOLERANCE || t - TRACE_TOLERANCE > 1.0)
+  {
+    return false;
+  }
+
+  const auto u = (v0mu0(0) * r(1) - v0mu0(1) * r(0)) / rxs;
+  if (0 < u + TRACE_TOLERANCE && u - TRACE_TOLERANCE <= 1.0)
+  {
+    intersection_point = start + r * t;
+    intersection_distance = t * length;
+    return true;
+  }
+
+  // Not parallel, but don't intersect
+  return false;
+}