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