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 7247 : 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 7247 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
36 : // Create a list of sidesets involving the element to be split
37 7247 : std::vector<std::vector<boundary_id_type>> elem_side_list;
38 7247 : elem_side_list.resize(6);
39 7247 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 6, elem_side_list);
40 :
41 7247 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
42 7247 : 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 7247 : 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 7247 : std::vector<std::vector<unsigned int>> opt_option;
48 7247 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
49 7247 : mesh.elem_ptr(elem_id)->node_ptr(1),
50 7247 : mesh.elem_ptr(elem_id)->node_ptr(2),
51 7247 : mesh.elem_ptr(elem_id)->node_ptr(3),
52 7247 : mesh.elem_ptr(elem_id)->node_ptr(4),
53 7247 : mesh.elem_ptr(elem_id)->node_ptr(5),
54 7247 : mesh.elem_ptr(elem_id)->node_ptr(6),
55 50729 : mesh.elem_ptr(elem_id)->node_ptr(7)};
56 7247 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
57 :
58 7247 : std::vector<std::vector<const Node *>> optimized_node_list;
59 7247 : hexNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
60 :
61 7247 : std::vector<Elem *> elems_Tet4;
62 50729 : for (const auto i : index_range(optimized_node_list))
63 : {
64 43482 : auto new_elem = std::make_unique<Tet4>();
65 43482 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
66 43482 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
67 43482 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
68 43482 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
69 43482 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
70 43482 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
71 43482 : converted_elems_ids.push_back(elems_Tet4.back()->id());
72 :
73 217410 : 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 173928 : if (rotated_tet_face_indices[i][j] < 6)
79 : {
80 97056 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
81 10092 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
82 : }
83 : }
84 43482 : }
85 :
86 : // Retain element extra integers
87 50729 : for (unsigned int i = 0; i < 6; i++)
88 43482 : 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 7247 : }
93 :
94 : void
95 783 : 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 783 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
102 : // Create a list of sidesets involving the element to be split
103 783 : std::vector<std::vector<boundary_id_type>> elem_side_list;
104 783 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 5, elem_side_list);
105 :
106 783 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
107 783 : 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 783 : 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 783 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
114 783 : mesh.elem_ptr(elem_id)->node_ptr(1),
115 783 : mesh.elem_ptr(elem_id)->node_ptr(2),
116 783 : mesh.elem_ptr(elem_id)->node_ptr(3),
117 783 : mesh.elem_ptr(elem_id)->node_ptr(4),
118 3915 : mesh.elem_ptr(elem_id)->node_ptr(5)};
119 783 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
120 783 : std::vector<std::vector<const Node *>> optimized_node_list;
121 783 : prismNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
122 :
123 783 : std::vector<Elem *> elems_Tet4;
124 3132 : for (const auto i : index_range(optimized_node_list))
125 : {
126 2349 : auto new_elem = std::make_unique<Tet4>();
127 2349 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
128 2349 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
129 2349 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
130 2349 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
131 2349 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
132 2349 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
133 2349 : converted_elems_ids.push_back(elems_Tet4.back()->id());
134 :
135 11745 : 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 9396 : if (rotated_tet_face_indices[i][j] < 5)
141 : {
142 7452 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
143 1188 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
144 : }
145 : }
146 2349 : }
147 :
148 : // Retain element extra integers
149 3132 : for (unsigned int i = 0; i < 3; i++)
150 2349 : 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 783 : }
155 :
156 : void
157 576 : 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 576 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
164 : // Create a list of sidesets involving the element to be split
165 576 : std::vector<std::vector<boundary_id_type>> elem_side_list;
166 576 : elementBoundaryInfoCollector(bdry_side_list, elem_id, 5, elem_side_list);
167 :
168 576 : const unsigned int n_elem_extra_ids = mesh.n_elem_integers();
169 576 : 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 576 : 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 576 : std::vector<const Node *> elem_node_list = {mesh.elem_ptr(elem_id)->node_ptr(0),
175 576 : mesh.elem_ptr(elem_id)->node_ptr(1),
176 576 : mesh.elem_ptr(elem_id)->node_ptr(2),
177 576 : mesh.elem_ptr(elem_id)->node_ptr(3),
178 2304 : mesh.elem_ptr(elem_id)->node_ptr(4)};
179 576 : std::vector<std::vector<unsigned int>> rotated_tet_face_indices;
180 576 : std::vector<std::vector<const Node *>> optimized_node_list;
181 576 : pyramidNodesToTetNodesDeterminer(elem_node_list, rotated_tet_face_indices, optimized_node_list);
182 :
183 576 : std::vector<Elem *> elems_Tet4;
184 1728 : for (const auto i : index_range(optimized_node_list))
185 : {
186 1152 : auto new_elem = std::make_unique<Tet4>();
187 1152 : new_elem->set_node(0, const_cast<Node *>(optimized_node_list[i][0]));
188 1152 : new_elem->set_node(1, const_cast<Node *>(optimized_node_list[i][1]));
189 1152 : new_elem->set_node(2, const_cast<Node *>(optimized_node_list[i][2]));
190 1152 : new_elem->set_node(3, const_cast<Node *>(optimized_node_list[i][3]));
191 1152 : new_elem->subdomain_id() = mesh.elem_ptr(elem_id)->subdomain_id();
192 1152 : elems_Tet4.push_back(mesh.add_elem(std::move(new_elem)));
193 1152 : converted_elems_ids.push_back(elems_Tet4.back()->id());
194 :
195 5760 : 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 4608 : if (rotated_tet_face_indices[i][j] < 5)
201 : {
202 4284 : for (const auto & side_info : elem_side_list[rotated_tet_face_indices[i][j]])
203 828 : boundary_info.add_side(elems_Tet4.back(), j, side_info);
204 : }
205 : }
206 1152 : }
207 :
208 : // Retain element extra integers
209 1728 : for (unsigned int i = 0; i < 2; i++)
210 1152 : 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 576 : }
215 :
216 : std::vector<unsigned int>
217 7247 : neighborNodeIndicesHEX8(unsigned int min_id_index)
218 : {
219 : const std::vector<std::vector<unsigned int>> preset_indices = {
220 65223 : {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 7247 : if (min_id_index > 7)
222 0 : mooseError("The input node index is out of range.");
223 : else
224 14494 : return preset_indices[min_id_index];
225 21741 : }
226 :
227 : void
228 7247 : 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 7247 : std::vector<dof_id_type> node_ids(8);
234 65223 : for (unsigned int i = 0; i < 8; i++)
235 57976 : node_ids[i] = hex_nodes[i]->id();
236 :
237 7247 : const unsigned int min_node_id_index = std::distance(
238 7247 : 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 7247 : const auto neighbor_node_indices = neighborNodeIndicesHEX8(min_node_id_index);
243 :
244 7247 : const auto neighbor_node_ids = {node_ids[neighbor_node_indices[0]],
245 7247 : node_ids[neighbor_node_indices[1]],
246 14494 : node_ids[neighbor_node_indices[2]]};
247 : const unsigned int sec_min_pos =
248 7247 : std::distance(std::begin(neighbor_node_ids),
249 7247 : 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 7247 : std::vector<unsigned int> face_rotation;
256 7247 : std::vector<unsigned int> rotated_indices;
257 7247 : nodeRotationHEX8(min_node_id_index, sec_min_pos, face_rotation, rotated_indices);
258 7247 : std::vector<const Node *> rotated_hex_nodes;
259 65223 : for (unsigned int i = 0; i < 8; i++)
260 57976 : rotated_hex_nodes.push_back(hex_nodes[rotated_indices[i]]);
261 :
262 : // Find the selection of each face's cutting direction
263 7247 : 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 7247 : std::vector<std::vector<unsigned int>> tet_face_indices;
268 7247 : const auto tet_nodes_set = tetNodesForHex(diagonal_directions, tet_face_indices);
269 50729 : for (const auto & tet_face_index : tet_face_indices)
270 : {
271 43482 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
272 217410 : for (const auto & face_index : tet_face_index)
273 : {
274 173928 : if (face_index < 6)
275 86964 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
276 : else
277 86964 : rotated_tet_face_indices.back().push_back(6);
278 : }
279 : }
280 :
281 50729 : for (const auto & tet_nodes : tet_nodes_set)
282 : {
283 43482 : tet_nodes_list.push_back(std::vector<const Node *>());
284 217410 : for (const auto & tet_node : tet_nodes)
285 173928 : tet_nodes_list.back().push_back(rotated_hex_nodes[tet_node]);
286 : }
287 7247 : }
288 :
289 : std::vector<bool>
290 7247 : 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 50729 : {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 7247 : std::vector<bool> diagonal_directions;
296 50729 : for (const auto & face_index : face_indices)
297 : {
298 43482 : std::vector<const Node *> quad_nodes = {hex_nodes[face_index[0]],
299 43482 : hex_nodes[face_index[1]],
300 43482 : hex_nodes[face_index[2]],
301 130446 : hex_nodes[face_index[3]]};
302 43482 : diagonal_directions.push_back(quadFaceDiagonalDirection(quad_nodes));
303 43482 : }
304 14494 : return diagonal_directions;
305 21741 : }
306 :
307 : bool
308 55695 : quadFaceDiagonalDirection(const std::vector<const Node *> & quad_nodes)
309 : {
310 : const std::vector<dof_id_type> node_ids = {
311 55695 : quad_nodes[0]->id(), quad_nodes[1]->id(), quad_nodes[2]->id(), quad_nodes[3]->id()};
312 55695 : const unsigned int min_id_index = std::distance(
313 55695 : std::begin(node_ids), std::min_element(std::begin(node_ids), std::end(node_ids)));
314 55695 : if (min_id_index == 0 || min_id_index == 2)
315 36616 : return true;
316 : else
317 19079 : return false;
318 55695 : }
319 :
320 : std::vector<std::vector<unsigned int>>
321 7247 : 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 57976 : {true, false, true, false, false, false}};
331 :
332 7247 : const unsigned int input_index = std::distance(
333 : std::begin(possible_inputs),
334 7247 : std::find(std::begin(possible_inputs), std::end(possible_inputs), diagonal_directions));
335 :
336 7247 : switch (input_index)
337 : {
338 338 : case 0:
339 2366 : tet_face_indices = {
340 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}};
341 2366 : 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 6909 : case 2:
347 48363 : tet_face_indices = {
348 48363 : {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 48363 : 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 36235 : }
369 :
370 : void
371 7247 : 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 239151 : {{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 239151 : {{5, 3, 2, 1, 4, 0}, {4, 5, 1, 0, 3, 2}, {3, 4, 0, 2, 5, 1}}};
399 :
400 7247 : 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 7247 : face_rotation = preset_face_indices[min_id_index][sec_min_pos];
406 7247 : node_rotation = preset_indices[min_id_index][sec_min_pos];
407 : }
408 376844 : }
409 :
410 : void
411 12213 : 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 85491 : {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 85491 : {4, 2, 1, 3, 0}};
428 :
429 12213 : 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 12213 : face_rotation = preset_face_indices[min_id_index];
435 12213 : node_rotation = preset_indices[min_id_index];
436 : }
437 48852 : }
438 :
439 : void
440 12213 : 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 12213 : std::vector<dof_id_type> node_ids(6);
446 85491 : for (unsigned int i = 0; i < 6; i++)
447 73278 : node_ids[i] = prism_nodes[i]->id();
448 :
449 12213 : const unsigned int min_node_id_index = std::distance(
450 12213 : 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 12213 : std::vector<unsigned int> face_rotation;
456 12213 : std::vector<unsigned int> rotated_indices;
457 12213 : nodeRotationPRISM6(min_node_id_index, face_rotation, rotated_indices);
458 12213 : std::vector<const Node *> rotated_prism_nodes;
459 85491 : for (unsigned int i = 0; i < 6; i++)
460 73278 : rotated_prism_nodes.push_back(prism_nodes[rotated_indices[i]]);
461 :
462 12213 : std::vector<const Node *> key_quad_nodes = {rotated_prism_nodes[1],
463 12213 : rotated_prism_nodes[2],
464 12213 : rotated_prism_nodes[5],
465 12213 : rotated_prism_nodes[4]};
466 :
467 : // Find the selection of each face's cutting direction
468 12213 : 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 12213 : std::vector<std::vector<unsigned int>> tet_face_indices;
473 12213 : const auto tet_nodes_set = tetNodesForPrism(diagonal_direction, tet_face_indices);
474 48852 : for (const auto & tet_face_index : tet_face_indices)
475 : {
476 36639 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
477 183195 : for (const auto & face_index : tet_face_index)
478 : {
479 146556 : if (face_index < 5)
480 97704 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
481 : else
482 48852 : rotated_tet_face_indices.back().push_back(5);
483 : }
484 : }
485 :
486 48852 : for (const auto & tet_nodes : tet_nodes_set)
487 : {
488 36639 : tet_nodes_list.push_back(std::vector<const Node *>());
489 183195 : for (const auto & tet_node : tet_nodes)
490 146556 : tet_nodes_list.back().push_back(rotated_prism_nodes[tet_node]);
491 : }
492 12213 : }
493 :
494 : std::vector<std::vector<unsigned int>>
495 12213 : tetNodesForPrism(const bool diagonal_direction,
496 : std::vector<std::vector<unsigned int>> & tet_face_indices)
497 : {
498 :
499 12213 : if (diagonal_direction)
500 : {
501 29160 : tet_face_indices = {{4, 3, 5, 1}, {2, 1, 5, 5}, {2, 5, 3, 0}};
502 29160 : return {{3, 5, 4, 0}, {1, 4, 5, 0}, {1, 5, 2, 0}};
503 : }
504 : else
505 : {
506 19692 : tet_face_indices = {{4, 3, 5, 1}, {2, 1, 5, 0}, {2, 5, 5, 3}};
507 19692 : return {{3, 5, 4, 0}, {1, 4, 2, 0}, {2, 4, 5, 0}};
508 : }
509 36639 : }
510 :
511 : void
512 810 : 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 4050 : {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 4050 : {0, 1, 2, 3, 4}, {1, 2, 3, 0, 4}, {2, 3, 0, 1, 4}, {3, 0, 1, 2, 4}};
521 :
522 810 : 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 810 : face_rotation = preset_face_indices[min_id_index];
528 810 : node_rotation = preset_indices[min_id_index];
529 : }
530 3240 : }
531 :
532 : void
533 810 : 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 810 : std::vector<dof_id_type> node_ids(4);
539 4050 : for (unsigned int i = 0; i < 4; i++)
540 3240 : node_ids[i] = pyramid_nodes[i]->id();
541 :
542 810 : const unsigned int min_node_id_index = std::distance(
543 810 : 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 810 : std::vector<unsigned int> face_rotation;
549 810 : std::vector<unsigned int> rotated_indices;
550 810 : nodeRotationPYRAMID5(min_node_id_index, face_rotation, rotated_indices);
551 810 : std::vector<const Node *> rotated_pyramid_nodes;
552 4860 : for (unsigned int i = 0; i < 5; i++)
553 4050 : 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 2430 : const std::vector<std::vector<unsigned int>> tet_nodes_set = {{0, 1, 2, 4}, {0, 2, 3, 4}};
557 2430 : 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 2430 : for (const auto & tet_face_index : tet_face_indices)
562 : {
563 1620 : rotated_tet_face_indices.push_back(std::vector<unsigned int>());
564 8100 : for (const auto & face_index : tet_face_index)
565 : {
566 6480 : if (face_index < 5)
567 4860 : rotated_tet_face_indices.back().push_back(face_rotation[face_index]);
568 : else
569 1620 : rotated_tet_face_indices.back().push_back(5);
570 : }
571 : }
572 :
573 2430 : for (const auto & tet_nodes : tet_nodes_set)
574 : {
575 1620 : tet_nodes_list.push_back(std::vector<const Node *>());
576 8100 : for (const auto & tet_node : tet_nodes)
577 6480 : tet_nodes_list.back().push_back(rotated_pyramid_nodes[tet_node]);
578 : }
579 3240 : }
580 :
581 : void
582 160 : 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 160 : std::vector<dof_id_type> converted_elems_ids_to_retain;
589 : // Build boundary information of the mesh
590 160 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
591 160 : const auto bdry_side_list = boundary_info.build_side_list();
592 10206 : for (const auto & elem_to_process : elems_to_process)
593 : {
594 10046 : switch (mesh.elem_ptr(elem_to_process.first)->type())
595 : {
596 7247 : case ElemType::HEX8:
597 7247 : hexElemSplitter(mesh,
598 : bdry_side_list,
599 7247 : elem_to_process.first,
600 7247 : elem_to_process.second ? converted_elems_ids_to_track
601 : : converted_elems_ids_to_retain);
602 7247 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
603 7247 : break;
604 576 : case ElemType::PYRAMID5:
605 576 : pyramidElemSplitter(mesh,
606 : bdry_side_list,
607 576 : elem_to_process.first,
608 576 : elem_to_process.second ? converted_elems_ids_to_track
609 : : converted_elems_ids_to_retain);
610 576 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
611 576 : break;
612 783 : case ElemType::PRISM6:
613 783 : prismElemSplitter(mesh,
614 : bdry_side_list,
615 783 : elem_to_process.first,
616 783 : elem_to_process.second ? converted_elems_ids_to_track
617 : : converted_elems_ids_to_retain);
618 783 : mesh.elem_ptr(elem_to_process.first)->subdomain_id() = block_id_to_remove;
619 783 : break;
620 1440 : case ElemType::TET4:
621 1440 : if (elem_to_process.second)
622 576 : converted_elems_ids_to_track.push_back(elem_to_process.first);
623 : else
624 864 : converted_elems_ids_to_retain.push_back(elem_to_process.first);
625 1440 : break;
626 0 : default:
627 0 : mooseError("Unexpected element type.");
628 : }
629 : }
630 :
631 160 : if (delete_block_to_remove)
632 : {
633 1494 : for (auto elem_it = mesh.active_subdomain_elements_begin(block_id_to_remove);
634 1494 : elem_it != mesh.active_subdomain_elements_end(block_id_to_remove);
635 : elem_it++)
636 1494 : mesh.delete_elem(*elem_it);
637 :
638 88 : mesh.contract();
639 88 : mesh.prepare_for_use();
640 : }
641 160 : }
642 :
643 : void
644 92 : convert3DMeshToAllTet4(ReplicatedMesh & mesh)
645 : {
646 : // Subdomain ID for new utility blocks must be new
647 92 : std::set<subdomain_id_type> subdomain_ids_set;
648 92 : mesh.subdomain_ids(subdomain_ids_set);
649 92 : const subdomain_id_type max_subdomain_id = *subdomain_ids_set.rbegin();
650 92 : const subdomain_id_type block_id_to_remove = max_subdomain_id + 1;
651 92 : std::vector<std::pair<dof_id_type, bool>> original_elems;
652 :
653 2362 : for (auto elem_it = mesh.active_elements_begin(); elem_it != mesh.active_elements_end();
654 : elem_it++)
655 : {
656 2274 : if ((*elem_it)->default_order() != Order::FIRST)
657 4 : mooseError("Only first order elements are supported for cutting.");
658 2270 : original_elems.push_back(std::make_pair((*elem_it)->id(), false));
659 88 : }
660 :
661 88 : std::vector<dof_id_type> converted_elems_ids_to_track;
662 :
663 88 : convert3DMeshToAllTet4(
664 : mesh, original_elems, converted_elems_ids_to_track, block_id_to_remove, true);
665 88 : }
666 :
667 : void
668 37586 : 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 37586 : elem_side_list.resize(n_elem_sides);
674 : const auto selected_bdry_side_list =
675 37586 : std::equal_range(bdry_side_list.begin(), bdry_side_list.end(), elem_id, BCTupleKeyComp{});
676 37586 : for (auto selected_bdry_side = selected_bdry_side_list.first;
677 50111 : selected_bdry_side != selected_bdry_side_list.second;
678 12525 : ++selected_bdry_side)
679 : {
680 12525 : elem_side_list[std::get<1>(*selected_bdry_side)].push_back(std::get<2>(*selected_bdry_side));
681 : }
682 37586 : }
683 : }
|
