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 : const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map =
492 75 : _mesh.nodeToElemMap();
493 : std::map<dof_id_type, std::set<dof_id_type>> crack_front_node_to_elem_map;
494 :
495 533 : for (const auto & node_id : nodes)
496 : {
497 : const auto & node_to_elem_pair = node_to_elem_map.find(node_id);
498 : mooseAssert(node_to_elem_pair != node_to_elem_map.end(),
499 : "Could not find crack front node " << node_id << " in the node to elem map");
500 :
501 : const std::vector<dof_id_type> & connected_elems = node_to_elem_pair->second;
502 2058 : for (std::size_t i = 0; i < connected_elems.size(); ++i)
503 1600 : crack_front_node_to_elem_map[node_id].insert(connected_elems[i]);
504 : }
505 :
506 : // Determine which nodes are connected to each other via elements, and construct line elements
507 : // to represent
508 : // those connections
509 : std::vector<std::vector<dof_id_type>> line_elems;
510 : std::map<dof_id_type, std::vector<dof_id_type>> node_to_line_elem_map;
511 :
512 75 : for (auto cfnemit = crack_front_node_to_elem_map.begin();
513 533 : cfnemit != crack_front_node_to_elem_map.end();
514 : ++cfnemit)
515 : {
516 : auto cfnemit2 = cfnemit;
517 3935 : for (++cfnemit2; cfnemit2 != crack_front_node_to_elem_map.end(); ++cfnemit2)
518 : {
519 :
520 : std::vector<dof_id_type> common_elements;
521 : std::set<dof_id_type> & elements_connected_to_node1 = cfnemit->second;
522 : std::set<dof_id_type> & elements_connected_to_node2 = cfnemit2->second;
523 3477 : std::set_intersection(elements_connected_to_node1.begin(),
524 : elements_connected_to_node1.end(),
525 : elements_connected_to_node2.begin(),
526 : elements_connected_to_node2.end(),
527 : std::inserter(common_elements, common_elements.end()));
528 :
529 3477 : if (common_elements.size() > 0)
530 : {
531 : std::vector<dof_id_type> my_line_elem;
532 390 : my_line_elem.push_back(cfnemit->first);
533 390 : my_line_elem.push_back(cfnemit2->first);
534 390 : node_to_line_elem_map[cfnemit->first].push_back(line_elems.size());
535 390 : node_to_line_elem_map[cfnemit2->first].push_back(line_elems.size());
536 390 : line_elems.push_back(my_line_elem);
537 390 : }
538 3477 : }
539 : }
540 :
541 : // Find nodes on ends of line (those connected to only one line element)
542 : std::vector<dof_id_type> end_nodes;
543 533 : for (auto nlemit = node_to_line_elem_map.begin(); nlemit != node_to_line_elem_map.end();
544 : ++nlemit)
545 : {
546 : std::size_t num_connected_elems = nlemit->second.size();
547 458 : if (num_connected_elems == 1)
548 136 : end_nodes.push_back(nlemit->first);
549 322 : else if (num_connected_elems != 2)
550 0 : mooseError(
551 0 : "Node ", nlemit->first, " is connected to >2 line segments in CrackFrontDefinition");
552 : }
553 :
554 : // For embedded crack with closed loop of crack front nodes, must pick the end nodes
555 75 : if (end_nodes.size() == 0) // Crack front is a loop. Pick nodes to be end nodes.
556 : {
557 7 : pickLoopCrackEndNodes(end_nodes, nodes, node_to_line_elem_map, line_elems);
558 7 : _closed_loop = true;
559 7 : if (_end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR ||
560 : _end_direction_method == END_DIRECTION_METHOD::END_CRACK_TANGENT_VECTOR)
561 0 : paramError("end_direction_method",
562 : "In CrackFrontDefinition, end_direction_method cannot be CrackDirectionVector "
563 : "or CrackTangentVector for a closed-loop crack");
564 7 : if (_intersecting_boundary_names.size() > 0)
565 0 : paramError("intersecting_boundary",
566 : "In CrackFrontDefinition, intersecting_boundary cannot be specified for a "
567 : "closed-loop crack");
568 : }
569 68 : else if (end_nodes.size() == 2) // Rearrange the order of the end nodes if needed
570 68 : orderEndNodes(end_nodes);
571 : else
572 0 : mooseError("In CrackFrontDefinition wrong number of end nodes. Number end nodes = ",
573 0 : end_nodes.size());
574 :
575 : // Create an ordered list of the nodes going along the line of the crack front
576 75 : _ordered_crack_front_nodes.push_back(end_nodes[0]);
577 :
578 75 : dof_id_type last_node = end_nodes[0];
579 : dof_id_type second_last_node = last_node;
580 458 : while (last_node != end_nodes[1])
581 : {
582 383 : std::vector<dof_id_type> & curr_node_line_elems = node_to_line_elem_map[last_node];
583 : bool found_new_node = false;
584 411 : for (std::size_t i = 0; i < curr_node_line_elems.size(); ++i)
585 : {
586 411 : std::vector<dof_id_type> curr_line_elem = line_elems[curr_node_line_elems[i]];
587 488 : for (std::size_t j = 0; j < curr_line_elem.size(); ++j)
588 : {
589 460 : dof_id_type line_elem_node = curr_line_elem[j];
590 460 : if (_closed_loop &&
591 238 : (last_node == end_nodes[0] &&
592 28 : line_elem_node == end_nodes[1])) // wrong direction around closed loop
593 7 : continue;
594 453 : if (line_elem_node != last_node && line_elem_node != second_last_node)
595 : {
596 383 : _ordered_crack_front_nodes.push_back(line_elem_node);
597 : found_new_node = true;
598 383 : break;
599 : }
600 : }
601 : if (found_new_node)
602 : break;
603 411 : }
604 383 : second_last_node = last_node;
605 383 : last_node = _ordered_crack_front_nodes.back();
606 : }
607 150 : }
608 179 : }
609 :
610 : void
611 68 : CrackFrontDefinition::orderEndNodes(std::vector<dof_id_type> & end_nodes)
612 : {
613 : // Choose the node to be the first node. Do that based on undeformed coordinates for
614 : // repeatability.
615 68 : Node & node0 = _mesh.nodeRef(end_nodes[0]);
616 68 : Node & node1 = _mesh.nodeRef(end_nodes[1]);
617 :
618 : std::size_t num_positive_coor0 = 0;
619 : std::size_t num_positive_coor1 = 0;
620 : Real dist_from_origin0 = 0.0;
621 : Real dist_from_origin1 = 0.0;
622 272 : for (std::size_t i = 0; i < 3; ++i)
623 : {
624 204 : dist_from_origin0 += node0(i) * node0(i);
625 204 : dist_from_origin1 += node1(i) * node1(i);
626 204 : if (MooseUtils::absoluteFuzzyGreaterThan(node0(i), 0.0, _tol))
627 18 : ++num_positive_coor0;
628 204 : if (MooseUtils::absoluteFuzzyGreaterThan(node1(i), 0.0, _tol))
629 68 : ++num_positive_coor1;
630 : }
631 68 : dist_from_origin0 = std::sqrt(dist_from_origin0);
632 68 : dist_from_origin1 = std::sqrt(dist_from_origin1);
633 :
634 : bool switch_ends = false;
635 68 : if (num_positive_coor1 > num_positive_coor0)
636 : {
637 : switch_ends = true;
638 : }
639 : else
640 : {
641 18 : if (!MooseUtils::absoluteFuzzyEqual(dist_from_origin1, dist_from_origin0, _tol))
642 : {
643 18 : if (dist_from_origin1 < dist_from_origin0)
644 : switch_ends = true;
645 : }
646 : else
647 : {
648 0 : if (end_nodes[1] < end_nodes[0])
649 : switch_ends = true;
650 : }
651 : }
652 : if (switch_ends)
653 : {
654 68 : std::size_t tmp_node = end_nodes[1];
655 68 : end_nodes[1] = end_nodes[0];
656 68 : end_nodes[0] = tmp_node;
657 : }
658 68 : }
659 :
660 : void
661 7 : CrackFrontDefinition::pickLoopCrackEndNodes(
662 : std::vector<dof_id_type> & end_nodes,
663 : std::set<dof_id_type> & nodes,
664 : std::map<dof_id_type, std::vector<dof_id_type>> & node_to_line_elem_map,
665 : std::vector<std::vector<dof_id_type>> & line_elems)
666 : {
667 : dof_id_type max_dist_node = 0;
668 : Real min_dist = std::numeric_limits<Real>::max();
669 : Real max_dist = -std::numeric_limits<Real>::max();
670 : // Pick the node farthest from the origin as the end node, or the one with
671 : // the greatest x coordinate if the nodes are equidistant from the origin
672 217 : for (auto nit = nodes.begin(); nit != nodes.end(); ++nit)
673 : {
674 210 : Node & node = _mesh.nodeRef(*nit);
675 210 : Real dist = node.norm();
676 210 : if (dist > max_dist)
677 : {
678 : max_dist = dist;
679 14 : max_dist_node = *nit;
680 : }
681 196 : else if (dist < min_dist)
682 : min_dist = dist;
683 : }
684 :
685 : dof_id_type end_node;
686 7 : if (MooseUtils::absoluteFuzzyGreaterThan(max_dist, min_dist, _tol))
687 0 : end_node = max_dist_node;
688 : else
689 : {
690 : std::vector<Node *> node_vec;
691 217 : for (auto nit = nodes.begin(); nit != nodes.end(); ++nit)
692 210 : node_vec.push_back(_mesh.nodePtr(*nit));
693 7 : end_node = maxNodeCoor(node_vec);
694 7 : }
695 :
696 7 : end_nodes.push_back(end_node);
697 :
698 : // Find the two nodes connected to the node identified as the end node, and pick one of those to
699 : // be the other end node
700 7 : auto end_node_line_elems = node_to_line_elem_map[end_node];
701 7 : if (end_node_line_elems.size() != 2)
702 0 : mooseError(
703 : "Crack front nodes are in a loop, but crack end node is only connected to one other node");
704 : std::vector<Node *> candidate_other_end_nodes;
705 :
706 21 : for (std::size_t i = 0; i < 2; ++i)
707 : {
708 14 : auto end_line_elem = line_elems[end_node_line_elems[i]];
709 42 : for (std::size_t j = 0; j < end_line_elem.size(); ++j)
710 : {
711 28 : auto line_elem_node = end_line_elem[j];
712 28 : if (line_elem_node != end_node)
713 14 : candidate_other_end_nodes.push_back(_mesh.nodePtr(line_elem_node));
714 : }
715 14 : }
716 7 : if (candidate_other_end_nodes.size() != 2)
717 0 : mooseError(
718 : "Crack front nodes are in a loop, but crack end node is not connected to two other nodes");
719 7 : end_nodes.push_back(maxNodeCoor(candidate_other_end_nodes, 1));
720 7 : }
721 :
722 : dof_id_type
723 14 : CrackFrontDefinition::maxNodeCoor(std::vector<Node *> & nodes, unsigned int dir0)
724 : {
725 : Real dirs[3];
726 : if (dir0 == 0)
727 : {
728 : dirs[0] = 0;
729 : dirs[1] = 1;
730 : dirs[2] = 2;
731 : }
732 : else if (dir0 == 1)
733 : {
734 : dirs[0] = 1;
735 : dirs[1] = 2;
736 : dirs[2] = 0;
737 : }
738 : else if (dir0 == 2)
739 : {
740 : dirs[0] = 2;
741 : dirs[1] = 0;
742 : dirs[2] = 1;
743 : }
744 : else
745 0 : mooseError("Invalid dir0 in CrackFrontDefinition::maxNodeCoor()");
746 :
747 : Real max_coor0 = -std::numeric_limits<Real>::max();
748 : std::vector<Node *> max_coor0_nodes;
749 238 : for (std::size_t i = 0; i < nodes.size(); ++i)
750 : {
751 224 : Real coor0 = (*nodes[i])(dirs[0]);
752 224 : if (coor0 > max_coor0)
753 : max_coor0 = coor0;
754 : }
755 238 : for (std::size_t i = 0; i < nodes.size(); ++i)
756 : {
757 224 : Real coor0 = (*nodes[i])(dirs[0]);
758 224 : if (MooseUtils::absoluteFuzzyEqual(coor0, max_coor0, _tol))
759 217 : max_coor0_nodes.push_back(nodes[i]);
760 : }
761 14 : if (max_coor0_nodes.size() > 1)
762 : {
763 : Real max_coor1 = -std::numeric_limits<Real>::max();
764 : std::vector<Node *> max_coor1_nodes;
765 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
766 : {
767 210 : Real coor1 = (*nodes[i])(dirs[1]);
768 210 : if (coor1 > max_coor1)
769 : max_coor1 = coor1;
770 : }
771 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
772 : {
773 210 : Real coor1 = (*nodes[i])(dirs[1]);
774 210 : if (MooseUtils::absoluteFuzzyEqual(coor1, max_coor1, _tol))
775 14 : max_coor1_nodes.push_back(nodes[i]);
776 : }
777 7 : if (max_coor1_nodes.size() > 1)
778 : {
779 : Real max_coor2 = -std::numeric_limits<Real>::max();
780 : std::vector<Node *> max_coor2_nodes;
781 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
782 : {
783 210 : Real coor2 = (*nodes[i])(dirs[2]);
784 210 : if (coor2 > max_coor2)
785 : max_coor2 = coor2;
786 : }
787 217 : for (std::size_t i = 0; i < nodes.size(); ++i)
788 : {
789 210 : Real coor2 = (*nodes[i])(dirs[2]);
790 210 : if (MooseUtils::absoluteFuzzyEqual(coor2, max_coor2, _tol))
791 7 : max_coor2_nodes.push_back(nodes[i]);
792 : }
793 7 : if (max_coor2_nodes.size() > 1)
794 0 : mooseError("Multiple nodes with same x,y,z coordinates within tolerance");
795 : else
796 7 : return max_coor2_nodes[0]->id();
797 7 : }
798 : else
799 0 : return max_coor1_nodes[0]->id();
800 7 : }
801 : else
802 7 : return max_coor0_nodes[0]->id();
803 14 : }
804 :
805 : void
806 195 : CrackFrontDefinition::updateCrackFrontGeometry()
807 : {
808 195 : computeCrackMouthNodes();
809 195 : _segment_lengths.clear();
810 195 : _tangent_directions.clear();
811 195 : _crack_directions.clear();
812 195 : _overall_length = 0.0;
813 195 : _rot_matrix.clear();
814 195 : _distances_along_front.clear();
815 195 : _angles_along_front.clear();
816 195 : _strain_along_front.clear();
817 195 : _crack_plane_normals.clear();
818 :
819 195 : if (_treat_as_2d)
820 : {
821 120 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
822 120 : _segment_lengths.reserve(num_crack_front_points);
823 120 : _tangent_directions.reserve(num_crack_front_points);
824 120 : _crack_directions.reserve(num_crack_front_points);
825 120 : if (_use_mesh_cutter)
826 : _crack_plane_normals =
827 0 : _crack_front_points_provider->getCrackPlaneNormals(num_crack_front_points);
828 :
829 252 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
830 : {
831 : RealVectorValue tangent_direction;
832 : RealVectorValue crack_direction;
833 132 : tangent_direction(_axis_2d) = 1.0;
834 132 : _tangent_directions.push_back(tangent_direction);
835 132 : const Point * crack_front_point = getCrackFrontPoint(i);
836 :
837 : crack_direction =
838 132 : calculateCrackFrontDirection(*crack_front_point, tangent_direction, MIDDLE_NODE, i);
839 132 : _crack_directions.push_back(crack_direction);
840 :
841 132 : RankTwoTensor rot_mat;
842 132 : rot_mat.fillRow(0, crack_direction);
843 132 : rot_mat(2, _axis_2d) = 1.0;
844 :
845 132 : if (!_use_mesh_cutter)
846 132 : _crack_plane_normals.push_back(tangent_direction.cross(crack_direction));
847 :
848 : mooseAssert(i <= _crack_plane_normals.size(), "_crack_plane_normals is the wrong size.");
849 132 : rot_mat.fillRow(1, _crack_plane_normals[i]);
850 132 : _rot_matrix.push_back(rot_mat);
851 :
852 132 : _segment_lengths.push_back(std::make_pair(0.0, 0.0));
853 132 : _distances_along_front.push_back(0.0);
854 132 : _angles_along_front.push_back(0.0);
855 : }
856 : }
857 : else
858 : {
859 : // filling crack_plane_normals from mesh_cutter
860 75 : if (_use_mesh_cutter)
861 : _crack_plane_normals =
862 0 : _crack_front_points_provider->getCrackPlaneNormals(_num_points_from_provider);
863 : // else filling crack_plane_normals from curvedCrackFront
864 75 : computeCurvedCrackFrontCrackPlaneNormals();
865 :
866 75 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
867 75 : _segment_lengths.reserve(num_crack_front_points);
868 75 : _tangent_directions.reserve(num_crack_front_points);
869 75 : _crack_directions.reserve(num_crack_front_points);
870 :
871 : RealVectorValue back_segment;
872 : Real back_segment_len = 0.0;
873 75 : if (_closed_loop)
874 : {
875 7 : back_segment = *getCrackFrontPoint(0) - *getCrackFrontPoint(num_crack_front_points - 1);
876 7 : back_segment_len = back_segment.norm();
877 : }
878 :
879 561 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
880 : {
881 : CRACK_NODE_TYPE ntype;
882 486 : if (_closed_loop)
883 : ntype = MIDDLE_NODE;
884 276 : else if (i == 0)
885 : ntype = END_1_NODE;
886 208 : else if (i == num_crack_front_points - 1)
887 : ntype = END_2_NODE;
888 : else
889 : ntype = MIDDLE_NODE;
890 :
891 : RealVectorValue forward_segment;
892 : Real forward_segment_len;
893 486 : if (ntype == END_2_NODE)
894 : forward_segment_len = 0.0;
895 418 : else if (_closed_loop && i == num_crack_front_points - 1)
896 : {
897 7 : forward_segment = *getCrackFrontPoint(0) - *getCrackFrontPoint(i);
898 7 : forward_segment_len = forward_segment.norm();
899 : }
900 : else
901 : {
902 411 : forward_segment = *getCrackFrontPoint(i + 1) - *getCrackFrontPoint(i);
903 411 : forward_segment_len = forward_segment.norm();
904 411 : _overall_length += forward_segment_len;
905 : }
906 :
907 486 : _segment_lengths.push_back(std::make_pair(back_segment_len, forward_segment_len));
908 486 : if (i == 0)
909 75 : _distances_along_front.push_back(0.0);
910 : else
911 411 : _distances_along_front.push_back(back_segment_len + _distances_along_front[i - 1]);
912 :
913 : RealVectorValue tangent_direction = back_segment + forward_segment;
914 486 : tangent_direction = tangent_direction / tangent_direction.norm();
915 :
916 : // If end tangent directions are given, correct the tangent at the end nodes
917 486 : if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT &&
918 : _end_direction_method == END_DIRECTION_METHOD::END_CRACK_TANGENT_VECTOR)
919 : {
920 0 : if (ntype == END_1_NODE)
921 0 : tangent_direction = _crack_tangent_vector_end_1;
922 0 : else if (ntype == END_2_NODE)
923 0 : tangent_direction = _crack_tangent_vector_end_2;
924 : }
925 :
926 486 : _tangent_directions.push_back(tangent_direction);
927 : _crack_directions.push_back(
928 486 : calculateCrackFrontDirection(*getCrackFrontPoint(i), tangent_direction, ntype, i));
929 :
930 : // correct tangent direction in the case of _use_mesh_cutter
931 486 : if (_use_mesh_cutter)
932 0 : _tangent_directions[i] = _crack_plane_normals[i].cross(_crack_directions[i]);
933 :
934 : // If the end directions are given by the user, correct also the tangent at the end nodes
935 486 : if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT &&
936 96 : _end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR &&
937 96 : (ntype == END_1_NODE || ntype == END_2_NODE))
938 : {
939 32 : _tangent_directions[i] = _crack_plane_normals[i].cross(_crack_directions[i]);
940 : }
941 :
942 486 : back_segment = forward_segment;
943 : back_segment_len = forward_segment_len;
944 : }
945 :
946 : // For CURVED_CRACK_FRONT, _crack_plane_normals gets computed in
947 : // computeCurvedCrackFrontCrackPlaneNormals
948 75 : if (_direction_method != DIRECTION_METHOD::CURVED_CRACK_FRONT && !_use_mesh_cutter)
949 : {
950 52 : std::size_t mid_id = (num_crack_front_points - 1) / 2;
951 52 : _crack_plane_normals.assign(num_crack_front_points,
952 52 : _tangent_directions[mid_id].cross(_crack_directions[mid_id]));
953 :
954 : // Make sure the normal vector is non-zero
955 : RealVectorValue zero_vec(0.0);
956 52 : if (_crack_plane_normals.front().absolute_fuzzy_equals(zero_vec, _tol))
957 0 : mooseError("Crack plane normal vector evaluates to zero");
958 : }
959 :
960 : // Calculate angles of each point along the crack front for an elliptical crack projected
961 : // to a circle.
962 75 : if (hasAngleAlongFront())
963 : {
964 14 : RealVectorValue origin_to_first_node = *getCrackFrontPoint(0) - _crack_mouth_coordinates;
965 14 : Real hyp = origin_to_first_node.norm();
966 : RealVectorValue norm_origin_to_first_node = origin_to_first_node / hyp;
967 : RealVectorValue tangent_to_first_node =
968 14 : -norm_origin_to_first_node.cross(_crack_plane_normals.front());
969 14 : tangent_to_first_node /= tangent_to_first_node.norm();
970 :
971 80 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
972 : {
973 66 : RealVectorValue origin_to_curr_node = *getCrackFrontPoint(i) - _crack_mouth_coordinates;
974 :
975 : Real adj = origin_to_curr_node * norm_origin_to_first_node;
976 : Real opp = origin_to_curr_node * tangent_to_first_node;
977 :
978 66 : Real angle = acos(adj / hyp) * 180.0 / libMesh::pi;
979 66 : if (opp < 0.0)
980 56 : angle = 360.0 - angle;
981 66 : _angles_along_front.push_back(angle);
982 : }
983 :
984 : // Correct angle on end nodes if they are 0 or 360 to be consistent with neighboring node
985 14 : if (num_crack_front_points > 1)
986 : {
987 14 : if (MooseUtils::absoluteFuzzyEqual(_angles_along_front[0], 0.0, _tol) &&
988 10 : _angles_along_front[1] > 180.0)
989 10 : _angles_along_front[0] = 360.0;
990 4 : else if (MooseUtils::absoluteFuzzyEqual(_angles_along_front[0], 360.0, _tol) &&
991 0 : _angles_along_front[1] < 180.0)
992 0 : _angles_along_front[0] = 0.0;
993 :
994 : if (MooseUtils::absoluteFuzzyEqual(
995 14 : _angles_along_front[num_crack_front_points - 1], 0.0, _tol) &&
996 0 : _angles_along_front[num_crack_front_points - 2] > 180.0)
997 0 : _angles_along_front[num_crack_front_points - 1] = 360.0;
998 : else if (MooseUtils::absoluteFuzzyEqual(
999 14 : _angles_along_front[num_crack_front_points - 1], 360.0, _tol) &&
1000 0 : _angles_along_front[num_crack_front_points - 2] < 180.0)
1001 0 : _angles_along_front[num_crack_front_points - 1] = 0.0;
1002 : }
1003 : }
1004 : else
1005 61 : _angles_along_front.resize(num_crack_front_points, 0.0);
1006 :
1007 : // Create rotation matrix
1008 561 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
1009 : {
1010 486 : RankTwoTensor rot_mat;
1011 486 : rot_mat.fillRow(0, _crack_directions[i]);
1012 486 : rot_mat.fillRow(1, _crack_plane_normals[i]);
1013 486 : rot_mat.fillRow(2, _tangent_directions[i]);
1014 486 : _rot_matrix.push_back(rot_mat);
1015 : }
1016 :
1017 75 : _console << "Summary of crack front geometry (used for fracture integrals):" << std::endl;
1018 : _console << "index node id x coord y coord z coord x dir y dir "
1019 75 : " z dir angle position seg length"
1020 75 : << std::endl;
1021 561 : for (std::size_t i = 0; i < num_crack_front_points; ++i)
1022 : {
1023 : std::size_t point_id;
1024 486 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1025 434 : point_id = _ordered_crack_front_nodes[i];
1026 : else
1027 52 : point_id = i;
1028 486 : _console << std::left << std::setw(8) << i + 1 << std::setw(10) << point_id << std::setw(14)
1029 486 : << (*getCrackFrontPoint(i))(0) << std::setw(14) << (*getCrackFrontPoint(i))(1)
1030 486 : << std::setw(14) << (*getCrackFrontPoint(i))(2) << std::setw(14)
1031 486 : << _crack_directions[i](0) << std::setw(14) << _crack_directions[i](1)
1032 486 : << std::setw(14) << _crack_directions[i](2);
1033 486 : if (hasAngleAlongFront())
1034 66 : _console << std::left << std::setw(14) << _angles_along_front[i];
1035 : else
1036 420 : _console << std::left << std::setw(14) << "--";
1037 486 : _console << std::left << std::setw(14) << _distances_along_front[i] << std::setw(14)
1038 486 : << (_segment_lengths[i].first + _segment_lengths[i].second) / 2.0 << std::endl;
1039 : }
1040 75 : _console << "overall length: " << _overall_length << std::endl;
1041 : }
1042 195 : }
1043 :
1044 : void
1045 195 : CrackFrontDefinition::computeCrackMouthNodes()
1046 : {
1047 195 : if (_crack_mouth_boundary_ids.size() > 0)
1048 : {
1049 : _crack_mouth_coordinates.zero();
1050 :
1051 : std::set<Node *> crack_mouth_nodes;
1052 18 : ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
1053 2968 : for (auto nd = bnd_nodes.begin(); nd != bnd_nodes.end(); ++nd)
1054 : {
1055 2950 : const BndNode * bnode = *nd;
1056 2950 : BoundaryID boundary_id = bnode->_bnd_id;
1057 :
1058 5866 : for (std::size_t ibid = 0; ibid < _crack_mouth_boundary_ids.size(); ++ibid)
1059 : {
1060 2950 : if (boundary_id == _crack_mouth_boundary_ids[ibid])
1061 : {
1062 34 : crack_mouth_nodes.insert(bnode->_node);
1063 : break;
1064 : }
1065 : }
1066 : }
1067 :
1068 52 : for (auto nit = crack_mouth_nodes.begin(); nit != crack_mouth_nodes.end(); ++nit)
1069 : {
1070 34 : _crack_mouth_coordinates += **nit;
1071 : }
1072 18 : _crack_mouth_coordinates /= static_cast<Real>(crack_mouth_nodes.size());
1073 :
1074 18 : if (_has_symmetry_plane)
1075 10 : _crack_mouth_coordinates(_symmetry_plane) = 0.0;
1076 : }
1077 195 : }
1078 :
1079 : void
1080 75 : CrackFrontDefinition::computeCurvedCrackFrontCrackPlaneNormals()
1081 : {
1082 75 : if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT && !_use_mesh_cutter)
1083 : {
1084 23 : _crack_plane_normals.clear();
1085 :
1086 : // Get 3 nodes on crack front
1087 23 : std::size_t num_points = getNumCrackFrontPoints();
1088 23 : if (num_points < 3)
1089 : {
1090 0 : mooseError("Crack front must contain at least 3 nodes to use CurvedCrackFront option");
1091 : }
1092 : std::size_t start_id;
1093 : std::size_t mid_id;
1094 : std::size_t end_id;
1095 :
1096 23 : if (_closed_loop)
1097 : {
1098 : start_id = 0;
1099 7 : mid_id = (num_points - 1) / 3;
1100 7 : end_id = 2 * mid_id;
1101 : }
1102 : else
1103 : {
1104 : start_id = 0;
1105 16 : mid_id = (num_points - 1) / 2;
1106 : end_id = num_points - 1;
1107 : }
1108 23 : const Point * start = getCrackFrontPoint(start_id);
1109 23 : const Point * mid = getCrackFrontPoint(mid_id);
1110 23 : const Point * end = getCrackFrontPoint(end_id);
1111 :
1112 : // Create two vectors connecting them
1113 : RealVectorValue v1 = *mid - *start;
1114 : RealVectorValue v2 = *end - *mid;
1115 :
1116 : // Take cross product to get normal
1117 23 : Point crack_plane_normal = v1.cross(v2);
1118 23 : crack_plane_normal = crack_plane_normal.unit();
1119 23 : _crack_plane_normals.assign(num_points, crack_plane_normal);
1120 :
1121 : // Make sure they're not collinear
1122 : RealVectorValue zero_vec(0.0);
1123 23 : if (_crack_plane_normals.front().absolute_fuzzy_equals(zero_vec, _tol))
1124 : {
1125 0 : mooseError("Nodes on crack front are too close to being collinear");
1126 : }
1127 : }
1128 75 : }
1129 :
1130 : RealVectorValue
1131 618 : CrackFrontDefinition::calculateCrackFrontDirection(const Point & crack_front_point,
1132 : const RealVectorValue & tangent_direction,
1133 : const CRACK_NODE_TYPE ntype,
1134 : const std::size_t crack_front_point_index) const
1135 : {
1136 : RealVectorValue crack_dir;
1137 : RealVectorValue zero_vec(0.0);
1138 :
1139 : bool calc_dir = true;
1140 618 : if (_end_direction_method == END_DIRECTION_METHOD::END_CRACK_DIRECTION_VECTOR)
1141 : {
1142 150 : if (ntype == END_1_NODE)
1143 : {
1144 26 : crack_dir = _crack_direction_vector_end_1;
1145 : calc_dir = false;
1146 : }
1147 124 : else if (ntype == END_2_NODE)
1148 : {
1149 26 : crack_dir = _crack_direction_vector_end_2;
1150 : calc_dir = false;
1151 : }
1152 : }
1153 :
1154 : if (calc_dir)
1155 : {
1156 566 : if (_direction_method == DIRECTION_METHOD::CRACK_DIRECTION_VECTOR)
1157 : {
1158 242 : crack_dir = _crack_direction_vector;
1159 : }
1160 324 : else if (_direction_method == DIRECTION_METHOD::CRACK_MOUTH)
1161 : {
1162 50 : if (_crack_mouth_coordinates.absolute_fuzzy_equals(crack_front_point, _tol))
1163 : {
1164 0 : mooseError("Crack mouth too close to crack front node");
1165 : }
1166 : RealVectorValue mouth_to_front = crack_front_point - _crack_mouth_coordinates;
1167 :
1168 50 : RealVectorValue crack_plane_normal = mouth_to_front.cross(tangent_direction);
1169 50 : if (crack_plane_normal.absolute_fuzzy_equals(zero_vec, _tol))
1170 : {
1171 0 : mooseError(
1172 : "Vector from crack mouth to crack front node is collinear with crack front segment");
1173 : }
1174 :
1175 50 : crack_dir = tangent_direction.cross(crack_plane_normal);
1176 : Real dotprod = crack_dir * mouth_to_front;
1177 50 : if (dotprod < 0)
1178 : {
1179 0 : crack_dir = -crack_dir;
1180 : }
1181 : }
1182 274 : else if (_direction_method == DIRECTION_METHOD::CURVED_CRACK_FRONT)
1183 : {
1184 274 : crack_dir = tangent_direction.cross(_crack_plane_normals[crack_front_point_index]);
1185 : }
1186 : }
1187 618 : crack_dir = crack_dir.unit();
1188 :
1189 618 : return crack_dir;
1190 : }
1191 :
1192 : void
1193 0 : CrackFrontDefinition::updateNumberOfCrackFrontPoints(std::size_t num_points)
1194 : {
1195 0 : _num_points_from_provider = num_points;
1196 0 : }
1197 :
1198 : const Node *
1199 3540944 : CrackFrontDefinition::getCrackFrontNodePtr(const std::size_t node_index) const
1200 : {
1201 : mooseAssert(node_index < _ordered_crack_front_nodes.size(), "node_index out of range");
1202 3540944 : const Node * crack_front_node = _mesh.nodePtr(_ordered_crack_front_nodes[node_index]);
1203 : mooseAssert(crack_front_node != nullptr, "invalid crack front node");
1204 3540944 : return crack_front_node;
1205 : }
1206 :
1207 : const Point *
1208 3998543 : CrackFrontDefinition::getCrackFrontPoint(const std::size_t point_index) const
1209 : {
1210 3998543 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1211 : {
1212 3539882 : return getCrackFrontNodePtr(point_index);
1213 : }
1214 : else
1215 : {
1216 : mooseAssert(point_index < _crack_front_points.size(), "point_index out of range");
1217 458661 : return &_crack_front_points[point_index];
1218 : }
1219 : }
1220 :
1221 : const RealVectorValue &
1222 3657323 : CrackFrontDefinition::getCrackFrontTangent(const std::size_t point_index) const
1223 : {
1224 : mooseAssert(point_index < _tangent_directions.size(), "point_index out of range");
1225 3657323 : return _tangent_directions[point_index];
1226 : }
1227 :
1228 : const RealVectorValue &
1229 49320 : CrackFrontDefinition::getCrackFrontNormal(const std::size_t point_index) const
1230 : {
1231 : mooseAssert(point_index < _crack_plane_normals.size(), "point_index out of range");
1232 49320 : return _crack_plane_normals[point_index];
1233 : }
1234 :
1235 : Real
1236 1662929 : CrackFrontDefinition::getCrackFrontForwardSegmentLength(const std::size_t point_index) const
1237 : {
1238 1662929 : return _segment_lengths[point_index].second;
1239 : }
1240 :
1241 : Real
1242 1662929 : CrackFrontDefinition::getCrackFrontBackwardSegmentLength(const std::size_t point_index) const
1243 : {
1244 1662929 : return _segment_lengths[point_index].first;
1245 : }
1246 :
1247 : const RealVectorValue &
1248 2025360 : CrackFrontDefinition::getCrackDirection(const std::size_t point_index) const
1249 : {
1250 2025360 : return _crack_directions[point_index];
1251 : }
1252 :
1253 : Real
1254 2684 : CrackFrontDefinition::getDistanceAlongFront(const std::size_t point_index) const
1255 : {
1256 2684 : return _distances_along_front[point_index];
1257 : }
1258 :
1259 : bool
1260 750 : CrackFrontDefinition::hasAngleAlongFront() const
1261 : {
1262 750 : return (_crack_mouth_boundary_names.size() > 0);
1263 : }
1264 :
1265 : Real
1266 189 : CrackFrontDefinition::getAngleAlongFront(const std::size_t point_index) const
1267 : {
1268 189 : if (!hasAngleAlongFront())
1269 0 : paramError(
1270 : "crack_mouth_boundary",
1271 : "In CrackFrontDefinition, Requested angle along crack front, but definition of crack mouth "
1272 : "boundary using 'crack_mouth_boundary' parameter is necessary to do that.");
1273 189 : return _angles_along_front[point_index];
1274 : }
1275 :
1276 : std::size_t
1277 82698 : CrackFrontDefinition::getNumCrackFrontPoints() const
1278 : {
1279 82698 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1280 75620 : return _ordered_crack_front_nodes.size();
1281 : else
1282 7078 : return _crack_front_points.size();
1283 : }
1284 :
1285 : RealVectorValue
1286 577992 : CrackFrontDefinition::rotateToCrackFrontCoords(const RealVectorValue vector,
1287 : const std::size_t point_index) const
1288 : {
1289 577992 : return _rot_matrix[point_index] * vector;
1290 : }
1291 :
1292 : RealVectorValue
1293 0 : CrackFrontDefinition::rotateFromCrackFrontCoordsToGlobal(const RealVectorValue vector,
1294 : const std::size_t point_index) const
1295 : {
1296 0 : RealVectorValue vec = _rot_matrix[point_index].transpose() * vector;
1297 0 : return vec;
1298 : }
1299 :
1300 : RankTwoTensor
1301 242160 : CrackFrontDefinition::rotateToCrackFrontCoords(const RankTwoTensor tensor,
1302 : const std::size_t point_index) const
1303 : {
1304 242160 : RankTwoTensor tmp_tensor(tensor);
1305 242160 : tmp_tensor.rotate(_rot_matrix[point_index]);
1306 242160 : return tmp_tensor;
1307 : }
1308 :
1309 : void
1310 80720 : CrackFrontDefinition::calculateRThetaToCrackFront(const Point qp,
1311 : const std::size_t point_index,
1312 : Real & r,
1313 : Real & theta) const
1314 : {
1315 80720 : std::size_t num_points = getNumCrackFrontPoints();
1316 : Point closest_point(0.0);
1317 : RealVectorValue closest_point_to_p;
1318 :
1319 80720 : const Point * crack_front_point = getCrackFrontPoint(point_index);
1320 80720 : RealVectorValue crack_front_point_rot = rotateToCrackFrontCoords(*crack_front_point, point_index);
1321 :
1322 : RealVectorValue crack_front_edge =
1323 80720 : rotateToCrackFrontCoords(_tangent_directions[point_index], point_index);
1324 :
1325 80720 : Point p_rot = rotateToCrackFrontCoords(qp, point_index);
1326 : p_rot = p_rot - crack_front_point_rot;
1327 :
1328 80720 : if (_treat_as_2d)
1329 : {
1330 : // In 2D, the closest node is the crack tip node and the position of the crack tip node is
1331 : // (0,0,0) in the crack front coordinate system
1332 : // In case this is a 3D mesh treated as 2D, project point onto same plane as crack front node.
1333 : // Note: In the crack front coordinate system, z is always in the tangent direction to the crack
1334 : // front
1335 : p_rot(2) = closest_point(2);
1336 : closest_point_to_p = p_rot;
1337 :
1338 : // Find r, the distance between the qp and the crack front
1339 : RealVectorValue r_vec = p_rot;
1340 36368 : r = r_vec.norm();
1341 : }
1342 : else
1343 : {
1344 : // Loop over crack front points to find the one closest to the point qp
1345 : Real min_dist = std::numeric_limits<Real>::max();
1346 248576 : for (std::size_t pit = 0; pit != num_points; ++pit)
1347 : {
1348 204224 : const Point * crack_front_point = getCrackFrontPoint(pit);
1349 : RealVectorValue crack_point_to_current_point = qp - *crack_front_point;
1350 204224 : Real dist = crack_point_to_current_point.norm();
1351 :
1352 204224 : if (dist < min_dist)
1353 : {
1354 : min_dist = dist;
1355 : closest_point = *crack_front_point;
1356 : }
1357 : }
1358 :
1359 : // Rotate coordinates to crack front coordinate system
1360 44352 : closest_point = rotateToCrackFrontCoords(closest_point, point_index);
1361 : closest_point = closest_point - crack_front_point_rot;
1362 :
1363 : // Find r, the distance between the qp and the crack front
1364 : Real edge_length_sq = crack_front_edge.norm_sq();
1365 : closest_point_to_p = p_rot - closest_point;
1366 : Real perp = crack_front_edge * closest_point_to_p;
1367 44352 : Real dist_along_edge = perp / edge_length_sq;
1368 : RealVectorValue point_on_edge = closest_point + crack_front_edge * dist_along_edge;
1369 : RealVectorValue r_vec = p_rot - point_on_edge;
1370 44352 : r = r_vec.norm();
1371 : }
1372 :
1373 : // Find theta, the angle between r and the crack front plane
1374 : RealVectorValue crack_plane_normal =
1375 80720 : rotateToCrackFrontCoords(_crack_plane_normals[point_index], point_index);
1376 : Real p_to_plane_dist = std::abs(closest_point_to_p * crack_plane_normal);
1377 :
1378 : // Determine if qp is above or below the crack plane
1379 80720 : Real y_local = p_rot(1) - closest_point(1);
1380 :
1381 : // Determine if qp is in front of or behind the crack front
1382 : RealVectorValue p2(p_rot);
1383 : p2(1) = 0;
1384 : RealVectorValue p2_vec = p2 - closest_point;
1385 80720 : Real ahead = crack_front_edge(2) * p2_vec(0) - crack_front_edge(0) * p2_vec(2);
1386 :
1387 : Real x_local(0);
1388 80720 : if (ahead >= 0)
1389 : x_local = 1;
1390 : else
1391 : x_local = -1;
1392 :
1393 : // Calculate theta based on in which quadrant in the crack front coordinate
1394 : // system the qp is located
1395 80720 : if (r > 0)
1396 : {
1397 : Real theta_quadrant1(0.0);
1398 80720 : if (MooseUtils::absoluteFuzzyEqual(r, p_to_plane_dist, _tol))
1399 : theta_quadrant1 = 0.5 * libMesh::pi;
1400 80720 : else if (p_to_plane_dist > r)
1401 0 : mooseError(
1402 : "Invalid distance p_to_plane_dist in CrackFrontDefinition::calculateRThetaToCrackFront");
1403 : else
1404 80720 : theta_quadrant1 = std::asin(p_to_plane_dist / r);
1405 :
1406 80720 : if (x_local >= 0 && y_local >= 0)
1407 18536 : theta = theta_quadrant1;
1408 :
1409 62184 : else if (x_local < 0 && y_local >= 0)
1410 17832 : theta = libMesh::pi - theta_quadrant1;
1411 :
1412 17168 : else if (x_local < 0 && y_local < 0)
1413 17168 : theta = -(libMesh::pi - theta_quadrant1);
1414 :
1415 27184 : else if (x_local >= 0 && y_local < 0)
1416 27184 : theta = -theta_quadrant1;
1417 : }
1418 0 : else if (r == 0)
1419 0 : theta = 0;
1420 : else
1421 0 : mooseError("Invalid distance r in CrackFrontDefinition::calculateRThetaToCrackFront");
1422 80720 : }
1423 :
1424 : std::size_t
1425 0 : CrackFrontDefinition::calculateRThetaToCrackFront(const Point qp, Real & r, Real & theta) const
1426 : {
1427 0 : std::size_t num_points = getNumCrackFrontPoints();
1428 :
1429 : // Loop over crack front points to find the one closest to the point qp
1430 : Real min_dist = std::numeric_limits<Real>::max();
1431 : std::size_t point_index = 0;
1432 0 : for (std::size_t pit = 0; pit != num_points; ++pit)
1433 : {
1434 0 : const Point * crack_front_point = getCrackFrontPoint(pit);
1435 : RealVectorValue crack_point_to_current_point = qp - *crack_front_point;
1436 0 : Real dist = crack_point_to_current_point.norm();
1437 :
1438 0 : if (dist < min_dist)
1439 : {
1440 : min_dist = dist;
1441 : point_index = pit;
1442 : }
1443 : }
1444 :
1445 0 : calculateRThetaToCrackFront(qp, point_index, r, theta);
1446 :
1447 0 : return point_index;
1448 : }
1449 :
1450 : bool
1451 94315 : CrackFrontDefinition::isNodeOnIntersectingBoundary(const Node * const node) const
1452 : {
1453 : bool is_on_boundary = false;
1454 : mooseAssert(node, "Invalid node");
1455 : dof_id_type node_id = node->id();
1456 129440 : for (std::size_t i = 0; i < _intersecting_boundary_ids.size(); ++i)
1457 : {
1458 38204 : if (_mesh.isBoundaryNode(node_id, _intersecting_boundary_ids[i]))
1459 : {
1460 : is_on_boundary = true;
1461 : break;
1462 : }
1463 : }
1464 94315 : return is_on_boundary;
1465 : }
1466 :
1467 : bool
1468 736 : CrackFrontDefinition::isPointWithIndexOnIntersectingBoundary(const std::size_t point_index) const
1469 : {
1470 : bool is_on_boundary = false;
1471 736 : if (_geom_definition_method == CRACK_GEOM_DEFINITION::CRACK_FRONT_NODES)
1472 : {
1473 668 : const Node * crack_front_node = getCrackFrontNodePtr(point_index);
1474 668 : is_on_boundary = isNodeOnIntersectingBoundary(crack_front_node);
1475 : }
1476 : else
1477 : {
1478 : // If the intersecting boundary option is used with crack front points, the
1479 : // first and last points are assumed to be on the intersecting boundaries.
1480 68 : std::size_t num_crack_front_points = getNumCrackFrontPoints();
1481 68 : if (point_index == 0 || point_index == num_crack_front_points - 1)
1482 : is_on_boundary = true;
1483 : }
1484 736 : return is_on_boundary;
1485 : }
1486 :
1487 : void
1488 18 : CrackFrontDefinition::calculateTangentialStrainAlongFront()
1489 : {
1490 : RealVectorValue disp_current_node;
1491 : RealVectorValue disp_previous_node;
1492 : RealVectorValue disp_next_node;
1493 :
1494 : RealVectorValue forward_segment0;
1495 : RealVectorValue forward_segment1;
1496 : Real forward_segment0_len;
1497 : Real forward_segment1_len;
1498 : RealVectorValue back_segment0;
1499 : RealVectorValue back_segment1;
1500 : Real back_segment0_len;
1501 : Real back_segment1_len;
1502 :
1503 : std::size_t num_crack_front_nodes = _ordered_crack_front_nodes.size();
1504 : const Node * current_node;
1505 : const Node * previous_node;
1506 : const Node * next_node;
1507 :
1508 18 : _strain_along_front.reserve(num_crack_front_nodes);
1509 :
1510 : // In finalize(), gatherMax builds and distributes the complete strain vector on all processors
1511 : // -> reset the vector every time
1512 144 : for (std::size_t i = 0; i < num_crack_front_nodes; ++i)
1513 126 : _strain_along_front[i] = -std::numeric_limits<Real>::max();
1514 :
1515 18 : MooseVariable & disp_x_var = _subproblem.getStandardVariable(_tid, _disp_x_var_name);
1516 18 : MooseVariable & disp_y_var = _subproblem.getStandardVariable(_tid, _disp_y_var_name);
1517 18 : MooseVariable & disp_z_var = _subproblem.getStandardVariable(_tid, _disp_z_var_name);
1518 :
1519 18 : current_node = getCrackFrontNodePtr(0);
1520 18 : if (current_node->processor_id() == processor_id())
1521 : {
1522 14 : disp_current_node(0) = disp_x_var.getNodalValue(*current_node);
1523 14 : disp_current_node(1) = disp_y_var.getNodalValue(*current_node);
1524 14 : disp_current_node(2) = disp_z_var.getNodalValue(*current_node);
1525 :
1526 14 : next_node = getCrackFrontNodePtr(1);
1527 14 : disp_next_node(0) = disp_x_var.getNodalValue(*next_node);
1528 14 : disp_next_node(1) = disp_y_var.getNodalValue(*next_node);
1529 14 : disp_next_node(2) = disp_z_var.getNodalValue(*next_node);
1530 :
1531 : forward_segment0 = *next_node - *current_node;
1532 14 : forward_segment0 = (forward_segment0 * _tangent_directions[0]) * _tangent_directions[0];
1533 14 : forward_segment0_len = forward_segment0.norm();
1534 :
1535 : forward_segment1 = (*next_node + disp_next_node) - (*current_node + disp_current_node);
1536 14 : forward_segment1 = (forward_segment1 * _tangent_directions[0]) * _tangent_directions[0];
1537 14 : forward_segment1_len = forward_segment1.norm();
1538 :
1539 14 : _strain_along_front[0] = (forward_segment1_len - forward_segment0_len) / forward_segment0_len;
1540 : }
1541 :
1542 108 : for (std::size_t i = 1; i < num_crack_front_nodes - 1; ++i)
1543 : {
1544 90 : current_node = getCrackFrontNodePtr(i);
1545 90 : if (current_node->processor_id() == processor_id())
1546 : {
1547 70 : disp_current_node(0) = disp_x_var.getNodalValue(*current_node);
1548 70 : disp_current_node(1) = disp_y_var.getNodalValue(*current_node);
1549 70 : disp_current_node(2) = disp_z_var.getNodalValue(*current_node);
1550 :
1551 70 : previous_node = getCrackFrontNodePtr(i - 1);
1552 70 : disp_previous_node(0) = disp_x_var.getNodalValue(*previous_node);
1553 70 : disp_previous_node(1) = disp_y_var.getNodalValue(*previous_node);
1554 70 : disp_previous_node(2) = disp_z_var.getNodalValue(*previous_node);
1555 :
1556 70 : next_node = getCrackFrontNodePtr(i + 1);
1557 70 : disp_next_node(0) = disp_x_var.getNodalValue(*next_node);
1558 70 : disp_next_node(1) = disp_y_var.getNodalValue(*next_node);
1559 70 : disp_next_node(2) = disp_z_var.getNodalValue(*next_node);
1560 :
1561 : back_segment0 = *current_node - *previous_node;
1562 70 : back_segment0 = (back_segment0 * _tangent_directions[i]) * _tangent_directions[i];
1563 70 : back_segment0_len = back_segment0.norm();
1564 :
1565 : back_segment1 = (*current_node + disp_current_node) - (*previous_node + disp_previous_node);
1566 70 : back_segment1 = (back_segment1 * _tangent_directions[i]) * _tangent_directions[i];
1567 70 : back_segment1_len = back_segment1.norm();
1568 :
1569 : forward_segment0 = *next_node - *current_node;
1570 70 : forward_segment0 = (forward_segment0 * _tangent_directions[i]) * _tangent_directions[i];
1571 70 : forward_segment0_len = forward_segment0.norm();
1572 :
1573 : forward_segment1 = (*next_node + disp_next_node) - (*current_node + disp_current_node);
1574 70 : forward_segment1 = (forward_segment1 * _tangent_directions[i]) * _tangent_directions[i];
1575 70 : forward_segment1_len = forward_segment1.norm();
1576 :
1577 70 : _strain_along_front[i] =
1578 70 : 0.5 * ((back_segment1_len - back_segment0_len) / back_segment0_len +
1579 70 : (forward_segment1_len - forward_segment0_len) / forward_segment0_len);
1580 : }
1581 : }
1582 :
1583 18 : current_node = getCrackFrontNodePtr(num_crack_front_nodes - 1);
1584 18 : if (current_node->processor_id() == processor_id())
1585 : {
1586 14 : disp_current_node(0) = disp_x_var.getNodalValue(*current_node);
1587 14 : disp_current_node(1) = disp_y_var.getNodalValue(*current_node);
1588 14 : disp_current_node(2) = disp_z_var.getNodalValue(*current_node);
1589 :
1590 14 : previous_node = getCrackFrontNodePtr(num_crack_front_nodes - 2);
1591 14 : disp_previous_node(0) = disp_x_var.getNodalValue(*previous_node);
1592 14 : disp_previous_node(1) = disp_y_var.getNodalValue(*previous_node);
1593 14 : disp_previous_node(2) = disp_z_var.getNodalValue(*previous_node);
1594 :
1595 : back_segment0 = *current_node - *previous_node;
1596 14 : back_segment0 = (back_segment0 * _tangent_directions[num_crack_front_nodes - 1]) *
1597 : _tangent_directions[num_crack_front_nodes - 1];
1598 14 : back_segment0_len = back_segment0.norm();
1599 :
1600 : back_segment1 = (*current_node + disp_current_node) - (*previous_node + disp_previous_node);
1601 14 : back_segment1 = (back_segment1 * _tangent_directions[num_crack_front_nodes - 1]) *
1602 : _tangent_directions[num_crack_front_nodes - 1];
1603 14 : back_segment1_len = back_segment1.norm();
1604 :
1605 14 : _strain_along_front[num_crack_front_nodes - 1] =
1606 14 : (back_segment1_len - back_segment0_len) / back_segment0_len;
1607 : }
1608 18 : }
1609 :
1610 : Real
1611 210 : CrackFrontDefinition::getCrackFrontTangentialStrain(const std::size_t node_index) const
1612 : {
1613 : Real strain;
1614 210 : if (_t_stress)
1615 : {
1616 210 : strain = _strain_along_front[node_index];
1617 : mooseAssert(strain > -std::numeric_limits<Real>::max(),
1618 : "Failure in parallel communication of crack tangential strain");
1619 : }
1620 : else
1621 0 : mooseError("In CrackFrontDefinition, tangential strain not available");
1622 :
1623 210 : return strain;
1624 : }
1625 :
1626 : void
1627 20 : CrackFrontDefinition::createQFunctionRings()
1628 : {
1629 : // In the variable names, "cfn" = crack front node
1630 20 : if (_treat_as_2d && _use_mesh_cutter == false) // 2D: the q-function defines an integral domain
1631 : // that is constant along the crack front
1632 : {
1633 : std::vector<std::vector<const Elem *>> nodes_to_elem_map;
1634 12 : MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), nodes_to_elem_map);
1635 :
1636 : std::set<dof_id_type> nodes_prev_ring;
1637 12 : nodes_prev_ring.insert(_ordered_crack_front_nodes.begin(), _ordered_crack_front_nodes.end());
1638 :
1639 : std::set<dof_id_type> connected_nodes_this_cfn;
1640 12 : connected_nodes_this_cfn.insert(_ordered_crack_front_nodes.begin(),
1641 : _ordered_crack_front_nodes.end());
1642 :
1643 : std::set<dof_id_type> old_ring_nodes_this_cfn = connected_nodes_this_cfn;
1644 :
1645 : // The first ring contains only the crack front node(s)
1646 : std::pair<dof_id_type, std::size_t> node_ring_index =
1647 : std::make_pair(_ordered_crack_front_nodes[0], 1);
1648 12 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1649 : connected_nodes_this_cfn.end());
1650 :
1651 : // Build rings of nodes around the crack front node
1652 44 : for (std::size_t ring = 2; ring <= _last_ring; ++ring)
1653 : {
1654 :
1655 : // Find nodes connected to the nodes of the previous ring
1656 : std::set<dof_id_type> new_ring_nodes_this_cfn;
1657 136 : for (auto nit = old_ring_nodes_this_cfn.begin(); nit != old_ring_nodes_this_cfn.end(); ++nit)
1658 : {
1659 : std::vector<const Node *> neighbors;
1660 104 : MeshTools::find_nodal_neighbors(
1661 104 : _mesh.getMesh(), _mesh.nodeRef(*nit), nodes_to_elem_map, neighbors);
1662 488 : for (std::size_t inei = 0; inei < neighbors.size(); ++inei)
1663 : {
1664 384 : auto thisit = connected_nodes_this_cfn.find(neighbors[inei]->id());
1665 :
1666 : // Add only nodes that are not already present in any of the rings
1667 384 : if (thisit == connected_nodes_this_cfn.end())
1668 248 : new_ring_nodes_this_cfn.insert(neighbors[inei]->id());
1669 : }
1670 104 : }
1671 :
1672 : // Add new nodes to rings
1673 32 : connected_nodes_this_cfn.insert(new_ring_nodes_this_cfn.begin(),
1674 : new_ring_nodes_this_cfn.end());
1675 : old_ring_nodes_this_cfn = new_ring_nodes_this_cfn;
1676 :
1677 : std::pair<dof_id_type, std::size_t> node_ring_index =
1678 32 : std::make_pair(_ordered_crack_front_nodes[0], ring);
1679 64 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1680 : connected_nodes_this_cfn.end());
1681 : }
1682 12 : }
1683 : else // The q-function defines one integral domain around each crack front node
1684 : {
1685 : std::size_t num_crack_front_points = _ordered_crack_front_nodes.size();
1686 : std::vector<std::vector<const Elem *>> nodes_to_elem_map;
1687 8 : MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), nodes_to_elem_map);
1688 32 : for (std::size_t icfn = 0; icfn < num_crack_front_points; ++icfn)
1689 : {
1690 : std::set<dof_id_type> nodes_prev_ring;
1691 24 : nodes_prev_ring.insert(_ordered_crack_front_nodes[icfn]);
1692 :
1693 : std::set<dof_id_type> connected_nodes_prev_cfn;
1694 : std::set<dof_id_type> connected_nodes_this_cfn;
1695 : std::set<dof_id_type> connected_nodes_next_cfn;
1696 :
1697 24 : connected_nodes_this_cfn.insert(_ordered_crack_front_nodes[icfn]);
1698 :
1699 24 : if (_closed_loop && icfn == 0)
1700 : {
1701 0 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[num_crack_front_points - 1]);
1702 0 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[icfn + 1]);
1703 : }
1704 24 : else if (_closed_loop && icfn == num_crack_front_points - 1)
1705 : {
1706 0 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[icfn - 1]);
1707 0 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[0]);
1708 : }
1709 24 : else if (icfn == 0)
1710 : {
1711 8 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[icfn + 1]);
1712 : }
1713 16 : else if (icfn == num_crack_front_points - 1)
1714 : {
1715 8 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[icfn - 1]);
1716 : }
1717 : else
1718 : {
1719 8 : connected_nodes_prev_cfn.insert(_ordered_crack_front_nodes[icfn - 1]);
1720 8 : connected_nodes_next_cfn.insert(_ordered_crack_front_nodes[icfn + 1]);
1721 : }
1722 :
1723 : std::set<dof_id_type> old_ring_nodes_prev_cfn = connected_nodes_prev_cfn;
1724 : std::set<dof_id_type> old_ring_nodes_this_cfn = connected_nodes_this_cfn;
1725 : std::set<dof_id_type> old_ring_nodes_next_cfn = connected_nodes_next_cfn;
1726 :
1727 : // The first ring contains only the crack front node
1728 : std::pair<dof_id_type, std::size_t> node_ring_index =
1729 : std::make_pair(_ordered_crack_front_nodes[icfn], 1);
1730 24 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1731 : connected_nodes_this_cfn.end());
1732 :
1733 : // Build rings of nodes around the crack front node
1734 72 : for (std::size_t ring = 2; ring <= _last_ring; ++ring)
1735 : {
1736 :
1737 : // Find nodes connected to the nodes of the previous ring, but exclude nodes in rings of
1738 : // neighboring crack front nodes
1739 : std::set<dof_id_type> new_ring_nodes_this_cfn;
1740 48 : addNodesToQFunctionRing(new_ring_nodes_this_cfn,
1741 : old_ring_nodes_this_cfn,
1742 : connected_nodes_this_cfn,
1743 : connected_nodes_prev_cfn,
1744 : connected_nodes_next_cfn,
1745 : nodes_to_elem_map);
1746 :
1747 : std::set<dof_id_type> new_ring_nodes_prev_cfn;
1748 48 : addNodesToQFunctionRing(new_ring_nodes_prev_cfn,
1749 : old_ring_nodes_prev_cfn,
1750 : connected_nodes_prev_cfn,
1751 : connected_nodes_this_cfn,
1752 : connected_nodes_next_cfn,
1753 : nodes_to_elem_map);
1754 :
1755 : std::set<dof_id_type> new_ring_nodes_next_cfn;
1756 48 : addNodesToQFunctionRing(new_ring_nodes_next_cfn,
1757 : old_ring_nodes_next_cfn,
1758 : connected_nodes_next_cfn,
1759 : connected_nodes_prev_cfn,
1760 : connected_nodes_this_cfn,
1761 : nodes_to_elem_map);
1762 :
1763 : // Add new nodes to the three sets of nodes
1764 48 : connected_nodes_prev_cfn.insert(new_ring_nodes_prev_cfn.begin(),
1765 : new_ring_nodes_prev_cfn.end());
1766 48 : connected_nodes_this_cfn.insert(new_ring_nodes_this_cfn.begin(),
1767 : new_ring_nodes_this_cfn.end());
1768 48 : connected_nodes_next_cfn.insert(new_ring_nodes_next_cfn.begin(),
1769 : new_ring_nodes_next_cfn.end());
1770 : old_ring_nodes_prev_cfn = new_ring_nodes_prev_cfn;
1771 : old_ring_nodes_this_cfn = new_ring_nodes_this_cfn;
1772 : old_ring_nodes_next_cfn = new_ring_nodes_next_cfn;
1773 :
1774 : std::pair<dof_id_type, std::size_t> node_ring_index =
1775 48 : std::make_pair(_ordered_crack_front_nodes[icfn], ring);
1776 96 : _crack_front_node_to_node_map[node_ring_index].insert(connected_nodes_this_cfn.begin(),
1777 : connected_nodes_this_cfn.end());
1778 : }
1779 : }
1780 8 : }
1781 20 : }
1782 :
1783 : void
1784 144 : CrackFrontDefinition::addNodesToQFunctionRing(
1785 : std::set<dof_id_type> & nodes_new_ring,
1786 : const std::set<dof_id_type> & nodes_old_ring,
1787 : const std::set<dof_id_type> & nodes_all_rings,
1788 : const std::set<dof_id_type> & nodes_neighbor1,
1789 : const std::set<dof_id_type> & nodes_neighbor2,
1790 : std::vector<std::vector<const Elem *>> & nodes_to_elem_map)
1791 : {
1792 384 : for (auto nit = nodes_old_ring.begin(); nit != nodes_old_ring.end(); ++nit)
1793 : {
1794 : std::vector<const Node *> neighbors;
1795 240 : MeshTools::find_nodal_neighbors(
1796 240 : _mesh.getMesh(), _mesh.nodeRef(*nit), nodes_to_elem_map, neighbors);
1797 1352 : for (std::size_t inei = 0; inei < neighbors.size(); ++inei)
1798 : {
1799 1112 : auto previt = nodes_all_rings.find(neighbors[inei]->id());
1800 1112 : auto thisit = nodes_neighbor1.find(neighbors[inei]->id());
1801 1112 : auto nextit = nodes_neighbor2.find(neighbors[inei]->id());
1802 :
1803 : // Add only nodes that are not already present in any of the three sets of nodes
1804 1112 : if (previt == nodes_all_rings.end() && thisit == nodes_neighbor1.end() &&
1805 : nextit == nodes_neighbor2.end())
1806 672 : nodes_new_ring.insert(neighbors[inei]->id());
1807 : }
1808 240 : }
1809 144 : }
1810 :
1811 : bool
1812 57377 : CrackFrontDefinition::isNodeInRing(const std::size_t ring_index,
1813 : const dof_id_type connected_node_id,
1814 : const std::size_t node_index) const
1815 : {
1816 : bool is_node_in_ring = false;
1817 : std::pair<dof_id_type, std::size_t> node_ring_key =
1818 57377 : std::make_pair(_ordered_crack_front_nodes[node_index], ring_index);
1819 : auto nnmit = _crack_front_node_to_node_map.find(node_ring_key);
1820 :
1821 57377 : if (nnmit == _crack_front_node_to_node_map.end())
1822 0 : mooseError("Could not find crack front node ",
1823 : _ordered_crack_front_nodes[node_index],
1824 : " in the crack front node to q-function ring-node map for ring ",
1825 : ring_index);
1826 :
1827 : std::set<dof_id_type> q_func_nodes = nnmit->second;
1828 57377 : if (q_func_nodes.find(connected_node_id) != q_func_nodes.end())
1829 : is_node_in_ring = true;
1830 :
1831 57377 : return is_node_in_ring;
1832 : }
1833 :
1834 : Real
1835 3606784 : CrackFrontDefinition::DomainIntegralQFunction(std::size_t crack_front_point_index,
1836 : std::size_t ring_index,
1837 : const Node * const current_node) const
1838 : {
1839 : Real dist_to_crack_front;
1840 : Real dist_along_tangent;
1841 3606784 : projectToFrontAtPoint(
1842 : dist_to_crack_front, dist_along_tangent, crack_front_point_index, current_node);
1843 :
1844 : Real q = 1.0;
1845 3606784 : if (dist_to_crack_front > _j_integral_radius_inner[ring_index] &&
1846 3539282 : dist_to_crack_front < _j_integral_radius_outer[ring_index])
1847 21339 : q = (_j_integral_radius_outer[ring_index] - dist_to_crack_front) /
1848 21339 : (_j_integral_radius_outer[ring_index] - _j_integral_radius_inner[ring_index]);
1849 3585445 : else if (dist_to_crack_front >= _j_integral_radius_outer[ring_index])
1850 : q = 0.0;
1851 :
1852 21339 : if (q > 0.0)
1853 : {
1854 88809 : Real tangent_multiplier = 1.0;
1855 88809 : if (!_treat_as_2d)
1856 : {
1857 : const Real forward_segment_length =
1858 40577 : getCrackFrontForwardSegmentLength(crack_front_point_index);
1859 : const Real backward_segment_length =
1860 40577 : getCrackFrontBackwardSegmentLength(crack_front_point_index);
1861 :
1862 40577 : if (dist_along_tangent >= 0.0)
1863 : {
1864 26395 : if (forward_segment_length > 0.0)
1865 24280 : tangent_multiplier = 1.0 - dist_along_tangent / forward_segment_length;
1866 : }
1867 : else
1868 : {
1869 14182 : if (backward_segment_length > 0.0)
1870 14024 : tangent_multiplier = 1.0 + dist_along_tangent / backward_segment_length;
1871 : }
1872 : }
1873 :
1874 88809 : tangent_multiplier = std::max(tangent_multiplier, 0.0);
1875 88809 : tangent_multiplier = std::min(tangent_multiplier, 1.0);
1876 :
1877 : // Set to zero if a node is on a designated free surface and its crack front node is not.
1878 88809 : if (isNodeOnIntersectingBoundary(current_node) &&
1879 : !_is_point_on_intersecting_boundary[crack_front_point_index])
1880 2180 : tangent_multiplier = 0.0;
1881 :
1882 88809 : q *= tangent_multiplier;
1883 : }
1884 :
1885 3606784 : return q;
1886 : }
1887 :
1888 : Real
1889 46464 : CrackFrontDefinition::DomainIntegralTopologicalQFunction(std::size_t crack_front_point_index,
1890 : std::size_t ring_index,
1891 : const Node * const current_node) const
1892 : {
1893 : Real q = 0;
1894 46464 : bool is_node_in_ring = isNodeInRing(ring_index, current_node->id(), crack_front_point_index);
1895 46464 : if (is_node_in_ring)
1896 : q = 1;
1897 :
1898 46464 : return q;
1899 : }
1900 :
1901 : void
1902 3606784 : CrackFrontDefinition::projectToFrontAtPoint(Real & dist_to_front,
1903 : Real & dist_along_tangent,
1904 : std::size_t crack_front_point_index,
1905 : const Node * const current_node) const
1906 : {
1907 3606784 : const Point * crack_front_point = getCrackFrontPoint(crack_front_point_index);
1908 :
1909 3606784 : Point p = *current_node;
1910 3606784 : const RealVectorValue & crack_front_tangent = getCrackFrontTangent(crack_front_point_index);
1911 :
1912 : RealVectorValue crack_node_to_current_node = p - *crack_front_point;
1913 3606784 : dist_along_tangent = crack_node_to_current_node * crack_front_tangent;
1914 : RealVectorValue projection_point = *crack_front_point + dist_along_tangent * crack_front_tangent;
1915 : RealVectorValue axis_to_current_node = p - projection_point;
1916 3606784 : dist_to_front = axis_to_current_node.norm();
1917 3606784 : }
1918 :
1919 : void
1920 0 : CrackFrontDefinition::isCutterModified(const bool is_cutter_modified)
1921 : {
1922 0 : _is_cutter_modified = is_cutter_modified;
1923 0 : }
|
