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