Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CrackFrontDefinition.h"
11 :
12 : // MOOSE includes
13 : #include "MooseError.h"
14 : #include "MooseMesh.h"
15 : #include "MooseVariable.h"
16 : #include "RankTwoTensor.h"
17 :
18 : #include "libmesh/mesh_tools.h"
19 : #include "libmesh/string_to_enum.h"
20 : #include "libmesh/quadrature.h"
21 :
22 : registerMooseObject("SolidMechanicsApp", CrackFrontDefinition);
23 :
24 : InputParameters
25 608 : CrackFrontDefinition::validParams()
26 : {
27 608 : InputParameters params = GeneralUserObject::validParams();
28 608 : params.addClassDescription("Used to describe geometric characteristics of the crack front for "
29 : "fracture integral calculations");
30 608 : params += BoundaryRestrictable::validParams();
31 608 : CrackFrontDefinition::includeCrackFrontDefinitionParams(params);
32 608 : params.set<bool>("use_displaced_mesh") = false;
33 :
34 1216 : params.addRelationshipManager("ElementSideNeighborLayers",
35 : Moose::RelationshipManagerType::ALGEBRAIC,
36 237 : [](const InputParameters &, InputParameters & rm_params)
37 237 : { rm_params.set<unsigned short>("layers") = 2; });
38 608 : return params;
39 0 : }
40 :
41 : void
42 1746 : CrackFrontDefinition::includeCrackFrontDefinitionParams(InputParameters & params)
43 : {
44 3492 : MooseEnum direction_method("CrackDirectionVector CrackMouth CurvedCrackFront");
45 : MooseEnum end_direction_method("NoSpecialTreatment CrackDirectionVector CrackTangentVector",
46 3492 : "NoSpecialTreatment");
47 3492 : params.addParam<std::vector<Point>>("crack_front_points", "Set of points to define crack front");
48 3492 : params.addParam<bool>("closed_loop", false, "Set of points forms forms a closed loop");
49 3492 : params.addRequiredParam<MooseEnum>(
50 : "crack_direction_method",
51 : direction_method,
52 1746 : "Method to determine direction of crack propagation. Choices are: " +
53 1746 : direction_method.getRawNames());
54 3492 : params.addParam<MooseEnum>(
55 : "crack_end_direction_method",
56 : end_direction_method,
57 1746 : "Method to determine direction of crack propagation at ends of crack. Choices are: " +
58 1746 : end_direction_method.getRawNames());
59 3492 : params.addParam<RealVectorValue>("crack_direction_vector", "Direction of crack propagation");
60 3492 : params.addParam<RealVectorValue>(
61 : "crack_direction_vector_end_1",
62 : "Direction of crack propagation for the node at end 1 of the crack");
63 3492 : params.addParam<RealVectorValue>(
64 : "crack_direction_vector_end_2",
65 : "Direction of crack propagation for the node at end 2 of the crack");
66 3492 : params.addParam<RealVectorValue>("crack_tangent_vector_end_1",
67 : "Direction of crack tangent for the node at end 1 of the crack");
68 3492 : params.addParam<RealVectorValue>("crack_tangent_vector_end_2",
69 : "Direction of crack tangent for the node at end 2 of the crack");
70 3492 : params.addParam<std::vector<BoundaryName>>(
71 : "crack_mouth_boundary", "Boundaries whose average coordinate defines the crack mouth");
72 3492 : params.addParam<std::vector<BoundaryName>>("intersecting_boundary",
73 : "Boundaries intersected by ends of crack");
74 3492 : params.addParam<bool>("2d", false, "Treat body as two-dimensional");
75 5238 : params.addRangeCheckedParam<unsigned int>(
76 : "axis_2d",
77 3492 : 2,
78 : "axis_2d>=0 & axis_2d<=2",
79 : "Out of plane axis for models treated as two-dimensional (0=x, 1=y, 2=z)");
80 3492 : params.addParam<unsigned int>("symmetry_plane",
81 : "Account for a symmetry plane passing through "
82 : "the plane of the crack, normal to the specified "
83 : "axis (0=x, 1=y, 2=z)");
84 3492 : params.addParam<bool>("t_stress", false, "Calculate T-stress");
85 3492 : params.addParam<bool>("q_function_rings", false, "Generate rings of nodes for q-function");
86 3492 : params.addParam<unsigned int>("last_ring", "The number of rings of nodes to generate");
87 3492 : params.addParam<unsigned int>("first_ring", "The number of rings of nodes to generate");
88 3492 : params.addParam<unsigned int>("nrings", "The number of rings of nodes to generate");
89 3492 : params.addParam<VariableName>("disp_x", "Variable containing the x displacement");
90 3492 : params.addParam<VariableName>("disp_y", "Variable containing the y displacement");
91 3492 : params.addParam<VariableName>("disp_z", "Variable containing the z displacement");
92 3492 : params.addParam<std::vector<Real>>(
93 : "j_integral_radius_inner", {}, "Radius for J-Integral calculation");
94 3492 : params.addParam<std::vector<Real>>(
95 : "j_integral_radius_outer", {}, "Radius for J-Integral calculation");
96 3492 : MooseEnum q_function_type("Geometry Topology", "Geometry");
97 3492 : params.addParam<MooseEnum>("q_function_type",
98 : q_function_type,
99 1746 : "The method used to define the integration domain. Options are: " +
100 1746 : q_function_type.getRawNames());
101 3492 : params.addParam<UserObjectName>(
102 : "crack_front_points_provider",
103 : "The UserObject provides the crack front points from XFEM GeometricCutObject");
104 :
105 3492 : params.addParam<unsigned int>("number_points_from_provider",
106 : "The number of crack front points, only needed for "
107 : "crack_front_points_provider that do not use a cut mesh.");
108 3492 : params.addParam<Real>(
109 : "crack_front_node_tolerance",
110 3492 : 1e-10,
111 : "General tolerance for locating nodes on the crack front based on xyz coordinates.");
112 1746 : }
113 :
114 196 : CrackFrontDefinition::CrackFrontDefinition(const InputParameters & parameters)
115 : : GeneralUserObject(parameters),
116 : BoundaryRestrictable(this, true), // false means nodesets
117 196 : _direction_method(getParam<MooseEnum>("crack_direction_method").getEnum<DIRECTION_METHOD>()),
118 196 : _end_direction_method(
119 196 : getParam<MooseEnum>("crack_end_direction_method").getEnum<END_DIRECTION_METHOD>()),
120 196 : _aux(_fe_problem.getAuxiliarySystem()),
121 196 : _mesh(_subproblem.mesh()),
122 392 : _treat_as_2d(getParam<bool>("2d")),
123 196 : _use_mesh_cutter(false),
124 196 : _is_cutter_modified(false),
125 392 : _closed_loop(getParam<bool>("closed_loop")),
126 392 : _axis_2d(getParam<unsigned int>("axis_2d")),
127 392 : _has_symmetry_plane(isParamValid("symmetry_plane")),
128 466 : _symmetry_plane(_has_symmetry_plane ? getParam<unsigned int>("symmetry_plane")
129 : : std::numeric_limits<unsigned int>::max()),
130 392 : _t_stress(getParam<bool>("t_stress")),
131 392 : _q_function_rings(getParam<bool>("q_function_rings")),
132 588 : _q_function_type(getParam<MooseEnum>("q_function_type")),
133 196 : _crack_front_points_provider(nullptr),
134 588 : _tol(getParam<Real>("crack_front_node_tolerance"))
135 : {
136 392 : auto boundary = isParamValid("boundary") ? getParam<std::vector<BoundaryName>>("boundary")
137 556 : : std::vector<BoundaryName>{};
138 392 : if (isParamValid("crack_front_points"))
139 : {
140 16 : if (boundary.size())
141 0 : paramError("crack_front_points",
142 : "CrackFrontDefinition error: since boundary is defined, crack_front_points should "
143 : "not be added.");
144 32 : if (isParamValid("crack_front_points_provider"))
145 0 : paramError("crack_front_points_provider",
146 : "As crack_front_points have been provided, the crack_front_points_provider will "
147 : "not be used and needs to be removed.");
148 32 : _crack_front_points = getParam<std::vector<Point>>("crack_front_points");
149 16 : _geom_definition_method = CRACK_GEOM_DEFINITION::CRACK_FRONT_POINTS;
150 16 : if (_t_stress)
151 0 : paramError("t_stress", "t_stress not yet supported with crack_front_points");
152 16 : if (_q_function_rings)
153 0 : paramError("q_function_rings", "q_function_rings not supported with crack_front_points");
154 : }
155 360 : else if (isParamValid("crack_front_points_provider"))
156 : {
157 0 : if (boundary.size())
158 0 : paramError("crack_front_points_provider",
159 : "CrackFrontDefinition error: since boundary is defined, "
160 : "crack_front_points_provider should not be added.");
161 :
162 0 : _geom_definition_method = CRACK_GEOM_DEFINITION::CRACK_FRONT_POINTS;
163 : }
164 360 : else if (isParamValid("number_points_from_provider"))
165 0 : paramError("number_points_from_provider",
166 : "CrackFrontDefinition error: number_points_from_provider is provided but "
167 : "crack_front_points_provider cannot be found.");
168 180 : else if (boundary.size())
169 : {
170 180 : _geom_definition_method = CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES;
171 360 : if (parameters.isParamSetByUser("closed_loop"))
172 0 : paramError("closed_loop",
173 : "In CrackFrontDefinition, if 'boundary' is defined, 'closed_loop' should not be "
174 : "set by user because closed loops are detected automatically");
175 : }
176 : else
177 0 : mooseError("In CrackFrontDefinition, must define one of 'boundary', 'crack_front_points' "
178 : "and 'crack_front_points_provider'");
179 :
180 392 : if (isParamValid("crack_mouth_boundary"))
181 54 : _crack_mouth_boundary_names = getParam<std::vector<BoundaryName>>("crack_mouth_boundary");
182 :
183 196 : if (_has_symmetry_plane)
184 135 : if (_symmetry_plane > 2)
185 0 : paramError("symmetry_plane",
186 : "symmetry_plane out of bounds: ",
187 : _symmetry_plane,
188 : " Must be >=0 and <=2.");
189 :
190 196 : switch (_direction_method)
191 : {
192 155 : case DIRECTION_METHOD::CRACK_DIRECTION_VECTOR:
193 310 : if (!isParamValid("crack_direction_vector"))
194 0 : paramError("crack_direction_vector",
195 : "crack_direction_vector must be specified if crack_direction_method = "
196 : "CrackDirectionVector");
197 310 : _crack_direction_vector = getParam<RealVectorValue>("crack_direction_vector");
198 155 : break;
199 18 : case DIRECTION_METHOD::CRACK_MOUTH:
200 36 : if (isParamValid("crack_direction_vector"))
201 0 : paramError("crack_direction_vector",
202 : "crack_direction_vector must not be specified if crack_direction_method = "
203 : "CrackMouthNodes");
204 18 : if (_crack_mouth_boundary_names.size() == 0)
205 0 : paramError(
206 : "crack_mouth_boundary",
207 : "crack_mouth_boundary must be specified if crack_direction_method = CrackMouthNodes");
208 : break;
209 23 : case DIRECTION_METHOD::CURVED_CRACK_FRONT:
210 46 : if (isParamValid("crack_direction_vector"))
211 0 : paramError("crack_direction_vector",
212 : "crack_direction_vector must not be specified if crack_direction_method = "
213 : "CurvedCrackFront");
214 : break;
215 0 : default:
216 0 : paramError("crack_direction_method", "Invalid direction_method");
217 : }
218 :
219 392 : if (isParamValid("intersecting_boundary"))
220 66 : _intersecting_boundary_names = getParam<std::vector<BoundaryName>>("intersecting_boundary");
221 :
222 196 : if (_end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR)
223 : {
224 52 : if (!isParamValid("crack_direction_vector_end_1"))
225 0 : paramError("crack_direction_vector_end_1",
226 : "crack_direction_vector_end_1 must be specified if crack_end_direction_method = "
227 : "CrackDirectionVector");
228 52 : if (!isParamValid("crack_direction_vector_end_2"))
229 0 : paramError("crack_direction_vector_end_2",
230 : "crack_direction_vector_end_2 must be specified if crack_end_direction_method = "
231 : "CrackDirectionVector");
232 52 : _crack_direction_vector_end_1 = getParam<RealVectorValue>("crack_direction_vector_end_1");
233 52 : _crack_direction_vector_end_2 = getParam<RealVectorValue>("crack_direction_vector_end_2");
234 : }
235 :
236 196 : if (_end_direction_method == END_DIRECTION_METHOD::END_CRACK_TANGENT_VECTOR)
237 : {
238 0 : if (!isParamValid("crack_tangent_vector_end_1"))
239 0 : paramError("crack_tangent_vector_end_1",
240 : "crack_tangent_vector_end_1 must be specified if crack_end_tangent_method = "
241 : "CrackTangentVector");
242 0 : if (!isParamValid("crack_tangent_vector_end_2"))
243 0 : paramError("crack_tangent_vector_end_2",
244 : "crack_tangent_vector_end_2 must be specified if crack_end_tangent_method = "
245 : "CrackTangentVector");
246 0 : _crack_tangent_vector_end_1 = getParam<RealVectorValue>("crack_tangent_vector_end_1");
247 0 : _crack_tangent_vector_end_2 = getParam<RealVectorValue>("crack_tangent_vector_end_2");
248 : }
249 :
250 440 : if (isParamValid("disp_x") && isParamValid("disp_y") && isParamValid("disp_z"))
251 : {
252 16 : _disp_x_var_name = getParam<VariableName>("disp_x");
253 16 : _disp_y_var_name = getParam<VariableName>("disp_y");
254 24 : _disp_z_var_name = getParam<VariableName>("disp_z");
255 : }
256 188 : else if (_t_stress == true && _treat_as_2d == false)
257 0 : paramError("displacements", "Displacement variables must be provided for T-stress calculation");
258 :
259 196 : if (_q_function_rings)
260 : {
261 40 : if (!isParamValid("last_ring"))
262 0 : paramError("last_ring",
263 : "The max number of rings of nodes to generate must be provided if "
264 : "q_function_rings = true");
265 40 : _last_ring = getParam<unsigned int>("last_ring");
266 40 : _first_ring = getParam<unsigned int>("first_ring");
267 : }
268 : else
269 : {
270 352 : _j_integral_radius_inner = getParam<std::vector<Real>>("j_integral_radius_inner");
271 528 : _j_integral_radius_outer = getParam<std::vector<Real>>("j_integral_radius_outer");
272 : }
273 196 : }
274 :
275 : void
276 439 : CrackFrontDefinition::execute()
277 : {
278 : // Because J-Integral is based on original geometry, the crack front geometry
279 : // is never updated, so everything that needs to happen is done in initialSetup()
280 439 : if (_t_stress == true && _treat_as_2d == false)
281 18 : calculateTangentialStrainAlongFront();
282 439 : }
283 :
284 : void
285 195 : CrackFrontDefinition::initialSetup()
286 : {
287 390 : if (isParamValid("crack_front_points_provider"))
288 : {
289 0 : _crack_front_points_provider = &getUserObjectByName<CrackFrontPointsProvider>(
290 : getParam<UserObjectName>("crack_front_points_provider"));
291 0 : if (_crack_front_points_provider->usesMesh())
292 : {
293 0 : _use_mesh_cutter = true;
294 :
295 : // Automatically get number of crack front points from mesh-based provider
296 0 : _num_points_from_provider = _crack_front_points_provider->getNumberOfCrackFrontPoints();
297 0 : mooseInfo("CrackFrontDefinition: Automatically detected ",
298 0 : _num_points_from_provider,
299 : " crack front points from mesh-based provider");
300 :
301 0 : if (_direction_method != DIRECTION_METHOD::CURVED_CRACK_FRONT)
302 0 : paramError("crack_direction_method",
303 : "Using a `crack_front_points_provider` that uses an XFEM cutter mesh also "
304 : "requires setting 'crack_direction_method = CurvedCrackFront'");
305 0 : if (isParamValid("crack_mouth_boundary"))
306 0 : paramError("crack_mouth_boundary",
307 : "'crack_mouth_boundary' cannot be set when using a "
308 : "'crack_front_points_provider' that uses an XFEM cutter mesh");
309 : }
310 0 : else if (isParamValid("number_points_from_provider"))
311 : {
312 0 : _num_points_from_provider = getParam<unsigned int>("number_points_from_provider");
313 : }
314 0 : else if (_num_points_from_provider == 0)
315 : {
316 0 : paramError(
317 : "number_points_from_provider",
318 : "CrackFrontDefinition error: When using a non-mesh-based crack_front_points_provider, "
319 : "the number_points_from_provider parameter must be provided.");
320 : }
321 : }
322 195 : if (_crack_front_points_provider != nullptr)
323 : {
324 : _crack_front_points =
325 0 : _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider);
326 : }
327 :
328 195 : _crack_mouth_boundary_ids = _mesh.getBoundaryIDs(_crack_mouth_boundary_names, true);
329 195 : _intersecting_boundary_ids = _mesh.getBoundaryIDs(_intersecting_boundary_names, true);
330 :
331 195 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
332 : {
333 : std::set<dof_id_type> nodes;
334 179 : getCrackFrontNodes(nodes);
335 179 : orderCrackFrontNodes(nodes);
336 : }
337 :
338 195 : if (_closed_loop && _intersecting_boundary_names.size() > 0)
339 0 : paramError("intersecting_boundary", "Cannot use intersecting_boundary with closed-loop cracks");
340 :
341 195 : updateCrackFrontGeometry();
342 :
343 195 : if (_q_function_rings)
344 20 : createQFunctionRings();
345 :
346 195 : if (_t_stress)
347 : {
348 : std::size_t num_crack_front_nodes = _ordered_crack_front_nodes.size();
349 72 : for (std::size_t i = 0; i < num_crack_front_nodes; ++i)
350 60 : _strain_along_front.push_back(-std::numeric_limits<Real>::max());
351 : }
352 :
353 195 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
354 195 : if (_q_function_type == "GEOMETRY")
355 : {
356 175 : if (!_treat_as_2d)
357 67 : if (num_crack_front_points < 1)
358 0 : mooseError("num_crack_front_points is not > 0");
359 753 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
360 : {
361 578 : bool is_point_on_intersecting_boundary = isPointWithIndexOnIntersectingBoundary(i);
362 578 : _is_point_on_intersecting_boundary.push_back(is_point_on_intersecting_boundary);
363 : }
364 : }
365 195 : }
366 :
367 : void
368 439 : CrackFrontDefinition::initialize()
369 : {
370 : // Update the crack front for fracture integral calculations
371 : // This is only useful for growing cracks which are currently described by the mesh
372 : // cutter
373 439 : if (_use_mesh_cutter && _is_cutter_modified)
374 : {
375 : // Automatically update the number of crack front points as crack grows
376 0 : _num_points_from_provider = _crack_front_points_provider->getNumberOfCrackFrontPoints();
377 :
378 : _crack_front_points =
379 0 : _crack_front_points_provider->getCrackFrontPoints(_num_points_from_provider);
380 0 : updateCrackFrontGeometry();
381 0 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
382 0 : if (_q_function_type == "GEOMETRY")
383 0 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
384 : {
385 0 : bool is_point_on_intersecting_boundary = isPointWithIndexOnIntersectingBoundary(i);
386 0 : _is_point_on_intersecting_boundary.push_back(is_point_on_intersecting_boundary);
387 : }
388 : }
389 439 : }
390 :
391 : void
392 439 : CrackFrontDefinition::finalize()
393 : {
394 439 : if (_t_stress)
395 27 : _communicator.max(_strain_along_front);
396 439 : }
397 :
398 : void
399 179 : CrackFrontDefinition::getCrackFrontNodes(std::set<dof_id_type> & nodes)
400 : {
401 179 : ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
402 18301 : for (auto nd = bnd_nodes.begin(); nd != bnd_nodes.end(); ++nd)
403 : {
404 18122 : const BndNode * bnode = *nd;
405 18122 : BoundaryID boundary_id = bnode->_bnd_id;
406 :
407 18122 : if (hasBoundary(boundary_id))
408 562 : nodes.insert(bnode->_node->id());
409 : }
410 :
411 179 : if (_treat_as_2d && _use_mesh_cutter == false)
412 : {
413 116 : if (nodes.size() > 1)
414 : {
415 : // Check that the nodes are collinear in the axis normal to the 2d plane
416 : unsigned int axis0;
417 : unsigned int axis1;
418 :
419 12 : switch (_axis_2d)
420 : {
421 : case 0:
422 : axis0 = 1;
423 : axis1 = 2;
424 : break;
425 : case 1:
426 : axis0 = 0;
427 : axis1 = 2;
428 : break;
429 : case 2:
430 : axis0 = 0;
431 : axis1 = 1;
432 : break;
433 0 : default:
434 0 : mooseError("Invalid axis.");
435 : }
436 :
437 : Real node0coor0 = 0;
438 : Real node0coor1 = 0;
439 :
440 36 : for (auto sit = nodes.begin(); sit != nodes.end(); ++sit)
441 : {
442 24 : Node & curr_node = _mesh.nodeRef(*sit);
443 24 : if (sit == nodes.begin())
444 : {
445 12 : node0coor0 = curr_node(axis0);
446 12 : node0coor1 = curr_node(axis1);
447 : }
448 : else
449 : {
450 12 : if (!MooseUtils::absoluteFuzzyEqual(curr_node(axis0), node0coor0, _tol) ||
451 : !MooseUtils::absoluteFuzzyEqual(curr_node(axis1), node0coor1, _tol))
452 0 : mooseError("Boundary provided in CrackFrontDefinition contains ",
453 0 : nodes.size(),
454 : " nodes, which are not collinear in the ",
455 0 : _axis_2d,
456 : " axis. Must contain either 1 node or collinear nodes to treat as 2D.");
457 : }
458 : }
459 : }
460 : }
461 179 : }
462 :
463 : void
464 179 : CrackFrontDefinition::orderCrackFrontNodes(std::set<dof_id_type> & nodes)
465 : {
466 179 : _ordered_crack_front_nodes.clear();
467 179 : if (nodes.size() < 1)
468 0 : mooseError("No crack front nodes");
469 179 : else if (nodes.size() == 1)
470 : {
471 104 : _ordered_crack_front_nodes.push_back(*nodes.begin());
472 104 : if (!_treat_as_2d)
473 0 : mooseError("Boundary provided in CrackFrontDefinition contains 1 node, but model is not "
474 : "treated as 2d");
475 : }
476 75 : else if (_treat_as_2d && _use_mesh_cutter)
477 : {
478 : // This is for the 2D case that uses a mesh cutter object so every node is its own crack front
479 : // and is not connected to other crack front nodes. Copying the order here makes it the same
480 : // as that given in the MeshCut2DFractureUserObject
481 0 : std::copy(nodes.begin(), nodes.end(), _ordered_crack_front_nodes.begin());
482 : }
483 : else
484 : {
485 : // Loop through the set of crack front nodes, and create a node to element map for just the
486 : // crack front nodes
487 : // The main reason for creating a second map is that we need to do a sort prior to the
488 : // set_intersection.
489 : // The original map contains vectors, and we can't sort them, so we create sets in the local
490 : // map.
491 75 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
492 : std::map<dof_id_type, std::set<dof_id_type>> crack_front_node_to_elem_map;
493 :
494 533 : for (const auto & node_id : nodes)
495 : {
496 : const auto & node_to_elem_pair = node_to_elem_map.find(node_id);
497 : mooseAssert(node_to_elem_pair != node_to_elem_map.end(),
498 : "Could not find crack front node " << node_id << " in the node to elem map");
499 :
500 : const std::vector<dof_id_type> & connected_elems = node_to_elem_pair->second;
501 2058 : for (std::size_t i = 0; i < connected_elems.size(); ++i)
502 1600 : crack_front_node_to_elem_map[node_id].insert(connected_elems[i]);
503 : }
504 :
505 : // Determine which nodes are connected to each other via elements, and construct line elements
506 : // to represent
507 : // those connections
508 : std::vector<std::vector<dof_id_type>> line_elems;
509 : std::map<dof_id_type, std::vector<dof_id_type>> node_to_line_elem_map;
510 :
511 75 : for (auto cfnemit = crack_front_node_to_elem_map.begin();
512 533 : cfnemit != crack_front_node_to_elem_map.end();
513 : ++cfnemit)
514 : {
515 : auto cfnemit2 = cfnemit;
516 3935 : for (++cfnemit2; cfnemit2 != crack_front_node_to_elem_map.end(); ++cfnemit2)
517 : {
518 :
519 : std::vector<dof_id_type> common_elements;
520 : std::set<dof_id_type> & elements_connected_to_node1 = cfnemit->second;
521 : std::set<dof_id_type> & elements_connected_to_node2 = cfnemit2->second;
522 3477 : std::set_intersection(elements_connected_to_node1.begin(),
523 : elements_connected_to_node1.end(),
524 : elements_connected_to_node2.begin(),
525 : elements_connected_to_node2.end(),
526 : std::inserter(common_elements, common_elements.end()));
527 :
528 3477 : if (common_elements.size() > 0)
529 : {
530 : std::vector<dof_id_type> my_line_elem;
531 390 : my_line_elem.push_back(cfnemit->first);
532 390 : my_line_elem.push_back(cfnemit2->first);
533 390 : node_to_line_elem_map[cfnemit->first].push_back(line_elems.size());
534 390 : node_to_line_elem_map[cfnemit2->first].push_back(line_elems.size());
535 390 : line_elems.push_back(my_line_elem);
536 390 : }
537 3477 : }
538 : }
539 :
540 : // Find nodes on ends of line (those connected to only one line element)
541 : std::vector<dof_id_type> end_nodes;
542 533 : for (auto nlemit = node_to_line_elem_map.begin(); nlemit != node_to_line_elem_map.end();
543 : ++nlemit)
544 : {
545 : std::size_t num_connected_elems = nlemit->second.size();
546 458 : if (num_connected_elems == 1)
547 136 : end_nodes.push_back(nlemit->first);
548 322 : else if (num_connected_elems != 2)
549 0 : mooseError(
550 0 : "Node ", nlemit->first, " is connected to >2 line segments in CrackFrontDefinition");
551 : }
552 :
553 : // For embedded crack with closed loop of crack front nodes, must pick the end nodes
554 75 : if (end_nodes.size() == 0) // Crack front is a loop. Pick nodes to be end nodes.
555 : {
556 7 : pickLoopCrackEndNodes(end_nodes, nodes, node_to_line_elem_map, line_elems);
557 7 : _closed_loop = true;
558 7 : if (_end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR ||
559 : _end_direction_method == END_DIRECTION_METHOD::END_CRACK_TANGENT_VECTOR)
560 0 : paramError("end_direction_method",
561 : "In CrackFrontDefinition, end_direction_method cannot be CrackDirectionVector "
562 : "or CrackTangentVector for a closed-loop crack");
563 7 : if (_intersecting_boundary_names.size() > 0)
564 0 : paramError("intersecting_boundary",
565 : "In CrackFrontDefinition, intersecting_boundary cannot be specified for a "
566 : "closed-loop crack");
567 : }
568 68 : else if (end_nodes.size() == 2) // Rearrange the order of the end nodes if needed
569 68 : orderEndNodes(end_nodes);
570 : else
571 0 : mooseError("In CrackFrontDefinition wrong number of end nodes. Number end nodes = ",
572 0 : end_nodes.size());
573 :
574 : // Create an ordered list of the nodes going along the line of the crack front
575 75 : _ordered_crack_front_nodes.push_back(end_nodes[0]);
576 :
577 75 : dof_id_type last_node = end_nodes[0];
578 : dof_id_type second_last_node = last_node;
579 458 : while (last_node != end_nodes[1])
580 : {
581 383 : std::vector<dof_id_type> & curr_node_line_elems = node_to_line_elem_map[last_node];
582 : bool found_new_node = false;
583 411 : for (std::size_t i = 0; i < curr_node_line_elems.size(); ++i)
584 : {
585 411 : std::vector<dof_id_type> curr_line_elem = line_elems[curr_node_line_elems[i]];
586 488 : for (std::size_t j = 0; j < curr_line_elem.size(); ++j)
587 : {
588 460 : dof_id_type line_elem_node = curr_line_elem[j];
589 460 : if (_closed_loop &&
590 238 : (last_node == end_nodes[0] &&
591 28 : line_elem_node == end_nodes[1])) // wrong direction around closed loop
592 7 : continue;
593 453 : if (line_elem_node != last_node && line_elem_node != second_last_node)
594 : {
595 383 : _ordered_crack_front_nodes.push_back(line_elem_node);
596 : found_new_node = true;
597 383 : break;
598 : }
599 : }
600 : if (found_new_node)
601 : break;
602 411 : }
603 383 : second_last_node = last_node;
604 383 : last_node = _ordered_crack_front_nodes.back();
605 : }
606 150 : }
607 179 : }
608 :
609 : void
610 68 : CrackFrontDefinition::orderEndNodes(std::vector<dof_id_type> & end_nodes)
611 : {
612 : // Choose the node to be the first node. Do that based on undeformed coordinates for
613 : // repeatability.
614 68 : Node & node0 = _mesh.nodeRef(end_nodes[0]);
615 68 : Node & node1 = _mesh.nodeRef(end_nodes[1]);
616 :
617 : std::size_t num_positive_coor0 = 0;
618 : std::size_t num_positive_coor1 = 0;
619 : Real dist_from_origin0 = 0.0;
620 : Real dist_from_origin1 = 0.0;
621 272 : for (std::size_t i = 0; i < 3; ++i)
622 : {
623 204 : dist_from_origin0 += node0(i) * node0(i);
624 204 : dist_from_origin1 += node1(i) * node1(i);
625 204 : if (MooseUtils::absoluteFuzzyGreaterThan(node0(i), 0.0, _tol))
626 18 : ++num_positive_coor0;
627 204 : if (MooseUtils::absoluteFuzzyGreaterThan(node1(i), 0.0, _tol))
628 68 : ++num_positive_coor1;
629 : }
630 68 : dist_from_origin0 = std::sqrt(dist_from_origin0);
631 68 : dist_from_origin1 = std::sqrt(dist_from_origin1);
632 :
633 : bool switch_ends = false;
634 68 : if (num_positive_coor1 > num_positive_coor0)
635 : {
636 : switch_ends = true;
637 : }
638 : else
639 : {
640 18 : if (!MooseUtils::absoluteFuzzyEqual(dist_from_origin1, dist_from_origin0, _tol))
641 : {
642 18 : if (dist_from_origin1 < dist_from_origin0)
643 : switch_ends = true;
644 : }
645 : else
646 : {
647 0 : if (end_nodes[1] < end_nodes[0])
648 : switch_ends = true;
649 : }
650 : }
651 : if (switch_ends)
652 : {
653 68 : std::size_t tmp_node = end_nodes[1];
654 68 : end_nodes[1] = end_nodes[0];
655 68 : end_nodes[0] = tmp_node;
656 : }
657 68 : }
658 :
659 : void
660 7 : CrackFrontDefinition::pickLoopCrackEndNodes(
661 : std::vector<dof_id_type> & end_nodes,
662 : std::set<dof_id_type> & nodes,
663 : std::map<dof_id_type, std::vector<dof_id_type>> & node_to_line_elem_map,
664 : std::vector<std::vector<dof_id_type>> & line_elems)
665 : {
666 : dof_id_type max_dist_node = 0;
667 : Real min_dist = std::numeric_limits<Real>::max();
668 : Real max_dist = -std::numeric_limits<Real>::max();
669 : // Pick the node farthest from the origin as the end node, or the one with
670 : // the greatest x coordinate if the nodes are equidistant from the origin
671 217 : for (auto nit = nodes.begin(); nit != nodes.end(); ++nit)
672 : {
673 210 : Node & node = _mesh.nodeRef(*nit);
674 210 : Real dist = node.norm();
675 210 : if (dist > max_dist)
676 : {
677 : max_dist = dist;
678 14 : max_dist_node = *nit;
679 : }
680 196 : else if (dist < min_dist)
681 : min_dist = dist;
682 : }
683 :
684 : dof_id_type end_node;
685 7 : if (MooseUtils::absoluteFuzzyGreaterThan(max_dist, min_dist, _tol))
686 0 : end_node = max_dist_node;
687 : else
688 : {
689 : std::vector<Node *> node_vec;
690 217 : for (auto nit = nodes.begin(); nit != nodes.end(); ++nit)
691 210 : node_vec.push_back(_mesh.nodePtr(*nit));
692 7 : end_node = maxNodeCoor(node_vec);
693 7 : }
694 :
695 7 : end_nodes.push_back(end_node);
696 :
697 : // Find the two nodes connected to the node identified as the end node, and pick one of those to
698 : // be the other end node
699 7 : auto end_node_line_elems = node_to_line_elem_map[end_node];
700 7 : if (end_node_line_elems.size() != 2)
701 0 : mooseError(
702 : "Crack front nodes are in a loop, but crack end node is only connected to one other node");
703 : std::vector<Node *> candidate_other_end_nodes;
704 :
705 21 : for (std::size_t i = 0; i < 2; ++i)
706 : {
707 14 : auto end_line_elem = line_elems[end_node_line_elems[i]];
708 42 : for (std::size_t j = 0; j < end_line_elem.size(); ++j)
709 : {
710 28 : auto line_elem_node = end_line_elem[j];
711 28 : if (line_elem_node != end_node)
712 14 : candidate_other_end_nodes.push_back(_mesh.nodePtr(line_elem_node));
713 : }
714 14 : }
715 7 : if (candidate_other_end_nodes.size() != 2)
716 0 : mooseError(
717 : "Crack front nodes are in a loop, but crack end node is not connected to two other nodes");
718 7 : end_nodes.push_back(maxNodeCoor(candidate_other_end_nodes, 1));
719 7 : }
720 :
721 : dof_id_type
722 14 : CrackFrontDefinition::maxNodeCoor(std::vector<Node *> & nodes, unsigned int dir0)
723 : {
724 : Real dirs[3];
725 : if (dir0 == 0)
726 : {
727 : dirs[0] = 0;
728 : dirs[1] = 1;
729 : dirs[2] = 2;
730 : }
731 : else if (dir0 == 1)
732 : {
733 : dirs[0] = 1;
734 : dirs[1] = 2;
735 : dirs[2] = 0;
736 : }
737 : else if (dir0 == 2)
738 : {
739 : dirs[0] = 2;
740 : dirs[1] = 0;
741 : dirs[2] = 1;
742 : }
743 : else
744 0 : mooseError("Invalid dir0 in CrackFrontDefinition::maxNodeCoor()");
745 :
746 : Real max_coor0 = -std::numeric_limits<Real>::max();
747 : std::vector<Node *> max_coor0_nodes;
748 238 : for (std::size_t i = 0; i < nodes.size(); ++i)
749 : {
750 224 : Real coor0 = (*nodes[i])(dirs[0]);
751 224 : if (coor0 > max_coor0)
752 : max_coor0 = coor0;
753 : }
754 238 : for (std::size_t i = 0; i < nodes.size(); ++i)
755 : {
756 224 : Real coor0 = (*nodes[i])(dirs[0]);
757 224 : if (MooseUtils::absoluteFuzzyEqual(coor0, max_coor0, _tol))
758 217 : max_coor0_nodes.push_back(nodes[i]);
759 : }
760 14 : if (max_coor0_nodes.size() > 1)
761 : {
762 : Real max_coor1 = -std::numeric_limits<Real>::max();
763 : std::vector<Node *> max_coor1_nodes;
764 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
765 : {
766 210 : Real coor1 = (*nodes[i])(dirs[1]);
767 210 : if (coor1 > max_coor1)
768 : max_coor1 = coor1;
769 : }
770 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
771 : {
772 210 : Real coor1 = (*nodes[i])(dirs[1]);
773 210 : if (MooseUtils::absoluteFuzzyEqual(coor1, max_coor1, _tol))
774 14 : max_coor1_nodes.push_back(nodes[i]);
775 : }
776 7 : if (max_coor1_nodes.size() > 1)
777 : {
778 : Real max_coor2 = -std::numeric_limits<Real>::max();
779 : std::vector<Node *> max_coor2_nodes;
780 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
781 : {
782 210 : Real coor2 = (*nodes[i])(dirs[2]);
783 210 : if (coor2 > max_coor2)
784 : max_coor2 = coor2;
785 : }
786 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
787 : {
788 210 : Real coor2 = (*nodes[i])(dirs[2]);
789 210 : if (MooseUtils::absoluteFuzzyEqual(coor2, max_coor2, _tol))
790 7 : max_coor2_nodes.push_back(nodes[i]);
791 : }
792 7 : if (max_coor2_nodes.size() > 1)
793 0 : mooseError("Multiple nodes with same x,y,z coordinates within tolerance");
794 : else
795 7 : return max_coor2_nodes[0]->id();
796 7 : }
797 : else
798 0 : return max_coor1_nodes[0]->id();
799 7 : }
800 : else
801 7 : return max_coor0_nodes[0]->id();
802 14 : }
803 :
804 : void
805 195 : CrackFrontDefinition::updateCrackFrontGeometry()
806 : {
807 195 : computeCrackMouthNodes();
808 195 : _segment_lengths.clear();
809 195 : _tangent_directions.clear();
810 195 : _crack_directions.clear();
811 195 : _overall_length = 0.0;
812 195 : _rot_matrix.clear();
813 195 : _distances_along_front.clear();
814 195 : _angles_along_front.clear();
815 195 : _strain_along_front.clear();
816 195 : _crack_plane_normals.clear();
817 :
818 195 : if (_treat_as_2d)
819 : {
820 120 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
821 120 : _segment_lengths.reserve(num_crack_front_points);
822 120 : _tangent_directions.reserve(num_crack_front_points);
823 120 : _crack_directions.reserve(num_crack_front_points);
824 120 : if (_use_mesh_cutter)
825 : _crack_plane_normals =
826 0 : _crack_front_points_provider->getCrackPlaneNormals(num_crack_front_points);
827 :
828 252 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
829 : {
830 : RealVectorValue tangent_direction;
831 : RealVectorValue crack_direction;
832 132 : tangent_direction(_axis_2d) = 1.0;
833 132 : _tangent_directions.push_back(tangent_direction);
834 132 : const Point * crack_front_point = getCrackFrontPoint(i);
835 :
836 : crack_direction =
837 132 : calculateCrackFrontDirection(*crack_front_point, tangent_direction, MIDDLE_NODE, i);
838 132 : _crack_directions.push_back(crack_direction);
839 :
840 132 : RankTwoTensor rot_mat;
841 132 : rot_mat.fillRow(0, crack_direction);
842 132 : rot_mat(2, _axis_2d) = 1.0;
843 :
844 132 : if (!_use_mesh_cutter)
845 132 : _crack_plane_normals.push_back(tangent_direction.cross(crack_direction));
846 :
847 : mooseAssert(i <= _crack_plane_normals.size(), "_crack_plane_normals is the wrong size.");
848 132 : rot_mat.fillRow(1, _crack_plane_normals[i]);
849 132 : _rot_matrix.push_back(rot_mat);
850 :
851 132 : _segment_lengths.push_back(std::make_pair(0.0, 0.0));
852 132 : _distances_along_front.push_back(0.0);
853 132 : _angles_along_front.push_back(0.0);
854 : }
855 : }
856 : else
857 : {
858 : // filling crack_plane_normals from mesh_cutter
859 75 : if (_use_mesh_cutter)
860 : _crack_plane_normals =
861 0 : _crack_front_points_provider->getCrackPlaneNormals(_num_points_from_provider);
862 : // else filling crack_plane_normals from curvedCrackFront
863 75 : computeCurvedCrackFrontCrackPlaneNormals();
864 :
865 75 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
866 75 : _segment_lengths.reserve(num_crack_front_points);
867 75 : _tangent_directions.reserve(num_crack_front_points);
868 75 : _crack_directions.reserve(num_crack_front_points);
869 :
870 : RealVectorValue back_segment;
871 : Real back_segment_len = 0.0;
872 75 : if (_closed_loop)
873 : {
874 7 : back_segment = *getCrackFrontPoint(0) - *getCrackFrontPoint(num_crack_front_points - 1);
875 7 : back_segment_len = back_segment.norm();
876 : }
877 :
878 561 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
879 : {
880 : CRACK_NODE_TYPE ntype;
881 486 : if (_closed_loop)
882 : ntype = MIDDLE_NODE;
883 276 : else if (i == 0)
884 : ntype = END_1_NODE;
885 208 : else if (i == num_crack_front_points - 1)
886 : ntype = END_2_NODE;
887 : else
888 : ntype = MIDDLE_NODE;
889 :
890 : RealVectorValue forward_segment;
891 : Real forward_segment_len;
892 486 : if (ntype == END_2_NODE)
893 : forward_segment_len = 0.0;
894 418 : else if (_closed_loop && i == num_crack_front_points - 1)
895 : {
896 7 : forward_segment = *getCrackFrontPoint(0) - *getCrackFrontPoint(i);
897 7 : forward_segment_len = forward_segment.norm();
898 : }
899 : else
900 : {
901 411 : forward_segment = *getCrackFrontPoint(i + 1) - *getCrackFrontPoint(i);
902 411 : forward_segment_len = forward_segment.norm();
903 411 : _overall_length += forward_segment_len;
904 : }
905 :
906 486 : _segment_lengths.push_back(std::make_pair(back_segment_len, forward_segment_len));
907 486 : if (i == 0)
908 75 : _distances_along_front.push_back(0.0);
909 : else
910 411 : _distances_along_front.push_back(back_segment_len + _distances_along_front[i - 1]);
911 :
912 : RealVectorValue tangent_direction = back_segment + forward_segment;
913 486 : tangent_direction = tangent_direction / tangent_direction.norm();
914 :
915 : // If end tangent directions are given, correct the tangent at the end nodes
916 486 : if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT &&
917 : _end_direction_method == END_DIRECTION_METHOD::END_CRACK_TANGENT_VECTOR)
918 : {
919 0 : if (ntype == END_1_NODE)
920 0 : tangent_direction = _crack_tangent_vector_end_1;
921 0 : else if (ntype == END_2_NODE)
922 0 : tangent_direction = _crack_tangent_vector_end_2;
923 : }
924 :
925 486 : _tangent_directions.push_back(tangent_direction);
926 : _crack_directions.push_back(
927 486 : calculateCrackFrontDirection(*getCrackFrontPoint(i), tangent_direction, ntype, i));
928 :
929 : // correct tangent direction in the case of _use_mesh_cutter
930 486 : if (_use_mesh_cutter)
931 0 : _tangent_directions[i] = _crack_plane_normals[i].cross(_crack_directions[i]);
932 :
933 : // If the end directions are given by the user, correct also the tangent at the end nodes
934 486 : if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT &&
935 96 : _end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR &&
936 96 : (ntype == END_1_NODE || ntype == END_2_NODE))
937 : {
938 32 : _tangent_directions[i] = _crack_plane_normals[i].cross(_crack_directions[i]);
939 : }
940 :
941 486 : back_segment = forward_segment;
942 : back_segment_len = forward_segment_len;
943 : }
944 :
945 : // For CURVED_CRACK_FRONT, _crack_plane_normals gets computed in
946 : // computeCurvedCrackFrontCrackPlaneNormals
947 75 : if (_direction_method != DIRECTION_METHOD::CURVED_CRACK_FRONT && !_use_mesh_cutter)
948 : {
949 52 : std::size_t mid_id = (num_crack_front_points - 1) / 2;
950 52 : _crack_plane_normals.assign(num_crack_front_points,
951 52 : _tangent_directions[mid_id].cross(_crack_directions[mid_id]));
952 :
953 : // Make sure the normal vector is non-zero
954 : RealVectorValue zero_vec(0.0);
955 52 : if (_crack_plane_normals.front().absolute_fuzzy_equals(zero_vec, _tol))
956 0 : mooseError("Crack plane normal vector evaluates to zero");
957 : }
958 :
959 : // Calculate angles of each point along the crack front for an elliptical crack projected
960 : // to a circle.
961 75 : if (hasAngleAlongFront())
962 : {
963 14 : RealVectorValue origin_to_first_node = *getCrackFrontPoint(0) - _crack_mouth_coordinates;
964 14 : Real hyp = origin_to_first_node.norm();
965 : RealVectorValue norm_origin_to_first_node = origin_to_first_node / hyp;
966 : RealVectorValue tangent_to_first_node =
967 14 : -norm_origin_to_first_node.cross(_crack_plane_normals.front());
968 14 : tangent_to_first_node /= tangent_to_first_node.norm();
969 :
970 80 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
971 : {
972 66 : RealVectorValue origin_to_curr_node = *getCrackFrontPoint(i) - _crack_mouth_coordinates;
973 :
974 : Real adj = origin_to_curr_node * norm_origin_to_first_node;
975 : Real opp = origin_to_curr_node * tangent_to_first_node;
976 :
977 66 : Real angle = acos(adj / hyp) * 180.0 / libMesh::pi;
978 66 : if (opp < 0.0)
979 56 : angle = 360.0 - angle;
980 66 : _angles_along_front.push_back(angle);
981 : }
982 :
983 : // Correct angle on end nodes if they are 0 or 360 to be consistent with neighboring node
984 14 : if (num_crack_front_points > 1)
985 : {
986 14 : if (MooseUtils::absoluteFuzzyEqual(_angles_along_front[0], 0.0, _tol) &&
987 10 : _angles_along_front[1] > 180.0)
988 10 : _angles_along_front[0] = 360.0;
989 4 : else if (MooseUtils::absoluteFuzzyEqual(_angles_along_front[0], 360.0, _tol) &&
990 0 : _angles_along_front[1] < 180.0)
991 0 : _angles_along_front[0] = 0.0;
992 :
993 : if (MooseUtils::absoluteFuzzyEqual(
994 14 : _angles_along_front[num_crack_front_points - 1], 0.0, _tol) &&
995 0 : _angles_along_front[num_crack_front_points - 2] > 180.0)
996 0 : _angles_along_front[num_crack_front_points - 1] = 360.0;
997 : else if (MooseUtils::absoluteFuzzyEqual(
998 14 : _angles_along_front[num_crack_front_points - 1], 360.0, _tol) &&
999 0 : _angles_along_front[num_crack_front_points - 2] < 180.0)
1000 0 : _angles_along_front[num_crack_front_points - 1] = 0.0;
1001 : }
1002 : }
1003 : else
1004 61 : _angles_along_front.resize(num_crack_front_points, 0.0);
1005 :
1006 : // Create rotation matrix
1007 561 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
1008 : {
1009 486 : RankTwoTensor rot_mat;
1010 486 : rot_mat.fillRow(0, _crack_directions[i]);
1011 486 : rot_mat.fillRow(1, _crack_plane_normals[i]);
1012 486 : rot_mat.fillRow(2, _tangent_directions[i]);
1013 486 : _rot_matrix.push_back(rot_mat);
1014 : }
1015 :
1016 75 : _console << "Summary of crack front geometry (used for fracture integrals):" << std::endl;
1017 : _console << "index node id x coord y coord z coord x dir y dir "
1018 75 : " z dir angle position seg length"
1019 75 : << std::endl;
1020 561 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
1021 : {
1022 : std::size_t point_id;
1023 486 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1024 434 : point_id = _ordered_crack_front_nodes[i];
1025 : else
1026 52 : point_id = i;
1027 486 : _console << std::left << std::setw(8) << i + 1 << std::setw(10) << point_id << std::setw(14)
1028 486 : << (*getCrackFrontPoint(i))(0) << std::setw(14) << (*getCrackFrontPoint(i))(1)
1029 486 : << std::setw(14) << (*getCrackFrontPoint(i))(2) << std::setw(14)
1030 486 : << _crack_directions[i](0) << std::setw(14) << _crack_directions[i](1)
1031 486 : << std::setw(14) << _crack_directions[i](2);
1032 486 : if (hasAngleAlongFront())
1033 66 : _console << std::left << std::setw(14) << _angles_along_front[i];
1034 : else
1035 420 : _console << std::left << std::setw(14) << "--";
1036 486 : _console << std::left << std::setw(14) << _distances_along_front[i] << std::setw(14)
1037 486 : << (_segment_lengths[i].first + _segment_lengths[i].second) / 2.0 << std::endl;
1038 : }
1039 75 : _console << "overall length: " << _overall_length << std::endl;
1040 : }
1041 195 : }
1042 :
1043 : void
1044 195 : CrackFrontDefinition::computeCrackMouthNodes()
1045 : {
1046 195 : if (_crack_mouth_boundary_ids.size() > 0)
1047 : {
1048 : _crack_mouth_coordinates.zero();
1049 :
1050 : std::set<Node *> crack_mouth_nodes;
1051 18 : ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
1052 2968 : for (auto nd = bnd_nodes.begin(); nd != bnd_nodes.end(); ++nd)
1053 : {
1054 2950 : const BndNode * bnode = *nd;
1055 2950 : BoundaryID boundary_id = bnode->_bnd_id;
1056 :
1057 5866 : for (std::size_t ibid = 0; ibid < _crack_mouth_boundary_ids.size(); ++ibid)
1058 : {
1059 2950 : if (boundary_id == _crack_mouth_boundary_ids[ibid])
1060 : {
1061 34 : crack_mouth_nodes.insert(bnode->_node);
1062 : break;
1063 : }
1064 : }
1065 : }
1066 :
1067 52 : for (auto nit = crack_mouth_nodes.begin(); nit != crack_mouth_nodes.end(); ++nit)
1068 : {
1069 34 : _crack_mouth_coordinates += **nit;
1070 : }
1071 18 : _crack_mouth_coordinates /= static_cast<Real>(crack_mouth_nodes.size());
1072 :
1073 18 : if (_has_symmetry_plane)
1074 10 : _crack_mouth_coordinates(_symmetry_plane) = 0.0;
1075 : }
1076 195 : }
1077 :
1078 : void
1079 75 : CrackFrontDefinition::computeCurvedCrackFrontCrackPlaneNormals()
1080 : {
1081 75 : if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT && !_use_mesh_cutter)
1082 : {
1083 23 : _crack_plane_normals.clear();
1084 :
1085 : // Get 3 nodes on crack front
1086 23 : std::size_t num_points = getNumCrackFrontPoints();
1087 23 : if (num_points < 3)
1088 : {
1089 0 : mooseError("Crack front must contain at least 3 nodes to use CurvedCrackFront option");
1090 : }
1091 : std::size_t start_id;
1092 : std::size_t mid_id;
1093 : std::size_t end_id;
1094 :
1095 23 : if (_closed_loop)
1096 : {
1097 : start_id = 0;
1098 7 : mid_id = (num_points - 1) / 3;
1099 7 : end_id = 2 * mid_id;
1100 : }
1101 : else
1102 : {
1103 : start_id = 0;
1104 16 : mid_id = (num_points - 1) / 2;
1105 : end_id = num_points - 1;
1106 : }
1107 23 : const Point * start = getCrackFrontPoint(start_id);
1108 23 : const Point * mid = getCrackFrontPoint(mid_id);
1109 23 : const Point * end = getCrackFrontPoint(end_id);
1110 :
1111 : // Create two vectors connecting them
1112 : RealVectorValue v1 = *mid - *start;
1113 : RealVectorValue v2 = *end - *mid;
1114 :
1115 : // Take cross product to get normal
1116 23 : Point crack_plane_normal = v1.cross(v2);
1117 23 : crack_plane_normal = crack_plane_normal.unit();
1118 23 : _crack_plane_normals.assign(num_points, crack_plane_normal);
1119 :
1120 : // Make sure they're not collinear
1121 : RealVectorValue zero_vec(0.0);
1122 23 : if (_crack_plane_normals.front().absolute_fuzzy_equals(zero_vec, _tol))
1123 : {
1124 0 : mooseError("Nodes on crack front are too close to being collinear");
1125 : }
1126 : }
1127 75 : }
1128 :
1129 : RealVectorValue
1130 618 : CrackFrontDefinition::calculateCrackFrontDirection(const Point & crack_front_point,
1131 : const RealVectorValue & tangent_direction,
1132 : const CRACK_NODE_TYPE ntype,
1133 : const std::size_t crack_front_point_index) const
1134 : {
1135 : RealVectorValue crack_dir;
1136 : RealVectorValue zero_vec(0.0);
1137 :
1138 : bool calc_dir = true;
1139 618 : if (_end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR)
1140 : {
1141 150 : if (ntype == END_1_NODE)
1142 : {
1143 26 : crack_dir = _crack_direction_vector_end_1;
1144 : calc_dir = false;
1145 : }
1146 124 : else if (ntype == END_2_NODE)
1147 : {
1148 26 : crack_dir = _crack_direction_vector_end_2;
1149 : calc_dir = false;
1150 : }
1151 : }
1152 :
1153 : if (calc_dir)
1154 : {
1155 566 : if (_direction_method == DIRECTION_METHOD::CRACK_DIRECTION_VECTOR)
1156 : {
1157 242 : crack_dir = _crack_direction_vector;
1158 : }
1159 324 : else if (_direction_method == DIRECTION_METHOD::CRACK_MOUTH)
1160 : {
1161 50 : if (_crack_mouth_coordinates.absolute_fuzzy_equals(crack_front_point, _tol))
1162 : {
1163 0 : mooseError("Crack mouth too close to crack front node");
1164 : }
1165 : RealVectorValue mouth_to_front = crack_front_point - _crack_mouth_coordinates;
1166 :
1167 50 : RealVectorValue crack_plane_normal = mouth_to_front.cross(tangent_direction);
1168 50 : if (crack_plane_normal.absolute_fuzzy_equals(zero_vec, _tol))
1169 : {
1170 0 : mooseError(
1171 : "Vector from crack mouth to crack front node is collinear with crack front segment");
1172 : }
1173 :
1174 50 : crack_dir = tangent_direction.cross(crack_plane_normal);
1175 : Real dotprod = crack_dir * mouth_to_front;
1176 50 : if (dotprod < 0)
1177 : {
1178 0 : crack_dir = -crack_dir;
1179 : }
1180 : }
1181 274 : else if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT)
1182 : {
1183 274 : crack_dir = tangent_direction.cross(_crack_plane_normals[crack_front_point_index]);
1184 : }
1185 : }
1186 618 : crack_dir = crack_dir.unit();
1187 :
1188 618 : return crack_dir;
1189 : }
1190 :
1191 : void
1192 0 : CrackFrontDefinition::updateNumberOfCrackFrontPoints(std::size_t num_points)
1193 : {
1194 0 : _num_points_from_provider = num_points;
1195 0 : }
1196 :
1197 : const Node *
1198 3540944 : CrackFrontDefinition::getCrackFrontNodePtr(const std::size_t node_index) const
1199 : {
1200 : mooseAssert(node_index < _ordered_crack_front_nodes.size(), "node_index out of range");
1201 3540944 : const Node * crack_front_node = _mesh.nodePtr(_ordered_crack_front_nodes[node_index]);
1202 : mooseAssert(crack_front_node != nullptr, "invalid crack front node");
1203 3540944 : return crack_front_node;
1204 : }
1205 :
1206 : const Point *
1207 3998543 : CrackFrontDefinition::getCrackFrontPoint(const std::size_t point_index) const
1208 : {
1209 3998543 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1210 : {
1211 3539882 : return getCrackFrontNodePtr(point_index);
1212 : }
1213 : else
1214 : {
1215 : mooseAssert(point_index < _crack_front_points.size(), "point_index out of range");
1216 458661 : return &_crack_front_points[point_index];
1217 : }
1218 : }
1219 :
1220 : const RealVectorValue &
1221 3657323 : CrackFrontDefinition::getCrackFrontTangent(const std::size_t point_index) const
1222 : {
1223 : mooseAssert(point_index < _tangent_directions.size(), "point_index out of range");
1224 3657323 : return _tangent_directions[point_index];
1225 : }
1226 :
1227 : const RealVectorValue &
1228 49320 : CrackFrontDefinition::getCrackFrontNormal(const std::size_t point_index) const
1229 : {
1230 : mooseAssert(point_index < _crack_plane_normals.size(), "point_index out of range");
1231 49320 : return _crack_plane_normals[point_index];
1232 : }
1233 :
1234 : Real
1235 1662929 : CrackFrontDefinition::getCrackFrontForwardSegmentLength(const std::size_t point_index) const
1236 : {
1237 1662929 : return _segment_lengths[point_index].second;
1238 : }
1239 :
1240 : Real
1241 1662929 : CrackFrontDefinition::getCrackFrontBackwardSegmentLength(const std::size_t point_index) const
1242 : {
1243 1662929 : return _segment_lengths[point_index].first;
1244 : }
1245 :
1246 : const RealVectorValue &
1247 2025360 : CrackFrontDefinition::getCrackDirection(const std::size_t point_index) const
1248 : {
1249 2025360 : return _crack_directions[point_index];
1250 : }
1251 :
1252 : Real
1253 2684 : CrackFrontDefinition::getDistanceAlongFront(const std::size_t point_index) const
1254 : {
1255 2684 : return _distances_along_front[point_index];
1256 : }
1257 :
1258 : bool
1259 750 : CrackFrontDefinition::hasAngleAlongFront() const
1260 : {
1261 750 : return (_crack_mouth_boundary_names.size() > 0);
1262 : }
1263 :
1264 : Real
1265 189 : CrackFrontDefinition::getAngleAlongFront(const std::size_t point_index) const
1266 : {
1267 189 : if (!hasAngleAlongFront())
1268 0 : paramError(
1269 : "crack_mouth_boundary",
1270 : "In CrackFrontDefinition, Requested angle along crack front, but definition of crack mouth "
1271 : "boundary using 'crack_mouth_boundary' parameter is necessary to do that.");
1272 189 : return _angles_along_front[point_index];
1273 : }
1274 :
1275 : std::size_t
1276 82698 : CrackFrontDefinition::getNumCrackFrontPoints() const
1277 : {
1278 82698 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1279 75620 : return _ordered_crack_front_nodes.size();
1280 : else
1281 7078 : return _crack_front_points.size();
1282 : }
1283 :
1284 : RealVectorValue
1285 577992 : CrackFrontDefinition::rotateToCrackFrontCoords(const RealVectorValue vector,
1286 : const std::size_t point_index) const
1287 : {
1288 577992 : return _rot_matrix[point_index] * vector;
1289 : }
1290 :
1291 : RealVectorValue
1292 0 : CrackFrontDefinition::rotateFromCrackFrontCoordsToGlobal(const RealVectorValue vector,
1293 : const std::size_t point_index) const
1294 : {
1295 0 : RealVectorValue vec = _rot_matrix[point_index].transpose() * vector;
1296 0 : return vec;
1297 : }
1298 :
1299 : RankTwoTensor
1300 242160 : CrackFrontDefinition::rotateToCrackFrontCoords(const RankTwoTensor tensor,
1301 : const std::size_t point_index) const
1302 : {
1303 242160 : RankTwoTensor tmp_tensor(tensor);
1304 242160 : tmp_tensor.rotate(_rot_matrix[point_index]);
1305 242160 : return tmp_tensor;
1306 : }
1307 :
1308 : void
1309 80720 : CrackFrontDefinition::calculateRThetaToCrackFront(const Point qp,
1310 : const std::size_t point_index,
1311 : Real & r,
1312 : Real & theta) const
1313 : {
1314 80720 : std::size_t num_points = getNumCrackFrontPoints();
1315 : Point closest_point(0.0);
1316 : RealVectorValue closest_point_to_p;
1317 :
1318 80720 : const Point * crack_front_point = getCrackFrontPoint(point_index);
1319 80720 : RealVectorValue crack_front_point_rot = rotateToCrackFrontCoords(*crack_front_point, point_index);
1320 :
1321 : RealVectorValue crack_front_edge =
1322 80720 : rotateToCrackFrontCoords(_tangent_directions[point_index], point_index);
1323 :
1324 80720 : Point p_rot = rotateToCrackFrontCoords(qp, point_index);
1325 : p_rot = p_rot - crack_front_point_rot;
1326 :
1327 80720 : if (_treat_as_2d)
1328 : {
1329 : // In 2D, the closest node is the crack tip node and the position of the crack tip node is
1330 : // (0,0,0) in the crack front coordinate system
1331 : // In case this is a 3D mesh treated as 2D, project point onto same plane as crack front node.
1332 : // Note: In the crack front coordinate system, z is always in the tangent direction to the crack
1333 : // front
1334 : p_rot(2) = closest_point(2);
1335 : closest_point_to_p = p_rot;
1336 :
1337 : // Find r, the distance between the qp and the crack front
1338 : RealVectorValue r_vec = p_rot;
1339 36368 : r = r_vec.norm();
1340 : }
1341 : else
1342 : {
1343 : // Loop over crack front points to find the one closest to the point qp
1344 : Real min_dist = std::numeric_limits<Real>::max();
1345 248576 : for (std::size_t pit = 0; pit != num_points; ++pit)
1346 : {
1347 204224 : const Point * crack_front_point = getCrackFrontPoint(pit);
1348 : RealVectorValue crack_point_to_current_point = qp - *crack_front_point;
1349 204224 : Real dist = crack_point_to_current_point.norm();
1350 :
1351 204224 : if (dist < min_dist)
1352 : {
1353 : min_dist = dist;
1354 : closest_point = *crack_front_point;
1355 : }
1356 : }
1357 :
1358 : // Rotate coordinates to crack front coordinate system
1359 44352 : closest_point = rotateToCrackFrontCoords(closest_point, point_index);
1360 : closest_point = closest_point - crack_front_point_rot;
1361 :
1362 : // Find r, the distance between the qp and the crack front
1363 : Real edge_length_sq = crack_front_edge.norm_sq();
1364 : closest_point_to_p = p_rot - closest_point;
1365 : Real perp = crack_front_edge * closest_point_to_p;
1366 44352 : Real dist_along_edge = perp / edge_length_sq;
1367 : RealVectorValue point_on_edge = closest_point + crack_front_edge * dist_along_edge;
1368 : RealVectorValue r_vec = p_rot - point_on_edge;
1369 44352 : r = r_vec.norm();
1370 : }
1371 :
1372 : // Find theta, the angle between r and the crack front plane
1373 : RealVectorValue crack_plane_normal =
1374 80720 : rotateToCrackFrontCoords(_crack_plane_normals[point_index], point_index);
1375 : Real p_to_plane_dist = std::abs(closest_point_to_p * crack_plane_normal);
1376 :
1377 : // Determine if qp is above or below the crack plane
1378 80720 : Real y_local = p_rot(1) - closest_point(1);
1379 :
1380 : // Determine if qp is in front of or behind the crack front
1381 : RealVectorValue p2(p_rot);
1382 : p2(1) = 0;
1383 : RealVectorValue p2_vec = p2 - closest_point;
1384 80720 : Real ahead = crack_front_edge(2) * p2_vec(0) - crack_front_edge(0) * p2_vec(2);
1385 :
1386 : Real x_local(0);
1387 80720 : if (ahead >= 0)
1388 : x_local = 1;
1389 : else
1390 : x_local = -1;
1391 :
1392 : // Calculate theta based on in which quadrant in the crack front coordinate
1393 : // system the qp is located
1394 80720 : if (r > 0)
1395 : {
1396 : Real theta_quadrant1(0.0);
1397 80720 : if (MooseUtils::absoluteFuzzyEqual(r, p_to_plane_dist, _tol))
1398 : theta_quadrant1 = 0.5 * libMesh::pi;
1399 80720 : else if (p_to_plane_dist > r)
1400 0 : mooseError(
1401 : "Invalid distance p_to_plane_dist in CrackFrontDefinition::calculateRThetaToCrackFront");
1402 : else
1403 80720 : theta_quadrant1 = std::asin(p_to_plane_dist / r);
1404 :
1405 80720 : if (x_local >= 0 && y_local >= 0)
1406 18536 : theta = theta_quadrant1;
1407 :
1408 62184 : else if (x_local < 0 && y_local >= 0)
1409 17832 : theta = libMesh::pi - theta_quadrant1;
1410 :
1411 17168 : else if (x_local < 0 && y_local < 0)
1412 17168 : theta = -(libMesh::pi - theta_quadrant1);
1413 :
1414 27184 : else if (x_local >= 0 && y_local < 0)
1415 27184 : theta = -theta_quadrant1;
1416 : }
1417 0 : else if (r == 0)
1418 0 : theta = 0;
1419 : else
1420 0 : mooseError("Invalid distance r in CrackFrontDefinition::calculateRThetaToCrackFront");
1421 80720 : }
1422 :
1423 : std::size_t
1424 0 : CrackFrontDefinition::calculateRThetaToCrackFront(const Point qp, Real & r, Real & theta) const
1425 : {
1426 0 : std::size_t num_points = getNumCrackFrontPoints();
1427 :
1428 : // Loop over crack front points to find the one closest to the point qp
1429 : Real min_dist = std::numeric_limits<Real>::max();
1430 : std::size_t point_index = 0;
1431 0 : for (std::size_t pit = 0; pit != num_points; ++pit)
1432 : {
1433 0 : const Point * crack_front_point = getCrackFrontPoint(pit);
1434 : RealVectorValue crack_point_to_current_point = qp - *crack_front_point;
1435 0 : Real dist = crack_point_to_current_point.norm();
1436 :
1437 0 : if (dist < min_dist)
1438 : {
1439 : min_dist = dist;
1440 : point_index = pit;
1441 : }
1442 : }
1443 :
1444 0 : calculateRThetaToCrackFront(qp, point_index, r, theta);
1445 :
1446 0 : return point_index;
1447 : }
1448 :
1449 : bool
1450 94315 : CrackFrontDefinition::isNodeOnIntersectingBoundary(const Node * const node) const
1451 : {
1452 : bool is_on_boundary = false;
1453 : mooseAssert(node, "Invalid node");
1454 : dof_id_type node_id = node->id();
1455 129440 : for (std::size_t i = 0; i < _intersecting_boundary_ids.size(); ++i)
1456 : {
1457 38204 : if (_mesh.isBoundaryNode(node_id, _intersecting_boundary_ids[i]))
1458 : {
1459 : is_on_boundary = true;
1460 : break;
1461 : }
1462 : }
1463 94315 : return is_on_boundary;
1464 : }
1465 :
1466 : bool
1467 736 : CrackFrontDefinition::isPointWithIndexOnIntersectingBoundary(const std::size_t point_index) const
1468 : {
1469 : bool is_on_boundary = false;
1470 736 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1471 : {
1472 668 : const Node * crack_front_node = getCrackFrontNodePtr(point_index);
1473 668 : is_on_boundary = isNodeOnIntersectingBoundary(crack_front_node);
1474 : }
1475 : else
1476 : {
1477 : // If the intersecting boundary option is used with crack front points, the
1478 : // first and last points are assumed to be on the intersecting boundaries.
1479 68 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
1480 68 : if (point_index == 0 || point_index == num_crack_front_points - 1)
1481 : is_on_boundary = true;
1482 : }
1483 736 : return is_on_boundary;
1484 : }
1485 :
1486 : void
1487 18 : CrackFrontDefinition::calculateTangentialStrainAlongFront()
1488 : {
1489 : RealVectorValue disp_current_node;
1490 : RealVectorValue disp_previous_node;
1491 : RealVectorValue disp_next_node;
1492 :
1493 : RealVectorValue forward_segment0;
1494 : RealVectorValue forward_segment1;
1495 : Real forward_segment0_len;
1496 : Real forward_segment1_len;
1497 : RealVectorValue back_segment0;
1498 : RealVectorValue back_segment1;
1499 : Real back_segment0_len;
1500 : Real back_segment1_len;
1501 :
1502 : std::size_t num_crack_front_nodes = _ordered_crack_front_nodes.size();
1503 : const Node * current_node;
1504 : const Node * previous_node;
1505 : const Node * next_node;
1506 :
1507 18 : _strain_along_front.reserve(num_crack_front_nodes);
1508 :
1509 : // In finalize(), gatherMax builds and distributes the complete strain vector on all processors
1510 : // -> reset the vector every time
1511 144 : for (std::size_t i = 0; i < num_crack_front_nodes; ++i)
1512 126 : _strain_along_front[i] = -std::numeric_limits<Real>::max();
1513 :
1514 18 : MooseVariable & disp_x_var = _subproblem.getStandardVariable(_tid, _disp_x_var_name);
1515 18 : MooseVariable & disp_y_var = _subproblem.getStandardVariable(_tid, _disp_y_var_name);
1516 18 : MooseVariable & disp_z_var = _subproblem.getStandardVariable(_tid, _disp_z_var_name);
1517 :
1518 18 : current_node = getCrackFrontNodePtr(0);
1519 18 : if (current_node->processor_id() == processor_id())
1520 : {
1521 14 : disp_current_node(0) = disp_x_var.getNodalValue(*current_node);
1522 14 : disp_current_node(1) = disp_y_var.getNodalValue(*current_node);
1523 14 : disp_current_node(2) = disp_z_var.getNodalValue(*current_node);
1524 :
1525 14 : next_node = getCrackFrontNodePtr(1);
1526 14 : disp_next_node(0) = disp_x_var.getNodalValue(*next_node);
1527 14 : disp_next_node(1) = disp_y_var.getNodalValue(*next_node);
1528 14 : disp_next_node(2) = disp_z_var.getNodalValue(*next_node);
1529 :
1530 : forward_segment0 = *next_node - *current_node;
1531 14 : forward_segment0 = (forward_segment0 * _tangent_directions[0]) * _tangent_directions[0];
1532 14 : forward_segment0_len = forward_segment0.norm();
1533 :
1534 : forward_segment1 = (*next_node + disp_next_node) - (*current_node + disp_current_node);
1535 14 : forward_segment1 = (forward_segment1 * _tangent_directions[0]) * _tangent_directions[0];
1536 14 : forward_segment1_len = forward_segment1.norm();
1537 :
1538 14 : _strain_along_front[0] = (forward_segment1_len - forward_segment0_len) / forward_segment0_len;
1539 : }
1540 :
1541 108 : for (std::size_t i = 1; i < num_crack_front_nodes - 1; ++i)
1542 : {
1543 90 : current_node = getCrackFrontNodePtr(i);
1544 90 : if (current_node->processor_id() == processor_id())
1545 : {
1546 70 : disp_current_node(0) = disp_x_var.getNodalValue(*current_node);
1547 70 : disp_current_node(1) = disp_y_var.getNodalValue(*current_node);
1548 70 : disp_current_node(2) = disp_z_var.getNodalValue(*current_node);
1549 :
1550 70 : previous_node = getCrackFrontNodePtr(i - 1);
1551 70 : disp_previous_node(0) = disp_x_var.getNodalValue(*previous_node);
1552 70 : disp_previous_node(1) = disp_y_var.getNodalValue(*previous_node);
1553 70 : disp_previous_node(2) = disp_z_var.getNodalValue(*previous_node);
1554 :
1555 70 : next_node = getCrackFrontNodePtr(i + 1);
1556 70 : disp_next_node(0) = disp_x_var.getNodalValue(*next_node);
1557 70 : disp_next_node(1) = disp_y_var.getNodalValue(*next_node);
1558 70 : disp_next_node(2) = disp_z_var.getNodalValue(*next_node);
1559 :
1560 : back_segment0 = *current_node - *previous_node;
1561 70 : back_segment0 = (back_segment0 * _tangent_directions[i]) * _tangent_directions[i];
1562 70 : back_segment0_len = back_segment0.norm();
1563 :
1564 : back_segment1 = (*current_node + disp_current_node) - (*previous_node + disp_previous_node);
1565 70 : back_segment1 = (back_segment1 * _tangent_directions[i]) * _tangent_directions[i];
1566 70 : back_segment1_len = back_segment1.norm();
1567 :
1568 : forward_segment0 = *next_node - *current_node;
1569 70 : forward_segment0 = (forward_segment0 * _tangent_directions[i]) * _tangent_directions[i];
1570 70 : forward_segment0_len = forward_segment0.norm();
1571 :
1572 : forward_segment1 = (*next_node + disp_next_node) - (*current_node + disp_current_node);
1573 70 : forward_segment1 = (forward_segment1 * _tangent_directions[i]) * _tangent_directions[i];
1574 70 : forward_segment1_len = forward_segment1.norm();
1575 :
1576 70 : _strain_along_front[i] =
1577 70 : 0.5 * ((back_segment1_len - back_segment0_len) / back_segment0_len +
1578 70 : (forward_segment1_len - forward_segment0_len) / forward_segment0_len);
1579 : }
1580 : }
1581 :
1582 18 : current_node = getCrackFrontNodePtr(num_crack_front_nodes - 1);
1583 18 : if (current_node->processor_id() == processor_id())
1584 : {
1585 14 : disp_current_node(0) = disp_x_var.getNodalValue(*current_node);
1586 14 : disp_current_node(1) = disp_y_var.getNodalValue(*current_node);
1587 14 : disp_current_node(2) = disp_z_var.getNodalValue(*current_node);
1588 :
1589 14 : previous_node = getCrackFrontNodePtr(num_crack_front_nodes - 2);
1590 14 : disp_previous_node(0) = disp_x_var.getNodalValue(*previous_node);
1591 14 : disp_previous_node(1) = disp_y_var.getNodalValue(*previous_node);
1592 14 : disp_previous_node(2) = disp_z_var.getNodalValue(*previous_node);
1593 :
1594 : back_segment0 = *current_node - *previous_node;
1595 14 : back_segment0 = (back_segment0 * _tangent_directions[num_crack_front_nodes - 1]) *
1596 : _tangent_directions[num_crack_front_nodes - 1];
1597 14 : back_segment0_len = back_segment0.norm();
1598 :
1599 : back_segment1 = (*current_node + disp_current_node) - (*previous_node + disp_previous_node);
1600 14 : back_segment1 = (back_segment1 * _tangent_directions[num_crack_front_nodes - 1]) *
1601 : _tangent_directions[num_crack_front_nodes - 1];
1602 14 : back_segment1_len = back_segment1.norm();
1603 :
1604 14 : _strain_along_front[num_crack_front_nodes - 1] =
1605 14 : (back_segment1_len - back_segment0_len) / back_segment0_len;
1606 : }
1607 18 : }
1608 :
1609 : Real
1610 210 : CrackFrontDefinition::getCrackFrontTangentialStrain(const std::size_t node_index) const
1611 : {
1612 : Real strain;
1613 210 : if (_t_stress)
1614 : {
1615 210 : strain = _strain_along_front[node_index];
1616 : mooseAssert(strain > -std::numeric_limits<Real>::max(),
1617 : "Failure in parallel communication of crack tangential strain");
1618 : }
1619 : else
1620 0 : mooseError("In CrackFrontDefinition, tangential strain not available");
1621 :
1622 210 : return strain;
1623 : }
1624 :
1625 : void
1626 20 : CrackFrontDefinition::createQFunctionRings()
1627 : {
1628 : // In the variable names, "cfn" = crack front node
1629 20 : if (_treat_as_2d && _use_mesh_cutter == false) // 2D: the q-function defines an integral domain
1630 : // that is constant along the crack front
1631 : {
1632 : std::vector<std::vector<const Elem *>> nodes_to_elem_map;
1633 12 : MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), nodes_to_elem_map);
1634 :
1635 : std::set<dof_id_type> nodes_prev_ring;
1636 12 : nodes_prev_ring.insert(_ordered_crack_front_nodes.begin(), _ordered_crack_front_nodes.end());
1637 :
1638 : std::set<dof_id_type> connected_nodes_this_cfn;
1639 12 : connected_nodes_this_cfn.insert(_ordered_crack_front_nodes.begin(),
1640 : _ordered_crack_front_nodes.end());
1641 :
1642 : std::set<dof_id_type> old_ring_nodes_this_cfn = connected_nodes_this_cfn;
1643 :
1644 : // The first ring contains only the crack front node(s)
1645 : std::pair<dof_id_type, std::size_t> node_ring_index =
1646 : std::make_pair(_ordered_crack_front_nodes[0], 1);
1647 12 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1648 : connected_nodes_this_cfn.end());
1649 :
1650 : // Build rings of nodes around the crack front node
1651 44 : for (std::size_t ring = 2; ring <= _last_ring; ++ring)
1652 : {
1653 :
1654 : // Find nodes connected to the nodes of the previous ring
1655 : std::set<dof_id_type> new_ring_nodes_this_cfn;
1656 136 : for (auto nit = old_ring_nodes_this_cfn.begin(); nit != old_ring_nodes_this_cfn.end(); ++nit)
1657 : {
1658 : std::vector<const Node *> neighbors;
1659 104 : MeshTools::find_nodal_neighbors(
1660 104 : _mesh.getMesh(), _mesh.nodeRef(*nit), nodes_to_elem_map, neighbors);
1661 488 : for (std::size_t inei = 0; inei < neighbors.size(); ++inei)
1662 : {
1663 384 : auto thisit = connected_nodes_this_cfn.find(neighbors[inei]->id());
1664 :
1665 : // Add only nodes that are not already present in any of the rings
1666 384 : if (thisit == connected_nodes_this_cfn.end())
1667 248 : new_ring_nodes_this_cfn.insert(neighbors[inei]->id());
1668 : }
1669 104 : }
1670 :
1671 : // Add new nodes to rings
1672 32 : connected_nodes_this_cfn.insert(new_ring_nodes_this_cfn.begin(),
1673 : new_ring_nodes_this_cfn.end());
1674 : old_ring_nodes_this_cfn = new_ring_nodes_this_cfn;
1675 :
1676 : std::pair<dof_id_type, std::size_t> node_ring_index =
1677 32 : std::make_pair(_ordered_crack_front_nodes[0], ring);
1678 64 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1679 : connected_nodes_this_cfn.end());
1680 : }
1681 12 : }
1682 : else // The q-function defines one integral domain around each crack front node
1683 : {
1684 : std::size_t num_crack_front_points = _ordered_crack_front_nodes.size();
1685 : std::vector<std::vector<const Elem *>> nodes_to_elem_map;
1686 8 : MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), nodes_to_elem_map);
1687 32 : for (std::size_t icfn = 0; icfn < num_crack_front_points; ++icfn)
1688 : {
1689 : std::set<dof_id_type> nodes_prev_ring;
1690 24 : nodes_prev_ring.insert(_ordered_crack_front_nodes[icfn]);
1691 :
1692 : std::set<dof_id_type> connected_nodes_prev_cfn;
1693 : std::set<dof_id_type> connected_nodes_this_cfn;
1694 : std::set<dof_id_type> connected_nodes_next_cfn;
1695 :
1696 24 : connected_nodes_this_cfn.insert(_ordered_crack_front_nodes[icfn]);
1697 :
1698 24 : if (_closed_loop && icfn == 0)
1699 : {
1700 0 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[num_crack_front_points - 1]);
1701 0 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[icfn + 1]);
1702 : }
1703 24 : else if (_closed_loop && icfn == num_crack_front_points - 1)
1704 : {
1705 0 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[icfn - 1]);
1706 0 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[0]);
1707 : }
1708 24 : else if (icfn == 0)
1709 : {
1710 8 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[icfn + 1]);
1711 : }
1712 16 : else if (icfn == num_crack_front_points - 1)
1713 : {
1714 8 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[icfn - 1]);
1715 : }
1716 : else
1717 : {
1718 8 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[icfn - 1]);
1719 8 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[icfn + 1]);
1720 : }
1721 :
1722 : std::set<dof_id_type> old_ring_nodes_prev_cfn = connected_nodes_prev_cfn;
1723 : std::set<dof_id_type> old_ring_nodes_this_cfn = connected_nodes_this_cfn;
1724 : std::set<dof_id_type> old_ring_nodes_next_cfn = connected_nodes_next_cfn;
1725 :
1726 : // The first ring contains only the crack front node
1727 : std::pair<dof_id_type, std::size_t> node_ring_index =
1728 : std::make_pair(_ordered_crack_front_nodes[icfn], 1);
1729 24 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1730 : connected_nodes_this_cfn.end());
1731 :
1732 : // Build rings of nodes around the crack front node
1733 72 : for (std::size_t ring = 2; ring <= _last_ring; ++ring)
1734 : {
1735 :
1736 : // Find nodes connected to the nodes of the previous ring, but exclude nodes in rings of
1737 : // neighboring crack front nodes
1738 : std::set<dof_id_type> new_ring_nodes_this_cfn;
1739 48 : addNodesToQFunctionRing(new_ring_nodes_this_cfn,
1740 : old_ring_nodes_this_cfn,
1741 : connected_nodes_this_cfn,
1742 : connected_nodes_prev_cfn,
1743 : connected_nodes_next_cfn,
1744 : nodes_to_elem_map);
1745 :
1746 : std::set<dof_id_type> new_ring_nodes_prev_cfn;
1747 48 : addNodesToQFunctionRing(new_ring_nodes_prev_cfn,
1748 : old_ring_nodes_prev_cfn,
1749 : connected_nodes_prev_cfn,
1750 : connected_nodes_this_cfn,
1751 : connected_nodes_next_cfn,
1752 : nodes_to_elem_map);
1753 :
1754 : std::set<dof_id_type> new_ring_nodes_next_cfn;
1755 48 : addNodesToQFunctionRing(new_ring_nodes_next_cfn,
1756 : old_ring_nodes_next_cfn,
1757 : connected_nodes_next_cfn,
1758 : connected_nodes_prev_cfn,
1759 : connected_nodes_this_cfn,
1760 : nodes_to_elem_map);
1761 :
1762 : // Add new nodes to the three sets of nodes
1763 48 : connected_nodes_prev_cfn.insert(new_ring_nodes_prev_cfn.begin(),
1764 : new_ring_nodes_prev_cfn.end());
1765 48 : connected_nodes_this_cfn.insert(new_ring_nodes_this_cfn.begin(),
1766 : new_ring_nodes_this_cfn.end());
1767 48 : connected_nodes_next_cfn.insert(new_ring_nodes_next_cfn.begin(),
1768 : new_ring_nodes_next_cfn.end());
1769 : old_ring_nodes_prev_cfn = new_ring_nodes_prev_cfn;
1770 : old_ring_nodes_this_cfn = new_ring_nodes_this_cfn;
1771 : old_ring_nodes_next_cfn = new_ring_nodes_next_cfn;
1772 :
1773 : std::pair<dof_id_type, std::size_t> node_ring_index =
1774 48 : std::make_pair(_ordered_crack_front_nodes[icfn], ring);
1775 96 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1776 : connected_nodes_this_cfn.end());
1777 : }
1778 : }
1779 8 : }
1780 20 : }
1781 :
1782 : void
1783 144 : CrackFrontDefinition::addNodesToQFunctionRing(
1784 : std::set<dof_id_type> & nodes_new_ring,
1785 : const std::set<dof_id_type> & nodes_old_ring,
1786 : const std::set<dof_id_type> & nodes_all_rings,
1787 : const std::set<dof_id_type> & nodes_neighbor1,
1788 : const std::set<dof_id_type> & nodes_neighbor2,
1789 : std::vector<std::vector<const Elem *>> & nodes_to_elem_map)
1790 : {
1791 384 : for (auto nit = nodes_old_ring.begin(); nit != nodes_old_ring.end(); ++nit)
1792 : {
1793 : std::vector<const Node *> neighbors;
1794 240 : MeshTools::find_nodal_neighbors(
1795 240 : _mesh.getMesh(), _mesh.nodeRef(*nit), nodes_to_elem_map, neighbors);
1796 1352 : for (std::size_t inei = 0; inei < neighbors.size(); ++inei)
1797 : {
1798 1112 : auto previt = nodes_all_rings.find(neighbors[inei]->id());
1799 1112 : auto thisit = nodes_neighbor1.find(neighbors[inei]->id());
1800 1112 : auto nextit = nodes_neighbor2.find(neighbors[inei]->id());
1801 :
1802 : // Add only nodes that are not already present in any of the three sets of nodes
1803 1112 : if (previt == nodes_all_rings.end() && thisit == nodes_neighbor1.end() &&
1804 : nextit == nodes_neighbor2.end())
1805 672 : nodes_new_ring.insert(neighbors[inei]->id());
1806 : }
1807 240 : }
1808 144 : }
1809 :
1810 : bool
1811 57377 : CrackFrontDefinition::isNodeInRing(const std::size_t ring_index,
1812 : const dof_id_type connected_node_id,
1813 : const std::size_t node_index) const
1814 : {
1815 : bool is_node_in_ring = false;
1816 : std::pair<dof_id_type, std::size_t> node_ring_key =
1817 57377 : std::make_pair(_ordered_crack_front_nodes[node_index], ring_index);
1818 : auto nnmit = _crack_front_node_to_node_map.find(node_ring_key);
1819 :
1820 57377 : if (nnmit == _crack_front_node_to_node_map.end())
1821 0 : mooseError("Could not find crack front node ",
1822 : _ordered_crack_front_nodes[node_index],
1823 : " in the crack front node to q-function ring-node map for ring ",
1824 : ring_index);
1825 :
1826 : std::set<dof_id_type> q_func_nodes = nnmit->second;
1827 57377 : if (q_func_nodes.find(connected_node_id) != q_func_nodes.end())
1828 : is_node_in_ring = true;
1829 :
1830 57377 : return is_node_in_ring;
1831 : }
1832 :
1833 : Real
1834 3606784 : CrackFrontDefinition::DomainIntegralQFunction(std::size_t crack_front_point_index,
1835 : std::size_t ring_index,
1836 : const Node * const current_node) const
1837 : {
1838 : Real dist_to_crack_front;
1839 : Real dist_along_tangent;
1840 3606784 : projectToFrontAtPoint(
1841 : dist_to_crack_front, dist_along_tangent, crack_front_point_index, current_node);
1842 :
1843 : Real q = 1.0;
1844 3606784 : if (dist_to_crack_front > _j_integral_radius_inner[ring_index] &&
1845 3539282 : dist_to_crack_front < _j_integral_radius_outer[ring_index])
1846 21339 : q = (_j_integral_radius_outer[ring_index] - dist_to_crack_front) /
1847 21339 : (_j_integral_radius_outer[ring_index] - _j_integral_radius_inner[ring_index]);
1848 3585445 : else if (dist_to_crack_front >= _j_integral_radius_outer[ring_index])
1849 : q = 0.0;
1850 :
1851 21339 : if (q > 0.0)
1852 : {
1853 88809 : Real tangent_multiplier = 1.0;
1854 88809 : if (!_treat_as_2d)
1855 : {
1856 : const Real forward_segment_length =
1857 40577 : getCrackFrontForwardSegmentLength(crack_front_point_index);
1858 : const Real backward_segment_length =
1859 40577 : getCrackFrontBackwardSegmentLength(crack_front_point_index);
1860 :
1861 40577 : if (dist_along_tangent >= 0.0)
1862 : {
1863 26395 : if (forward_segment_length > 0.0)
1864 24280 : tangent_multiplier = 1.0 - dist_along_tangent / forward_segment_length;
1865 : }
1866 : else
1867 : {
1868 14182 : if (backward_segment_length > 0.0)
1869 14024 : tangent_multiplier = 1.0 + dist_along_tangent / backward_segment_length;
1870 : }
1871 : }
1872 :
1873 88809 : tangent_multiplier = std::max(tangent_multiplier, 0.0);
1874 88809 : tangent_multiplier = std::min(tangent_multiplier, 1.0);
1875 :
1876 : // Set to zero if a node is on a designated free surface and its crack front node is not.
1877 88809 : if (isNodeOnIntersectingBoundary(current_node) &&
1878 : !_is_point_on_intersecting_boundary[crack_front_point_index])
1879 2180 : tangent_multiplier = 0.0;
1880 :
1881 88809 : q *= tangent_multiplier;
1882 : }
1883 :
1884 3606784 : return q;
1885 : }
1886 :
1887 : Real
1888 46464 : CrackFrontDefinition::DomainIntegralTopologicalQFunction(std::size_t crack_front_point_index,
1889 : std::size_t ring_index,
1890 : const Node * const current_node) const
1891 : {
1892 : Real q = 0;
1893 46464 : bool is_node_in_ring = isNodeInRing(ring_index, current_node->id(), crack_front_point_index);
1894 46464 : if (is_node_in_ring)
1895 : q = 1;
1896 :
1897 46464 : return q;
1898 : }
1899 :
1900 : void
1901 3606784 : CrackFrontDefinition::projectToFrontAtPoint(Real & dist_to_front,
1902 : Real & dist_along_tangent,
1903 : std::size_t crack_front_point_index,
1904 : const Node * const current_node) const
1905 : {
1906 3606784 : const Point * crack_front_point = getCrackFrontPoint(crack_front_point_index);
1907 :
1908 3606784 : Point p = *current_node;
1909 3606784 : const RealVectorValue & crack_front_tangent = getCrackFrontTangent(crack_front_point_index);
1910 :
1911 : RealVectorValue crack_node_to_current_node = p - *crack_front_point;
1912 3606784 : dist_along_tangent = crack_node_to_current_node * crack_front_tangent;
1913 : RealVectorValue projection_point = *crack_front_point + dist_along_tangent * crack_front_tangent;
1914 : RealVectorValue axis_to_current_node = p - projection_point;
1915 3606784 : dist_to_front = axis_to_current_node.norm();
1916 3606784 : }
1917 :
1918 : void
1919 0 : CrackFrontDefinition::isCutterModified(const bool is_cutter_modified)
1920 : {
1921 0 : _is_cutter_modified = is_cutter_modified;
1922 0 : }
|
