Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // Local includes
20 : #include "libmesh/cell_pyramid18.h"
21 : #include "libmesh/edge_edge3.h"
22 : #include "libmesh/face_tri7.h"
23 : #include "libmesh/face_quad9.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 :
32 :
33 : // ------------------------------------------------------------
34 : // Pyramid18 class static member initializations
35 : const int Pyramid18::num_nodes;
36 : const int Pyramid18::nodes_per_side;
37 : const int Pyramid18::nodes_per_edge;
38 :
39 : const unsigned int Pyramid18::side_nodes_map[Pyramid18::num_sides][Pyramid18::nodes_per_side] =
40 : {
41 : {0, 1, 4, 5, 10, 9, 14, 99, 99}, // Side 0 (front)
42 : {1, 2, 4, 6, 11, 10, 15, 99, 99}, // Side 1 (right)
43 : {2, 3, 4, 7, 12, 11, 16, 99, 99}, // Side 2 (back)
44 : {3, 0, 4, 8, 9, 12, 17, 99, 99}, // Side 3 (left)
45 : {0, 3, 2, 1, 8, 7, 6, 5, 13} // Side 4 (base)
46 : };
47 :
48 : const unsigned int Pyramid18::edge_nodes_map[Pyramid18::num_edges][Pyramid18::nodes_per_edge] =
49 : {
50 : {0, 1, 5}, // Edge 0
51 : {1, 2, 6}, // Edge 1
52 : {2, 3, 7}, // Edge 2
53 : {0, 3, 8}, // Edge 3
54 : {0, 4, 9}, // Edge 4
55 : {1, 4, 10}, // Edge 5
56 : {2, 4, 11}, // Edge 6
57 : {3, 4, 12} // Edge 7
58 : };
59 :
60 : // ------------------------------------------------------------
61 : // Pyramid18 class member functions
62 :
63 151632 : bool Pyramid18::is_vertex(const unsigned int i) const
64 : {
65 151632 : if (i < 5)
66 42120 : return true;
67 6552 : return false;
68 : }
69 :
70 :
71 :
72 44928 : bool Pyramid18::is_edge(const unsigned int i) const
73 : {
74 44928 : if (i < 5)
75 0 : return false;
76 44928 : if (i > 12)
77 17280 : return false;
78 2304 : return true;
79 : }
80 :
81 :
82 :
83 110181 : bool Pyramid18::is_face(const unsigned int i) const
84 : {
85 110181 : if (i > 12)
86 22899 : return true;
87 7020 : return false;
88 : }
89 :
90 :
91 :
92 67131 : bool Pyramid18::is_node_on_side(const unsigned int n,
93 : const unsigned int s) const
94 : {
95 5556 : libmesh_assert_less (s, n_sides());
96 5556 : return std::find(std::begin(side_nodes_map[s]),
97 67131 : std::end(side_nodes_map[s]),
98 67131 : n) != std::end(side_nodes_map[s]);
99 : }
100 :
101 : std::vector<unsigned>
102 32832 : Pyramid18::nodes_on_side(const unsigned int s) const
103 : {
104 2736 : libmesh_assert_less(s, n_sides());
105 32832 : auto trim = (s == 4) ? 0 : 2;
106 32832 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
107 : }
108 :
109 : std::vector<unsigned>
110 27648 : Pyramid18::nodes_on_edge(const unsigned int e) const
111 : {
112 2304 : libmesh_assert_less(e, n_edges());
113 29952 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
114 : }
115 :
116 13824 : bool Pyramid18::is_node_on_edge(const unsigned int n,
117 : const unsigned int e) const
118 : {
119 1152 : libmesh_assert_less (e, n_edges());
120 1152 : return std::find(std::begin(edge_nodes_map[e]),
121 13824 : std::end(edge_nodes_map[e]),
122 13824 : n) != std::end(edge_nodes_map[e]);
123 : }
124 :
125 :
126 :
127 67618 : bool Pyramid18::has_affine_map() const
128 : {
129 : // TODO: If the base is a parallelogram and all the triangular faces are planar,
130 : // the map should be linear, but I need to test this theory...
131 67618 : return false;
132 : }
133 :
134 :
135 :
136 9922045 : Order Pyramid18::default_order() const
137 : {
138 9922045 : return THIRD;
139 : }
140 :
141 :
142 :
143 0 : dof_id_type Pyramid18::key (const unsigned int s) const
144 : {
145 0 : libmesh_assert_less (s, this->n_sides());
146 :
147 0 : switch (s)
148 : {
149 0 : case 0: // triangular face 1
150 : case 1: // triangular face 2
151 : case 2: // triangular face 3
152 : case 3: // triangular face 4
153 0 : return this->compute_key (this->node_id(s+14));
154 :
155 0 : case 4: // the quad face at z=0
156 0 : return this->compute_key (this->node_id(13));
157 :
158 0 : default:
159 0 : libmesh_error_msg("Invalid side s = " << s);
160 : }
161 : }
162 :
163 :
164 :
165 355 : unsigned int Pyramid18::local_side_node(unsigned int side,
166 : unsigned int side_node) const
167 : {
168 10 : libmesh_assert_less (side, this->n_sides());
169 :
170 : // Never more than 9 nodes per side.
171 10 : libmesh_assert_less(side_node, Pyramid18::nodes_per_side);
172 :
173 : // Some sides have 7 nodes.
174 10 : libmesh_assert(side == 4 || side_node < 7);
175 :
176 355 : return Pyramid18::side_nodes_map[side][side_node];
177 : }
178 :
179 :
180 :
181 120092 : unsigned int Pyramid18::local_edge_node(unsigned int edge,
182 : unsigned int edge_node) const
183 : {
184 9992 : libmesh_assert_less(edge, this->n_edges());
185 9992 : libmesh_assert_less(edge_node, Pyramid18::nodes_per_edge);
186 :
187 120092 : return Pyramid18::edge_nodes_map[edge][edge_node];
188 : }
189 :
190 :
191 :
192 280114 : std::unique_ptr<Elem> Pyramid18::build_side_ptr (const unsigned int i)
193 : {
194 9796 : libmesh_assert_less (i, this->n_sides());
195 :
196 280114 : std::unique_ptr<Elem> face;
197 :
198 280114 : switch (i)
199 : {
200 142649 : case 0: // triangular face 1
201 : case 1: // triangular face 2
202 : case 2: // triangular face 3
203 : case 3: // triangular face 4
204 : {
205 142649 : face = std::make_unique<Tri7>();
206 142649 : break;
207 : }
208 137465 : case 4: // the quad face at z=0
209 : {
210 137465 : face = std::make_unique<Quad9>();
211 137465 : break;
212 : }
213 0 : default:
214 0 : libmesh_error_msg("Invalid side i = " << i);
215 : }
216 :
217 : // Set the nodes
218 2515842 : for (auto n : face->node_index_range())
219 2313664 : face->set_node(n, this->node_ptr(Pyramid18::side_nodes_map[i][n]));
220 :
221 280114 : face->set_interior_parent(this);
222 270318 : face->inherit_data_from(*this);
223 :
224 280114 : return face;
225 0 : }
226 :
227 :
228 :
229 549000 : void Pyramid18::build_side_ptr (std::unique_ptr<Elem> & side,
230 : const unsigned int i)
231 : {
232 18672 : libmesh_assert_less (i, this->n_sides());
233 :
234 549000 : switch (i)
235 : {
236 14040 : case 0: // triangular face 1
237 : case 1: // triangular face 2
238 : case 2: // triangular face 3
239 : case 3: // triangular face 4
240 : {
241 412182 : if (!side.get() || side->type() != TRI7)
242 : {
243 1194 : side = this->build_side_ptr(i);
244 647 : return;
245 : }
246 13990 : break;
247 : }
248 4632 : case 4: // the quad face at z=0
249 : {
250 136818 : if (!side.get() || side->type() != QUAD9)
251 : {
252 262398 : side = this->build_side_ptr(i);
253 135737 : return;
254 : }
255 94 : break;
256 : }
257 0 : default:
258 0 : libmesh_error_msg("Invalid side i = " << i);
259 : }
260 :
261 398532 : side->inherit_data_from(*this);
262 :
263 : // Set the nodes
264 3303090 : for (auto n : side->node_index_range())
265 2989250 : side->set_node(n, this->node_ptr(Pyramid18::side_nodes_map[i][n]));
266 : }
267 :
268 :
269 :
270 0 : std::unique_ptr<Elem> Pyramid18::build_edge_ptr (const unsigned int i)
271 : {
272 0 : return this->simple_build_edge_ptr<Edge3,Pyramid18>(i);
273 : }
274 :
275 :
276 :
277 0 : void Pyramid18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
278 : {
279 0 : this->simple_build_edge_ptr<Pyramid18>(edge, i, EDGE3);
280 0 : }
281 :
282 :
283 :
284 0 : void Pyramid18::connectivity(const unsigned int libmesh_dbg_var(sc),
285 : const IOPackage iop,
286 : std::vector<dof_id_type> & /*conn*/) const
287 : {
288 0 : libmesh_assert(_nodes);
289 0 : libmesh_assert_less (sc, this->n_sub_elem());
290 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
291 :
292 0 : switch (iop)
293 : {
294 0 : case TECPLOT:
295 : {
296 : // TODO
297 0 : libmesh_not_implemented();
298 : }
299 :
300 0 : case VTK:
301 : {
302 : // TODO
303 0 : libmesh_not_implemented();
304 : }
305 :
306 0 : default:
307 0 : libmesh_error_msg("Unsupported IO package " << iop);
308 : }
309 : }
310 :
311 :
312 :
313 891618 : unsigned int Pyramid18::n_second_order_adjacent_vertices (const unsigned int n) const
314 : {
315 25116 : switch (n)
316 : {
317 15456 : case 5:
318 : case 6:
319 : case 7:
320 : case 8:
321 : case 9:
322 : case 10:
323 : case 11:
324 : case 12:
325 15456 : return 2;
326 :
327 1932 : case 13:
328 1932 : return 4;
329 :
330 7728 : case 14:
331 : case 15:
332 : case 16:
333 : case 17:
334 7728 : return 3;
335 :
336 0 : default:
337 0 : libmesh_error_msg("Invalid node n = " << n);
338 : }
339 : }
340 :
341 :
342 2194752 : unsigned short int Pyramid18::second_order_adjacent_vertex (const unsigned int n,
343 : const unsigned int v) const
344 : {
345 61824 : libmesh_assert_greater_equal (n, this->n_vertices());
346 61824 : libmesh_assert_less (n, this->n_nodes());
347 :
348 2194752 : switch (n)
349 : {
350 1097376 : case 5:
351 : case 6:
352 : case 7:
353 : case 8:
354 : case 9:
355 : case 10:
356 : case 11:
357 : case 12:
358 : {
359 30912 : libmesh_assert_less (v, 2);
360 :
361 : // This is the analog of the static, const arrays
362 : // {Hex,Prism,Tet10}::_second_order_adjacent_vertices
363 : // defined in the respective source files... possibly treat
364 : // this similarly once the Pyramid13 has been added?
365 1097376 : constexpr unsigned short node_list[8][2] =
366 : {
367 : {0,1},
368 : {1,2},
369 : {2,3},
370 : {0,3},
371 : {0,4},
372 : {1,4},
373 : {2,4},
374 : {3,4}
375 : };
376 :
377 1097376 : return node_list[n-5][v];
378 : }
379 :
380 : // mid-face node on bottom
381 7728 : case 13:
382 : {
383 7728 : libmesh_assert_less (v, 4);
384 :
385 : // The vertex nodes surrounding node 13 are 0, 1, 2, and 3.
386 : // Thus, the v'th node is simply = v.
387 274344 : return cast_int<unsigned short>(v);
388 : }
389 :
390 : // mid-face nodes on triangles
391 823032 : case 14:
392 : case 15:
393 : case 16:
394 : case 17:
395 : {
396 23184 : libmesh_assert_less (v, 3);
397 :
398 823032 : constexpr unsigned short node_list[4][3] =
399 : {
400 : {0,1,4},
401 : {1,2,4},
402 : {2,3,4},
403 : {0,3,4}
404 : };
405 :
406 823032 : return node_list[n-14][v];
407 : }
408 :
409 0 : default:
410 0 : libmesh_error_msg("Invalid n = " << n);
411 :
412 : }
413 : }
414 :
415 :
416 :
417 4788 : void Pyramid18::permute(unsigned int perm_num)
418 : {
419 300 : libmesh_assert_less (perm_num, 4);
420 :
421 11556 : for (unsigned int i = 0; i != perm_num; ++i)
422 : {
423 6768 : swap4nodes(0,1,2,3);
424 6768 : swap4nodes(5,6,7,8);
425 6768 : swap4nodes(9,10,11,12);
426 6768 : swap4nodes(14,15,16,17);
427 6336 : swap4neighbors(0,1,2,3);
428 : }
429 4788 : }
430 :
431 :
432 1728 : void Pyramid18::flip(BoundaryInfo * boundary_info)
433 : {
434 144 : libmesh_assert(boundary_info);
435 :
436 1728 : swap2nodes(0,1);
437 1728 : swap2nodes(2,3);
438 1728 : swap2nodes(6,8);
439 1728 : swap2nodes(9,10);
440 1728 : swap2nodes(11,12);
441 1728 : swap2nodes(15,17);
442 144 : swap2neighbors(1,3);
443 1728 : swap2boundarysides(1,3,boundary_info);
444 1728 : swap2boundaryedges(1,3,boundary_info);
445 1728 : swap2boundaryedges(4,5,boundary_info);
446 1728 : swap2boundaryedges(6,7,boundary_info);
447 1728 : }
448 :
449 :
450 2880 : unsigned int Pyramid18::center_node_on_side(const unsigned short side) const
451 : {
452 240 : libmesh_assert_less (side, Pyramid18::num_sides);
453 2880 : return side == 4 ? 13 : side+14;
454 : }
455 :
456 :
457 8640 : ElemType Pyramid18::side_type (const unsigned int s) const
458 : {
459 720 : libmesh_assert_less (s, 5);
460 8640 : if (s < 4)
461 6912 : return TRI7;
462 144 : return QUAD9;
463 : }
464 :
465 :
466 : } // namespace libMesh
|
