diff --git a/modules/xfem/include/efa/EFAElement3D.h b/modules/xfem/include/efa/EFAElement3D.h index 68df51ca63c1..8d5dbd54c264 100644 --- a/modules/xfem/include/efa/EFAElement3D.h +++ b/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 & xi_2d, diff --git a/modules/xfem/include/userobjects/CrackMeshCut3DUserObject.h b/modules/xfem/include/userobjects/CrackMeshCut3DUserObject.h index 25680bc4dbb5..4fcb700496a7 100644 --- a/modules/xfem/include/userobjects/CrackMeshCut3DUserObject.h +++ b/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 & _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 & 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 & 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 diff --git a/modules/xfem/src/efa/EFAElement3D.C b/modules/xfem/src/efa/EFAElement3D.C index 1fce96aea187..3cf84cb0534f 100644 --- a/modules/xfem/src/efa/EFAElement3D.C +++ b/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 diff --git a/modules/xfem/src/userobjects/CrackMeshCut3DUserObject.C b/modules/xfem/src/userobjects/CrackMeshCut3DUserObject.C index baebf446acb9..cbecb2860198 100644 --- a/modules/xfem/src/userobjects/CrackMeshCut3DUserObject.C +++ b/modules/xfem/src/userobjects/CrackMeshCut3DUserObject.C @@ -11,6 +11,7 @@ #include "XFEMFuncs.h" #include "MooseError.h" +#include "libmesh/libmesh_common.h" #include "libmesh/string_to_enum.h" #include "MooseMesh.h" #include "MooseEnum.h" @@ -51,19 +52,21 @@ CrackMeshCut3DUserObject::validParams() params.addParam( "size_control", 0, "Criterion for refining elements while growing the crack"); params.addParam("n_step_growth", 0, "Number of steps for crack growth"); + params.addParam( + "min_elem_area", 1e-6, "Growth elements smaller than min_elem_area will not be added."); params.addClassDescription("Creates a UserObject for a mesh cutter in 3D problems"); return params; } -// This code does not allow predefined crack growth as a function of time -// all inital cracks are defined at t_start = t_end = 0 CrackMeshCut3DUserObject::CrackMeshCut3DUserObject(const InputParameters & parameters) : MeshCutUserObjectBase(parameters), _mesh(_subproblem.mesh()), _growth_dir_method(getParam("growth_dir_method").getEnum()), _growth_increment_method( getParam("growth_increment_method").getEnum()), + _size_control(getParam("size_control")), _n_step_growth(getParam("n_step_growth")), + _min_elem_area(getParam("min_elem_area")), _is_mesh_modified(false), _func_x(parameters.isParamValid("growth_direction_x") ? &getFunction("growth_direction_x") : nullptr), @@ -91,11 +94,8 @@ CrackMeshCut3DUserObject::CrackMeshCut3DUserObject(const InputParameters & param if (_grow) { - if (!isParamValid("size_control")) + if (!isParamSetByUser("size_control")) paramError("size_control", "Crack growth needs size control."); - - _size_control = getParam("size_control"); - if (_growth_dir_method == GrowthDirectionEnum::FUNCTION && (_func_x == nullptr || _func_y == nullptr || _func_z == nullptr)) mooseError("function is not specified for the function method that defines growth direction"); @@ -158,8 +158,6 @@ CrackMeshCut3DUserObject::initialize() _is_mesh_modified = true; growFront(); sortFrontNodes(); - if (_inactive_boundary_pos.size() != 0) - findFrontIntersection(); refineFront(); triangulation(); joinBoundary(); @@ -185,7 +183,6 @@ CrackMeshCut3DUserObject::cutElementByGeometry(const Elem * elem, std::vector & cut_faces) const // With the crack defined by a planar mesh, this method cuts a solid element by all elements in the // planar mesh -// TODO: Time evolving cuts not yet supported in 3D (hence the lack of use of the time variable) { bool elem_cut = false; @@ -349,6 +346,15 @@ CrackMeshCut3DUserObject::getRelativePosition(const Point & p1, return len_p1_p / full_len; } +bool +CrackMeshCut3DUserObject::isTriAreaAboveTol(const Point & p1, + const Point & p2, + const Point & p3) const +{ + Real area = 0.5 * ((p2 - p1).cross(p3 - p1)).norm(); + return area >= _min_elem_area; +} + bool CrackMeshCut3DUserObject::isInsideCutPlane(const std::vector & vertices, const Point & p) const @@ -590,60 +596,6 @@ CrackMeshCut3DUserObject::findDistance(dof_id_type node1, dof_id_type node2) return distance; } -void -CrackMeshCut3DUserObject::refineBoundary() -{ - std::vector new_boundary_order(_boundary.begin(), _boundary.end()); - - mooseAssert(_boundary.size() >= 2, "Boundary must be at least two nodes"); - - for (unsigned int i = _boundary.size() - 1; i >= 1; --i) - { - dof_id_type node1 = _boundary[i - 1]; - dof_id_type node2 = _boundary[i]; - - Real distance = findDistance(node1, node2); - - if (distance > _size_control) - { - unsigned int n = static_cast(distance / _size_control); - std::array x1; - std::array x2; - - Node * n1 = _cutter_mesh->node_ptr(node1); - mooseAssert(n1 != nullptr, "Node is NULL"); - Point & p1 = *n1; - Node * n2 = _cutter_mesh->node_ptr(node2); - mooseAssert(n2 != nullptr, "Node is NULL"); - Point & p2 = *n2; - - for (unsigned int j = 0; j < 3; ++j) - { - x1[j] = p1(j); - x2[j] = p2(j); - } - - for (unsigned int j = 0; j < n; ++j) - { - Point x; - for (unsigned int k = 0; k < 3; ++k) - x(k) = x2[k] - (x2[k] - x1[k]) * (j + 1) / (n + 1); - - Node * this_node = Node::build(x, _cutter_mesh->n_nodes()).release(); - _cutter_mesh->add_node(this_node); - - dof_id_type id = _cutter_mesh->n_nodes() - 1; - auto it = new_boundary_order.begin(); - new_boundary_order.insert(it + i, id); - } - } - } - - _boundary = new_boundary_order; - mooseAssert(_boundary.size() > 0, "Boundary should not have zero size"); - _boundary.pop_back(); -} - void CrackMeshCut3DUserObject::findActiveBoundaryNodes() { @@ -697,7 +649,9 @@ CrackMeshCut3DUserObject::findActiveBoundaryNodes() std::vector temp; for (unsigned int j = _inactive_boundary_pos[n_inactive_boundary - 1]; j < n_boundary; ++j) temp.push_back(_boundary[j]); - for (unsigned int j = 0; j <= _inactive_boundary_pos[0]; ++j) + // _boundary is closed (_boundary[n_boundary - 1] == _boundary[0]); skip j=0 to avoid + // pushing the same node twice in a row at the wrap-around + for (unsigned int j = 1; j <= _inactive_boundary_pos[0]; ++j) temp.push_back(_boundary[j]); _active_boundary.push_back(temp); } @@ -718,13 +672,7 @@ CrackMeshCut3DUserObject::findActiveBoundaryDirection() std::vector temp; Point dir; - if (_inactive_boundary_pos.size() != 0) - { - for (unsigned int j = 0; j < 3; ++j) - dir(j) = 0; - temp.push_back(dir); - } - + // Compute direction for all interior active nodes unsigned int i1 = 1; unsigned int i2 = _active_boundary[i].size() - 1; if (_inactive_boundary_pos.size() == 0) @@ -740,9 +688,9 @@ CrackMeshCut3DUserObject::findActiveBoundaryDirection() Node * this_node = _cutter_mesh->node_ptr(_active_boundary[i][j]); mooseAssert(this_node, "Node is NULL"); Point & this_point = *this_node; - dir(0) = _func_x->value(0, this_point); - dir(1) = _func_y->value(0, this_point); - dir(2) = _func_z->value(0, this_point); + dir(0) = _func_x->value(_t, this_point); + dir(1) = _func_y->value(_t, this_point); + dir(2) = _func_z->value(_t, this_point); temp.push_back(dir); } @@ -810,11 +758,11 @@ CrackMeshCut3DUserObject::findActiveBoundaryDirection() else mooseError("This growth_dir_method is not pre-defined!"); - if (_inactive_boundary_pos.size() != 0) + // Boundary nodes copy direction from their connected active neighbor + if (_inactive_boundary_pos.size() != 0 && temp.size() > 0) { - for (unsigned int j = 0; j < 3; ++j) - dir(j) = 0; - temp.push_back(dir); + temp.insert(temp.begin(), temp.front()); + temp.push_back(temp.back()); } _active_direction.push_back(temp); @@ -842,203 +790,250 @@ CrackMeshCut3DUserObject::growFront() { _front.clear(); + mooseAssert(!(_cfd && _active_boundary.size() != 1), + "crack-front-definition using the cutter mesh only supports one active crack front " + "segment for now"); + + const std::vector front_point_index = _cfd ? getFrontPointsIndex() : std::vector(); + + // Track already-grown boundary nodes to avoid duplicates when a boundary node + // appears in two active segments (as the last node of one and the first of the next) + std::map grown_node_map; + for (unsigned int i = 0; i < _active_boundary.size(); ++i) { std::vector temp; - unsigned int i1 = 1; - unsigned int i2 = _active_boundary[i].size() - 1; - if (_inactive_boundary_pos.size() == 0) + unsigned int i2 = _active_boundary[i].size(); + for (unsigned int j = 0; j < i2; ++j) { - i1 = 0; - i2 = _active_boundary[i].size(); - } + dof_id_type orig_id = _active_boundary[i][j]; - std::vector index = getFrontPointsIndex(); - for (unsigned int j = i1; j < i2; ++j) - { - Node * this_node = _cutter_mesh->node_ptr(_active_boundary[i][j]); + // If this node was already grown in a previous segment, reuse the grown node + auto git = grown_node_map.find(orig_id); + if (git != grown_node_map.end()) + { + temp.push_back(git->second); + continue; + } + + Node * this_node = _cutter_mesh->node_ptr(orig_id); mooseAssert(this_node, "Node is NULL"); Point & this_point = *this_node; Point dir = _active_direction[i][j]; + const Real growth_increment = computeGrowthIncrement(j, i2, front_point_index); + Point x; - Real growth_increment = 0; - switch (_growth_increment_method) - { - case GrowthRateEnum::FUNCTION: - { - growth_increment = _func_v->value(0, Point(0, 0, 0)); - break; - } - case GrowthRateEnum::REPORTER: - { - int ind = index[j]; - if (index[j] == -1) - growth_increment = 0; - else - growth_increment = _growth_inc_reporter->at(ind); - break; - } - default: - { - mooseError("This growth_increment_method is not pre-defined!"); - break; - } - } for (unsigned int k = 0; k < 3; ++k) x(k) = this_point(k) + dir(k) * growth_increment; - this_node = Node::build(x, _cutter_mesh->n_nodes()).release(); - _cutter_mesh->add_node(this_node); + x = projectInteriorInactiveEndpoint(i, j, i2, x); - dof_id_type id = _cutter_mesh->n_nodes() - 1; - temp.push_back(id); - - if (_cfd) - { - auto it = std::find(_tracked_crack_front_points.begin(), - _tracked_crack_front_points.end(), - _active_boundary[0][j]); - if (it != _tracked_crack_front_points.end()) - { - unsigned int pos = std::distance(_tracked_crack_front_points.begin(), it); - _tracked_crack_front_points[pos] = id; - } - } + const dof_id_type id = + appendAdvancedFrontNodeCheckingDegenerateTriangles(i, j, orig_id, x, temp); + grown_node_map[orig_id] = id; + updateTrackedCrackFrontPoint(orig_id, id); } _front.push_back(temp); } } -void -CrackMeshCut3DUserObject::sortFrontNodes() -// TBD; it is not needed for current problems but will be useful for fracture growth +bool +CrackMeshCut3DUserObject::isInactiveEndpoint(unsigned int front_node_index, + unsigned int front_size) const { + return _inactive_boundary_pos.size() != 0 && + (front_node_index == 0 || front_node_index + 1 == front_size); } -void -CrackMeshCut3DUserObject::findFrontIntersection() +Real +CrackMeshCut3DUserObject::computeGrowthIncrement(unsigned int front_node_index, + unsigned int front_size, + const std::vector & front_point_index) const { - ConstBndElemRange & range = *_mesh.getBoundaryElementRange(); - - for (unsigned int i = 0; i < _front.size(); ++i) + switch (_growth_increment_method) { - if (_front[i].size() >= 2) + case GrowthRateEnum::FUNCTION: + return _func_v->value(_t, Point(0, 0, 0)); + + case GrowthRateEnum::REPORTER: { - std::vector pint1; - std::vector pint2; - std::vector length1; - std::vector length2; + mooseAssert(!front_point_index.empty(), + "Crack-front indices must be available for reporter-based growth."); - Real node_id = _front[i][0]; - Node * this_node = _cutter_mesh->node_ptr(node_id); - mooseAssert(this_node, "Node is NULL"); - Point & p2 = *this_node; + if (isInactiveEndpoint(front_node_index, front_size)) + { + const bool has_adjacent_front_point = front_point_index.size() > 1; + const bool is_first_inactive_endpoint = + front_node_index == 0 && has_adjacent_front_point && front_point_index[1] != -1; + const bool is_last_inactive_endpoint = front_node_index + 1 == front_size && + has_adjacent_front_point && + front_point_index[front_size - 2] != -1; + // Boundary node: use increment from nearest active neighbor + if (is_first_inactive_endpoint) + return _growth_inc_reporter->at(front_point_index[1]); + if (is_last_inactive_endpoint) + return _growth_inc_reporter->at(front_point_index[front_size - 2]); + + mooseError("Inactive crackfront node is not connected to a neighbor."); + } - if (_front[i].size() >= 4) - node_id = _front[i][2]; - else - node_id = _front[i][1]; + return _growth_inc_reporter->at(front_point_index[front_node_index]); + } + } - this_node = _cutter_mesh->node_ptr(node_id); - mooseAssert(this_node, "Node is NULL"); - Point & p1 = *this_node; + mooseError("This growth_increment_method is not pre-defined!"); +} - node_id = _front[i].back(); - this_node = _cutter_mesh->node_ptr(node_id); - mooseAssert(this_node, "Node is NULL"); - Point & p4 = *this_node; +Point +CrackMeshCut3DUserObject::projectInteriorInactiveEndpoint(unsigned int segment_index, + unsigned int front_node_index, + unsigned int front_size, + const Point & candidate_point) const +{ + if (!isInactiveEndpoint(front_node_index, front_size)) + return candidate_point; - if (_front[i].size() >= 4) - node_id = _front[i][_front[i].size() - 3]; - else - node_id = _front[i][_front[i].size() - 2]; + std::unique_ptr pl = _mesh.getPointLocator(); + pl->enable_out_of_mesh_mode(); + if ((*pl)(candidate_point) == nullptr) + return candidate_point; - this_node = _cutter_mesh->node_ptr(node_id); - mooseAssert(this_node, "Node is NULL"); - Point & p3 = *this_node; + const bool at_start_of_crack_front = (front_node_index == 0 && front_size > 1); + const bool at_end_of_crack_front = (front_node_index + 1 == front_size && front_size > 1); + if (!(at_start_of_crack_front || at_end_of_crack_front)) + return candidate_point; - bool do_inter1 = 1; - bool do_inter2 = 1; + const unsigned int active_neighbor_index = at_start_of_crack_front ? 1 : (front_size - 2); + Node * active_nd = _cutter_mesh->node_ptr(_active_boundary[segment_index][active_neighbor_index]); + mooseAssert(active_nd, "Active neighbor node is NULL"); - std::unique_ptr pl = _mesh.getPointLocator(); - pl->enable_out_of_mesh_mode(); - const Elem * elem = (*pl)(p1); - if (elem == nullptr) - do_inter1 = 0; - elem = (*pl)(p4); - if (elem == nullptr) - do_inter2 = 0; + Point projected_point = candidate_point; + const Point A = *active_nd; + Point ray_dir = projected_point - A; + const Real ray_len = ray_dir.norm(); + if (ray_len <= libMesh::TOLERANCE) + return projected_point; + ray_dir /= ray_len; - for (const auto & belem : range) - { - Point pt; - std::vector vertices; + ConstBndElemRange & bnd_range = *_mesh.getBoundaryElementRange(); - elem = belem->_elem; - std::unique_ptr curr_side = elem->side_ptr(belem->_side); - for (unsigned int j = 0; j < curr_side->n_nodes(); ++j) - { - const Node * node = curr_side->node_ptr(j); - const Point & this_point = *node; - vertices.push_back(this_point); - } - - if (findIntersection(p1, p2, vertices, pt)) - { - pint1.push_back(pt); - length1.push_back((pt - p1) * (pt - p1)); - } - if (findIntersection(p3, p4, vertices, pt)) - { - pint2.push_back(pt); - length2.push_back((pt - p3) * (pt - p3)); - } - } + for (const auto & belem : bnd_range) + { + const Elem * elem = belem->_elem; + std::unique_ptr curr_side = elem->side_ptr(belem->_side); + std::vector vertices; + for (unsigned int i = 0; i < curr_side->n_nodes(); ++i) + vertices.push_back(*(curr_side->node_ptr(i))); + + Point pt; + if (findIntersection(A, projected_point, vertices, pt)) + return pt + ray_dir * libMesh::TOLERANCE; + } - if (length1.size() != 0 && do_inter1) - { - auto it1 = std::min_element(length1.begin(), length1.end()); - Point inter1 = pint1[std::distance(length1.begin(), it1)]; - inter1 += (inter1 - p1) * _const_intersection; + return projected_point; +} - Node * this_node = Node::build(inter1, _cutter_mesh->n_nodes()).release(); - _cutter_mesh->add_node(this_node); +dof_id_type +CrackMeshCut3DUserObject::appendAdvancedFrontNodeCheckingDegenerateTriangles( + unsigned int segment_index, + unsigned int front_node_index, + dof_id_type orig_id, + const Point & candidate_point, + std::vector & front_nodes) +{ + Node * this_node = _cutter_mesh->node_ptr(orig_id); + mooseAssert(this_node, "Node is NULL"); + const Point & this_point = *this_node; + + const unsigned int front_size = _active_boundary[segment_index].size(); + const dof_id_type prev_id = (front_node_index > 0) + ? _active_boundary[segment_index][front_node_index - 1] + : DofObject::invalid_id; + const dof_id_type next_id = (front_node_index + 1 < front_size) + ? _active_boundary[segment_index][front_node_index + 1] + : DofObject::invalid_id; + + bool prev_degenerate = false; + unsigned int n_relevant_neighbors = 0; + unsigned int n_degenerate = 0; + + if (prev_id != DofObject::invalid_id) + { + ++n_relevant_neighbors; + Point prev_pt = *_cutter_mesh->node_ptr(prev_id); + if (!isTriAreaAboveTol(this_point, prev_pt, candidate_point)) + { + ++n_degenerate; + prev_degenerate = true; + } + } - mooseAssert(_cutter_mesh->n_nodes() - 1 > 0, "The cut mesh must be at least one element."); - unsigned int n = _cutter_mesh->n_nodes() - 1; + if (next_id != DofObject::invalid_id) + { + ++n_relevant_neighbors; + Point next_pt = *_cutter_mesh->node_ptr(next_id); + if (!isTriAreaAboveTol(this_point, next_pt, candidate_point)) + ++n_degenerate; + } - auto it = _front[i].begin(); - _front[i].insert(it, n); + const bool all_degenerate = (n_relevant_neighbors > 0 && n_degenerate == n_relevant_neighbors); + const bool some_degenerate = (n_degenerate > 0 && !all_degenerate); - if (_cfd) - _tracked_crack_front_points[_tracked_crack_front_points.size() - 1] = n; - } + if (all_degenerate) + { + front_nodes.push_back(orig_id); + return orig_id; + } - if (length2.size() != 0 && do_inter2) - { - auto it2 = std::min_element(length2.begin(), length2.end()); - Point inter2 = pint2[std::distance(length2.begin(), it2)]; - inter2 += (inter2 - p2) * _const_intersection; + Node * new_node = Node::build(candidate_point, _cutter_mesh->n_nodes()).release(); + _cutter_mesh->add_node(new_node); + const dof_id_type new_id = _cutter_mesh->n_nodes() - 1; - Node * this_node = Node::build(inter2, _cutter_mesh->n_nodes()).release(); - _cutter_mesh->add_node(this_node); + if (some_degenerate) + { + // Place orig adjacent to the degenerate (stuck) neighbor, new adjacent to the non-degenerate + // neighbor, so triangulation pairs the right nodes on each edge. + if (prev_degenerate) + { + front_nodes.push_back(orig_id); + front_nodes.push_back(new_id); + } + else + { + front_nodes.push_back(new_id); + front_nodes.push_back(orig_id); + } + } + else + front_nodes.push_back(new_id); - dof_id_type n = _cutter_mesh->n_nodes() - 1; + return new_id; +} - auto it = _front[i].begin(); - unsigned int m = _front[i].size(); - _front[i].insert(it + m, n); +void +CrackMeshCut3DUserObject::updateTrackedCrackFrontPoint(dof_id_type orig_id, dof_id_type new_id) +{ + if (!_cfd) + return; - if (_cfd) - _tracked_crack_front_points[0] = n; - } - } + auto it = + std::find(_tracked_crack_front_points.begin(), _tracked_crack_front_points.end(), orig_id); + if (it != _tracked_crack_front_points.end()) + { + const unsigned int pos = std::distance(_tracked_crack_front_points.begin(), it); + _tracked_crack_front_points[pos] = new_id; } } +void +CrackMeshCut3DUserObject::sortFrontNodes() +// TBD; it is not needed for current problems but will be useful for fracture growth +{ +} + void CrackMeshCut3DUserObject::refineFront() { @@ -1050,6 +1045,14 @@ CrackMeshCut3DUserObject::refineFront() if (_inactive_boundary_pos.size() == 0) i1 = _front[ifront].size(); + // When there are inactive boundary nodes, determine which segment indices are + // adjacent to the boundary endpoints (first and last in the segment). + // Refinement nodes on those segments may land outside the FEM mesh and be + // incorrectly treated as inactive endpoints on the next step. + bool has_inactive = (_inactive_boundary_pos.size() != 0); + unsigned int first_boundary_i = 1; // segment between pos 0 and 1 + unsigned int last_boundary_i = _front[ifront].size() - 1; // segment between pos size-2 and size-1 + for (unsigned int i = i1; i >= 1; --i) { unsigned int i2 = i; @@ -1079,12 +1082,31 @@ CrackMeshCut3DUserObject::refineFront() x2[j] = p2(j); } + // Get the corresponding old boundary node for area check + unsigned int ab_idx = (i - 1 < _active_boundary[ifront].size()) ? i - 1 : 0; + Point ab_pt = *_cutter_mesh->node_ptr(_active_boundary[ifront][ab_idx]); + for (unsigned int j = 0; j < n; ++j) { Point x; for (unsigned int k = 0; k < 3; ++k) x(k) = x2[k] - (x2[k] - x1[k]) * (j + 1) / (n + 1); + // Skip if both triangles with old boundary node would be degenerate + if (!isTriAreaAboveTol(x, p1, ab_pt) && !isTriAreaAboveTol(x, p2, ab_pt)) + continue; + + // For segments adjacent to inactive boundary endpoints, skip refinement nodes + // that are outside the FEM mesh. A node outside would be treated as an + // inactive endpoint on the next step, corrupting the active boundary structure. + if (has_inactive && (i == first_boundary_i || i == last_boundary_i)) + { + std::unique_ptr pl = _mesh.getPointLocator(); + pl->enable_out_of_mesh_mode(); + if ((*pl)(x) == nullptr) // outside the FEM mesh + continue; + } + Node * this_node = Node::build(x, _cutter_mesh->n_nodes()).release(); _cutter_mesh->add_node(this_node); @@ -1197,6 +1219,12 @@ CrackMeshCut3DUserObject::triangulation() i1++; } + // Skip degenerate elements + if (!isTriAreaAboveTol(Point(*_cutter_mesh->node_ptr(elem[0])), + Point(*_cutter_mesh->node_ptr(elem[1])), + Point(*_cutter_mesh->node_ptr(elem[2])))) + continue; + Elem * new_elem = Elem::build(TRI3).release(); for (unsigned int i = 0; i < _cut_elem_nnode; ++i) diff --git a/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_domain.i b/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_domain.i index b1bd322c139e..9564143258f0 100644 --- a/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_domain.i +++ b/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_domain.i @@ -119,8 +119,4 @@ [Outputs] execute_on = 'timestep_end' exodus = true - [console] - type = Console - output_linear = true - [] [] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_fatigue.i b/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_fatigue.i index 29940d537906..7a035dbf6197 100644 --- a/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_fatigue.i +++ b/modules/xfem/test/tests/solid_mechanics_basic/edge_crack_3d_fatigue.i @@ -4,7 +4,7 @@ growth_increment_name = "growth_increment" cycles_to_max_growth_increment_name = "fatigue" crackMeshCut3DUserObject_name = cut_mesh - max_growth_increment = 0.1 + max_growth_increment = 0.1001 paris_law_c = 1e-13 paris_law_m = 2.5 [] @@ -68,5 +68,4 @@ [Outputs] file_base = edge_crack_3d_fatigue_out - json = true [] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_domain.i b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_domain.i new file mode 100644 index 000000000000..8400a4edf4ce --- /dev/null +++ b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_domain.i @@ -0,0 +1,165 @@ +[GlobalParams] + displacements = 'disp_x disp_y disp_z' + volumetric_locking_correction = true +[] + +rad=0.1 +offset = 0 +spin = 10 +[Mesh] + #---- CUTTER MESH + [cicle_outline] + type = ParsedCurveGenerator + x_formula = '${rad}*cos(t)' + y_formula = '${rad}*sin(t)' + section_bounding_t_values = '0 ${fparse 2*pi}' + nums_segments = '10' + is_closed_loop = true + [] + [circle_surface] + type = XYDelaunayGenerator + boundary = 'cicle_outline' + desired_area = 0.001 + output_subdomain_id = 1 + [] + [circle_spin] + type = TransformGenerator + input = circle_surface + transform = ROTATE + vector_value = '0 0 ${spin}' + [] + [circle_rotate] + type = TransformGenerator + input = circle_spin + transform = ROTATE + vector_value = '0 90 0' + [] + [circle_move] + type = TransformGenerator + input = circle_rotate + transform = TRANSLATE + vector_value = '${offset} 0 -0.01' + save_with_name = mesh_cutter + [] + #---- FEM MESH + [FEM_mesh] + type = GeneratedMeshGenerator + dim = 3 + nx = 14 + ny = 3 + nz = 8 + xmin = ${fparse -2*rad} + xmax = ${fparse 2*rad} + ymin = ${fparse -rad} + ymax = ${fparse rad} + zmin = 0 + zmax = ${fparse 2*rad} + elem_type = HEX8 + [] + [pin] + type = ExtraNodesetGenerator + input = FEM_mesh + new_boundary = 'pin' + coord = '${fparse 2*rad} ${fparse -rad} ${fparse rad}' + use_closest_node = true + [] + [center] + type = ExtraNodesetGenerator + input = pin + new_boundary = 'center' + coord = '0 ${fparse rad} 0' + use_closest_node = true + [] + final_generator = center +[] + +[Physics/SolidMechanics/QuasiStatic] + [all] + strain = FINITE + add_variables = true + generate_output = 'stress_xx stress_yy stress_zz vonmises_stress' + [] +[] + +[BCs] + [left_x] + type = DirichletBC + boundary = left + variable = disp_x + value = 0.0 + [] + [right_x] + type = DirichletBC + boundary = right + variable = disp_x + value = 0.0 + [] + # Ramp traction in z so KI varies along the crack front: + # near z=0 (surface, inactive nodes): below k_low -> no growth + # mid-front: transition zone -> moderate growth + # deepest nodes (z~0.09): above k_high -> max growth + [top_y] + type = FunctionNeumannBC + boundary = top + variable = disp_y + function = '20 + 500 * z' + [] + [bottom_y] + type = DirichletBC + boundary = bottom + variable = disp_y + value = 0 + [] + [pin_z] + type = DirichletBC + boundary = pin + variable = disp_z + value = 0.0 + [] +[] + +[Materials] + [elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 207000 + poissons_ratio = 0.3 + block = 0 + [] + [stress] + type = ComputeFiniteStrainElasticStress + block = 0 + [] +[] + +[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 = 'none' + + [Predictor] + type = SimplePredictor + scale = 1.0 + [] + + # controls for linear iterations + l_max_its = 100 + l_tol = 1e-2 + + # controls for nonlinear iterations + nl_max_its = 15 + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-10 + + # time control + start_time = 0.0 + dt = 1.0 + end_time = 4 + max_xfem_update = 1 +[] + +[Outputs] +[] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_function.i b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_function.i new file mode 100644 index 000000000000..3431b6a32ecb --- /dev/null +++ b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_function.i @@ -0,0 +1,48 @@ + +[XFEM] + geometric_cut_userobjects = 'cut_mesh' + qrule = volfrac + output_cut_plane = true +[] + +[Functions] + [growth_func_x] + type = ParsedFunction + expression = '0.2*(x-${offset})/abs(${rad})' + [] + [growth_func_y] + type = ParsedFunction + expression = 0.0 + [] + [growth_func_z] + type = ParsedFunction + expression = 'if(t<3,0,0.1)' + [] + [growth_func_v] + type = ParsedFunction + expression = 'if(t<3,0.02,0.015)' + [] +[] + +[UserObjects] + [cut_mesh] + type = CrackMeshCut3DUserObject + mesh_generator_name = mesh_cutter + growth_dir_method = FUNCTION + growth_direction_x = growth_func_x + growth_direction_y = growth_func_y + growth_direction_z = growth_func_z + size_control = .05 + n_step_growth = 1 + growth_rate = growth_func_v + [] +[] + +[Outputs] + file_base = face_crack_function_off_${offset}_spin_${spin} + [xfemcutter] + type = XFEMCutMeshOutput + xfem_cutter_uo = cut_mesh + [] + execute_on=final +[] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_hoopStress.i b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_hoopStress.i new file mode 100644 index 000000000000..d98aa0ee00af --- /dev/null +++ b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_hoopStress.i @@ -0,0 +1,39 @@ + +[XFEM] + geometric_cut_userobjects = 'cut_mesh' + qrule = volfrac + output_cut_plane = true +[] + +[UserObjects] + [cut_mesh] + type = CrackMeshCut3DUserObject + mesh_generator_name = mesh_cutter + growth_dir_method = MAX_HOOP_STRESS + size_control = .05 + n_step_growth = 1 + growth_rate = 0.02 + [] +[] + +[DomainIntegral] + integrals = 'Jintegral InteractionIntegralKI InteractionIntegralKII' + displacements = 'disp_x disp_y disp_z' + crack_front_points_provider = cut_mesh + crack_direction_method = CurvedCrackFront + radius_inner = '0.025' + radius_outer = '0.1' + poissons_ratio = 0.3 + youngs_modulus = 207000 + block = 0 + incremental = true +[] + +[Outputs] + file_base = face_crack_hoopStress_off_${offset}_spin_${spin} + [xfemcutter] + type = XFEMCutMeshOutput + xfem_cutter_uo = cut_mesh + [] + execute_on=final +[] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_radialFunc.i b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_radialFunc.i new file mode 100644 index 000000000000..9e30b2fa76fc --- /dev/null +++ b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_radialFunc.i @@ -0,0 +1,48 @@ + +[XFEM] + geometric_cut_userobjects = 'cut_mesh' + qrule = volfrac + output_cut_plane = true +[] + +[Functions] + [growth_func_x] + type = ParsedFunction + value = '(x-${offset})/sqrt((x-${offset})^2 + z^2)' + [] + [growth_func_y] + type = ParsedFunction + expression = 0.0 + [] + [growth_func_z] + type = ParsedFunction + value = 'z/sqrt((x-${offset})^2 + z^2)' + [] + [growth_func_v] + type = ParsedFunction + expression = '.02' + [] +[] + +[UserObjects] + [cut_mesh] + type = CrackMeshCut3DUserObject + mesh_generator_name = mesh_cutter + growth_dir_method = FUNCTION + growth_direction_x = growth_func_x + growth_direction_y = growth_func_y + growth_direction_z = growth_func_z + size_control = .05 + n_step_growth = 1 + growth_rate = growth_func_v + [] +[] + +[Outputs] + file_base = face_crack_radialFunc_off_${offset}_spin_${spin} + [xfemcutter] + type = XFEMCutMeshOutput + xfem_cutter_uo = cut_mesh + [] + execute_on=final +[] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_scc.i b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_scc.i new file mode 100644 index 000000000000..e3996c25d197 --- /dev/null +++ b/modules/xfem/test/tests/solid_mechanics_basic/face_crack_3d_xfem_scc.i @@ -0,0 +1,54 @@ + +[XFEM] + geometric_cut_userobjects = 'cut_mesh' + qrule = volfrac + output_cut_plane = true +[] + +[Reporters] + [scc_crack_growth] + type = StressCorrosionCrackingExponential + growth_increment_name = "crack_growth" + time_to_max_growth_increment_name = "max_growth_timestep" + crackMeshCut3DUserObject_name = cut_mesh + max_growth_increment = 0.02 + k_low = 10 + k_high = 20 + growth_rate_mid_multiplier = 0.00075 + growth_rate_mid_exp_factor = 1 + [] +[] + +[UserObjects] + [cut_mesh] + type = CrackMeshCut3DUserObject + mesh_generator_name = mesh_cutter + growth_dir_method = MAX_HOOP_STRESS + size_control = .05 + n_step_growth = 1 + growth_increment_method = REPORTER + growth_reporter = "scc_crack_growth/crack_growth" + [] +[] + +[DomainIntegral] + integrals = 'Jintegral InteractionIntegralKI InteractionIntegralKII' + displacements = 'disp_x disp_y disp_z' + crack_front_points_provider = cut_mesh + crack_direction_method = CurvedCrackFront + radius_inner = '0.025' + radius_outer = '0.1' + poissons_ratio = 0.3 + youngs_modulus = 207000 + block = 0 + incremental = true +[] + +[Outputs] + file_base = face_crack_scc_off_${offset}_spin_${spin} + [xfemcutter] + type = XFEMCutMeshOutput + xfem_cutter_uo = cut_mesh + [] + execute_on=final +[] diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/double_cantilever_crack_out_II_KI_1_0009.csv b/modules/xfem/test/tests/solid_mechanics_basic/gold/double_cantilever_crack_out_II_KI_1_0009.csv index c5b75b20d7ef..c6748e9d9e82 100644 --- a/modules/xfem/test/tests/solid_mechanics_basic/gold/double_cantilever_crack_out_II_KI_1_0009.csv +++ b/modules/xfem/test/tests/solid_mechanics_basic/gold/double_cantilever_crack_out_II_KI_1_0009.csv @@ -1,8 +1,8 @@ II_KI_1,id,x,y,z -4176.7831338884,0,5.4376712625919,0.6902300599503,2.0334835 -2254.7722712718,0.6851334999989,5.4376712625897,0.69023005995375,1.3483500000011 -2291.750183827,1.3634834999986,5.4376712625892,0.69023005995503,0.67000000000138 -4490.5871620418,2.0334834999981,5.4376712625853,0.6902300599606,1.879670737393e-12 -2291.7501838092,2.703483499999,5.4376712625833,0.69023005996411,-0.66999999999897 -2234.7381816512,3.3818334999994,5.4376712625831,0.69023005996468,-1.3483499999994 -4346.1658455183,4.0534835,5.437671262582,0.69023005996671,-2.02 +4423.5634730749,0,5.4376712625889,0.69023005995495,2.01 +2245.2935810301,0.66999999999988,5.4376712625889,0.69023005995495,1.3400000000001 +2278.8054255322,1.3399999999989,5.4376712625886,0.69023005995556,0.67000000000112 +4490.587162043,2.0099999999982,5.4376712625852,0.69023005996058,1.7661076625498e-12 +2278.8054255201,2.6799999999992,5.4376712625838,0.69023005996362,-0.66999999999919 +2245.2935810234,3.35,5.4376712625837,0.69023005996365,-1.34 +4423.563473058,4.02,5.4376712625837,0.69023005996365,-2.01 diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e deleted file mode 100644 index 8056d299c9d0..000000000000 Binary files a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e and /dev/null differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e-s002 b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e-s002 index e3756b40d28a..f211c6bfc084 100644 Binary files a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e-s002 and b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_fatigue_out.e-s002 differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out.e-s003 b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out.e-s003 index 837585bbbb4e..ad4fb1892e6a 100644 Binary files a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out.e-s003 and b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out.e-s003 differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out_XFEMCutMeshOutput.e-s0004 b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out_XFEMCutMeshOutput.e-s0004 index 9fcaa0da1b3a..59bf352fd1d7 100644 Binary files a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out_XFEMCutMeshOutput.e-s0004 and b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_mhs_out_XFEMCutMeshOutput.e-s0004 differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_crit_out.json b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_crit_out.json index 0b9f7fe9c230..ce09bc1e6090 100644 --- a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_crit_out.json +++ b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_crit_out.json @@ -30,33 +30,33 @@ "crack_growth": [ 0.0, 0.1, - 0.09999843381232325, + 0.09999951666015743, 0.0 ], "id": [ 0.0, - 0.06766966700463402, - 0.13500336700775803, - 0.2026733370094032 + 0.067333, + 0.13466670000312403, + 0.20200000000312401 ], - "max_growth_timestep": 11.121705655963058, + "max_growth_timestep": 11.150025746130249, "x": [ - 0.5006265231632837, - 0.5006271726089853, - 0.5006278188303126, - 0.500628468278922 + 0.5006271726095989, + 0.5006271726085989, + 0.5006278188302707, + 0.5006278188292708 ], "y": [ - 0.5086195006503782, - 0.5086195565940692, - 0.5086196122600105, - 0.5086196682039521 + 0.508619556594082, + 0.508619556594082, + 0.508619612260494, + 0.508619612260494 ], "z": [ - 0.20133666700001482, + 0.2009999999985204, 0.13366699999852039, 0.06633329999852038, - -0.001336669999985204 + -0.0010000000014796257 ] }, "time": 2.0, diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_out_II_KI_1_0012.csv b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_out_II_KI_1_0012.csv index 034bc697547e..99be54ae0eaf 100644 --- a/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_out_II_KI_1_0012.csv +++ b/modules/xfem/test/tests/solid_mechanics_basic/gold/edge_crack_3d_scc_out_II_KI_1_0012.csv @@ -1,5 +1,5 @@ II_KI_1,id,x,y,z -24.107023438153,0,0.49354320067,0.50460660459144,0.20133666440821 -12.06471099363,0.067669405231478,0.49354359840248,0.50460659434588,0.1336672591779 -12.064767682571,0.13500310523157,0.49354399416184,0.50460658415114,0.066333559178974 -24.106823970657,0.20267333700349,0.49354439189918,0.50460657390546,-0.001336672591779 +24.350155518531,0,0.4933438029025,0.50459573950793,0.20100014533238 +12.035755324636,0.067333,0.4933438029015,0.50459573950793,0.13366714533238 +12.035797701633,0.13466669999972,0.49334422545983,0.50459574549825,0.066333445333988 +24.350152017561,0.20199999999972,0.49334422545883,0.50459574549825,-0.00099985466601235 diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_function_off_0_spin_0_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_function_off_0_spin_0_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..599acfe5632b Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_function_off_0_spin_0_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_hoopStress_off_0_spin_0_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_hoopStress_off_0_spin_0_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..1a9984c1cfae Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_hoopStress_off_0_spin_0_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_radialFunc_off_0_spin_4_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_radialFunc_off_0_spin_4_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..25cec5ef0304 Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_radialFunc_off_0_spin_4_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0.21_spin_0_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0.21_spin_0_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..5d5b74c15a41 Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0.21_spin_0_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_0_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_0_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..20db33093ba8 Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_0_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_15_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_15_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..ce33cd1a9e74 Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_15_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_25_XFEMCutMeshOutput.e b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_25_XFEMCutMeshOutput.e new file mode 100644 index 000000000000..de6f42b2a922 Binary files /dev/null and b/modules/xfem/test/tests/solid_mechanics_basic/gold/face_crack_scc_off_0_spin_25_XFEMCutMeshOutput.e differ diff --git a/modules/xfem/test/tests/solid_mechanics_basic/tests b/modules/xfem/test/tests/solid_mechanics_basic/tests index 3872ec2e6e9a..328a19002613 100644 --- a/modules/xfem/test/tests/solid_mechanics_basic/tests +++ b/modules/xfem/test/tests/solid_mechanics_basic/tests @@ -1,12 +1,16 @@ [Tests] - issues = '#6320 #23572 #32418' + issues = '#6320 #23572 #32418 #32779' + # XFEM exodiff tests usually require the following options + # newly created cut elements can lie on top of one another so geometrical matching will not work + # map = false + # XFEM requires --enable-unique-ids in libmesh + # capabilities = 'unique_id' [crack_propagation_ave] type = Exodiff input = crack_propagation_2d.i exodiff = 'crack_propagation_2d_out.e crack_propagation_2d_out.e-s002' abs_zero = 1e-8 map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' restep = false # Issue #31054 requirement = 'The XFEM module shall represent a propagating crack in a 2D mechanics problem ' @@ -23,7 +27,6 @@ exodiff = 'crack_propagation_2d_out.e crack_propagation_2d_out.e-s002' abs_zero = 1e-8 map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' requirement = 'The XFEM module shall represent a propagating crack in a 2D mechanics problem in ' 'which crack growth occurs when the average stress in the element at the crack tip ' @@ -38,7 +41,6 @@ exodiff = 'crack_propagation_2d_single_point_out.e crack_propagation_2d_single_point_out.e-s002' abs_zero = 1e-8 map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' restep = false # Issue #31054 requirement = 'The XFEM module shall represent a propagating crack in a 2D mechanics problem in ' @@ -52,7 +54,6 @@ exodiff = 'edge_crack_3d_out.e' abs_zero = 1e-8 map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' requirement = 'The XFEM module shall permit definition of a stationary crack in a 3D mechanics ' 'model with XFEM, where the crack is defined using a rectangular cutting plane by RectangleCutUserObject' @@ -66,10 +67,9 @@ [fatigue] type = Exodiff input = 'edge_crack_3d_domain.i edge_crack_3d_fatigue.i' - exodiff = 'edge_crack_3d_fatigue_out.e edge_crack_3d_fatigue_out.e-s002' + exodiff = 'edge_crack_3d_fatigue_out.e-s002' abs_zero = 1e-8 map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' growth speeds determined by the fatigue cracking Paris law, or' [] @@ -78,7 +78,6 @@ cli_args = 'Executioner/end_time=11' input = 'edge_crack_3d_domain.i edge_crack_3d_scc.i' csvdiff = 'edge_crack_3d_scc_out_II_KI_1_0012.csv' - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' growth speeds determined by the stress corrosion cracking constant timestep, or' [] @@ -87,7 +86,6 @@ cli_args = 'Executioner/end_time=2' input = 'edge_crack_3d_domain.i edge_crack_3d_scc_crit.i' jsondiff = 'edge_crack_3d_scc_crit_out.json' - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' growth speeds determined by the stress corrosion cracking maximum growth length, or' [] @@ -98,7 +96,6 @@ abs_zero = 1e-4 #crack face is normal to z so this gets rid of exodiffs caused by z-direction variables map = false partial=true - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' growth directions determined by computing the maximum hoop stress law using interaction integrals, or' [] @@ -108,7 +105,6 @@ exodiff = 'edge_crack_3d_propagation_out.e edge_crack_3d_propagation_out.e-s002' abs_zero = 1e-8 map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' growth directions determined by an parsed function, or' [] @@ -117,9 +113,71 @@ input = 'edge_crack_3d_domain.i edge_crack_3d_mhs.i' cli_args = 'DomainIntegral/number_points_from_provider=3' expect_err = 'This must match the number of points provided by the XFEM mesh cutter object.' - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' - detail = ' error when the number of points on the cutter mesh do not match the number of points in the DomainIntegralAction.' + detail = ' error when the number of points on the cutter mesh do not match the number of points in the DomainIntegralAction, or' + [] + [face_cut_mid_function] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_function.i' + cli_args = 'offset=0 spin=0' + exodiff = 'face_crack_function_off_0_spin_0_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' without producing degenerate elements and properly moving crack front nodes outside the FEM volume using a prescribed function growth direction, or' + [] + [face_cut_mid_hoopStress] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_hoopStress.i' + cli_args = 'offset=0 spin=0' + exodiff = 'face_crack_hoopStress_off_0_spin_0_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' properly moving crack front nodes using the maximum hoop stress growth direction, or' + [] + [face_cut_mid_scc] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_scc.i' + cli_args = 'offset=0 spin=0' + exodiff = 'face_crack_scc_off_0_spin_0_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' properly moving crack front nodes using SCC-driven non-uniform growth, or' + [] + [face_cut_corner_scc] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_scc.i' + cli_args = 'offset=0.21 spin=0' + exodiff = 'face_crack_scc_off_0.21_spin_0_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' properly moving crack front nodes for a corner crack using SCC-driven non-uniform growth, or' + [] + [inactive_refine] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_scc.i' + cli_args = 'offset=0 spin=25 Executioner/end_time=2' + exodiff = 'face_crack_scc_off_0_spin_25_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' properly refining surface cutting crack front if refinement occurs inside volume, or' + [] + [inactive_no_refine] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_scc.i' + cli_args = 'offset=0 spin=15 Executioner/end_time=2' + exodiff = 'face_crack_scc_off_0_spin_15_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' not refining surface cutting crack front if refinement occurs outside volume, or' + [] + [double_cut_edge] + type = Exodiff + input = 'face_crack_3d_domain.i face_crack_3d_xfem_radialFunc.i' + cli_args = 'offset=0 spin=4' + exodiff = 'face_crack_radialFunc_off_0_spin_4_XFEMCutMeshOutput.e' + map = true # only matching cutter mesh, not xfem volume mesh + capabilities = 'unique_id' + detail = ' growing the crack front and double cutting a volume element edge.' [] [] [penny_crack] @@ -132,7 +190,6 @@ input = penny_crack.i exodiff = 'penny_crack_out.e' map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' explicitly specified.' [] @@ -141,7 +198,6 @@ input = penny_crack_cfp.i exodiff = 'penny_crack_cfp_out.e' map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' detail = ' are provided by a class that derives from CrackFrontPointsProvider.' [] @@ -151,7 +207,6 @@ input = elliptical_crack.i exodiff = 'elliptical_crack_out.e' map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' requirement = 'The XFEM module shall permit definition of a stationary crack in a 3D mechanics ' 'model with XFEM, where the crack is defined using an elliptical cutting plane by EllipseCutUserObject' @@ -162,7 +217,6 @@ input = square_branch_quad_2d.i exodiff = 'square_branch_quad_2d_out.e square_branch_quad_2d_out.e-s002 square_branch_quad_2d_out.e-s003' map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' requirement = 'The XFEM system shall permit branched cracks to be represented in 2D by sequentially cutting a ' '4-noded quadrilateral element by two prescribed evolving cutting planes' @@ -173,7 +227,6 @@ input = square_branch_tri_2d.i exodiff = 'square_branch_tri_2d_out.e square_branch_tri_2d_out.e-s002 square_branch_tri_2d_out.e-s003' map = false - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' requirement = 'The XFEM system shall permit branched cracks to be represented in 2D by sequentially cutting a ' '3-noded triangle element by two prescribed evolving cutting planes' @@ -183,7 +236,6 @@ type = CSVDiff input = crack_counter.i csvdiff = 'crack_counter_out.csv' - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id' requirement = 'The XFEM system shall provide an accessor function to the crack_tip_origin_direction_map' allow_test_objects = true @@ -195,7 +247,6 @@ type = CSVDiff input = 'double_cantilever_crack.i' csvdiff = 'double_cantilever_crack_out_II_KI_1_0009.csv' - # XFEM requires --enable-unique-ids in libmesh capabilities = 'unique_id & method=opt' requirement = 'The XFEM module shall reproduce published results for a moving crack in a 3D double cantilever beam.' [] diff --git a/modules/xfem/unit/src/ElementFragmentAlgorithmTest.C b/modules/xfem/unit/src/ElementFragmentAlgorithmTest.C index 477860966856..4522fc3b3634 100644 --- a/modules/xfem/unit/src/ElementFragmentAlgorithmTest.C +++ b/modules/xfem/unit/src/ElementFragmentAlgorithmTest.C @@ -9,8 +9,14 @@ #include "gtest/gtest.h" +#include + #include "MooseUtils.h" +#include "MooseException.h" #include "ElementFragmentAlgorithm.h" +#include "EFAElement3D.h" +#include "EFAFace.h" +#include "EFAEdge.h" // set this global to true to enable a lot of mesh data output to the console const bool debug_print_mesh = false; @@ -3369,3 +3375,52 @@ TEST(ElementFragmentAlgorithm, test7a) std::vector en_gold = {0, 1, 2, 3, 4, 5}; CheckNodes(embedded_nodes, en_gold); } + +TEST(ElementFragmentAlgorithm, duplicateEmbeddedNodeReuseOnSharedEdgeCurrentFailure) +{ + ElementFragmentAlgorithm mesh(Moose::out); + + // Single HEX8 with the standard local coordinates: + // node 0 = (0,0,0), node 1 = (1,0,0), node 2 = (1,1,0), node 3 = (0,1,0) + // node 4 = (0,0,1), node 5 = (1,0,1), node 6 = (1,1,1), node 7 = (0,1,1) + // + // This unit test does not reproduce the full application propagation history across multiple + // structural elements. It reproduces the local condition at the exact throw site in + // EFAElement3D::addFaceEdgeCut(): + // two different faces of the same HEX8 reference the same physical edge point, and the second + // path arrives with a different embedded node pointer. + // + // In the application failure, face 5 and face 4 both touch the shared physical edge and a later + // propagation path reaches that same edge position with a different embedded node id. Here: + // 1. face 5 edge 3 seeds the shared physical edge and mirrors it onto face 4 edge 2 through + // add_to_adjacent=true + // 2. face 4 edge 2 is then revisited directly with a different embedded node + // This is narrower than the full application history, but it matches the same "different faces, + // same physical edge, different embedded node" failure condition. + std::vector elem = {0, 1, 2, 3, 4, 5, 6, 7}; + mesh.add3DElement(elem, 0); + + auto * base_elem = dynamic_cast(mesh.getElemByID(0)); + ASSERT_NE(base_elem, nullptr); + + std::map embedded_nodes; + EFANode embedded_a(1000, EFANode::N_CATEGORY_EMBEDDED); + EFANode embedded_b(1001, EFANode::N_CATEGORY_EMBEDDED); + auto * face5_edge3 = base_elem->getFace(5)->getEdge(3); + auto * face4_edge2 = base_elem->getFace(4)->getEdge(2); + + // Face 5 edge 3 and face 4 edge 2 are the same physical HEX8 edge with opposite face-local + // orientation. Seed that shared edge from face 5 first, allowing adjacent-face propagation to + // populate face 4 with the same embedded node. + base_elem->addFaceEdgeCut(5, 3, 0.5, &embedded_a, embedded_nodes, false, true); + ASSERT_EQ(face5_edge3->numEmbeddedNodes(), 1); + ASSERT_EQ(face4_edge2->numEmbeddedNodes(), 1); + + // Revisit the same physical edge through the other face with a different embedded node. + // Current code throws here. The intended behavior after the fix is to reuse the existing + // embedded node on that edge/position instead of erroring out, so this test is expected to + // fail until that logic is changed. + base_elem->addFaceEdgeCut(4, 2, 0.5, &embedded_b, embedded_nodes, false, true); + ASSERT_EQ(face5_edge3->numEmbeddedNodes(), 1); + ASSERT_EQ(face4_edge2->numEmbeddedNodes(), 1); +}