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