Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CrackMeshCut3DUserObject.h"
11 :
12 : #include "XFEMFuncs.h"
13 : #include "MooseError.h"
14 : #include "libmesh/string_to_enum.h"
15 : #include "MooseMesh.h"
16 : #include "MooseEnum.h"
17 : #include "libmesh/face_tri3.h"
18 : #include "libmesh/edge_edge2.h"
19 : #include "libmesh/serial_mesh.h"
20 : #include "libmesh/plane.h"
21 : #include "libmesh/mesh_tools.h"
22 : #include "Function.h"
23 :
24 : registerMooseObject("XFEMApp", CrackMeshCut3DUserObject);
25 :
26 : InputParameters
27 24 : CrackMeshCut3DUserObject::validParams()
28 : {
29 24 : InputParameters params = GeometricCutUserObject::validParams();
30 48 : params.addRequiredParam<MeshFileName>(
31 : "mesh_file",
32 : "Mesh file for the XFEM geometric cut; currently only the xda type is supported");
33 48 : MooseEnum growthDirection("MAX_HOOP_STRESS FUNCTION", "FUNCTION");
34 48 : params.addParam<MooseEnum>(
35 : "growth_dir_method", growthDirection, "choose from FUNCTION, MAX_HOOP_STRESS");
36 48 : MooseEnum growthRate("FATIGUE FUNCTION", "FUNCTION");
37 48 : params.addParam<MooseEnum>("growth_rate_method", growthRate, "choose from FUNCTION, FATIGUE");
38 48 : params.addParam<FunctionName>("growth_direction_x",
39 : "Function defining x-component of crack growth direction");
40 48 : params.addParam<FunctionName>("growth_direction_y",
41 : "Function defining y-component of crack growth direction");
42 48 : params.addParam<FunctionName>("growth_direction_z",
43 : "Function defining z-component of crack growth direction");
44 :
45 48 : params.addParam<FunctionName>("growth_rate", "Function defining crack growth rate");
46 48 : params.addParam<Real>(
47 48 : "size_control", 0, "Criterion for refining elements while growing the crack");
48 48 : params.addParam<unsigned int>("n_step_growth", 0, "Number of steps for crack growth");
49 48 : params.addParam<std::vector<dof_id_type>>("crack_front_nodes",
50 : "Set of nodes to define crack front");
51 24 : params.addClassDescription("Creates a UserObject for a mesh cutter in 3D problems");
52 24 : return params;
53 24 : }
54 :
55 : // This code does not allow predefined crack growth as a function of time
56 : // all inital cracks are defined at t_start = t_end = 0
57 12 : CrackMeshCut3DUserObject::CrackMeshCut3DUserObject(const InputParameters & parameters)
58 : : GeometricCutUserObject(parameters, true),
59 12 : _mesh(_subproblem.mesh()),
60 24 : _growth_dir_method(getParam<MooseEnum>("growth_dir_method").getEnum<GrowthDirectionEnum>()),
61 24 : _growth_rate_method(getParam<MooseEnum>("growth_rate_method").getEnum<GrowthRateEnum>()),
62 12 : _n_step_growth(getParam<unsigned int>("n_step_growth")),
63 12 : _is_mesh_modified(false),
64 20 : _func_x(parameters.isParamValid("growth_direction_x") ? &getFunction("growth_direction_x")
65 : : nullptr),
66 20 : _func_y(parameters.isParamValid("growth_direction_y") ? &getFunction("growth_direction_y")
67 : : nullptr),
68 20 : _func_z(parameters.isParamValid("growth_direction_z") ? &getFunction("growth_direction_z")
69 : : nullptr),
70 48 : _func_v(parameters.isParamValid("growth_rate") ? &getFunction("growth_rate") : nullptr)
71 : {
72 12 : _grow = (_n_step_growth == 0 ? 0 : 1);
73 :
74 12 : if (_grow)
75 : {
76 24 : if (!isParamValid("size_control"))
77 0 : mooseError("Crack growth needs size control");
78 :
79 24 : _size_control = getParam<Real>("size_control");
80 :
81 12 : if (_growth_dir_method == GrowthDirectionEnum::FUNCTION &&
82 8 : (_func_x == nullptr || _func_y == nullptr || _func_z == nullptr))
83 0 : mooseError("function is not specified for the function method that defines growth direction");
84 :
85 12 : if (_growth_dir_method == GrowthDirectionEnum::FUNCTION && _func_v == nullptr)
86 0 : mooseError("function is not specified for the function method that defines growth rate");
87 :
88 12 : if (_growth_dir_method == GrowthDirectionEnum::FUNCTION && _func_v == nullptr)
89 0 : mooseError("function with a variable is not specified for the fatigue method that defines "
90 : "growth rate");
91 :
92 24 : if (isParamValid("crack_front_nodes"))
93 : {
94 24 : _tracked_crack_front_points = getParam<std::vector<dof_id_type>>("crack_front_nodes");
95 8 : _num_crack_front_points = _tracked_crack_front_points.size();
96 8 : _cfd = true;
97 : }
98 : else
99 4 : _cfd = false;
100 : }
101 :
102 12 : if ((_growth_dir_method == GrowthDirectionEnum::MAX_HOOP_STRESS ||
103 8 : _growth_rate_method == GrowthRateEnum::FATIGUE) &&
104 8 : !_cfd)
105 0 : mooseError("'crack_front_nodes' is not specified to use crack growth criteria!");
106 :
107 : // only the xda type is currently supported
108 12 : MeshFileName xfem_cut_mesh_file = getParam<MeshFileName>("mesh_file");
109 12 : _cut_mesh = std::make_unique<ReplicatedMesh>(_communicator);
110 12 : _cut_mesh->read(xfem_cut_mesh_file);
111 :
112 : // test element type; only tri3 elements are allowed
113 168 : for (const auto & cut_elem : _cut_mesh->element_ptr_range())
114 : {
115 72 : if (cut_elem->n_nodes() != _cut_elem_nnode)
116 0 : mooseError("The input cut mesh should include tri elements only!");
117 72 : if (cut_elem->dim() != _cut_elem_dim)
118 0 : mooseError("The input cut mesh should have 2D elements only!");
119 12 : }
120 12 : }
121 :
122 : void
123 12 : CrackMeshCut3DUserObject::initialSetup()
124 : {
125 12 : if (_cfd)
126 : {
127 8 : _crack_front_definition =
128 8 : &_fe_problem.getUserObject<CrackFrontDefinition>("crackFrontDefinition");
129 8 : _crack_front_points = _tracked_crack_front_points;
130 : }
131 :
132 12 : if (_grow)
133 : {
134 12 : findBoundaryNodes();
135 12 : findBoundaryEdges();
136 12 : sortBoundaryNodes();
137 : }
138 :
139 12 : if (_growth_rate_method == GrowthRateEnum::FATIGUE)
140 : {
141 4 : _dn.clear();
142 4 : _n.clear();
143 : }
144 12 : }
145 :
146 : void
147 59 : CrackMeshCut3DUserObject::initialize()
148 : {
149 59 : _is_mesh_modified = false;
150 :
151 59 : if (_grow)
152 : {
153 59 : if (_t_step == 1)
154 12 : _last_step_initialized = 1;
155 :
156 59 : _stop = 0;
157 :
158 59 : if (_t_step > 1 && _t_step != _last_step_initialized)
159 : {
160 32 : _last_step_initialized = _t_step;
161 :
162 64 : for (unsigned int i = 0; i < _n_step_growth; ++i)
163 : {
164 32 : if (_stop != 1)
165 : {
166 32 : findActiveBoundaryNodes();
167 32 : findActiveBoundaryDirection();
168 32 : _is_mesh_modified = true;
169 32 : growFront();
170 32 : sortFrontNodes();
171 32 : if (_inactive_boundary_pos.size() != 0)
172 32 : findFrontIntersection();
173 32 : refineFront();
174 32 : triangulation();
175 32 : joinBoundary();
176 : }
177 : }
178 : }
179 : }
180 59 : if (_cfd)
181 42 : _crack_front_definition->isCutterModified(_is_mesh_modified);
182 59 : }
183 :
184 : bool
185 0 : CrackMeshCut3DUserObject::cutElementByGeometry(const Elem * /*elem*/,
186 : std::vector<Xfem::CutEdge> & /*cut_edges*/,
187 : std::vector<Xfem::CutNode> & /*cut_nodes*/) const
188 : {
189 0 : mooseError("invalid method for 3D mesh cutting");
190 : return false;
191 : }
192 :
193 : bool
194 2406 : CrackMeshCut3DUserObject::cutElementByGeometry(const Elem * elem,
195 : std::vector<Xfem::CutFace> & cut_faces) const
196 : // With the crack defined by a planar mesh, this method cuts a solid element by all elements in the
197 : // planar mesh
198 : // TODO: Time evolving cuts not yet supported in 3D (hence the lack of use of the time variable)
199 : {
200 : bool elem_cut = false;
201 :
202 2406 : if (elem->dim() != _elem_dim)
203 0 : mooseError("The structural mesh to be cut by a surface mesh must be 3D!");
204 :
205 16842 : for (unsigned int i = 0; i < elem->n_sides(); ++i)
206 : {
207 : // This returns the lowest-order type of side.
208 14436 : std::unique_ptr<const Elem> curr_side = elem->side_ptr(i);
209 14436 : if (curr_side->dim() != 2)
210 0 : mooseError("In cutElementByGeometry dimension of side must be 2, but it is ",
211 0 : curr_side->dim());
212 14436 : unsigned int n_edges = curr_side->n_sides();
213 :
214 : std::vector<unsigned int> cut_edges;
215 : std::vector<Real> cut_pos;
216 :
217 72180 : for (unsigned int j = 0; j < n_edges; j++)
218 : {
219 : // This returns the lowest-order type of side.
220 57744 : std::unique_ptr<const Elem> curr_edge = curr_side->side_ptr(j);
221 57744 : if (curr_edge->type() != EDGE2)
222 0 : mooseError("In cutElementByGeometry face edge must be EDGE2, but type is: ",
223 0 : libMesh::Utility::enum_to_string(curr_edge->type()),
224 : " base element type is: ",
225 0 : libMesh::Utility::enum_to_string(elem->type()));
226 : const Node * node1 = curr_edge->node_ptr(0);
227 : const Node * node2 = curr_edge->node_ptr(1);
228 :
229 1565280 : for (const auto & cut_elem : _cut_mesh->element_ptr_range())
230 : {
231 : std::vector<Point> vertices;
232 :
233 2899584 : for (auto & node : cut_elem->node_ref_range())
234 : {
235 2174688 : Point & this_point = node;
236 2174688 : vertices.push_back(this_point);
237 : }
238 :
239 : Point intersection;
240 724896 : if (intersectWithEdge(*node1, *node2, vertices, intersection))
241 : {
242 3336 : cut_edges.push_back(j);
243 3336 : cut_pos.emplace_back(getRelativePosition(*node1, *node2, intersection));
244 : }
245 782640 : }
246 57744 : }
247 :
248 : // if two edges of an element are cut, it is considered as an element being cut
249 14436 : if (cut_edges.size() == 2)
250 : {
251 : elem_cut = true;
252 : Xfem::CutFace mycut;
253 1578 : mycut._face_id = i;
254 1578 : mycut._face_edge.push_back(cut_edges[0]);
255 1578 : mycut._face_edge.push_back(cut_edges[1]);
256 1578 : mycut._position.push_back(cut_pos[0]);
257 1578 : mycut._position.push_back(cut_pos[1]);
258 1578 : cut_faces.push_back(mycut);
259 : }
260 14436 : }
261 2406 : return elem_cut;
262 : }
263 :
264 : bool
265 0 : CrackMeshCut3DUserObject::cutFragmentByGeometry(std::vector<std::vector<Point>> & /*frag_edges*/,
266 : std::vector<Xfem::CutEdge> & /*cut_edges*/) const
267 : {
268 0 : mooseError("invalid method for 3D mesh cutting");
269 : return false;
270 : }
271 :
272 : bool
273 0 : CrackMeshCut3DUserObject::cutFragmentByGeometry(std::vector<std::vector<Point>> & /*frag_faces*/,
274 : std::vector<Xfem::CutFace> & /*cut_faces*/) const
275 : {
276 : // TODO: Need this for branching in 3D
277 0 : mooseError("cutFragmentByGeometry not yet implemented for 3D mesh cutting");
278 : return false;
279 : }
280 :
281 : bool
282 724896 : CrackMeshCut3DUserObject::intersectWithEdge(const Point & p1,
283 : const Point & p2,
284 : const std::vector<Point> & vertices,
285 : Point & pint) const
286 : {
287 : bool has_intersection = false;
288 :
289 724896 : Plane elem_plane(vertices[0], vertices[1], vertices[2]);
290 724896 : Point point = vertices[0];
291 724896 : Point normal = elem_plane.unit_normal(point);
292 :
293 724896 : std::array<Real, 3> plane_point = {{point(0), point(1), point(2)}};
294 724896 : std::array<Real, 3> planenormal = {{normal(0), normal(1), normal(2)}};
295 724896 : std::array<Real, 3> edge_point1 = {{p1(0), p1(1), p1(2)}};
296 724896 : std::array<Real, 3> edge_point2 = {{p2(0), p2(1), p2(2)}};
297 724896 : std::array<Real, 3> cut_point = {{0.0, 0.0, 0.0}};
298 :
299 724896 : if (Xfem::plane_normal_line_exp_int_3d(
300 : &plane_point[0], &planenormal[0], &edge_point1[0], &edge_point2[0], &cut_point[0]) == 1)
301 : {
302 328752 : Point temp_p(cut_point[0], cut_point[1], cut_point[2]);
303 328752 : if (isInsideCutPlane(vertices, temp_p) && isInsideEdge(p1, p2, temp_p))
304 : {
305 3336 : pint = temp_p;
306 : has_intersection = true;
307 : }
308 : }
309 724896 : return has_intersection;
310 724896 : }
311 :
312 : bool
313 6176 : CrackMeshCut3DUserObject::findIntersection(const Point & p1,
314 : const Point & p2,
315 : const std::vector<Point> & vertices,
316 : Point & pint) const
317 : {
318 : bool has_intersection = false;
319 :
320 6176 : Plane elem_plane(vertices[0], vertices[1], vertices[2]);
321 6176 : Point point = vertices[0];
322 6176 : Point normal = elem_plane.unit_normal(point);
323 :
324 6176 : std::array<Real, 3> plane_point = {{point(0), point(1), point(2)}};
325 6176 : std::array<Real, 3> planenormal = {{normal(0), normal(1), normal(2)}};
326 6176 : std::array<Real, 3> p_begin = {{p1(0), p1(1), p1(2)}};
327 6176 : std::array<Real, 3> p_end = {{p2(0), p2(1), p2(2)}};
328 6176 : std::array<Real, 3> cut_point = {{0.0, 0.0, 0.0}};
329 :
330 6176 : if (Xfem::plane_normal_line_exp_int_3d(
331 : &plane_point[0], &planenormal[0], &p_begin[0], &p_end[0], &cut_point[0]) == 1)
332 : {
333 5024 : Point p(cut_point[0], cut_point[1], cut_point[2]);
334 5024 : Real dotp = ((p - p1) * (p2 - p1)) / ((p2 - p1) * (p2 - p1));
335 5024 : if (isInsideCutPlane(vertices, p) && dotp > 1)
336 : {
337 64 : pint = p;
338 : has_intersection = true;
339 : }
340 : }
341 6176 : return has_intersection;
342 6176 : }
343 :
344 : bool
345 11688 : CrackMeshCut3DUserObject::isInsideEdge(const Point & p1, const Point & p2, const Point & p) const
346 : {
347 : Real dotp1 = (p1 - p) * (p2 - p1);
348 : Real dotp2 = (p2 - p) * (p2 - p1);
349 11688 : return (dotp1 * dotp2 <= 0.0);
350 : }
351 :
352 : Real
353 3336 : CrackMeshCut3DUserObject::getRelativePosition(const Point & p1,
354 : const Point & p2,
355 : const Point & p) const
356 : {
357 3336 : Real full_len = (p2 - p1).norm();
358 3336 : Real len_p1_p = (p - p1).norm();
359 3336 : return len_p1_p / full_len;
360 : }
361 :
362 : bool
363 333776 : CrackMeshCut3DUserObject::isInsideCutPlane(const std::vector<Point> & vertices,
364 : const Point & p) const
365 : {
366 333776 : unsigned int n_node = vertices.size();
367 :
368 333776 : Plane elem_plane(vertices[0], vertices[1], vertices[2]);
369 333776 : Point normal = elem_plane.unit_normal(vertices[0]);
370 :
371 : bool inside = false;
372 : unsigned int counter = 0;
373 :
374 1340128 : for (unsigned int i = 0; i < n_node; ++i)
375 : {
376 1006352 : unsigned int iplus1 = (i < n_node - 1 ? i + 1 : 0);
377 1006352 : Point middle2p = p - 0.5 * (vertices[i] + vertices[iplus1]);
378 : const Point side_tang = vertices[iplus1] - vertices[i];
379 1006352 : Point side_norm = side_tang.cross(normal);
380 1006352 : Xfem::normalizePoint(middle2p);
381 1006352 : Xfem::normalizePoint(side_norm);
382 1006352 : if (middle2p * side_norm <= 0.0)
383 589980 : counter += 1;
384 : }
385 333776 : if (counter == n_node)
386 : inside = true;
387 333776 : return inside;
388 333776 : }
389 :
390 : void
391 12 : CrackMeshCut3DUserObject::findBoundaryNodes()
392 : {
393 12 : auto boundary_node_ids = MeshTools::find_boundary_nodes(*_cut_mesh);
394 108 : for (auto it = boundary_node_ids.cbegin(); it != boundary_node_ids.cend(); it++)
395 : {
396 96 : dof_id_type id = *it;
397 : std::vector<dof_id_type> neighbors;
398 96 : _boundary_map[id] = neighbors;
399 96 : }
400 12 : }
401 :
402 : void
403 12 : CrackMeshCut3DUserObject::findBoundaryEdges()
404 : {
405 : _boundary_edges.clear();
406 :
407 : std::vector<dof_id_type> corner_elem_id;
408 : unsigned int counter = 0;
409 :
410 12 : std::vector<dof_id_type> node_id(_cut_elem_nnode);
411 12 : std::vector<bool> is_node_on_boundary(_cut_elem_nnode);
412 :
413 96 : for (const auto & cut_elem : _cut_mesh->element_ptr_range())
414 : {
415 288 : for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
416 : {
417 216 : node_id[i] = cut_elem->node_ptr(i)->id();
418 216 : is_node_on_boundary[i] = (_boundary_map.find(node_id[i]) != _boundary_map.end());
419 : }
420 :
421 72 : if (is_node_on_boundary[0] && is_node_on_boundary[1] && is_node_on_boundary[2])
422 : {
423 : // this is an element at the corner; all nodes are on the boundary but not all edges are on
424 : // the boundary
425 72 : corner_elem_id.push_back(counter);
426 : }
427 : else
428 : {
429 : // for other elements, find and store boundary edges
430 0 : for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
431 : {
432 : // if both nodes on an edge are on the boundary, it is a boundary edge.
433 0 : if (is_node_on_boundary[i] && is_node_on_boundary[(i + 1 <= 2) ? i + 1 : 0])
434 : {
435 0 : dof_id_type node1 = node_id[i];
436 0 : dof_id_type node2 = node_id[(i + 1 <= 2) ? i + 1 : 0];
437 0 : if (node1 > node2)
438 : std::swap(node1, node2);
439 :
440 : Xfem::CutEdge ce;
441 :
442 0 : if (node1 > node2)
443 : std::swap(node1, node2);
444 0 : ce._id1 = node1;
445 0 : ce._id2 = node2;
446 :
447 0 : _boundary_edges.insert(ce);
448 : }
449 : }
450 : }
451 72 : ++counter;
452 12 : }
453 :
454 : // loop over edges in corner elements
455 : // if an edge is shared by two elements, it is not an boundary edge (is_edge_inside = 1)
456 84 : for (unsigned int i = 0; i < corner_elem_id.size(); ++i)
457 : {
458 72 : auto elem_it = _cut_mesh->elements_begin();
459 :
460 252 : for (dof_id_type j = 0; j < corner_elem_id[i]; ++j)
461 180 : ++elem_it;
462 72 : Elem * cut_elem = *elem_it;
463 :
464 288 : for (unsigned int j = 0; j < _cut_elem_nnode; ++j)
465 : {
466 : bool is_edge_inside = 0;
467 :
468 : dof_id_type node1 = cut_elem->node_ptr(j)->id();
469 216 : dof_id_type node2 = cut_elem->node_ptr((j + 1 <= 2) ? j + 1 : 0)->id();
470 216 : if (node1 > node2)
471 : std::swap(node1, node2);
472 :
473 : unsigned int counter = 0;
474 1332 : for (const auto & cut_elem2 : _cut_mesh->element_ptr_range())
475 : {
476 1020 : if (counter != corner_elem_id[i])
477 : {
478 3228 : for (unsigned int k = 0; k < _cut_elem_nnode; ++k)
479 : {
480 2484 : dof_id_type node3 = cut_elem2->node_ptr(k)->id();
481 2484 : dof_id_type node4 = cut_elem2->node_ptr((k + 1 <= 2) ? k + 1 : 0)->id();
482 2484 : if (node3 > node4)
483 : std::swap(node3, node4);
484 :
485 2484 : if (node1 == node3 && node2 == node4)
486 : {
487 : is_edge_inside = 1;
488 120 : goto endloop;
489 : }
490 : }
491 : }
492 900 : ++counter;
493 216 : }
494 : endloop:
495 : if (is_edge_inside == 0)
496 : {
497 : // store boundary edges
498 : Xfem::CutEdge ce;
499 :
500 96 : if (node1 > node2)
501 : std::swap(node1, node2);
502 96 : ce._id1 = node1;
503 96 : ce._id2 = node2;
504 :
505 96 : _boundary_edges.insert(ce);
506 : }
507 : else
508 : {
509 : // this is not a boundary edge; remove it from existing edge list
510 684 : for (auto it = _boundary_edges.begin(); it != _boundary_edges.end();)
511 : {
512 564 : if ((*it)._id1 == node1 && (*it)._id2 == node2)
513 0 : it = _boundary_edges.erase(it);
514 : else
515 : ++it;
516 : }
517 : }
518 : }
519 : }
520 12 : }
521 :
522 : void
523 12 : CrackMeshCut3DUserObject::sortBoundaryNodes()
524 : {
525 12 : _boundary.clear();
526 :
527 108 : for (auto it = _boundary_edges.begin(); it != _boundary_edges.end(); ++it)
528 : {
529 96 : dof_id_type node1 = (*it)._id1;
530 96 : dof_id_type node2 = (*it)._id2;
531 :
532 : mooseAssert(_boundary_map.find(node1) != _boundary_map.end(),
533 : "_boundary_map does not have this key");
534 : mooseAssert(_boundary_map.find(node2) != _boundary_map.end(),
535 : "_boundary_map does not have this key");
536 :
537 96 : _boundary_map.find(node1)->second.push_back(node2);
538 96 : _boundary_map.find(node2)->second.push_back(node1);
539 : }
540 :
541 : auto it = _boundary_map.begin();
542 108 : while (it != _boundary_map.end())
543 : {
544 96 : if (it->second.size() != 2)
545 0 : mooseError(
546 : "Boundary nodes in the cutter mesh must have exactly two neighbors; this one has: ",
547 0 : it->second.size());
548 : ++it;
549 : }
550 :
551 : auto it2 = _boundary_edges.begin();
552 12 : dof_id_type node1 = (*it2)._id1;
553 12 : dof_id_type node2 = (*it2)._id2;
554 12 : _boundary.push_back(node1);
555 12 : _boundary.push_back(node2);
556 :
557 96 : for (unsigned int i = 0; i < _boundary_edges.size() - 1; ++i)
558 : {
559 : mooseAssert(_boundary_map.find(node2) != _boundary_map.end(),
560 : "_boundary_map does not have this key");
561 :
562 84 : dof_id_type node3 = _boundary_map.find(node2)->second[0];
563 84 : dof_id_type node4 = _boundary_map.find(node2)->second[1];
564 :
565 84 : if (node3 == node1)
566 : {
567 48 : _boundary.push_back(node4);
568 48 : node1 = node2;
569 48 : node2 = node4;
570 : }
571 36 : else if (node4 == node1)
572 : {
573 36 : _boundary.push_back(node3);
574 36 : node1 = node2;
575 36 : node2 = node3;
576 : }
577 : else
578 0 : mooseError("Discontinuity in cutter boundary");
579 : }
580 12 : }
581 :
582 : Real
583 416 : CrackMeshCut3DUserObject::findDistance(dof_id_type node1, dof_id_type node2)
584 : {
585 416 : Node * n1 = _cut_mesh->node_ptr(node1);
586 : mooseAssert(n1 != nullptr, "Node is NULL");
587 416 : Node * n2 = _cut_mesh->node_ptr(node2);
588 : mooseAssert(n2 != nullptr, "Node is NULL");
589 416 : Real distance = (*n1 - *n2).norm();
590 416 : return distance;
591 : }
592 :
593 : void
594 0 : CrackMeshCut3DUserObject::refineBoundary()
595 : {
596 0 : std::vector<dof_id_type> new_boundary_order(_boundary.begin(), _boundary.end());
597 :
598 : mooseAssert(_boundary.size() >= 2, "Boundary should have at least two nodes");
599 :
600 0 : for (unsigned int i = _boundary.size() - 1; i >= 1; --i)
601 : {
602 0 : dof_id_type node1 = _boundary[i - 1];
603 0 : dof_id_type node2 = _boundary[i];
604 :
605 0 : Real distance = findDistance(node1, node2);
606 :
607 0 : if (distance > _size_control)
608 : {
609 0 : unsigned int n = static_cast<unsigned int>(distance / _size_control);
610 : std::array<Real, 3> x1;
611 : std::array<Real, 3> x2;
612 :
613 0 : Node * n1 = _cut_mesh->node_ptr(node1);
614 : mooseAssert(n1 != nullptr, "Node is NULL");
615 : Point & p1 = *n1;
616 0 : Node * n2 = _cut_mesh->node_ptr(node2);
617 : mooseAssert(n2 != nullptr, "Node is NULL");
618 : Point & p2 = *n2;
619 :
620 0 : for (unsigned int j = 0; j < 3; ++j)
621 : {
622 0 : x1[j] = p1(j);
623 0 : x2[j] = p2(j);
624 : }
625 :
626 0 : for (unsigned int j = 0; j < n; ++j)
627 : {
628 : Point x;
629 0 : for (unsigned int k = 0; k < 3; ++k)
630 0 : x(k) = x2[k] - (x2[k] - x1[k]) * (j + 1) / (n + 1);
631 :
632 0 : Node * this_node = Node::build(x, _cut_mesh->n_nodes()).release();
633 0 : _cut_mesh->add_node(this_node);
634 :
635 0 : dof_id_type id = _cut_mesh->n_nodes() - 1;
636 : auto it = new_boundary_order.begin();
637 0 : new_boundary_order.insert(it + i, id);
638 : }
639 : }
640 : }
641 :
642 0 : _boundary = new_boundary_order;
643 : mooseAssert(_boundary.size() > 0, "Boundary should not have zero size");
644 : _boundary.pop_back();
645 0 : }
646 :
647 : void
648 49 : CrackMeshCut3DUserObject::findActiveBoundaryNodes()
649 : {
650 49 : _active_boundary.clear();
651 49 : _inactive_boundary_pos.clear();
652 :
653 49 : std::unique_ptr<PointLocatorBase> pl = _mesh.getPointLocator();
654 49 : pl->enable_out_of_mesh_mode();
655 :
656 49 : unsigned int n_boundary = _boundary.size();
657 :
658 : // if the node is outside of the structural model, store its position in _boundary to
659 : // _inactive_boundary_pos
660 596 : for (unsigned int j = 0; j < n_boundary; ++j)
661 : {
662 547 : Node * this_node = _cut_mesh->node_ptr(_boundary[j]);
663 : mooseAssert(this_node, "Node is NULL");
664 547 : Point & this_point = *this_node;
665 :
666 547 : const Elem * elem = (*pl)(this_point);
667 547 : if (elem == nullptr)
668 449 : _inactive_boundary_pos.push_back(j);
669 : }
670 :
671 49 : unsigned int n_inactive_boundary = _inactive_boundary_pos.size();
672 :
673 : // all nodes are inactive, stop
674 49 : if (n_inactive_boundary == n_boundary)
675 0 : _stop = 1;
676 :
677 : // find and store active boundary segments in "_active_boundary"
678 49 : if (n_inactive_boundary == 0)
679 0 : _active_boundary.push_back(_boundary);
680 : else
681 : {
682 449 : for (unsigned int i = 0; i < n_inactive_boundary - 1; ++i)
683 : {
684 400 : if (_inactive_boundary_pos[i + 1] - _inactive_boundary_pos[i] != 1)
685 : {
686 : std::vector<dof_id_type> temp;
687 245 : for (unsigned int j = _inactive_boundary_pos[i]; j <= _inactive_boundary_pos[i + 1]; ++j)
688 : {
689 196 : temp.push_back(_boundary[j]);
690 : }
691 49 : _active_boundary.push_back(temp);
692 49 : }
693 : }
694 49 : if (_inactive_boundary_pos[n_inactive_boundary - 1] - _inactive_boundary_pos[0] <
695 49 : n_boundary - 1)
696 : {
697 : std::vector<dof_id_type> temp;
698 0 : for (unsigned int j = _inactive_boundary_pos[n_inactive_boundary - 1]; j < n_boundary; ++j)
699 0 : temp.push_back(_boundary[j]);
700 0 : for (unsigned int j = 0; j <= _inactive_boundary_pos[0]; ++j)
701 0 : temp.push_back(_boundary[j]);
702 0 : _active_boundary.push_back(temp);
703 0 : }
704 : }
705 49 : }
706 :
707 : void
708 32 : CrackMeshCut3DUserObject::findActiveBoundaryDirection()
709 : {
710 : mooseAssert(!(_cfd && _active_boundary.size() != 1),
711 : "crack-front-definition using the cutter mesh only supports one active crack front "
712 : "segment for now");
713 :
714 32 : _active_direction.clear();
715 :
716 64 : for (unsigned int i = 0; i < _active_boundary.size(); ++i)
717 : {
718 : std::vector<Point> temp;
719 : Point dir;
720 :
721 32 : if (_inactive_boundary_pos.size() != 0)
722 : {
723 128 : for (unsigned int j = 0; j < 3; ++j)
724 96 : dir(j) = 0;
725 32 : temp.push_back(dir);
726 : }
727 :
728 : unsigned int i1 = 1;
729 32 : unsigned int i2 = _active_boundary[i].size() - 1;
730 32 : if (_inactive_boundary_pos.size() == 0)
731 : {
732 : i1 = 0;
733 : i2 = _active_boundary[i].size();
734 : }
735 :
736 32 : if (_growth_dir_method == GrowthDirectionEnum::FUNCTION)
737 : // loop over active front points
738 60 : for (unsigned int j = i1; j < i2; ++j)
739 : {
740 40 : Node * this_node = _cut_mesh->node_ptr(_active_boundary[i][j]);
741 : mooseAssert(this_node, "Node is NULL");
742 40 : Point & this_point = *this_node;
743 40 : dir(0) = _func_x->value(0, this_point);
744 40 : dir(1) = _func_y->value(0, this_point);
745 40 : dir(2) = _func_z->value(0, this_point);
746 :
747 40 : temp.push_back(dir);
748 : }
749 : // determine growth direction based on KI and KII at the crack front
750 12 : else if (_growth_dir_method == GrowthDirectionEnum::MAX_HOOP_STRESS)
751 : {
752 36 : const VectorPostprocessorValue & k1 = getVectorPostprocessorValueByName("II_KI_1", "II_KI_1");
753 : const VectorPostprocessorValue & k2 =
754 36 : getVectorPostprocessorValueByName("II_KII_1", "II_KII_1");
755 : mooseAssert(k1.size() == k2.size(), "KI and KII VPPs should have the same size");
756 : mooseAssert(k1.size() == _active_boundary[0].size(),
757 : "the number of crack front nodes in the self-similar method should equal to the "
758 : "size of VPP defined at the crack front");
759 : mooseAssert(_crack_front_points.size() == _active_boundary[0].size(),
760 : "the number of crack front nodes should be the same in _crack_front_points and "
761 : "_active_boundary[0]");
762 :
763 : // the node order in _active_boundary[0] and _crack_front_points may be the same or opposite,
764 : // their correspondence is needed
765 12 : std::vector<int> index = getFrontPointsIndex();
766 :
767 36 : for (unsigned int j = i1; j < i2; ++j)
768 : {
769 24 : int ind = index[j];
770 24 : Real theta = 2 * std::atan((k1[ind] - std::sqrt(k1[ind] * k1[ind] + k2[ind] * k2[ind])) /
771 24 : (4 * k2[ind]));
772 : RealVectorValue dir_cfc; // growth direction in crack front coord (cfc) system based on the
773 : // max hoop stress criterion
774 : RealVectorValue
775 : dir; // growth direction in global coord system based on the max hoop stress criterion
776 24 : dir_cfc(0) = std::cos(theta);
777 24 : dir_cfc(1) = std::sin(theta);
778 : dir_cfc(2) = 0;
779 24 : dir = _crack_front_definition->rotateFromCrackFrontCoordsToGlobal(dir_cfc, ind);
780 :
781 24 : temp.push_back(dir);
782 : }
783 12 : }
784 : else
785 0 : mooseError("This growth_dir_method is not pre-defined!");
786 :
787 32 : if (_inactive_boundary_pos.size() != 0)
788 : {
789 128 : for (unsigned int j = 0; j < 3; ++j)
790 96 : dir(j) = 0;
791 32 : temp.push_back(dir);
792 : }
793 :
794 32 : _active_direction.push_back(temp);
795 32 : }
796 :
797 : // normalize the directional vector
798 : Real maxl = 0;
799 :
800 64 : for (unsigned int i = 0; i < _active_direction.size(); ++i)
801 160 : for (unsigned int j = 0; j < _active_direction[i].size(); ++j)
802 : {
803 128 : Point pt = _active_direction[i][j];
804 128 : Real length = std::sqrt(pt * pt);
805 128 : if (length > maxl)
806 : maxl = length;
807 : }
808 :
809 64 : for (unsigned int i = 0; i < _active_direction.size(); ++i)
810 160 : for (unsigned int j = 0; j < _active_direction[i].size(); ++j)
811 : _active_direction[i][j] /= maxl;
812 32 : }
813 :
814 : void
815 32 : CrackMeshCut3DUserObject::growFront()
816 : {
817 32 : _front.clear();
818 :
819 64 : for (unsigned int i = 0; i < _active_boundary.size(); ++i)
820 : {
821 : std::vector<dof_id_type> temp;
822 :
823 : unsigned int i1 = 1;
824 32 : unsigned int i2 = _active_boundary[i].size() - 1;
825 32 : if (_inactive_boundary_pos.size() == 0)
826 : {
827 : i1 = 0;
828 : i2 = _active_boundary[i].size();
829 : }
830 :
831 96 : for (unsigned int j = i1; j < i2; ++j)
832 : {
833 64 : Node * this_node = _cut_mesh->node_ptr(_active_boundary[i][j]);
834 : mooseAssert(this_node, "Node is NULL");
835 : Point & this_point = *this_node;
836 64 : Point dir = _active_direction[i][j];
837 :
838 : Point x;
839 :
840 64 : if (_growth_rate_method == GrowthRateEnum::FUNCTION)
841 160 : for (unsigned int k = 0; k < 3; ++k)
842 : {
843 120 : Real velo = _func_v->value(0, Point(0, 0, 0));
844 120 : x(k) = this_point(k) + dir(k) * velo;
845 : }
846 24 : else if (_growth_rate_method == GrowthRateEnum::FATIGUE)
847 : {
848 : // get the number of loading cycles for this growth increament
849 24 : if (j == i1)
850 : {
851 12 : unsigned long int dn = (unsigned long int)_func_v->value(0, Point(0, 0, 0));
852 12 : _dn.push_back(dn);
853 12 : _n.push_back(_n.size() == 0 ? dn : dn + _n[_n.size() - 1]);
854 : }
855 :
856 24 : Real growth_size = _growth_size[j];
857 96 : for (unsigned int k = 0; k < 3; ++k)
858 72 : x(k) = this_point(k) + dir(k) * growth_size;
859 : }
860 : else
861 0 : mooseError("This growth_rate_method is not pre-defined!");
862 :
863 128 : this_node = Node::build(x, _cut_mesh->n_nodes()).release();
864 64 : _cut_mesh->add_node(this_node);
865 :
866 64 : dof_id_type id = _cut_mesh->n_nodes() - 1;
867 64 : temp.push_back(id);
868 :
869 64 : if (_cfd)
870 : {
871 48 : auto it = std::find(_tracked_crack_front_points.begin(),
872 : _tracked_crack_front_points.end(),
873 : _active_boundary[0][j]);
874 48 : if (it != _tracked_crack_front_points.end())
875 : {
876 : unsigned int pos = std::distance(_tracked_crack_front_points.begin(), it);
877 48 : _tracked_crack_front_points[pos] = id;
878 : }
879 : }
880 : }
881 :
882 32 : _front.push_back(temp);
883 32 : }
884 32 : }
885 :
886 : void
887 32 : CrackMeshCut3DUserObject::sortFrontNodes()
888 : // TBD; it is not needed for current problems but will be useful for fracture growth
889 : {
890 32 : }
891 :
892 : void
893 32 : CrackMeshCut3DUserObject::findFrontIntersection()
894 : {
895 32 : ConstBndElemRange & range = *_mesh.getBoundaryElementRange();
896 :
897 64 : for (unsigned int i = 0; i < _front.size(); ++i)
898 : {
899 32 : if (_front[i].size() >= 2)
900 : {
901 : std::vector<Point> pint1;
902 : std::vector<Point> pint2;
903 : std::vector<Real> length1;
904 : std::vector<Real> length2;
905 :
906 32 : Real node_id = _front[i][0];
907 32 : Node * this_node = _cut_mesh->node_ptr(node_id);
908 : mooseAssert(this_node, "Node is NULL");
909 32 : Point & p2 = *this_node;
910 :
911 32 : if (_front[i].size() >= 4)
912 0 : node_id = _front[i][2];
913 : else
914 32 : node_id = _front[i][1];
915 :
916 32 : this_node = _cut_mesh->node_ptr(node_id);
917 : mooseAssert(this_node, "Node is NULL");
918 32 : Point & p1 = *this_node;
919 :
920 32 : node_id = _front[i].back();
921 32 : this_node = _cut_mesh->node_ptr(node_id);
922 : mooseAssert(this_node, "Node is NULL");
923 32 : Point & p4 = *this_node;
924 :
925 32 : if (_front[i].size() >= 4)
926 0 : node_id = _front[i][_front[i].size() - 3];
927 : else
928 32 : node_id = _front[i][_front[i].size() - 2];
929 :
930 32 : this_node = _cut_mesh->node_ptr(node_id);
931 : mooseAssert(this_node, "Node is NULL");
932 32 : Point & p3 = *this_node;
933 :
934 : bool do_inter1 = 1;
935 : bool do_inter2 = 1;
936 :
937 32 : std::unique_ptr<PointLocatorBase> pl = _mesh.getPointLocator();
938 32 : pl->enable_out_of_mesh_mode();
939 32 : const Elem * elem = (*pl)(p1);
940 32 : if (elem == nullptr)
941 : do_inter1 = 0;
942 32 : elem = (*pl)(p4);
943 32 : if (elem == nullptr)
944 : do_inter2 = 0;
945 :
946 3120 : for (const auto & belem : range)
947 : {
948 : Point pt;
949 : std::vector<Point> vertices;
950 :
951 3088 : elem = belem->_elem;
952 3088 : std::unique_ptr<const Elem> curr_side = elem->side_ptr(belem->_side);
953 15440 : for (unsigned int j = 0; j < curr_side->n_nodes(); ++j)
954 : {
955 : const Node * node = curr_side->node_ptr(j);
956 12352 : const Point & this_point = *node;
957 12352 : vertices.push_back(this_point);
958 : }
959 :
960 3088 : if (findIntersection(p1, p2, vertices, pt))
961 : {
962 32 : pint1.push_back(pt);
963 32 : length1.push_back((pt - p1) * (pt - p1));
964 : }
965 3088 : if (findIntersection(p3, p4, vertices, pt))
966 : {
967 32 : pint2.push_back(pt);
968 32 : length2.push_back((pt - p3) * (pt - p3));
969 : }
970 3088 : }
971 :
972 32 : if (length1.size() != 0 && do_inter1)
973 : {
974 : auto it1 = std::min_element(length1.begin(), length1.end());
975 32 : Point inter1 = pint1[std::distance(length1.begin(), it1)];
976 : inter1 += (inter1 - p1) * _const_intersection;
977 :
978 64 : Node * this_node = Node::build(inter1, _cut_mesh->n_nodes()).release();
979 32 : _cut_mesh->add_node(this_node);
980 :
981 : mooseAssert(_cut_mesh->n_nodes() - 1 > 0, "The cut mesh should have at least one element.");
982 32 : unsigned int n = _cut_mesh->n_nodes() - 1;
983 :
984 : auto it = _front[i].begin();
985 32 : _front[i].insert(it, n);
986 :
987 32 : if (_cfd)
988 24 : _tracked_crack_front_points[_tracked_crack_front_points.size() - 1] = n;
989 : }
990 :
991 32 : if (length2.size() != 0 && do_inter2)
992 : {
993 : auto it2 = std::min_element(length2.begin(), length2.end());
994 32 : Point inter2 = pint2[std::distance(length2.begin(), it2)];
995 : inter2 += (inter2 - p2) * _const_intersection;
996 :
997 64 : Node * this_node = Node::build(inter2, _cut_mesh->n_nodes()).release();
998 32 : _cut_mesh->add_node(this_node);
999 :
1000 32 : dof_id_type n = _cut_mesh->n_nodes() - 1;
1001 :
1002 : auto it = _front[i].begin();
1003 : unsigned int m = _front[i].size();
1004 32 : _front[i].insert(it + m, n);
1005 :
1006 32 : if (_cfd)
1007 24 : _tracked_crack_front_points[0] = n;
1008 : }
1009 32 : }
1010 : }
1011 32 : }
1012 :
1013 : void
1014 32 : CrackMeshCut3DUserObject::refineFront()
1015 : {
1016 32 : std::vector<std::vector<dof_id_type>> new_front(_front.begin(), _front.end());
1017 :
1018 64 : for (unsigned int ifront = 0; ifront < _front.size(); ++ifront)
1019 : {
1020 32 : unsigned int i1 = _front[ifront].size() - 1;
1021 32 : if (_inactive_boundary_pos.size() == 0)
1022 : i1 = _front[ifront].size();
1023 :
1024 128 : for (unsigned int i = i1; i >= 1; --i)
1025 : {
1026 : unsigned int i2 = i;
1027 96 : if (_inactive_boundary_pos.size() == 0)
1028 0 : i2 = (i <= _front[ifront].size() - 1 ? i : 0);
1029 :
1030 96 : dof_id_type node1 = _front[ifront][i - 1];
1031 96 : dof_id_type node2 = _front[ifront][i2];
1032 96 : Real distance = findDistance(node1, node2);
1033 :
1034 96 : if (distance > _size_control)
1035 : {
1036 0 : unsigned int n = static_cast<int>(distance / _size_control);
1037 : std::array<Real, 3> x1;
1038 : std::array<Real, 3> x2;
1039 :
1040 0 : Node * this_node = _cut_mesh->node_ptr(node1);
1041 : mooseAssert(this_node, "Node is NULL");
1042 : Point & p1 = *this_node;
1043 0 : this_node = _cut_mesh->node_ptr(node2);
1044 : mooseAssert(this_node, "Node is NULL");
1045 : Point & p2 = *this_node;
1046 :
1047 0 : for (unsigned int j = 0; j < 3; ++j)
1048 : {
1049 0 : x1[j] = p1(j);
1050 0 : x2[j] = p2(j);
1051 : }
1052 :
1053 0 : for (unsigned int j = 0; j < n; ++j)
1054 : {
1055 : Point x;
1056 0 : for (unsigned int k = 0; k < 3; ++k)
1057 0 : x(k) = x2[k] - (x2[k] - x1[k]) * (j + 1) / (n + 1);
1058 :
1059 0 : Node * this_node = Node::build(x, _cut_mesh->n_nodes()).release();
1060 0 : _cut_mesh->add_node(this_node);
1061 :
1062 0 : dof_id_type id = _cut_mesh->n_nodes() - 1;
1063 :
1064 : auto it = new_front[ifront].begin();
1065 0 : new_front[ifront].insert(it + i, id);
1066 : }
1067 : }
1068 : }
1069 : }
1070 :
1071 32 : _front = new_front;
1072 :
1073 32 : if (_cfd)
1074 : {
1075 24 : if (_front[0][0] == _tracked_crack_front_points[0] &&
1076 0 : _front[0].back() == _tracked_crack_front_points.back())
1077 0 : _crack_front_points = _front[0];
1078 24 : else if (_front[0][0] == _tracked_crack_front_points.back() &&
1079 24 : _front[0].back() == _tracked_crack_front_points[0])
1080 : {
1081 24 : _crack_front_points = _front[0];
1082 24 : std::reverse(_crack_front_points.begin(), _crack_front_points.end());
1083 : }
1084 : else
1085 0 : mooseError("the crack front and the tracked crack front definition must match in terms of "
1086 : "their end nodes");
1087 :
1088 24 : _num_crack_front_points = _crack_front_points.size();
1089 24 : _crack_front_definition->updateNumberOfCrackFrontPoints(_num_crack_front_points);
1090 : }
1091 32 : }
1092 :
1093 : void
1094 32 : CrackMeshCut3DUserObject::triangulation()
1095 : {
1096 :
1097 : mooseAssert(_active_boundary.size() == _front.size(),
1098 : "_active_boundary and _front should have the same size!");
1099 :
1100 32 : if (_inactive_boundary_pos.size() == 0)
1101 : {
1102 0 : _active_boundary[0].push_back(_active_boundary[0][0]);
1103 0 : _front[0].push_back(_front[0][0]);
1104 : }
1105 :
1106 : // loop over active segments
1107 64 : for (unsigned int k = 0; k < _front.size(); ++k)
1108 : {
1109 32 : unsigned int n1 = _active_boundary[k].size();
1110 32 : unsigned int n2 = _front[k].size();
1111 :
1112 : unsigned int i1 = 0;
1113 : unsigned int i2 = 0;
1114 :
1115 : // stop when all nodes are associated with an element
1116 224 : while (!(i1 == n1 - 1 && i2 == n2 - 1))
1117 : {
1118 : std::vector<dof_id_type> elem;
1119 :
1120 192 : dof_id_type p1 = _active_boundary[k][i1]; // node in the old front
1121 192 : dof_id_type p2 = _front[k][i2]; // node in the new front
1122 :
1123 192 : if (i1 != n1 - 1 && i2 != n2 - 1)
1124 : {
1125 160 : dof_id_type p3 = _active_boundary[k][i1 + 1]; // next node in the old front
1126 160 : dof_id_type p4 = _front[k][i2 + 1]; // next node in the new front
1127 :
1128 160 : elem.push_back(p1);
1129 160 : elem.push_back(p2);
1130 :
1131 160 : Real d1 = findDistance(p1, p4);
1132 160 : Real d2 = findDistance(p3, p2);
1133 :
1134 160 : if (d1 < d2)
1135 : {
1136 76 : elem.push_back(p4);
1137 : i2++;
1138 : }
1139 :
1140 : else
1141 : {
1142 84 : elem.push_back(p3);
1143 : i1++;
1144 : }
1145 160 : }
1146 :
1147 32 : else if (i1 == n1 - 1)
1148 : {
1149 20 : dof_id_type p4 = _front[k][i2 + 1]; // next node in the new front
1150 :
1151 20 : elem.push_back(p1);
1152 20 : elem.push_back(p2);
1153 20 : elem.push_back(p4);
1154 : i2++;
1155 : }
1156 :
1157 12 : else if (i2 == n2 - 1)
1158 : {
1159 12 : dof_id_type p3 = _active_boundary[k][i1 + 1]; // next node in the old front
1160 :
1161 12 : elem.push_back(p1);
1162 12 : elem.push_back(p2);
1163 12 : elem.push_back(p3);
1164 : i1++;
1165 : }
1166 :
1167 192 : Elem * new_elem = Elem::build(TRI3).release();
1168 :
1169 768 : for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
1170 : {
1171 : mooseAssert(_cut_mesh->node_ptr(elem[i]) != nullptr, "Node is NULL");
1172 576 : new_elem->set_node(i, _cut_mesh->node_ptr(elem[i]));
1173 : }
1174 :
1175 192 : _cut_mesh->add_elem(new_elem);
1176 192 : }
1177 : }
1178 32 : }
1179 :
1180 : void
1181 32 : CrackMeshCut3DUserObject::joinBoundary()
1182 : {
1183 32 : if (_inactive_boundary_pos.size() == 0)
1184 : {
1185 0 : _boundary = _front[0];
1186 : _boundary.pop_back();
1187 0 : return;
1188 : }
1189 :
1190 : std::vector<dof_id_type> full_front;
1191 :
1192 32 : unsigned int size1 = _active_boundary.size();
1193 :
1194 64 : for (unsigned int i = 0; i < size1; ++i)
1195 : {
1196 32 : unsigned int size2 = _active_boundary[i].size();
1197 :
1198 32 : dof_id_type bd1 = _active_boundary[i][size2 - 1];
1199 32 : dof_id_type bd2 = _active_boundary[i + 1 < size1 ? i + 1 : 0][0];
1200 :
1201 32 : full_front.insert(full_front.end(), _front[i].begin(), _front[i].end());
1202 :
1203 32 : auto it1 = std::find(_boundary.begin(), _boundary.end(), bd1);
1204 32 : unsigned int pos1 = std::distance(_boundary.begin(), it1);
1205 32 : auto it2 = std::find(_boundary.begin(), _boundary.end(), bd2);
1206 32 : unsigned int pos2 = std::distance(_boundary.begin(), it2);
1207 :
1208 32 : if (pos1 <= pos2)
1209 0 : full_front.insert(full_front.end(), _boundary.begin() + pos1, _boundary.begin() + pos2 + 1);
1210 : else
1211 : {
1212 32 : full_front.insert(full_front.end(), _boundary.begin() + pos1, _boundary.end());
1213 32 : full_front.insert(full_front.end(), _boundary.begin(), _boundary.begin() + pos2 + 1);
1214 : }
1215 : }
1216 :
1217 32 : _boundary = full_front;
1218 32 : }
1219 :
1220 : const std::vector<Point>
1221 162 : CrackMeshCut3DUserObject::getCrackFrontPoints(unsigned int number_crack_front_points) const
1222 : {
1223 162 : std::vector<Point> crack_front_points(number_crack_front_points);
1224 : // number_crack_front_points is updated via
1225 : // _crack_front_definition->updateNumberOfCrackFrontPoints(_crack_front_points.size())
1226 162 : if (number_crack_front_points != _crack_front_points.size())
1227 0 : mooseError("number_points_from_provider does not match the number of nodes given in "
1228 : "crack_front_nodes");
1229 810 : for (unsigned int i = 0; i < number_crack_front_points; ++i)
1230 : {
1231 648 : dof_id_type id = _crack_front_points[i];
1232 648 : Node * this_node = _cut_mesh->node_ptr(id);
1233 : mooseAssert(this_node, "Node is NULL");
1234 : Point & this_point = *this_node;
1235 648 : crack_front_points[i] = this_point;
1236 : }
1237 162 : return crack_front_points;
1238 0 : }
1239 :
1240 : const std::vector<RealVectorValue>
1241 162 : CrackMeshCut3DUserObject::getCrackPlaneNormals(unsigned int number_crack_front_points) const
1242 : {
1243 162 : std::vector<RealVectorValue> crack_plane_normals(number_crack_front_points);
1244 :
1245 : // build the node-to-elems map
1246 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> node_to_elems_map;
1247 : node_to_elems_map.clear();
1248 6132 : for (const auto & elem : _cut_mesh->element_ptr_range())
1249 11616 : for (auto & node : elem->node_ref_range())
1250 8874 : node_to_elems_map[node.id()].push_back(elem->id());
1251 :
1252 : // build the elem-to-normal map
1253 : std::unordered_map<dof_id_type, RealVectorValue> elem_to_normal_map;
1254 : elem_to_normal_map.clear();
1255 6132 : for (const auto & elem : _cut_mesh->element_ptr_range())
1256 : {
1257 2904 : Point & p1 = *elem->node_ptr(0);
1258 2904 : Point & p2 = *elem->node_ptr(1);
1259 2904 : Point & p3 = *elem->node_ptr(2);
1260 2904 : Plane elem_plane(p3, p2, p1); // to match the current normal of 0,0,-1;
1261 2904 : RealVectorValue normal = elem_plane.unit_normal(p1);
1262 2904 : elem_to_normal_map[elem->id()] = normal;
1263 3066 : }
1264 :
1265 : // for any front node, the normal is averaged based on the normals of all elements sharing this
1266 : // node this code may fail when the front node has no element connected to it, e.g. refinement at
1267 : // step 1 has to be disabled
1268 810 : for (unsigned int i = 0; i < number_crack_front_points; ++i)
1269 : {
1270 648 : dof_id_type id = _crack_front_points[i];
1271 648 : std::vector<dof_id_type> elems = node_to_elems_map[id];
1272 648 : unsigned int n_elem = elems.size();
1273 :
1274 : RealVectorValue normal_avr = 0;
1275 2090 : for (unsigned int j = 0; j < n_elem; ++j)
1276 1442 : normal_avr += elem_to_normal_map[elems[j]];
1277 : normal_avr = normal_avr / n_elem;
1278 648 : crack_plane_normals[i] = normal_avr;
1279 648 : }
1280 162 : return crack_plane_normals;
1281 0 : }
1282 :
1283 : std::vector<int>
1284 29 : CrackMeshCut3DUserObject::getFrontPointsIndex()
1285 : {
1286 : // Crack front definition using the cutter mesh currently only supports one active crack front
1287 : // segment
1288 : unsigned int ibnd = 0;
1289 29 : unsigned int size_this_segment = _active_boundary[ibnd].size();
1290 29 : unsigned int n_inactive_nodes = _inactive_boundary_pos.size();
1291 :
1292 29 : std::vector<int> index(size_this_segment, -1);
1293 :
1294 29 : unsigned int i1 = n_inactive_nodes == 0 ? 0 : 1;
1295 29 : unsigned int i2 = n_inactive_nodes == 0 ? size_this_segment : size_this_segment - 1;
1296 :
1297 : // loop over active front points
1298 87 : for (unsigned int j = i1; j < i2; ++j)
1299 : {
1300 58 : dof_id_type id = _active_boundary[ibnd][j];
1301 58 : auto it = std::find(_crack_front_points.begin(), _crack_front_points.end(), id);
1302 58 : index[j] = std::distance(_crack_front_points.begin(), it);
1303 : }
1304 :
1305 29 : return index;
1306 : }
1307 :
1308 : void
1309 17 : CrackMeshCut3DUserObject::setSubCriticalGrowthSize(std::vector<Real> & growth_size)
1310 : {
1311 17 : _growth_size = growth_size;
1312 17 : }
1313 :
1314 : unsigned int
1315 0 : CrackMeshCut3DUserObject::getNumberOfCrackFrontPoints() const
1316 : {
1317 0 : return _num_crack_front_points;
1318 : }
|
