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