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