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 :
14 : #include "libmesh/elem.h"
15 : #include "libmesh/enum_order.h"
16 : #include "libmesh/boundary_info.h"
17 : #include "libmesh/mesh_base.h"
18 : #include "libmesh/parallel.h"
19 : #include "libmesh/parallel_algebra.h"
20 : #include "libmesh/utility.h"
21 : #include "libmesh/cell_tet4.h"
22 : #include "libmesh/face_tri3.h"
23 :
24 : using namespace libMesh;
25 :
26 : namespace MooseMeshElementConversionUtils
27 : {
28 : void
29 6456 : hexElemSplitter(ReplicatedMesh & mesh,
30 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
31 : const dof_id_type elem_id,
32 : std::vector<dof_id_type> & converted_elems_ids)
33 : {
34 : // Build boundary information of the mesh
35 6456 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
36 : // Create a list of sidesets involving the element to be split
37 6456 : std::vector<std::vector<boundary_id_type>> elem_side_list;
38 6456 : elem_side_list.resize(6);
39 6456 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 6, elem_side_list);
40 :
41 6456 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
42 6456 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
43 : // Record all the element extra integers of the original quad element
44 6456 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
45 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
46 :
47 6456 : std::vector<std::vector<unsigned int>> opt_option;
48 6456 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
49 6456 : mesh.elem_ptr(elem_id)->node_ptr(1),
50 6456 : mesh.elem_ptr(elem_id)->node_ptr(2),
51 6456 : mesh.elem_ptr(elem_id)->node_ptr(3),
52 6456 : mesh.elem_ptr(elem_id)->node_ptr(4),
53 6456 : mesh.elem_ptr(elem_id)->node_ptr(5),
54 6456 : mesh.elem_ptr(elem_id)->node_ptr(6),
55 45192 : mesh.elem_ptr(elem_id)->node_ptr(7)};
56 6456 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
57 :
58 6456 : std::vector<std::vector<const Node *>> optimized_node_list;
59 6456 : hexNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
60 :
61 6456 : std::vector<Elem *> elems_Tet4;
62 45192 : for (const auto i : index_range(optimized_node_list))
63 : {
64 38736 : auto new_elem = std::make_unique<Tet4>();
65 38736 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
66 38736 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
67 38736 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
68 38736 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
69 38736 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
70 38736 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
71 38736 : converted_elems_ids.push_back(elems_Tet4.back()->id());
72 :
73 193680 : for (unsigned int j = 0; j < 4; j++)
74 : {
75 : // A hex element has 6 faces indexed from 0 to 5
76 : // a <6 value indicates that the face of the tet element corresponds to the face of the
77 : // original hex; a =6 value means the face of the tet is an interior face of the hex
78 154944 : if (rotated_tet_face_indices[i][j] < 6)
79 : {
80 86528 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
81 9056 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
82 : }
83 : }
84 38736 : }
85 :
86 : // Retain element extra integers
87 45192 : for (unsigned int i = 0; i < 6; i++)
88 38736 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
89 : {
90 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
91 : }
92 6456 : }
93 :
94 : void
95 696 : prismElemSplitter(ReplicatedMesh & mesh,
96 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
97 : const dof_id_type elem_id,
98 : std::vector<dof_id_type> & converted_elems_ids)
99 : {
100 : // Build boundary information of the mesh
101 696 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
102 : // Create a list of sidesets involving the element to be split
103 696 : std::vector<std::vector<boundary_id_type>> elem_side_list;
104 696 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 5, elem_side_list);
105 :
106 696 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
107 696 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
108 :
109 : // Record all the element extra integers of the original quad element
110 696 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
111 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
112 :
113 696 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
114 696 : mesh.elem_ptr(elem_id)->node_ptr(1),
115 696 : mesh.elem_ptr(elem_id)->node_ptr(2),
116 696 : mesh.elem_ptr(elem_id)->node_ptr(3),
117 696 : mesh.elem_ptr(elem_id)->node_ptr(4),
118 3480 : mesh.elem_ptr(elem_id)->node_ptr(5)};
119 696 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
120 696 : std::vector<std::vector<const Node *>> optimized_node_list;
121 696 : prismNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
122 :
123 696 : std::vector<Elem *> elems_Tet4;
124 2784 : for (const auto i : index_range(optimized_node_list))
125 : {
126 2088 : auto new_elem = std::make_unique<Tet4>();
127 2088 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
128 2088 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
129 2088 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
130 2088 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
131 2088 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
132 2088 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
133 2088 : converted_elems_ids.push_back(elems_Tet4.back()->id());
134 :
135 10440 : for (unsigned int j = 0; j < 4; j++)
136 : {
137 : // A prism element has 5 faces indexed from 0 to 4
138 : // a <4 value indicates that the face of the tet element corresponds to the face of the
139 : // original prism; a =5 value means the face of the tet is an interior face of the prism
140 8352 : if (rotated_tet_face_indices[i][j] < 5)
141 : {
142 6624 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
143 1056 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
144 : }
145 : }
146 2088 : }
147 :
148 : // Retain element extra integers
149 2784 : for (unsigned int i = 0; i < 3; i++)
150 2088 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
151 : {
152 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
153 : }
154 696 : }
155 :
156 : void
157 512 : pyramidElemSplitter(ReplicatedMesh & mesh,
158 : const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
159 : const dof_id_type elem_id,
160 : std::vector<dof_id_type> & converted_elems_ids)
161 : {
162 : // Build boundary information of the mesh
163 512 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
164 : // Create a list of sidesets involving the element to be split
165 512 : std::vector<std::vector<boundary_id_type>> elem_side_list;
166 512 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 5, elem_side_list);
167 :
168 512 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
169 512 : std::vector<dof_id_type> exist_extra_ids(n_elem_extra_ids);
170 : // Record all the element extra integers of the original quad element
171 512 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
172 0 : exist_extra_ids[j] = mesh.elem_ptr(elem_id)->get_extra_integer(j);
173 :
174 512 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
175 512 : mesh.elem_ptr(elem_id)->node_ptr(1),
176 512 : mesh.elem_ptr(elem_id)->node_ptr(2),
177 512 : mesh.elem_ptr(elem_id)->node_ptr(3),
178 2048 : mesh.elem_ptr(elem_id)->node_ptr(4)};
179 512 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
180 512 : std::vector<std::vector<const Node *>> optimized_node_list;
181 512 : pyramidNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
182 :
183 512 : std::vector<Elem *> elems_Tet4;
184 1536 : for (const auto i : index_range(optimized_node_list))
185 : {
186 1024 : auto new_elem = std::make_unique<Tet4>();
187 1024 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
188 1024 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
189 1024 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
190 1024 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
191 1024 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
192 1024 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
193 1024 : converted_elems_ids.push_back(elems_Tet4.back()->id());
194 :
195 5120 : for (unsigned int j = 0; j < 4; j++)
196 : {
197 : // A pyramid element has 5 faces indexed from 0 to 4
198 : // a <4 value indicates that the face of the tet element corresponds to the face of the
199 : // original pyramid; a =5 value means the face of the tet is an interior face of the pyramid
200 4096 : if (rotated_tet_face_indices[i][j] < 5)
201 : {
202 3808 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
203 736 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
204 : }
205 : }
206 1024 : }
207 :
208 : // Retain element extra integers
209 1536 : for (unsigned int i = 0; i < 2; i++)
210 1024 : for (unsigned int j = 0; j < n_elem_extra_ids; j++)
211 : {
212 0 : elems_Tet4[i]->set_extra_integer(j, exist_extra_ids[j]);
213 : }
214 512 : }
215 :
216 : std::vector<unsigned int>
217 6456 : neighborNodeIndicesHEX8(unsigned int min_id_index)
218 : {
219 : const std::vector<std::vector<unsigned int>> preset_indices = {
220 58104 : {1, 3, 4}, {0, 2, 5}, {3, 1, 6}, {2, 0, 7}, {5, 7, 0}, {4, 6, 1}, {7, 5, 2}, {6, 4, 3}};
221 6456 : if (min_id_index > 7)
222 0 : mooseError("The input node index is out of range.");
223 : else
224 12912 : return preset_indices[min_id_index];
225 19368 : }
226 :
227 : void
228 6456 : hexNodesToTetNodesDeterminer(std::vector<const Node *> & hex_nodes,
229 : std::vector<std::vector<unsigned int>> & rotated_tet_face_indices,
230 : std::vector<std::vector<const Node *>> & tet_nodes_list)
231 : {
232 : // Find the node with the minimum id
233 6456 : std::vector<dof_id_type> node_ids(8);
234 58104 : for (unsigned int i = 0; i < 8; i++)
235 51648 : node_ids[i] = hex_nodes[i]->id();
236 :
237 6456 : const unsigned int min_node_id_index = std::distance(
238 6456 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
239 : // Get the index of the three neighbor nodes of the minimum node
240 : // The order is consistent with the description in nodeRotationHEX8()
241 : // Then determine the index of the second minimum node
242 6456 : const auto neighbor_node_indices = neighborNodeIndicesHEX8(min_node_id_index);
243 :
244 6456 : const auto neighbor_node_ids = {node_ids[neighbor_node_indices[0]],
245 6456 : node_ids[neighbor_node_indices[1]],
246 12912 : node_ids[neighbor_node_indices[2]]};
247 : const unsigned int sec_min_pos =
248 6456 : std::distance(std::begin(neighbor_node_ids),
249 6456 : std::min_element(std::begin(neighbor_node_ids), std::end(neighbor_node_ids)));
250 :
251 : // Rotate the node and face indices based on the identified minimum and second minimum nodes
252 : // After the rotation, we guarantee that the minimum node is the first node (Node 0)
253 : // And the second node (Node 1) has the minium global id among the three neighbor nodes of Node 0
254 : // This makes the splitting process simpler
255 6456 : std::vector<unsigned int> face_rotation;
256 6456 : std::vector<unsigned int> rotated_indices;
257 6456 : nodeRotationHEX8(min_node_id_index, sec_min_pos, face_rotation, rotated_indices);
258 6456 : std::vector<const Node *> rotated_hex_nodes;
259 58104 : for (unsigned int i = 0; i < 8; i++)
260 51648 : rotated_hex_nodes.push_back(hex_nodes[rotated_indices[i]]);
261 :
262 : // Find the selection of each face's cutting direction
263 6456 : const auto diagonal_directions = quadFaceDiagonalDirectionsHex(rotated_hex_nodes);
264 :
265 : // Based on the determined splitting directions of all the faces, determine the nodes of each
266 : // resulting TET4 elements after the splitting.
267 6456 : std::vector<std::vector<unsigned int>> tet_face_indices;
268 6456 : const auto tet_nodes_set = tetNodesForHex(diagonal_directions, tet_face_indices);
269 45192 : for (const auto & tet_face_index : tet_face_indices)
270 : {
271 38736 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
272 193680 : for (const auto & face_index : tet_face_index)
273 : {
274 154944 : if (face_index < 6)
275 77472 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
276 : else
277 77472 : rotated_tet_face_indices.back().push_back(6);
278 : }
279 : }
280 :
281 45192 : for (const auto & tet_nodes : tet_nodes_set)
282 : {
283 38736 : tet_nodes_list.push_back(std::vector<const Node *>());
284 193680 : for (const auto & tet_node : tet_nodes)
285 154944 : tet_nodes_list.back().push_back(rotated_hex_nodes[tet_node]);
286 : }
287 6456 : }
288 :
289 : std::vector<bool>
290 6456 : quadFaceDiagonalDirectionsHex(const std::vector<const Node *> & hex_nodes)
291 : {
292 : // Bottom/Top; Front/Back; Right/Left
293 : const std::vector<std::vector<unsigned int>> face_indices = {
294 45192 : {0, 1, 2, 3}, {4, 5, 6, 7}, {0, 1, 5, 4}, {2, 3, 7, 6}, {1, 2, 6, 5}, {3, 0, 4, 7}};
295 6456 : std::vector<bool> diagonal_directions;
296 45192 : for (const auto & face_index : face_indices)
297 : {
298 38736 : std::vector<const Node *> quad_nodes = {hex_nodes[face_index[0]],
299 38736 : hex_nodes[face_index[1]],
300 38736 : hex_nodes[face_index[2]],
301 116208 : hex_nodes[face_index[3]]};
302 38736 : diagonal_directions.push_back(quadFaceDiagonalDirection(quad_nodes));
303 38736 : }
304 12912 : return diagonal_directions;
305 19368 : }
306 :
307 : bool
308 49592 : quadFaceDiagonalDirection(const std::vector<const Node *> & quad_nodes)
309 : {
310 : const std::vector<dof_id_type> node_ids = {
311 49592 : quad_nodes[0]->id(), quad_nodes[1]->id(), quad_nodes[2]->id(), quad_nodes[3]->id()};
312 49592 : const unsigned int min_id_index = std::distance(
313 49592 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
314 49592 : if (min_id_index == 0 || min_id_index == 2)
315 32608 : return true;
316 : else
317 16984 : return false;
318 49592 : }
319 :
320 : std::vector<std::vector<unsigned int>>
321 6456 : tetNodesForHex(const std::vector<bool> diagonal_directions,
322 : std::vector<std::vector<unsigned int>> & tet_face_indices)
323 : {
324 : const std::vector<std::vector<bool>> possible_inputs = {{true, true, true, true, true, false},
325 : {true, true, true, true, false, false},
326 : {true, true, true, false, true, false},
327 : {true, false, true, true, true, false},
328 : {true, false, true, true, false, false},
329 : {true, false, true, false, true, false},
330 51648 : {true, false, true, false, false, false}};
331 :
332 6456 : const unsigned int input_index = std::distance(
333 : std::begin(possible_inputs),
334 6456 : std::find(std::begin(possible_inputs), std::end(possible_inputs), diagonal_directions));
335 :
336 6456 : switch (input_index)
337 : {
338 304 : case 0:
339 2128 : tet_face_indices = {
340 2128 : {0, 6, 2, 6}, {1, 6, 2, 6}, {1, 6, 5, 6}, {0, 6, 3, 4}, {6, 6, 3, 6}, {6, 4, 5, 6}};
341 2128 : 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}};
342 0 : case 1:
343 0 : tet_face_indices = {
344 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}};
345 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}};
346 6152 : case 2:
347 43064 : tet_face_indices = {
348 43064 : {0, 6, 2, 6}, {1, 6, 2, 6}, {1, 6, 5, 6}, {4, 6, 5, 6}, {4, 6, 3, 6}, {0, 6, 3, 6}};
349 43064 : 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}};
350 0 : case 3:
351 0 : tet_face_indices = {
352 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}};
353 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}};
354 0 : case 4:
355 0 : tet_face_indices = {{0, 1, 2, 6}, {0, 6, 3, 4}, {5, 4, 6, 1}, {5, 6, 3, 2}, {6, 6, 6, 6}};
356 0 : return {{0, 1, 2, 5}, {0, 2, 3, 7}, {4, 7, 5, 0}, {5, 7, 6, 2}, {0, 2, 7, 5}};
357 0 : case 5:
358 0 : tet_face_indices = {
359 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}};
360 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}};
361 0 : case 6:
362 0 : tet_face_indices = {
363 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}};
364 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}};
365 0 : default:
366 0 : mooseError("Unexpected input.");
367 : }
368 32280 : }
369 :
370 : void
371 6456 : nodeRotationHEX8(const unsigned int min_id_index,
372 : const unsigned int sec_min_pos,
373 : std::vector<unsigned int> & face_rotation,
374 : std::vector<unsigned int> & node_rotation)
375 : {
376 : // Assuming the original hex element is a cube, the vectors formed by nodes 0-1, 0-3, and 0-4 are
377 : // overlapped with the x, y, and z axes, respectively. sec_min_pos = 0 means the second minimum
378 : // node is in the x direction, sec_min_pos = 1 means the second minimum node is in the y
379 : // direction, and sec_min_pos = 2 means the second minimum node is in the z direction.
380 : const std::vector<std::vector<std::vector<unsigned int>>> preset_indices = {
381 : {{0, 1, 2, 3, 4, 5, 6, 7}, {0, 3, 7, 4, 1, 2, 6, 5}, {0, 4, 5, 1, 3, 7, 6, 2}},
382 : {{1, 0, 4, 5, 2, 3, 7, 6}, {1, 2, 3, 0, 5, 6, 7, 4}, {1, 5, 6, 2, 0, 4, 7, 3}},
383 : {{2, 3, 0, 1, 6, 7, 4, 5}, {2, 1, 5, 6, 3, 0, 4, 7}, {2, 6, 7, 3, 1, 5, 4, 0}},
384 : {{3, 2, 6, 7, 0, 1, 5, 4}, {3, 0, 1, 2, 7, 4, 5, 6}, {3, 7, 4, 0, 2, 6, 5, 1}},
385 : {{4, 5, 1, 0, 7, 6, 2, 3}, {4, 7, 6, 5, 0, 3, 2, 1}, {4, 0, 3, 7, 5, 1, 2, 6}},
386 : {{5, 4, 7, 6, 1, 0, 3, 2}, {5, 6, 2, 1, 4, 7, 3, 0}, {5, 1, 0, 4, 6, 2, 3, 7}},
387 : {{6, 7, 3, 2, 5, 4, 0, 1}, {6, 5, 4, 7, 2, 1, 0, 3}, {6, 2, 1, 5, 7, 3, 0, 4}},
388 213048 : {{7, 6, 5, 4, 3, 2, 1, 0}, {7, 4, 0, 3, 6, 5, 1, 2}, {7, 3, 2, 6, 4, 0, 1, 5}}};
389 :
390 : const std::vector<std::vector<std::vector<unsigned int>>> preset_face_indices = {
391 : {{0, 1, 2, 3, 4, 5}, {4, 0, 3, 5, 1, 2}, {1, 4, 5, 2, 0, 3}},
392 : {{1, 0, 4, 5, 2, 3}, {0, 2, 3, 4, 1, 5}, {2, 1, 5, 3, 0, 4}},
393 : {{0, 3, 4, 1, 2, 5}, {2, 0, 1, 5, 3, 4}, {3, 2, 5, 4, 0, 1}},
394 : {{3, 0, 2, 5, 4, 1}, {0, 4, 1, 2, 3, 5}, {4, 3, 5, 1, 0, 2}},
395 : {{1, 5, 2, 0, 4, 3}, {5, 4, 3, 2, 1, 0}, {4, 1, 0, 3, 5, 2}},
396 : {{5, 1, 4, 3, 2, 0}, {2, 5, 3, 0, 1, 4}, {1, 2, 0, 4, 5, 3}},
397 : {{3, 5, 4, 0, 2, 1}, {5, 2, 1, 4, 3, 0}, {2, 3, 0, 1, 5, 4}},
398 213048 : {{5, 3, 2, 1, 4, 0}, {4, 5, 1, 0, 3, 2}, {3, 4, 0, 2, 5, 1}}};
399 :
400 6456 : if (min_id_index > 7 || sec_min_pos > 2)
401 0 : mooseError("The input node index is out of range.");
402 : else
403 : {
404 : // index: new face index; value: old face index
405 6456 : face_rotation = preset_face_indices[min_id_index][sec_min_pos];
406 6456 : node_rotation = preset_indices[min_id_index][sec_min_pos];
407 : }
408 335712 : }
409 :
410 : void
411 10856 : nodeRotationPRISM6(unsigned int min_id_index,
412 : std::vector<unsigned int> & face_rotation,
413 : std::vector<unsigned int> & node_rotation)
414 : {
415 : const std::vector<std::vector<unsigned int>> preset_indices = {{0, 1, 2, 3, 4, 5},
416 : {1, 2, 0, 4, 5, 3},
417 : {2, 0, 1, 5, 3, 4},
418 : {3, 5, 4, 0, 2, 1},
419 : {4, 3, 5, 1, 0, 2},
420 75992 : {5, 4, 3, 2, 1, 0}};
421 :
422 : const std::vector<std::vector<unsigned int>> preset_face_indices = {{0, 1, 2, 3, 4},
423 : {0, 2, 3, 1, 4},
424 : {0, 3, 1, 2, 4},
425 : {4, 3, 2, 1, 0},
426 : {4, 1, 3, 2, 0},
427 75992 : {4, 2, 1, 3, 0}};
428 :
429 10856 : if (min_id_index > 5)
430 0 : mooseError("The input node index is out of range.");
431 : else
432 : {
433 : // index: new face index; value: old face index
434 10856 : face_rotation = preset_face_indices[min_id_index];
435 10856 : node_rotation = preset_indices[min_id_index];
436 : }
437 43424 : }
438 :
439 : void
440 10856 : prismNodesToTetNodesDeterminer(std::vector<const Node *> & prism_nodes,
441 : std::vector<std::vector<unsigned int>> & rotated_tet_face_indices,
442 : std::vector<std::vector<const Node *>> & tet_nodes_list)
443 : {
444 : // Find the node with the minimum id
445 10856 : std::vector<dof_id_type> node_ids(6);
446 75992 : for (unsigned int i = 0; i < 6; i++)
447 65136 : node_ids[i] = prism_nodes[i]->id();
448 :
449 10856 : const unsigned int min_node_id_index = std::distance(
450 10856 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
451 :
452 : // Rotate the node and face indices based on the identified minimum node
453 : // After the rotation, we guarantee that the minimum node is the first node (Node 0)
454 : // This makes the splitting process simpler
455 10856 : std::vector<unsigned int> face_rotation;
456 10856 : std::vector<unsigned int> rotated_indices;
457 10856 : nodeRotationPRISM6(min_node_id_index, face_rotation, rotated_indices);
458 10856 : std::vector<const Node *> rotated_prism_nodes;
459 75992 : for (unsigned int i = 0; i < 6; i++)
460 65136 : rotated_prism_nodes.push_back(prism_nodes[rotated_indices[i]]);
461 :
462 10856 : std::vector<const Node *> key_quad_nodes = {rotated_prism_nodes[1],
463 10856 : rotated_prism_nodes[2],
464 10856 : rotated_prism_nodes[5],
465 10856 : rotated_prism_nodes[4]};
466 :
467 : // Find the selection of each face's cutting direction
468 10856 : const bool diagonal_direction = quadFaceDiagonalDirection(key_quad_nodes);
469 :
470 : // Based on the determined splitting directions of all the faces, determine the nodes of each
471 : // resulting TET4 elements after the splitting.
472 10856 : std::vector<std::vector<unsigned int>> tet_face_indices;
473 10856 : const auto tet_nodes_set = tetNodesForPrism(diagonal_direction, tet_face_indices);
474 43424 : for (const auto & tet_face_index : tet_face_indices)
475 : {
476 32568 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
477 162840 : for (const auto & face_index : tet_face_index)
478 : {
479 130272 : if (face_index < 5)
480 86848 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
481 : else
482 43424 : rotated_tet_face_indices.back().push_back(5);
483 : }
484 : }
485 :
486 43424 : for (const auto & tet_nodes : tet_nodes_set)
487 : {
488 32568 : tet_nodes_list.push_back(std::vector<const Node *>());
489 162840 : for (const auto & tet_node : tet_nodes)
490 130272 : tet_nodes_list.back().push_back(rotated_prism_nodes[tet_node]);
491 : }
492 10856 : }
493 :
494 : std::vector<std::vector<unsigned int>>
495 10856 : tetNodesForPrism(const bool diagonal_direction,
496 : std::vector<std::vector<unsigned int>> & tet_face_indices)
497 : {
498 :
499 10856 : if (diagonal_direction)
500 : {
501 25920 : tet_face_indices = {{4, 3, 5, 1}, {2, 1, 5, 5}, {2, 5, 3, 0}};
502 25920 : return {{3, 5, 4, 0}, {1, 4, 5, 0}, {1, 5, 2, 0}};
503 : }
504 : else
505 : {
506 17504 : tet_face_indices = {{4, 3, 5, 1}, {2, 1, 5, 0}, {2, 5, 5, 3}};
507 17504 : return {{3, 5, 4, 0}, {1, 4, 2, 0}, {2, 4, 5, 0}};
508 : }
509 32568 : }
510 :
511 : void
512 720 : nodeRotationPYRAMID5(unsigned int min_id_index,
513 : std::vector<unsigned int> & face_rotation,
514 : std::vector<unsigned int> & node_rotation)
515 : {
516 : const std::vector<std::vector<unsigned int>> preset_indices = {
517 3600 : {0, 1, 2, 3, 4}, {1, 2, 3, 0, 4}, {2, 3, 0, 1, 4}, {3, 0, 1, 2, 4}};
518 :
519 : const std::vector<std::vector<unsigned int>> preset_face_indices = {
520 3600 : {0, 1, 2, 3, 4}, {1, 2, 3, 0, 4}, {2, 3, 0, 1, 4}, {3, 0, 1, 2, 4}};
521 :
522 720 : if (min_id_index > 3)
523 0 : mooseError("The input node index is out of range.");
524 : else
525 : {
526 : // index: new face index; value: old face index
527 720 : face_rotation = preset_face_indices[min_id_index];
528 720 : node_rotation = preset_indices[min_id_index];
529 : }
530 2880 : }
531 :
532 : void
533 720 : pyramidNodesToTetNodesDeterminer(std::vector<const Node *> & pyramid_nodes,
534 : std::vector<std::vector<unsigned int>> & rotated_tet_face_indices,
535 : std::vector<std::vector<const Node *>> & tet_nodes_list)
536 : {
537 : // Find the node with the minimum id, ignoring the top node
538 720 : std::vector<dof_id_type> node_ids(4);
539 3600 : for (unsigned int i = 0; i < 4; i++)
540 2880 : node_ids[i] = pyramid_nodes[i]->id();
541 :
542 720 : const unsigned int min_node_id_index = std::distance(
543 720 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
544 :
545 : // Rotate the node and face indices based on the identified minimum nodes
546 : // After the rotation, we guarantee that the minimum node is the first node (Node 0)
547 : // This makes the splitting process simpler
548 720 : std::vector<unsigned int> face_rotation;
549 720 : std::vector<unsigned int> rotated_indices;
550 720 : nodeRotationPYRAMID5(min_node_id_index, face_rotation, rotated_indices);
551 720 : std::vector<const Node *> rotated_pyramid_nodes;
552 4320 : for (unsigned int i = 0; i < 5; i++)
553 3600 : rotated_pyramid_nodes.push_back(pyramid_nodes[rotated_indices[i]]);
554 :
555 : // There is only one quad face in a pyramid element, so the splitting selection is binary
556 2160 : const std::vector<std::vector<unsigned int>> tet_nodes_set = {{0, 1, 2, 4}, {0, 2, 3, 4}};
557 2160 : const std::vector<std::vector<unsigned int>> tet_face_indices = {{4, 0, 1, 5}, {4, 5, 2, 3}};
558 :
559 : // Based on the determined splitting direction, determine the nodes of each resulting TET4
560 : // elements after the splitting.
561 2160 : for (const auto & tet_face_index : tet_face_indices)
562 : {
563 1440 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
564 7200 : for (const auto & face_index : tet_face_index)
565 : {
566 5760 : if (face_index < 5)
567 4320 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
568 : else
569 1440 : rotated_tet_face_indices.back().push_back(5);
570 : }
571 : }
572 :
573 2160 : for (const auto & tet_nodes : tet_nodes_set)
574 : {
575 1440 : tet_nodes_list.push_back(std::vector<const Node *>());
576 7200 : for (const auto & tet_node : tet_nodes)
577 5760 : tet_nodes_list.back().push_back(rotated_pyramid_nodes[tet_node]);
578 : }
579 2880 : }
580 :
581 : void
582 144 : convert3DMeshToAllTet4(ReplicatedMesh & mesh,
583 : const std::vector<std::pair<dof_id_type, bool>> & elems_to_process,
584 : std::vector<dof_id_type> & converted_elems_ids_to_track,
585 : const subdomain_id_type block_id_to_remove,
586 : const bool delete_block_to_remove)
587 : {
588 144 : std::vector<dof_id_type> converted_elems_ids_to_retain;
589 : // Build boundary information of the mesh
590 144 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
591 144 : const auto bdry_side_list = boundary_info.build_side_list();
592 9088 : for (const auto & elem_to_process : elems_to_process)
593 : {
594 8944 : switch (mesh.elem_ptr(elem_to_process.first)->type())
595 : {
596 6456 : case ElemType::HEX8:
597 6456 : hexElemSplitter(mesh,
598 : bdry_side_list,
599 6456 : elem_to_process.first,
600 6456 : elem_to_process.second ? converted_elems_ids_to_track
601 : : converted_elems_ids_to_retain);
602 6456 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
603 6456 : break;
604 512 : case ElemType::PYRAMID5:
605 512 : pyramidElemSplitter(mesh,
606 : bdry_side_list,
607 512 : elem_to_process.first,
608 512 : elem_to_process.second ? converted_elems_ids_to_track
609 : : converted_elems_ids_to_retain);
610 512 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
611 512 : break;
612 696 : case ElemType::PRISM6:
613 696 : prismElemSplitter(mesh,
614 : bdry_side_list,
615 696 : elem_to_process.first,
616 696 : elem_to_process.second ? converted_elems_ids_to_track
617 : : converted_elems_ids_to_retain);
618 696 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
619 696 : break;
620 1280 : case ElemType::TET4:
621 1280 : if (elem_to_process.second)
622 512 : converted_elems_ids_to_track.push_back(elem_to_process.first);
623 : else
624 768 : converted_elems_ids_to_retain.push_back(elem_to_process.first);
625 1280 : break;
626 0 : default:
627 0 : mooseError("Unexpected element type.");
628 : }
629 : }
630 :
631 144 : if (delete_block_to_remove)
632 : {
633 1344 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
634 1344 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
635 : elem_it++)
636 1344 : mesh.delete_elem(*elem_it);
637 :
638 80 : mesh.contract();
639 80 : mesh.prepare_for_use();
640 : }
641 144 : }
642 :
643 : void
644 84 : convert3DMeshToAllTet4(ReplicatedMesh & mesh)
645 : {
646 : // Subdomain ID for new utility blocks must be new
647 84 : std::set<subdomain_id_type> subdomain_ids_set;
648 84 : mesh.subdomain_ids(subdomain_ids_set);
649 84 : const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin();
650 84 : const subdomain_id_type block_id_to_remove = max_subdomain_id + 1;
651 84 : std::vector<std::pair<dof_id_type, bool>> original_elems;
652 :
653 2116 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
654 : elem_it++)
655 : {
656 2036 : if ((*elem_it)->default_order() != Order::FIRST)
657 4 : mooseError("Only first order elements are supported for cutting.");
658 2032 : original_elems.push_back(std::make_pair((*elem_it)->id(), false));
659 80 : }
660 :
661 80 : std::vector<dof_id_type> converted_elems_ids_to_track;
662 :
663 80 : convert3DMeshToAllTet4(
664 : mesh, original_elems, converted_elems_ids_to_track, block_id_to_remove, true);
665 80 : }
666 :
667 : void
668 33424 : elementBoundaryInfoCollector(const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
669 : const dof_id_type elem_id,
670 : const unsigned short n_elem_sides,
671 : std::vector<std::vector<boundary_id_type>> & elem_side_list)
672 : {
673 33424 : elem_side_list.resize(n_elem_sides);
674 : const auto selected_bdry_side_list =
675 33424 : std::equal_range(bdry_side_list.begin(), bdry_side_list.end(), elem_id, BCTupleKeyComp{});
676 33424 : for (auto selected_bdry_side = selected_bdry_side_list.first;
677 44600 : selected_bdry_side != selected_bdry_side_list.second;
678 11176 : ++selected_bdry_side)
679 : {
680 11176 : elem_side_list[std::get<1>(*selected_bdry_side)].push_back(std::get<2>(*selected_bdry_side));
681 : }
682 33424 : }
683 : }
|
