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