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 : // MOOSE includes
11 : #include "MooseMeshElementConversionUtils.h"
12 : #include "MooseError.h"
13 : #include "MathUtils.h"
14 : #include "MooseMeshUtils.h"
15 :
16 : #include "libmesh/elem.h"
17 : #include "libmesh/enum_order.h"
18 : #include "libmesh/boundary_info.h"
19 : #include "libmesh/mesh_base.h"
20 : #include "libmesh/parallel.h"
21 : #include "libmesh/parallel_algebra.h"
22 : #include "libmesh/utility.h"
23 : #include "libmesh/cell_tet4.h"
24 : #include "libmesh/face_tri3.h"
25 : #include "libmesh/cell_pyramid5.h"
26 : #include "libmesh/cell_c0polyhedron.h"
27 :
28 : using namespace libMesh;
29 :
30 : namespace MooseMeshElementConversionUtils
31 : {
32 : void
33 8198 : hexElemSplitter(MeshBase & mesh,
34 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
35 : const dof_id_type elem_id,
36 : std::vector<dof_id_type> & converted_elems_ids)
37 : {
38 : // Build boundary information of the mesh
39 8198 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
40 : // Create a list of sidesets involving the element to be split
41 8198 : std::vector<std::vector<boundary_id_type>> elem_side_list;
42 8198 : elem_side_list.resize(6);
43 8198 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 6, elem_side_list);
44 :
45 8198 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
46 8198 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
47 : // Record all the element extra integers of the original quad element
48 8198 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
49 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
50 :
51 8198 : std::vector<std::vector<unsigned int>> opt_option;
52 8198 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
53 8198 : mesh.elem_ptr(elem_id)->node_ptr(1),
54 8198 : mesh.elem_ptr(elem_id)->node_ptr(2),
55 8198 : mesh.elem_ptr(elem_id)->node_ptr(3),
56 8198 : mesh.elem_ptr(elem_id)->node_ptr(4),
57 8198 : mesh.elem_ptr(elem_id)->node_ptr(5),
58 8198 : mesh.elem_ptr(elem_id)->node_ptr(6),
59 65584 : mesh.elem_ptr(elem_id)->node_ptr(7)};
60 8198 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
61 :
62 8198 : std::vector<std::vector<const Node *>> optimized_node_list;
63 8198 : hexNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
64 :
65 8198 : std::vector<Elem *> elems_Tet4;
66 57386 : for (const auto i : index_range(optimized_node_list))
67 : {
68 49188 : auto new_elem = std::make_unique<Tet4>();
69 49188 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
70 49188 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
71 49188 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
72 49188 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
73 49188 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
74 49188 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
75 49188 : converted_elems_ids.push_back(elems_Tet4.back()->id());
76 :
77 245940 : for (unsigned int j = 0; j < 4; j++)
78 : {
79 : // A hex element has 6 faces indexed from 0 to 5
80 : // a <6 value indicates that the face of the tet element corresponds to the face of the
81 : // original hex; a =6 value means the face of the tet is an interior face of the hex
82 196752 : if (rotated_tet_face_indices[i][j] < 6)
83 : {
84 111032 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
85 12656 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
86 : }
87 : }
88 49188 : }
89 :
90 : // Retain element extra integers
91 57386 : for (unsigned int i = 0; i < 6; i++)
92 49188 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
93 : {
94 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
95 : }
96 8198 : }
97 :
98 : void
99 696 : prismElemSplitter(MeshBase & mesh,
100 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
101 : const dof_id_type elem_id,
102 : std::vector<dof_id_type> & converted_elems_ids)
103 : {
104 : // Build boundary information of the mesh
105 696 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
106 : // Create a list of sidesets involving the element to be split
107 696 : std::vector<std::vector<boundary_id_type>> elem_side_list;
108 696 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 5, elem_side_list);
109 :
110 696 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
111 696 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
112 :
113 : // Record all the element extra integers of the original quad element
114 696 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
115 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
116 :
117 696 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
118 696 : mesh.elem_ptr(elem_id)->node_ptr(1),
119 696 : mesh.elem_ptr(elem_id)->node_ptr(2),
120 696 : mesh.elem_ptr(elem_id)->node_ptr(3),
121 696 : mesh.elem_ptr(elem_id)->node_ptr(4),
122 4176 : mesh.elem_ptr(elem_id)->node_ptr(5)};
123 696 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
124 696 : std::vector<std::vector<const Node *>> optimized_node_list;
125 696 : prismNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
126 :
127 696 : std::vector<Elem *> elems_Tet4;
128 2784 : for (const auto i : index_range(optimized_node_list))
129 : {
130 2088 : auto new_elem = std::make_unique<Tet4>();
131 2088 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
132 2088 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
133 2088 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
134 2088 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
135 2088 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
136 2088 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
137 2088 : converted_elems_ids.push_back(elems_Tet4.back()->id());
138 :
139 10440 : for (unsigned int j = 0; j < 4; j++)
140 : {
141 : // A prism element has 5 faces indexed from 0 to 4
142 : // a <4 value indicates that the face of the tet element corresponds to the face of the
143 : // original prism; a =5 value means the face of the tet is an interior face of the prism
144 8352 : if (rotated_tet_face_indices[i][j] < 5)
145 : {
146 6624 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
147 1056 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
148 : }
149 : }
150 2088 : }
151 :
152 : // Retain element extra integers
153 2784 : for (unsigned int i = 0; i < 3; i++)
154 2088 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
155 : {
156 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
157 : }
158 696 : }
159 :
160 : void
161 512 : pyramidElemSplitter(MeshBase & mesh,
162 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
163 : const dof_id_type elem_id,
164 : std::vector<dof_id_type> & converted_elems_ids)
165 : {
166 : // Build boundary information of the mesh
167 512 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
168 : // Create a list of sidesets involving the element to be split
169 512 : std::vector<std::vector<boundary_id_type>> elem_side_list;
170 512 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 5, elem_side_list);
171 :
172 512 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
173 512 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
174 : // Record all the element extra integers of the original quad element
175 512 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
176 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
177 :
178 512 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
179 512 : mesh.elem_ptr(elem_id)->node_ptr(1),
180 512 : mesh.elem_ptr(elem_id)->node_ptr(2),
181 512 : mesh.elem_ptr(elem_id)->node_ptr(3),
182 2560 : mesh.elem_ptr(elem_id)->node_ptr(4)};
183 512 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
184 512 : std::vector<std::vector<const Node *>> optimized_node_list;
185 512 : pyramidNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
186 :
187 512 : std::vector<Elem *> elems_Tet4;
188 1536 : for (const auto i : index_range(optimized_node_list))
189 : {
190 1024 : auto new_elem = std::make_unique<Tet4>();
191 1024 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
192 1024 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
193 1024 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
194 1024 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
195 1024 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
196 1024 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
197 1024 : converted_elems_ids.push_back(elems_Tet4.back()->id());
198 :
199 5120 : for (unsigned int j = 0; j < 4; j++)
200 : {
201 : // A pyramid element has 5 faces indexed from 0 to 4
202 : // a <4 value indicates that the face of the tet element corresponds to the face of the
203 : // original pyramid; a =5 value means the face of the tet is an interior face of the pyramid
204 4096 : if (rotated_tet_face_indices[i][j] < 5)
205 : {
206 3808 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
207 736 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
208 : }
209 : }
210 1024 : }
211 :
212 : // Retain element extra integers
213 1536 : for (unsigned int i = 0; i < 2; i++)
214 1024 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
215 : {
216 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
217 : }
218 512 : }
219 :
220 : void
221 150 : polyhedronElemSplitter(MeshBase & mesh,
222 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
223 : const dof_id_type elem_id,
224 : std::vector<dof_id_type> & converted_elems_ids)
225 : {
226 150 : auto elem = mesh.elem_ptr(elem_id);
227 : // Build boundary information of the mesh
228 150 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
229 : // Create a list of sidesets involving the element to be split
230 150 : std::vector<std::vector<boundary_id_type>> elem_side_list;
231 150 : elementBoundaryInfoCollector(bdry_side_list, elem_id, elem->n_sides(), elem_side_list);
232 :
233 150 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
234 150 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
235 :
236 : // Record all the element extra integers of the original quad element
237 150 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
238 0 : exist_extra_ids[j] = elem->get_extra_integer(j);
239 :
240 : // Split the polygon using its tetrahedralization
241 150 : const auto poly = dynamic_cast<C0Polyhedron *>(elem);
242 150 : Node * v_avg_node = nullptr;
243 : // Currently the only extra node ever use to tetrahedralize
244 150 : if (dynamic_cast<C0Polyhedron *>(elem))
245 : v_avg_node =
246 150 : mesh.add_point(elem->vertex_average(), mesh.max_node_id() + elem_id, elem->processor_id());
247 :
248 150 : std::vector<Elem *> elems_Tet4;
249 3560 : for (const auto tri_i : make_range(poly->n_subelements()))
250 : {
251 3410 : const auto tri_indices = poly->subelement(tri_i);
252 :
253 3410 : auto new_elem = std::make_unique<Tet4>();
254 17050 : for (const auto i : make_range(4))
255 13640 : if (tri_indices[i] >= 0)
256 13640 : new_elem->set_node(i, const_cast<Node *>(elem->node_ptr(tri_indices[i])));
257 : else
258 0 : new_elem->set_node(i, v_avg_node);
259 :
260 3410 : new_elem->subdomain_id() = elem->subdomain_id();
261 3410 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
262 3410 : converted_elems_ids.push_back(elems_Tet4.back()->id());
263 :
264 : // Get subelement to side map
265 3410 : const auto tet_face_indices = poly->subelement_sides_to_poly_sides(tri_i);
266 :
267 : // Add the sides of the tets to the relevant boundaries
268 17050 : for (unsigned int j = 0; j < 4; j++)
269 13640 : if (tet_face_indices[j] < int(elem->n_sides()))
270 3480 : for (const auto & side_info : elem_side_list[tet_face_indices[j]])
271 0 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
272 3410 : }
273 :
274 : // Retain element extra integers
275 3560 : for (const auto i : make_range(poly->n_subelements()))
276 3410 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
277 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
278 150 : }
279 :
280 : std::vector<unsigned int>
281 8198 : neighborNodeIndicesHEX8(unsigned int min_id_index)
282 : {
283 : const std::vector<std::vector<unsigned int>> preset_indices = {
284 81980 : {1, 3, 4}, {0, 2, 5}, {3, 1, 6}, {2, 0, 7}, {5, 7, 0}, {4, 6, 1}, {7, 5, 2}, {6, 4, 3}};
285 8198 : if (min_id_index > 7)
286 0 : mooseError("The input node index is out of range.");
287 : else
288 16396 : return preset_indices[min_id_index];
289 32792 : }
290 :
291 : void
292 8198 : hexNodesToTetNodesDeterminer(std::vector<const Node *> & hex_nodes,
293 : std::vector<std::vector<unsigned int>> & rotated_tet_face_indices,
294 : std::vector<std::vector<const Node *>> & tet_nodes_list)
295 : {
296 : // Find the node with the minimum id
297 8198 : std::vector<dof_id_type> node_ids(8);
298 73782 : for (unsigned int i = 0; i < 8; i++)
299 65584 : node_ids[i] = hex_nodes[i]->id();
300 :
301 8198 : const unsigned int min_node_id_index = std::distance(
302 8198 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
303 : // Get the index of the three neighbor nodes of the minimum node
304 : // The order is consistent with the description in nodeRotationHEX8()
305 : // Then determine the index of the second minimum node
306 8198 : const auto neighbor_node_indices = neighborNodeIndicesHEX8(min_node_id_index);
307 :
308 8198 : const auto neighbor_node_ids = {node_ids[neighbor_node_indices[0]],
309 8198 : node_ids[neighbor_node_indices[1]],
310 16396 : node_ids[neighbor_node_indices[2]]};
311 : const unsigned int sec_min_pos =
312 8198 : std::distance(std::begin(neighbor_node_ids),
313 8198 : std::min_element(std::begin(neighbor_node_ids), std::end(neighbor_node_ids)));
314 :
315 : // Rotate the node and face indices based on the identified minimum and second minimum nodes
316 : // After the rotation, we guarantee that the minimum node is the first node (Node 0)
317 : // And the second node (Node 1) has the minium global id among the three neighbor nodes of Node 0
318 : // This makes the splitting process simpler
319 8198 : std::vector<unsigned int> face_rotation;
320 8198 : std::vector<unsigned int> rotated_indices;
321 8198 : nodeRotationHEX8(min_node_id_index, sec_min_pos, face_rotation, rotated_indices);
322 8198 : std::vector<const Node *> rotated_hex_nodes;
323 73782 : for (unsigned int i = 0; i < 8; i++)
324 65584 : rotated_hex_nodes.push_back(hex_nodes[rotated_indices[i]]);
325 :
326 : // Find the selection of each face's cutting direction
327 8198 : const auto diagonal_directions = quadFaceDiagonalDirectionsHex(rotated_hex_nodes);
328 :
329 : // Based on the determined splitting directions of all the faces, determine the nodes of each
330 : // resulting TET4 elements after the splitting.
331 8198 : std::vector<std::vector<unsigned int>> tet_face_indices;
332 8198 : const auto tet_nodes_set = tetNodesForHex(diagonal_directions, tet_face_indices);
333 57386 : for (const auto & tet_face_index : tet_face_indices)
334 : {
335 49188 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
336 245940 : for (const auto & face_index : tet_face_index)
337 : {
338 196752 : if (face_index < 6)
339 98376 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
340 : else
341 98376 : rotated_tet_face_indices.back().push_back(6);
342 : }
343 : }
344 :
345 57386 : for (const auto & tet_nodes : tet_nodes_set)
346 : {
347 49188 : tet_nodes_list.push_back(std::vector<const Node *>());
348 245940 : for (const auto & tet_node : tet_nodes)
349 196752 : tet_nodes_list.back().push_back(rotated_hex_nodes[tet_node]);
350 : }
351 8198 : }
352 :
353 : std::vector<bool>
354 8198 : quadFaceDiagonalDirectionsHex(const std::vector<const Node *> & hex_nodes)
355 : {
356 : // Bottom/Top; Front/Back; Right/Left
357 : const std::vector<std::vector<unsigned int>> face_indices = {
358 65584 : {0, 1, 2, 3}, {4, 5, 6, 7}, {0, 1, 5, 4}, {2, 3, 7, 6}, {1, 2, 6, 5}, {3, 0, 4, 7}};
359 8198 : std::vector<bool> diagonal_directions;
360 57386 : for (const auto & face_index : face_indices)
361 : {
362 49188 : std::vector<const Node *> quad_nodes = {hex_nodes[face_index[0]],
363 49188 : hex_nodes[face_index[1]],
364 49188 : hex_nodes[face_index[2]],
365 196752 : hex_nodes[face_index[3]]};
366 49188 : diagonal_directions.push_back(quadFaceDiagonalDirection(quad_nodes));
367 49188 : }
368 16396 : return diagonal_directions;
369 32792 : }
370 :
371 : bool
372 64010 : quadFaceDiagonalDirection(const std::vector<const Node *> & quad_nodes)
373 : {
374 : const std::vector<dof_id_type> node_ids = {
375 192030 : quad_nodes[0]->id(), quad_nodes[1]->id(), quad_nodes[2]->id(), quad_nodes[3]->id()};
376 64010 : const unsigned int min_id_index = std::distance(
377 64010 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
378 64010 : if (min_id_index == 0 || min_id_index == 2)
379 41730 : return true;
380 : else
381 22280 : return false;
382 64010 : }
383 :
384 : std::vector<std::vector<unsigned int>>
385 8198 : tetNodesForHex(const std::vector<bool> diagonal_directions,
386 : std::vector<std::vector<unsigned int>> & tet_face_indices)
387 : {
388 : const std::vector<std::vector<bool>> possible_inputs = {{true, true, true, true, true, false},
389 : {true, true, true, true, false, false},
390 : {true, true, true, false, true, false},
391 : {true, false, true, true, true, false},
392 : {true, false, true, true, false, false},
393 : {true, false, true, false, true, false},
394 73782 : {true, false, true, false, false, false}};
395 :
396 8198 : const unsigned int input_index = std::distance(
397 : std::begin(possible_inputs),
398 8198 : std::find(std::begin(possible_inputs), std::end(possible_inputs), diagonal_directions));
399 :
400 8198 : switch (input_index)
401 : {
402 300 : case 0:
403 2100 : tet_face_indices = {
404 2100 : {0, 6, 2, 6}, {1, 6, 2, 6}, {1, 6, 5, 6}, {0, 6, 3, 4}, {6, 6, 3, 6}, {6, 4, 5, 6}};
405 2400 : return {{0, 1, 2, 6}, {0, 5, 1, 6}, {0, 4, 5, 6}, {0, 2, 3, 7}, {0, 6, 2, 7}, {0, 4, 6, 7}};
406 0 : case 1:
407 0 : tet_face_indices = {
408 0 : {0, 1, 2, 6}, {6, 6, 2, 6}, {6, 6, 5, 1}, {0, 6, 3, 4}, {6, 6, 3, 6}, {6, 4, 5, 6}};
409 0 : return {{0, 1, 2, 5}, {0, 2, 6, 5}, {0, 6, 4, 5}, {0, 2, 3, 7}, {0, 6, 2, 7}, {0, 4, 6, 7}};
410 7898 : case 2:
411 55286 : tet_face_indices = {
412 55286 : {0, 6, 2, 6}, {1, 6, 2, 6}, {1, 6, 5, 6}, {4, 6, 5, 6}, {4, 6, 3, 6}, {0, 6, 3, 6}};
413 63184 : return {{0, 1, 2, 6}, {0, 5, 1, 6}, {0, 4, 5, 6}, {0, 7, 4, 6}, {0, 3, 7, 6}, {0, 2, 3, 6}};
414 0 : case 3:
415 0 : tet_face_indices = {
416 0 : {4, 6, 5, 1}, {6, 6, 5, 6}, {6, 1, 2, 6}, {4, 0, 3, 6}, {6, 6, 3, 6}, {6, 6, 2, 0}};
417 0 : return {{0, 7, 4, 5}, {0, 6, 7, 5}, {0, 1, 6, 5}, {0, 3, 7, 2}, {0, 7, 6, 2}, {0, 6, 1, 2}};
418 0 : case 4:
419 0 : tet_face_indices = {{0, 1, 2, 6}, {0, 6, 3, 4}, {5, 4, 6, 1}, {5, 6, 3, 2}, {6, 6, 6, 6}};
420 0 : return {{0, 1, 2, 5}, {0, 2, 3, 7}, {4, 7, 5, 0}, {5, 7, 6, 2}, {0, 2, 7, 5}};
421 0 : case 5:
422 0 : tet_face_indices = {
423 0 : {4, 6, 5, 1}, {6, 6, 5, 6}, {6, 1, 2, 6}, {2, 6, 6, 0}, {3, 6, 6, 0}, {3, 6, 6, 4}};
424 0 : return {{0, 7, 4, 5}, {0, 6, 7, 5}, {0, 1, 6, 5}, {1, 6, 2, 0}, {2, 6, 3, 0}, {3, 6, 7, 0}};
425 0 : case 6:
426 0 : tet_face_indices = {
427 0 : {1, 4, 5, 6}, {6, 6, 5, 6}, {6, 6, 3, 4}, {1, 6, 2, 0}, {6, 6, 2, 6}, {6, 0, 3, 6}};
428 0 : return {{0, 4, 5, 7}, {0, 5, 6, 7}, {0, 6, 3, 7}, {0, 5, 1, 2}, {0, 6, 5, 2}, {0, 3, 6, 2}};
429 0 : default:
430 0 : mooseError("Unexpected input.");
431 : }
432 65584 : }
433 :
434 : void
435 8198 : nodeRotationHEX8(const unsigned int min_id_index,
436 : const unsigned int sec_min_pos,
437 : std::vector<unsigned int> & face_rotation,
438 : std::vector<unsigned int> & node_rotation)
439 : {
440 : // Assuming the original hex element is a cube, the vectors formed by nodes 0-1, 0-3, and 0-4 are
441 : // overlapped with the x, y, and z axes, respectively. sec_min_pos = 0 means the second minimum
442 : // node is in the x direction, sec_min_pos = 1 means the second minimum node is in the y
443 : // direction, and sec_min_pos = 2 means the second minimum node is in the z direction.
444 : const std::vector<std::vector<std::vector<unsigned int>>> preset_indices = {
445 : {{0, 1, 2, 3, 4, 5, 6, 7}, {0, 3, 7, 4, 1, 2, 6, 5}, {0, 4, 5, 1, 3, 7, 6, 2}},
446 : {{1, 0, 4, 5, 2, 3, 7, 6}, {1, 2, 3, 0, 5, 6, 7, 4}, {1, 5, 6, 2, 0, 4, 7, 3}},
447 : {{2, 3, 0, 1, 6, 7, 4, 5}, {2, 1, 5, 6, 3, 0, 4, 7}, {2, 6, 7, 3, 1, 5, 4, 0}},
448 : {{3, 2, 6, 7, 0, 1, 5, 4}, {3, 0, 1, 2, 7, 4, 5, 6}, {3, 7, 4, 0, 2, 6, 5, 1}},
449 : {{4, 5, 1, 0, 7, 6, 2, 3}, {4, 7, 6, 5, 0, 3, 2, 1}, {4, 0, 3, 7, 5, 1, 2, 6}},
450 : {{5, 4, 7, 6, 1, 0, 3, 2}, {5, 6, 2, 1, 4, 7, 3, 0}, {5, 1, 0, 4, 6, 2, 3, 7}},
451 : {{6, 7, 3, 2, 5, 4, 0, 1}, {6, 5, 4, 7, 2, 1, 0, 3}, {6, 2, 1, 5, 7, 3, 0, 4}},
452 278732 : {{7, 6, 5, 4, 3, 2, 1, 0}, {7, 4, 0, 3, 6, 5, 1, 2}, {7, 3, 2, 6, 4, 0, 1, 5}}};
453 :
454 : const std::vector<std::vector<std::vector<unsigned int>>> preset_face_indices = {
455 : {{0, 1, 2, 3, 4, 5}, {4, 0, 3, 5, 1, 2}, {1, 4, 5, 2, 0, 3}},
456 : {{1, 0, 4, 5, 2, 3}, {0, 2, 3, 4, 1, 5}, {2, 1, 5, 3, 0, 4}},
457 : {{0, 3, 4, 1, 2, 5}, {2, 0, 1, 5, 3, 4}, {3, 2, 5, 4, 0, 1}},
458 : {{3, 0, 2, 5, 4, 1}, {0, 4, 1, 2, 3, 5}, {4, 3, 5, 1, 0, 2}},
459 : {{1, 5, 2, 0, 4, 3}, {5, 4, 3, 2, 1, 0}, {4, 1, 0, 3, 5, 2}},
460 : {{5, 1, 4, 3, 2, 0}, {2, 5, 3, 0, 1, 4}, {1, 2, 0, 4, 5, 3}},
461 : {{3, 5, 4, 0, 2, 1}, {5, 2, 1, 4, 3, 0}, {2, 3, 0, 1, 5, 4}},
462 278732 : {{5, 3, 2, 1, 4, 0}, {4, 5, 1, 0, 3, 2}, {3, 4, 0, 2, 5, 1}}};
463 :
464 8198 : if (min_id_index > 7 || sec_min_pos > 2)
465 0 : mooseError("The input node index is out of range.");
466 : else
467 : {
468 : // index: new face index; value: old face index
469 8198 : face_rotation = preset_face_indices[min_id_index][sec_min_pos];
470 8198 : node_rotation = preset_indices[min_id_index][sec_min_pos];
471 : }
472 442692 : }
473 :
474 : void
475 14822 : nodeRotationPRISM6(unsigned int min_id_index,
476 : std::vector<unsigned int> & face_rotation,
477 : std::vector<unsigned int> & node_rotation)
478 : {
479 : const std::vector<std::vector<unsigned int>> preset_indices = {{0, 1, 2, 3, 4, 5},
480 : {1, 2, 0, 4, 5, 3},
481 : {2, 0, 1, 5, 3, 4},
482 : {3, 5, 4, 0, 2, 1},
483 : {4, 3, 5, 1, 0, 2},
484 118576 : {5, 4, 3, 2, 1, 0}};
485 :
486 : const std::vector<std::vector<unsigned int>> preset_face_indices = {{0, 1, 2, 3, 4},
487 : {0, 2, 3, 1, 4},
488 : {0, 3, 1, 2, 4},
489 : {4, 3, 2, 1, 0},
490 : {4, 1, 3, 2, 0},
491 118576 : {4, 2, 1, 3, 0}};
492 :
493 14822 : if (min_id_index > 5)
494 0 : mooseError("The input node index is out of range.");
495 : else
496 : {
497 : // index: new face index; value: old face index
498 14822 : face_rotation = preset_face_indices[min_id_index];
499 14822 : node_rotation = preset_indices[min_id_index];
500 : }
501 88932 : }
502 :
503 : void
504 14822 : prismNodesToTetNodesDeterminer(std::vector<const Node *> & prism_nodes,
505 : std::vector<std::vector<unsigned int>> & rotated_tet_face_indices,
506 : std::vector<std::vector<const Node *>> & tet_nodes_list)
507 : {
508 : // Find the node with the minimum id
509 14822 : std::vector<dof_id_type> node_ids(6);
510 103754 : for (unsigned int i = 0; i < 6; i++)
511 88932 : node_ids[i] = prism_nodes[i]->id();
512 :
513 14822 : const unsigned int min_node_id_index = std::distance(
514 14822 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
515 :
516 : // Rotate the node and face indices based on the identified minimum node
517 : // After the rotation, we guarantee that the minimum node is the first node (Node 0)
518 : // This makes the splitting process simpler
519 14822 : std::vector<unsigned int> face_rotation;
520 14822 : std::vector<unsigned int> rotated_indices;
521 14822 : nodeRotationPRISM6(min_node_id_index, face_rotation, rotated_indices);
522 14822 : std::vector<const Node *> rotated_prism_nodes;
523 103754 : for (unsigned int i = 0; i < 6; i++)
524 88932 : rotated_prism_nodes.push_back(prism_nodes[rotated_indices[i]]);
525 :
526 14822 : std::vector<const Node *> key_quad_nodes = {rotated_prism_nodes[1],
527 14822 : rotated_prism_nodes[2],
528 14822 : rotated_prism_nodes[5],
529 29644 : rotated_prism_nodes[4]};
530 :
531 : // Find the selection of each face's cutting direction
532 14822 : const bool diagonal_direction = quadFaceDiagonalDirection(key_quad_nodes);
533 :
534 : // Based on the determined splitting directions of all the faces, determine the nodes of each
535 : // resulting TET4 elements after the splitting.
536 14822 : std::vector<std::vector<unsigned int>> tet_face_indices;
537 14822 : const auto tet_nodes_set = tetNodesForPrism(diagonal_direction, tet_face_indices);
538 59288 : for (const auto & tet_face_index : tet_face_indices)
539 : {
540 44466 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
541 222330 : for (const auto & face_index : tet_face_index)
542 : {
543 177864 : if (face_index < 5)
544 118576 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
545 : else
546 59288 : rotated_tet_face_indices.back().push_back(5);
547 : }
548 : }
549 :
550 59288 : for (const auto & tet_nodes : tet_nodes_set)
551 : {
552 44466 : tet_nodes_list.push_back(std::vector<const Node *>());
553 222330 : for (const auto & tet_node : tet_nodes)
554 177864 : tet_nodes_list.back().push_back(rotated_prism_nodes[tet_node]);
555 : }
556 14822 : }
557 :
558 : std::vector<std::vector<unsigned int>>
559 14822 : tetNodesForPrism(const bool diagonal_direction,
560 : std::vector<std::vector<unsigned int>> & tet_face_indices)
561 : {
562 :
563 14822 : if (diagonal_direction)
564 : {
565 34552 : tet_face_indices = {{4, 3, 5, 1}, {2, 1, 5, 5}, {2, 5, 3, 0}};
566 43190 : return {{3, 5, 4, 0}, {1, 4, 5, 0}, {1, 5, 2, 0}};
567 : }
568 : else
569 : {
570 24736 : tet_face_indices = {{4, 3, 5, 1}, {2, 1, 5, 0}, {2, 5, 5, 3}};
571 30920 : return {{3, 5, 4, 0}, {1, 4, 2, 0}, {2, 4, 5, 0}};
572 : }
573 74110 : }
574 :
575 : void
576 804 : nodeRotationPYRAMID5(unsigned int min_id_index,
577 : std::vector<unsigned int> & face_rotation,
578 : std::vector<unsigned int> & node_rotation)
579 : {
580 : const std::vector<std::vector<unsigned int>> preset_indices = {
581 4824 : {0, 1, 2, 3, 4}, {1, 2, 3, 0, 4}, {2, 3, 0, 1, 4}, {3, 0, 1, 2, 4}};
582 :
583 : const std::vector<std::vector<unsigned int>> preset_face_indices = {
584 4824 : {0, 1, 2, 3, 4}, {1, 2, 3, 0, 4}, {2, 3, 0, 1, 4}, {3, 0, 1, 2, 4}};
585 :
586 804 : if (min_id_index > 3)
587 0 : mooseError("The input node index is out of range.");
588 : else
589 : {
590 : // index: new face index; value: old face index
591 804 : face_rotation = preset_face_indices[min_id_index];
592 804 : node_rotation = preset_indices[min_id_index];
593 : }
594 4824 : }
595 :
596 : void
597 804 : pyramidNodesToTetNodesDeterminer(std::vector<const Node *> & pyramid_nodes,
598 : std::vector<std::vector<unsigned int>> & rotated_tet_face_indices,
599 : std::vector<std::vector<const Node *>> & tet_nodes_list)
600 : {
601 : // Find the node with the minimum id, ignoring the top node
602 804 : std::vector<dof_id_type> node_ids(4);
603 4020 : for (unsigned int i = 0; i < 4; i++)
604 3216 : node_ids[i] = pyramid_nodes[i]->id();
605 :
606 804 : const unsigned int min_node_id_index = std::distance(
607 804 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
608 :
609 : // Rotate the node and face indices based on the identified minimum nodes
610 : // After the rotation, we guarantee that the minimum node is the first node (Node 0)
611 : // This makes the splitting process simpler
612 804 : std::vector<unsigned int> face_rotation;
613 804 : std::vector<unsigned int> rotated_indices;
614 804 : nodeRotationPYRAMID5(min_node_id_index, face_rotation, rotated_indices);
615 804 : std::vector<const Node *> rotated_pyramid_nodes;
616 4824 : for (unsigned int i = 0; i < 5; i++)
617 4020 : rotated_pyramid_nodes.push_back(pyramid_nodes[rotated_indices[i]]);
618 :
619 : // There is only one quad face in a pyramid element, so the splitting selection is binary
620 3216 : const std::vector<std::vector<unsigned int>> tet_nodes_set = {{0, 1, 2, 4}, {0, 2, 3, 4}};
621 3216 : const std::vector<std::vector<unsigned int>> tet_face_indices = {{4, 0, 1, 5}, {4, 5, 2, 3}};
622 :
623 : // Based on the determined splitting direction, determine the nodes of each resulting TET4
624 : // elements after the splitting.
625 2412 : for (const auto & tet_face_index : tet_face_indices)
626 : {
627 1608 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
628 8040 : for (const auto & face_index : tet_face_index)
629 : {
630 6432 : if (face_index < 5)
631 4824 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
632 : else
633 1608 : rotated_tet_face_indices.back().push_back(5);
634 : }
635 : }
636 :
637 2412 : for (const auto & tet_nodes : tet_nodes_set)
638 : {
639 1608 : tet_nodes_list.push_back(std::vector<const Node *>());
640 8040 : for (const auto & tet_node : tet_nodes)
641 6432 : tet_nodes_list.back().push_back(rotated_pyramid_nodes[tet_node]);
642 : }
643 4824 : }
644 :
645 : void
646 268 : convert3DMeshToAllTet4(MeshBase & mesh,
647 : const std::vector<std::pair<dof_id_type, bool>> & elems_to_process,
648 : std::vector<dof_id_type> & converted_elems_ids_to_track,
649 : const subdomain_id_type block_id_to_remove,
650 : const bool delete_block_to_remove)
651 : {
652 : mooseAssert(mesh.is_serial(),
653 : "This method only supports serial meshes. If you are calling this method with a "
654 : "distributed mesh, please serialize it first.");
655 268 : std::vector<dof_id_type> converted_elems_ids_to_retain;
656 : // Build boundary information of the mesh
657 268 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
658 268 : const auto bdry_side_list = boundary_info.build_side_list();
659 11104 : for (const auto & elem_to_process : elems_to_process)
660 : {
661 10836 : switch (mesh.elem_ptr(elem_to_process.first)->type())
662 : {
663 8198 : case ElemType::HEX8:
664 8198 : hexElemSplitter(mesh,
665 : bdry_side_list,
666 8198 : elem_to_process.first,
667 8198 : elem_to_process.second ? converted_elems_ids_to_track
668 : : converted_elems_ids_to_retain);
669 8198 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
670 8198 : break;
671 512 : case ElemType::PYRAMID5:
672 512 : pyramidElemSplitter(mesh,
673 : bdry_side_list,
674 512 : elem_to_process.first,
675 512 : elem_to_process.second ? converted_elems_ids_to_track
676 : : converted_elems_ids_to_retain);
677 512 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
678 512 : break;
679 696 : case ElemType::PRISM6:
680 696 : prismElemSplitter(mesh,
681 : bdry_side_list,
682 696 : elem_to_process.first,
683 696 : elem_to_process.second ? converted_elems_ids_to_track
684 : : converted_elems_ids_to_retain);
685 696 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
686 696 : break;
687 150 : case ElemType::C0POLYHEDRON:
688 150 : polyhedronElemSplitter(mesh,
689 : bdry_side_list,
690 150 : elem_to_process.first,
691 150 : elem_to_process.second ? converted_elems_ids_to_track
692 : : converted_elems_ids_to_retain);
693 150 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
694 150 : break;
695 1280 : case ElemType::TET4:
696 1280 : if (elem_to_process.second)
697 512 : converted_elems_ids_to_track.push_back(elem_to_process.first);
698 : else
699 768 : converted_elems_ids_to_retain.push_back(elem_to_process.first);
700 1280 : break;
701 0 : default:
702 0 : mooseError("Unexpected element type.");
703 : }
704 : }
705 :
706 268 : if (delete_block_to_remove)
707 : {
708 1506 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
709 1506 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
710 : elem_it++)
711 1506 : mesh.delete_elem(*elem_it);
712 :
713 108 : mesh.contract();
714 108 : mesh.prepare_for_use();
715 : }
716 268 : }
717 :
718 : void
719 111 : convert3DMeshToAllTet4(MeshBase & mesh)
720 : {
721 : mooseAssert(mesh.is_serial(),
722 : "This method only supports serial meshes. If you are calling this method with a "
723 : "distributed mesh, please serialize it first.");
724 : // Subdomain ID for new utility blocks must be new
725 111 : std::set<subdomain_id_type> subdomain_ids_set;
726 111 : mesh.subdomain_ids(subdomain_ids_set);
727 111 : const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin();
728 111 : const subdomain_id_type block_id_to_remove = max_subdomain_id + 1;
729 111 : std::vector<std::pair<dof_id_type, bool>> original_elems;
730 :
731 2277 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
732 : elem_it++)
733 : {
734 2169 : if ((*elem_it)->default_order() != Order::FIRST)
735 3 : mooseError("Only first order elements are supported for cutting.");
736 2166 : original_elems.push_back(std::make_pair((*elem_it)->id(), false));
737 108 : }
738 :
739 108 : std::vector<dof_id_type> converted_elems_ids_to_track;
740 :
741 108 : convert3DMeshToAllTet4(
742 : mesh, original_elems, converted_elems_ids_to_track, block_id_to_remove, true);
743 108 : }
744 :
745 : void
746 45000 : elementBoundaryInfoCollector(const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
747 : const dof_id_type elem_id,
748 : const unsigned short n_elem_sides,
749 : std::vector<std::vector<boundary_id_type>> & elem_side_list)
750 : {
751 45000 : elem_side_list.resize(n_elem_sides);
752 : const auto selected_bdry_side_list =
753 45000 : std::equal_range(bdry_side_list.begin(), bdry_side_list.end(), elem_id, BCTupleKeyComp{});
754 45000 : for (auto selected_bdry_side = selected_bdry_side_list.first;
755 61000 : selected_bdry_side != selected_bdry_side_list.second;
756 16000 : ++selected_bdry_side)
757 : {
758 16000 : elem_side_list[std::get<1>(*selected_bdry_side)].push_back(std::get<2>(*selected_bdry_side));
759 : }
760 45000 : }
761 :
762 : void
763 1136 : convertElem(MeshBase & mesh,
764 : const dof_id_type & elem_id,
765 : const std::vector<unsigned int> & side_indices,
766 : const std::vector<std::vector<boundary_id_type>> & elem_side_info,
767 : const SubdomainID & subdomain_id_shift_base)
768 : {
769 1136 : const auto & elem_type = mesh.elem_ptr(elem_id)->type();
770 1136 : switch (elem_type)
771 : {
772 944 : case HEX8:
773 : // HEX8 to PYRAMID5 (+2*subdomain_id_shift_base)
774 : // HEX8 to TET4 (+subdomain_id_shift_base)
775 944 : convertHex8Elem(mesh, elem_id, side_indices, elem_side_info, subdomain_id_shift_base);
776 944 : break;
777 96 : case PRISM6:
778 : // PRISM6 to TET4 (+subdomain_id_shift_base)
779 : // PRISM6 to PYRAMID5 (+2*subdomain_id_shift_base)
780 96 : convertPrism6Elem(mesh, elem_id, side_indices, elem_side_info, subdomain_id_shift_base);
781 96 : break;
782 96 : case PYRAMID5:
783 : // PYRAMID5 to TET4 (+subdomain_id_shift_base)
784 96 : convertPyramid5Elem(mesh, elem_id, elem_side_info, subdomain_id_shift_base);
785 96 : break;
786 1136 : default:
787 : mooseAssert(false,
788 : "The provided element type '" + std::to_string(elem_type) +
789 : "' is not supported and is not supposed to be passed to this function. "
790 : "Only HEX8, PRISM6 and PYRAMID5 are supported.");
791 : }
792 1136 : }
793 :
794 : void
795 944 : convertHex8Elem(MeshBase & mesh,
796 : const dof_id_type & elem_id,
797 : const std::vector<unsigned int> & side_indices,
798 : const std::vector<std::vector<boundary_id_type>> & elem_side_info,
799 : const SubdomainID & subdomain_id_shift_base)
800 : {
801 : // We add a node at the centroid of the HEX8 element
802 : // With this node, the HEX8 can be converted into 6 PYRAMID5 elements
803 : // For the PYRAMID5 element at the 'side_indices', they can further be converted into 2 TET4
804 : // elements
805 944 : const Point elem_cent = mesh.elem_ptr(elem_id)->true_centroid();
806 944 : auto new_node = mesh.add_point(elem_cent);
807 6608 : for (const auto & i_side : make_range(mesh.elem_ptr(elem_id)->n_sides()))
808 : {
809 5664 : if (std::find(side_indices.begin(), side_indices.end(), i_side) != side_indices.end())
810 1792 : createUnitTet4FromHex8(
811 1792 : mesh, elem_id, i_side, new_node, elem_side_info[i_side], subdomain_id_shift_base);
812 : else
813 3872 : createUnitPyramid5FromHex8(
814 3872 : mesh, elem_id, i_side, new_node, elem_side_info[i_side], subdomain_id_shift_base);
815 : }
816 944 : }
817 :
818 : void
819 4064 : createUnitPyramid5FromHex8(MeshBase & mesh,
820 : const dof_id_type & elem_id,
821 : const unsigned int & side_index,
822 : const Node * new_node,
823 : const std::vector<boundary_id_type> & side_info,
824 : const SubdomainID & subdomain_id_shift_base)
825 : {
826 4064 : auto new_elem = std::make_unique<Pyramid5>();
827 4064 : new_elem->set_node(0, mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(3));
828 4064 : new_elem->set_node(1, mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(2));
829 4064 : new_elem->set_node(2, mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(1));
830 4064 : new_elem->set_node(3, mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(0));
831 4064 : new_elem->set_node(4, const_cast<Node *>(new_node));
832 4064 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base * 2;
833 4064 : auto new_elem_ptr = mesh.add_elem(std::move(new_elem));
834 4064 : retainEEID(mesh, elem_id, new_elem_ptr);
835 4786 : for (const auto & bid : side_info)
836 722 : mesh.get_boundary_info().add_side(new_elem_ptr, 4, bid);
837 4064 : }
838 :
839 : void
840 1792 : createUnitTet4FromHex8(MeshBase & mesh,
841 : const dof_id_type & elem_id,
842 : const unsigned int & side_index,
843 : const Node * new_node,
844 : const std::vector<boundary_id_type> & side_info,
845 : const SubdomainID & subdomain_id_shift_base)
846 : {
847 : // We want to make sure that the QUAD4 is divided by the diagonal that involves the node with
848 : // the lowest node id This may help maintain consistency for future applications
849 1792 : unsigned int lid_index = 0;
850 7168 : for (const auto & i : make_range(1, 4))
851 : {
852 10752 : if (mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(i)->id() <
853 5376 : mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(lid_index)->id())
854 708 : lid_index = i;
855 : }
856 :
857 1792 : auto new_elem_0 = std::make_unique<Tet4>();
858 5376 : new_elem_0->set_node(0,
859 1792 : mesh.elem_ptr(elem_id)
860 3584 : ->side_ptr(side_index)
861 1792 : ->node_ptr(MathUtils::euclideanMod(2 - lid_index % 2, 4)));
862 5376 : new_elem_0->set_node(1,
863 1792 : mesh.elem_ptr(elem_id)
864 3584 : ->side_ptr(side_index)
865 1792 : ->node_ptr(MathUtils::euclideanMod(1 - lid_index % 2, 4)));
866 5376 : new_elem_0->set_node(2,
867 1792 : mesh.elem_ptr(elem_id)
868 3584 : ->side_ptr(side_index)
869 1792 : ->node_ptr(MathUtils::euclideanMod(0 - lid_index % 2, 4)));
870 1792 : new_elem_0->set_node(3, const_cast<Node *>(new_node));
871 1792 : new_elem_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base;
872 1792 : auto new_elem_ptr_0 = mesh.add_elem(std::move(new_elem_0));
873 1792 : retainEEID(mesh, elem_id, new_elem_ptr_0);
874 :
875 1792 : auto new_elem_1 = std::make_unique<Tet4>();
876 5376 : new_elem_1->set_node(0,
877 1792 : mesh.elem_ptr(elem_id)
878 3584 : ->side_ptr(side_index)
879 1792 : ->node_ptr(MathUtils::euclideanMod(3 - lid_index % 2, 4)));
880 5376 : new_elem_1->set_node(1,
881 1792 : mesh.elem_ptr(elem_id)
882 3584 : ->side_ptr(side_index)
883 1792 : ->node_ptr(MathUtils::euclideanMod(2 - lid_index % 2, 4)));
884 5376 : new_elem_1->set_node(2,
885 1792 : mesh.elem_ptr(elem_id)
886 3584 : ->side_ptr(side_index)
887 1792 : ->node_ptr(MathUtils::euclideanMod(0 - lid_index % 2, 4)));
888 1792 : new_elem_1->set_node(3, const_cast<Node *>(new_node));
889 1792 : new_elem_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base;
890 1792 : auto new_elem_ptr_1 = mesh.add_elem(std::move(new_elem_1));
891 1792 : retainEEID(mesh, elem_id, new_elem_ptr_1);
892 :
893 3604 : for (const auto & bid : side_info)
894 : {
895 1812 : mesh.get_boundary_info().add_side(new_elem_ptr_0, 0, bid);
896 1812 : mesh.get_boundary_info().add_side(new_elem_ptr_1, 0, bid);
897 : }
898 1792 : }
899 :
900 : void
901 96 : convertPrism6Elem(MeshBase & mesh,
902 : const dof_id_type & elem_id,
903 : const std::vector<unsigned int> & side_indices,
904 : const std::vector<std::vector<boundary_id_type>> & elem_side_info,
905 : const SubdomainID & subdomain_id_shift_base)
906 : {
907 : // We add a node at the centroid of the PRISM6 element
908 : // With this node, the PRISM6 can be converted into 3 PYRAMID5 elements and 2 TET4 elements
909 : // For the PYRAMID5 element, it can further be converted into 2 TET3 elements
910 96 : const Point elem_cent = mesh.elem_ptr(elem_id)->true_centroid();
911 96 : auto new_node = mesh.add_point(elem_cent);
912 576 : for (const auto & i_side : make_range(mesh.elem_ptr(elem_id)->n_sides()))
913 : {
914 768 : if (i_side % 4 == 0 ||
915 768 : std::find(side_indices.begin(), side_indices.end(), i_side) != side_indices.end())
916 288 : createUnitTet4FromPrism6(
917 288 : mesh, elem_id, i_side, new_node, elem_side_info[i_side], subdomain_id_shift_base);
918 : else
919 192 : createUnitPyramid5FromPrism6(
920 192 : mesh, elem_id, i_side, new_node, elem_side_info[i_side], subdomain_id_shift_base);
921 : }
922 96 : }
923 :
924 : void
925 288 : createUnitTet4FromPrism6(MeshBase & mesh,
926 : const dof_id_type & elem_id,
927 : const unsigned int & side_index,
928 : const Node * new_node,
929 : const std::vector<boundary_id_type> & side_info,
930 : const SubdomainID & subdomain_id_shift_base)
931 : {
932 : // For side 1 and side 4, they are already TRI3, so only one TET is created
933 : // For side 0, 2, and 3, they are QUAD4, so we create 2 TETs
934 : // We want to make sure that the QUAD4 is divided by the diagonal that involves
935 : // the node with the lowest node id This may help maintain consistency for future applications
936 288 : bool is_side_quad = (side_index % 4 != 0);
937 288 : unsigned int lid_index = 0;
938 288 : if (is_side_quad)
939 384 : for (const auto & i : make_range(1, 4))
940 : {
941 576 : if (mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(i)->id() <
942 288 : mesh.elem_ptr(elem_id)->side_ptr(side_index)->node_ptr(lid_index)->id())
943 96 : lid_index = i;
944 : }
945 : // For a TRI3 side, lid_index is always 0, so the indices are always 2,1,0 here
946 288 : auto new_elem_0 = std::make_unique<Tet4>();
947 864 : new_elem_0->set_node(0,
948 288 : mesh.elem_ptr(elem_id)
949 576 : ->side_ptr(side_index)
950 288 : ->node_ptr(MathUtils::euclideanMod(2 - lid_index % 2, 4)));
951 864 : new_elem_0->set_node(1,
952 288 : mesh.elem_ptr(elem_id)
953 576 : ->side_ptr(side_index)
954 288 : ->node_ptr(MathUtils::euclideanMod(1 - lid_index % 2, 4)));
955 864 : new_elem_0->set_node(2,
956 288 : mesh.elem_ptr(elem_id)
957 576 : ->side_ptr(side_index)
958 288 : ->node_ptr(MathUtils::euclideanMod(0 - lid_index % 2, 4)));
959 288 : new_elem_0->set_node(3, const_cast<Node *>(new_node));
960 288 : new_elem_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base;
961 288 : auto new_elem_ptr_0 = mesh.add_elem(std::move(new_elem_0));
962 288 : retainEEID(mesh, elem_id, new_elem_ptr_0);
963 :
964 288 : Elem * new_elem_ptr_1 = nullptr;
965 288 : if (is_side_quad)
966 : {
967 96 : auto new_elem_1 = std::make_unique<Tet4>();
968 288 : new_elem_1->set_node(0,
969 96 : mesh.elem_ptr(elem_id)
970 192 : ->side_ptr(side_index)
971 96 : ->node_ptr(MathUtils::euclideanMod(3 - lid_index % 2, 4)));
972 288 : new_elem_1->set_node(1,
973 96 : mesh.elem_ptr(elem_id)
974 192 : ->side_ptr(side_index)
975 96 : ->node_ptr(MathUtils::euclideanMod(2 - lid_index % 2, 4)));
976 288 : new_elem_1->set_node(2,
977 96 : mesh.elem_ptr(elem_id)
978 192 : ->side_ptr(side_index)
979 96 : ->node_ptr(MathUtils::euclideanMod(0 - lid_index % 2, 4)));
980 96 : new_elem_1->set_node(3, const_cast<Node *>(new_node));
981 96 : new_elem_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base;
982 96 : new_elem_ptr_1 = mesh.add_elem(std::move(new_elem_1));
983 96 : retainEEID(mesh, elem_id, new_elem_ptr_1);
984 96 : }
985 :
986 576 : for (const auto & bid : side_info)
987 : {
988 288 : mesh.get_boundary_info().add_side(new_elem_ptr_0, 0, bid);
989 288 : if (new_elem_ptr_1)
990 192 : mesh.get_boundary_info().add_side(new_elem_ptr_1, 0, bid);
991 : }
992 288 : }
993 :
994 : void
995 192 : createUnitPyramid5FromPrism6(MeshBase & mesh,
996 : const dof_id_type & elem_id,
997 : const unsigned int & side_index,
998 : const Node * new_node,
999 : const std::vector<boundary_id_type> & side_info,
1000 : const SubdomainID & subdomain_id_shift_base)
1001 : {
1002 : // Same as Hex8
1003 192 : createUnitPyramid5FromHex8(
1004 : mesh, elem_id, side_index, new_node, side_info, subdomain_id_shift_base);
1005 192 : }
1006 :
1007 : void
1008 96 : convertPyramid5Elem(MeshBase & mesh,
1009 : const dof_id_type & elem_id,
1010 : const std::vector<std::vector<boundary_id_type>> & elem_side_info,
1011 : const SubdomainID & subdomain_id_shift_base)
1012 : {
1013 : // A Pyramid5 element has only one QUAD4 face, so we can convert it to 2 TET4 elements
1014 96 : unsigned int lid_index = 0;
1015 384 : for (const auto & i : make_range(1, 4))
1016 : {
1017 576 : if (mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(i)->id() <
1018 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(lid_index)->id())
1019 80 : lid_index = i;
1020 : }
1021 96 : auto new_elem_0 = std::make_unique<Tet4>();
1022 288 : new_elem_0->set_node(
1023 : 0,
1024 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(MathUtils::euclideanMod(2 - lid_index % 2, 4)));
1025 288 : new_elem_0->set_node(
1026 : 1,
1027 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(MathUtils::euclideanMod(1 - lid_index % 2, 4)));
1028 288 : new_elem_0->set_node(
1029 : 2,
1030 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(MathUtils::euclideanMod(0 - lid_index % 2, 4)));
1031 96 : new_elem_0->set_node(3, mesh.elem_ptr(elem_id)->node_ptr(4));
1032 96 : new_elem_0->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base;
1033 96 : auto new_elem_ptr_0 = mesh.add_elem(std::move(new_elem_0));
1034 96 : retainEEID(mesh, elem_id, new_elem_ptr_0);
1035 :
1036 96 : auto new_elem_1 = std::make_unique<Tet4>();
1037 288 : new_elem_1->set_node(
1038 : 0,
1039 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(MathUtils::euclideanMod(3 - lid_index % 2, 4)));
1040 288 : new_elem_1->set_node(
1041 : 1,
1042 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(MathUtils::euclideanMod(2 - lid_index % 2, 4)));
1043 288 : new_elem_1->set_node(
1044 : 2,
1045 288 : mesh.elem_ptr(elem_id)->side_ptr(4)->node_ptr(MathUtils::euclideanMod(0 - lid_index % 2, 4)));
1046 96 : new_elem_1->set_node(3, mesh.elem_ptr(elem_id)->node_ptr(4));
1047 96 : new_elem_1->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id() + subdomain_id_shift_base;
1048 96 : auto new_elem_ptr_1 = mesh.add_elem(std::move(new_elem_1));
1049 96 : retainEEID(mesh, elem_id, new_elem_ptr_1);
1050 :
1051 96 : for (const auto & bid : elem_side_info[0])
1052 0 : mesh.get_boundary_info().add_side(new_elem_ptr_0, 2 - lid_index % 2, bid);
1053 96 : for (const auto & bid : elem_side_info[1])
1054 0 : if (lid_index % 2)
1055 0 : mesh.get_boundary_info().add_side(new_elem_ptr_1, 2, bid);
1056 : else
1057 0 : mesh.get_boundary_info().add_side(new_elem_ptr_0, 1, bid);
1058 96 : for (const auto & bid : elem_side_info[2])
1059 0 : mesh.get_boundary_info().add_side(new_elem_ptr_1, 2 + lid_index % 2, bid);
1060 96 : for (const auto & bid : elem_side_info[3])
1061 0 : if (lid_index % 2)
1062 0 : mesh.get_boundary_info().add_side(new_elem_ptr_0, 1, bid);
1063 : else
1064 0 : mesh.get_boundary_info().add_side(new_elem_ptr_1, 3, bid);
1065 288 : for (const auto & bid : elem_side_info[4])
1066 : {
1067 192 : mesh.get_boundary_info().add_side(new_elem_ptr_0, 0, bid);
1068 192 : mesh.get_boundary_info().add_side(new_elem_ptr_1, 0, bid);
1069 : }
1070 96 : }
1071 :
1072 : void
1073 8224 : retainEEID(MeshBase & mesh, const dof_id_type & elem_id, Elem * new_elem_ptr)
1074 : {
1075 8224 : const unsigned int n_eeid = mesh.n_elem_integers();
1076 8224 : for (const auto & i : make_range(n_eeid))
1077 0 : new_elem_ptr->set_extra_integer(i, mesh.elem_ptr(elem_id)->get_extra_integer(i));
1078 8224 : }
1079 :
1080 : void
1081 81 : transitionLayerGenerator(MeshBase & mesh,
1082 : const std::vector<BoundaryName> & boundary_names,
1083 : const unsigned int conversion_element_layer_number,
1084 : const bool external_boundaries_checking)
1085 : {
1086 : mooseAssert(mesh.is_serial(),
1087 : "This method only supports serial meshes. If you are calling this method with a "
1088 : "distributed mesh, please serialize it first.");
1089 : // The base subdomain ID to shift the original elements because of the element type change
1090 81 : const auto sid_shift_base = MooseMeshUtils::getNextFreeSubdomainID(mesh);
1091 : // The maximum subdomain ID that would be involved is sid_shift_base * 3, we would like to make
1092 : // sure it would not overflow
1093 81 : if (sid_shift_base * 3 > std::numeric_limits<subdomain_id_type>::max())
1094 3 : throw MooseException("subdomain id overflow");
1095 :
1096 : // It would be convenient to have a single boundary id instead of a vector.
1097 78 : const auto uniform_tmp_bid = MooseMeshUtils::getNextFreeBoundaryID(mesh);
1098 :
1099 : // Check the boundaries and merge them
1100 78 : std::vector<boundary_id_type> boundary_ids;
1101 177 : for (const auto & sideset : boundary_names)
1102 : {
1103 102 : if (!MooseMeshUtils::hasBoundaryName(mesh, sideset))
1104 3 : throw MooseException("The provided boundary '", sideset, "' was not found within the mesh");
1105 99 : boundary_ids.push_back(MooseMeshUtils::getBoundaryID(sideset, mesh));
1106 99 : MooseMeshUtils::changeBoundaryId(mesh, boundary_ids.back(), uniform_tmp_bid, false);
1107 : }
1108 :
1109 75 : auto & sideset_map = mesh.get_boundary_info().get_sideset_map();
1110 75 : auto side_list = mesh.get_boundary_info().build_side_list();
1111 :
1112 75 : std::vector<std::pair<dof_id_type, std::vector<unsigned int>>> elems_list;
1113 75 : std::vector<std::set<dof_id_type>> layered_elems_list;
1114 75 : layered_elems_list.push_back(std::set<dof_id_type>());
1115 : // Need to collect the list of elements that need to be converted
1116 75 : if (!mesh.is_prepared())
1117 75 : mesh.find_neighbors();
1118 4229 : for (const auto & side_info : side_list)
1119 : {
1120 4160 : if (std::get<2>(side_info) == uniform_tmp_bid)
1121 : {
1122 : // Check if the involved side is TRI3 or QUAD4
1123 : // We do not limit the element type in the input mesh
1124 : // As long as the involved boundaries only consist of TRI3 and QUAD4 sides,
1125 : // this generator will work
1126 : // As the side element of a quadratic element is still a linear element in libMesh,
1127 : // we need to check the element's default_side_order() first
1128 918 : if (mesh.elem_ptr(std::get<0>(side_info))->default_side_order() != 1)
1129 : throw MooseException(
1130 3 : "The provided boundary set contains non-linear side elements, which is not supported.");
1131 : const auto side_type =
1132 915 : mesh.elem_ptr(std::get<0>(side_info))->side_ptr(std::get<1>(side_info))->type();
1133 915 : layered_elems_list.back().emplace(std::get<0>(side_info));
1134 : // If we enforce external boundary, then a non-null neighbor leads to an error
1135 : // Otherwise, we need to convert both sides, that means the neighbor information also needs to
1136 : // be added
1137 : const auto neighbor_ptr =
1138 915 : mesh.elem_ptr(std::get<0>(side_info))->neighbor_ptr(std::get<1>(side_info));
1139 915 : if (neighbor_ptr)
1140 : {
1141 75 : if (external_boundaries_checking)
1142 : throw MooseException(
1143 : "The provided boundary contains non-external sides, which is required when "
1144 3 : "external_boundaries_checking is enabled.");
1145 : else
1146 72 : layered_elems_list.back().emplace(neighbor_ptr->id());
1147 : }
1148 :
1149 912 : if (conversion_element_layer_number == 1)
1150 : {
1151 828 : if (side_type == TRI3)
1152 0 : continue; // Already TRI3, no need to convert
1153 828 : else if (side_type == QUAD4)
1154 : {
1155 828 : auto pit = std::find_if(elems_list.begin(),
1156 : elems_list.end(),
1157 828 : [elem_id = std::get<0>(side_info)](const auto & p)
1158 7828 : { return p.first == elem_id; });
1159 828 : if (elems_list.size() && pit != elems_list.end())
1160 : {
1161 296 : pit->second.push_back(std::get<1>(side_info));
1162 : }
1163 : else
1164 532 : elems_list.push_back(std::make_pair(
1165 1064 : std::get<0>(side_info), std::vector<unsigned int>({std::get<1>(side_info)})));
1166 828 : if (neighbor_ptr)
1167 : {
1168 0 : auto sit = std::find_if(elems_list.begin(),
1169 : elems_list.end(),
1170 0 : [elem_id = neighbor_ptr->id()](const auto & p)
1171 0 : { return p.first == elem_id; });
1172 0 : if (elems_list.size() && sit != elems_list.end())
1173 : {
1174 0 : sit->second.push_back(
1175 0 : neighbor_ptr->which_neighbor_am_i(mesh.elem_ptr(std::get<0>(side_info))));
1176 : }
1177 : else
1178 : {
1179 0 : elems_list.push_back(std::make_pair(
1180 0 : neighbor_ptr->id(),
1181 0 : std::vector<unsigned int>(
1182 0 : {neighbor_ptr->which_neighbor_am_i(mesh.elem_ptr(std::get<0>(side_info)))})));
1183 : }
1184 : }
1185 : }
1186 0 : else if (side_type == C0POLYGON)
1187 : throw MooseException("The provided boundary set contains C0POLYGON side elements, which "
1188 0 : "is not supported.");
1189 : else
1190 : mooseAssert(false,
1191 : "Impossible scenario: a linear non-polygon side element that is neither TRI3 "
1192 : "nor QUAD4.");
1193 : }
1194 : }
1195 : }
1196 :
1197 69 : if (conversion_element_layer_number > 1)
1198 : {
1199 11 : std::set<dof_id_type> total_elems_set(layered_elems_list.back());
1200 :
1201 47 : while (layered_elems_list.back().size() &&
1202 22 : layered_elems_list.size() < conversion_element_layer_number)
1203 : {
1204 14 : layered_elems_list.push_back(std::set<dof_id_type>());
1205 182 : for (const auto & elem_id : *(layered_elems_list.end() - 2))
1206 : {
1207 1176 : for (const auto & i_side : make_range(mesh.elem_ptr(elem_id)->n_sides()))
1208 : {
1209 1008 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i_side) != nullptr)
1210 : {
1211 744 : const auto & neighbor_id = mesh.elem_ptr(elem_id)->neighbor_ptr(i_side)->id();
1212 744 : if (total_elems_set.find(neighbor_id) == total_elems_set.end())
1213 : {
1214 156 : layered_elems_list.back().emplace(neighbor_id);
1215 156 : total_elems_set.emplace(neighbor_id);
1216 : }
1217 : }
1218 : }
1219 : }
1220 : }
1221 11 : }
1222 :
1223 : // Remove the last empty layer
1224 69 : if (layered_elems_list.back().empty())
1225 3 : layered_elems_list.pop_back();
1226 :
1227 69 : if (conversion_element_layer_number > layered_elems_list.size())
1228 : throw MooseException("There is fewer layers of elements in the input mesh than the requested "
1229 3 : "number of layers to convert.");
1230 :
1231 66 : std::vector<std::pair<dof_id_type, bool>> original_elems;
1232 : // construct a list of the element to convert to tet4
1233 : // Convert at most n_layer_conversion layers of elements
1234 66 : const unsigned int n_layer_conversion = layered_elems_list.size() - 1;
1235 74 : for (unsigned int i = 0; i < n_layer_conversion; ++i)
1236 152 : for (const auto & elem_id : layered_elems_list[i])
1237 : {
1238 : // As these elements will become TET4 elements, we need to shift the subdomain ID
1239 : // But we do not need to convert original TET4 elements
1240 144 : if (mesh.elem_ptr(elem_id)->type() != TET4)
1241 : {
1242 144 : original_elems.push_back(std::make_pair(elem_id, false));
1243 144 : mesh.elem_ptr(elem_id)->subdomain_id() += sid_shift_base;
1244 : }
1245 : }
1246 :
1247 66 : const subdomain_id_type block_id_to_remove = sid_shift_base * 3;
1248 :
1249 66 : std::vector<dof_id_type> converted_elems_ids_to_track;
1250 66 : MooseMeshElementConversionUtils::convert3DMeshToAllTet4(
1251 : mesh, original_elems, converted_elems_ids_to_track, block_id_to_remove, false);
1252 :
1253 : // Now we need to convert the elements on the transition layer
1254 : // First, we need to identify that the sides that are on the interface with previous layer (all
1255 : // TET layers)
1256 66 : if (n_layer_conversion)
1257 : {
1258 152 : for (const auto & elem_id : layered_elems_list[n_layer_conversion])
1259 : {
1260 1008 : for (const auto & i_side : make_range(mesh.elem_ptr(elem_id)->n_sides()))
1261 : {
1262 1536 : if (mesh.elem_ptr(elem_id)->neighbor_ptr(i_side) != nullptr &&
1263 1344 : layered_elems_list[n_layer_conversion - 1].count(
1264 1536 : mesh.elem_ptr(elem_id)->neighbor_ptr(i_side)->id()))
1265 : {
1266 144 : if (elems_list.size() && elems_list.back().first == elem_id)
1267 0 : elems_list.back().second.push_back(i_side);
1268 : else
1269 144 : elems_list.push_back(std::make_pair(elem_id, std::vector<unsigned int>({i_side})));
1270 : }
1271 : }
1272 : }
1273 : }
1274 :
1275 : // Now convert the elements
1276 742 : for (const auto & elem_info : elems_list)
1277 : {
1278 : // Find the involved sidesets of the element so that we can retain them
1279 : std::vector<std::vector<boundary_id_type>> elem_side_info(
1280 676 : mesh.elem_ptr(elem_info.first)->n_sides());
1281 676 : auto side_range = sideset_map.equal_range(mesh.elem_ptr(elem_info.first));
1282 3352 : for (auto i = side_range.first; i != side_range.second; ++i)
1283 2676 : elem_side_info[i->second.first].push_back(i->second.second);
1284 :
1285 676 : MooseMeshElementConversionUtils::convertElem(
1286 676 : mesh, elem_info.first, elem_info.second, elem_side_info, sid_shift_base);
1287 676 : }
1288 :
1289 : // delete the original elements that were converted
1290 742 : for (const auto & elem_info : elems_list)
1291 676 : mesh.elem_ptr(elem_info.first)->subdomain_id() = block_id_to_remove;
1292 886 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
1293 886 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
1294 : elem_it++)
1295 886 : mesh.delete_elem(*elem_it);
1296 : // delete temporary boundary id
1297 66 : mesh.get_boundary_info().remove_id(uniform_tmp_bid);
1298 :
1299 66 : mesh.contract();
1300 66 : mesh.unset_is_prepared();
1301 105 : }
1302 :
1303 : void
1304 86 : assignConvertedElementsSubdomainNameSuffix(
1305 : MeshBase & mesh,
1306 : const std::set<subdomain_id_type> & original_subdomain_ids,
1307 : const subdomain_id_type sid_shift_base,
1308 : const SubdomainName & tet_suffix,
1309 : const SubdomainName & pyramid_suffix)
1310 : {
1311 : // If we have an unprepared mesh, we at least need its element
1312 : // caches prepared for subdomain_ids
1313 86 : if (!mesh.preparation().has_cached_elem_data)
1314 56 : mesh.cache_elem_data();
1315 :
1316 190 : for (const auto & subdomain_id : original_subdomain_ids)
1317 : {
1318 110 : if (MooseMeshUtils::hasSubdomainID(mesh, subdomain_id + sid_shift_base))
1319 : {
1320 : const SubdomainName new_name =
1321 220 : (mesh.subdomain_name(subdomain_id).empty() ? std::to_string(subdomain_id)
1322 220 : : mesh.subdomain_name(subdomain_id)) +
1323 220 : '_' + tet_suffix;
1324 110 : if (MooseMeshUtils::hasSubdomainName(mesh, new_name))
1325 : throw MooseException(
1326 6 : "This suffix for converted TET4 elements results in a subdomain name, " + new_name +
1327 9 : ", that already exists in the mesh. Please choose a different suffix.");
1328 107 : mesh.subdomain_name(subdomain_id + sid_shift_base) = new_name;
1329 110 : }
1330 107 : if (MooseMeshUtils::hasSubdomainID(mesh, subdomain_id + 2 * sid_shift_base))
1331 : {
1332 : const SubdomainName new_name =
1333 182 : (mesh.subdomain_name(subdomain_id).empty() ? std::to_string(subdomain_id)
1334 182 : : mesh.subdomain_name(subdomain_id)) +
1335 182 : '_' + pyramid_suffix;
1336 91 : if (MooseMeshUtils::hasSubdomainName(mesh, new_name))
1337 : throw MooseException(
1338 6 : "This suffix for converted PYRAMID5 elements results in a subdomain name, " + new_name +
1339 9 : ", that already exists in the mesh. Please choose a different suffix.");
1340 88 : mesh.subdomain_name(subdomain_id + 2 * sid_shift_base) = new_name;
1341 91 : }
1342 : }
1343 80 : }
1344 : }
|
