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