fixes to xfem 3d cutter mesh

modules/xfem/include/efa/EFAElement3D.h

@@ -161,7 +161,6 @@ class EFAElement3D : public EFAElement
                             unsigned int edge_id,
                             double position,
                             EFANode * from_node,
-                            EFANode * embedded_node,
                             EFANode *& local_embedded);
   void mapParametricCoordinateFrom2DTo3D(unsigned int face_id,
                                          std::vector<double> & xi_2d,

modules/xfem/include/userobjects/CrackMeshCut3DUserObject.h

@@ -109,10 +109,13 @@ class CrackMeshCut3DUserObject : public MeshCutUserObjectBase
   const Real _const_intersection = 0.01;
 
   /// Used for cutter mesh refinement and front advancement
-  Real _size_control;
+  const Real _size_control;
 
   /// Number of steps to grow the mesh
-  unsigned int _n_step_growth;
+  const unsigned int _n_step_growth;
+
+  /// Minimum element area for crack growth elements
+  const Real _min_elem_area;
 
   /// Variables to help control the work flow
   bool _stop;
@@ -183,6 +186,16 @@ class CrackMeshCut3DUserObject : public MeshCutUserObjectBase
    */
   bool isInsideCutPlane(const std::vector<Point> & _vertices, const Point & p) const;
 
+  /**
+    Compute the area of a triangle defined by three points.
+    Returns true if the area is above min_elem_area.
+    @param p1 first vertex of the triangle
+    @param p2 second vertex of the triangle
+    @param p3 third vertex of the triangle
+    @return true if the triangle area exceeds the minimum threshold
+   */
+  bool isTriAreaAboveTol(const Point & p1, const Point & p2, const Point & p3) const;
+
   /**
     Find boundary nodes of the cutter mesh
     This is a simple algorithm simply based on the added angle = 360 degrees
@@ -206,11 +219,6 @@ class CrackMeshCut3DUserObject : public MeshCutUserObjectBase
    */
   Real findDistance(dof_id_type node1, dof_id_type node2);
 
-  /**
-    If boundary nodes are too sparse, add nodes in between
-   */
-  void refineBoundary();
-
   /**
     Find growth direction at each active node
    */
@@ -222,14 +230,66 @@ class CrackMeshCut3DUserObject : public MeshCutUserObjectBase
   void growFront();
 
   /**
-    Sort the front nodes
+    Determine whether a front node is an inactive endpoint of an active boundary segment.
+    @param front_node_index index of the node within the active boundary segment
+    @param front_size number of nodes in the active boundary segment
+    @return true if the node is an inactive endpoint
    */
-  void sortFrontNodes();
+  bool isInactiveEndpoint(unsigned int front_node_index, unsigned int front_size) const;
+
+  /**
+    Compute the crack growth increment for a front node.
+    @param front_node_index index of the node within the active boundary segment
+    @param front_size number of nodes in the active boundary segment
+    @param front_point_index mapping from active boundary nodes to crack-front-definition points
+    @return the growth increment applied to the node
+   */
+  Real computeGrowthIncrement(unsigned int front_node_index,
+                              unsigned int front_size,
+                              const std::vector<int> & front_point_index) const;
 
   /**
-    Find front-structure intersections
+    Project an inactive endpoint back outside the structural mesh if its proposed position lands
+    inside the volume.
+    @param segment_index index of the active boundary segment
+    @param front_node_index index of the node within the active boundary segment
+    @param front_size number of nodes in the active boundary segment
+    @param candidate_point proposed grown node position
+    @return the adjusted point position
    */
-  void findFrontIntersection();
+  Point projectInteriorInactiveEndpoint(unsigned int segment_index,
+                                        unsigned int front_node_index,
+                                        unsigned int front_size,
+                                        const Point & candidate_point) const;
+
+  /**
+    Add a grown front node while checking whether the newly formed crack-front triangles are
+    degenerate.
+    @param segment_index index of the active boundary segment
+    @param front_node_index index of the node within the active boundary segment
+    @param orig_id original boundary node id
+    @param candidate_point proposed grown node position
+    @param front_nodes output front-node ordering for the current active segment
+    @return the id that should represent this node on the grown front
+   */
+  dof_id_type
+  appendAdvancedFrontNodeCheckingDegenerateTriangles(unsigned int segment_index,
+                                                     unsigned int front_node_index,
+                                                     dof_id_type orig_id,
+                                                     const Point & candidate_point,
+                                                     std::vector<dof_id_type> & front_nodes);
+
+  /**
+    Update the tracked crack-front-definition node when a front node advances.
+    @param orig_id original crack front node id
+    @param new_id replacement node id after growth
+   */
+  void updateTrackedCrackFrontPoint(dof_id_type orig_id, dof_id_type new_id);
+
+  /**
+    Sort the front nodes
+   */
+  void sortFrontNodes();
 
   /**
     Refine the mesh at the front

modules/xfem/src/efa/EFAElement3D.C

@@ -2010,14 +2010,8 @@ EFAElement3D::addFaceEdgeCut(unsigned int face_id,
   {
     unsigned int emb_id = cut_edge->getEmbeddedNodeIndex(position, edge_node1);
     EFANode * old_emb = cut_edge->getEmbeddedNode(emb_id);
-    if (embedded_node && embedded_node != old_emb)
-      EFAError("Attempting to add edge intersection when one already exists with different node.",
-               " elem: ",
-               _id,
-               " edge: ",
-               edge_id,
-               " position: ",
-               position);
+    // If the same physical edge point is rediscovered through a different face/neighbor path,
+    // reuse the existing embedded node on that edge instead of treating it as a contradictory cut.
     local_embedded = old_emb;
     cut_exist = true;
   }
@@ -2026,9 +2020,8 @@ EFAElement3D::addFaceEdgeCut(unsigned int face_id,
       isPhysicalEdgeCut(face_id, edge_id, position))
   {
     // check if cut has already been added to the neighbor edges
-    checkNeighborFaceCut(face_id, edge_id, position, edge_node1, embedded_node, local_embedded);
-    checkNeighborFaceCut(
-        adj_face_id, adj_edge_id, position, edge_node1, embedded_node, local_embedded);
+    checkNeighborFaceCut(face_id, edge_id, position, edge_node1, local_embedded);
+    checkNeighborFaceCut(adj_face_id, adj_edge_id, position, edge_node1, local_embedded);
 
     if (!local_embedded) // need to create new embedded node
     {
@@ -2123,7 +2116,6 @@ EFAElement3D::checkNeighborFaceCut(unsigned int face_id,
                                    unsigned int edge_id,
                                    double position,
                                    EFANode * from_node,
-                                   EFANode * embedded_node,
                                    EFANode *& local_embedded)
 {
   // N.B. this is important. We are checking if the corresponding edge of the neighbor face or of
@@ -2142,12 +2134,8 @@ EFAElement3D::checkNeighborFaceCut(unsigned int face_id,
       unsigned int emb_id = neigh_edge->getEmbeddedNodeIndex(position, from_node);
       EFANode * old_emb = neigh_edge->getEmbeddedNode(emb_id);
 
-      if (embedded_node && embedded_node != old_emb)
-        EFAError(
-            "attempting to add edge intersection when one already exists with different node.");
-      if (local_embedded && local_embedded != old_emb)
-        EFAError("attempting to assign contradictory pointer to local_embedded.");
-
+      // Same edge and same position on a neighboring face should canonicalize to the existing
+      // embedded node, regardless of which propagation path discovered it first.
       local_embedded = old_emb;
     }
   } // en_iter