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