Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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/fe_reference_element_traits.h"
22 : #include "libmesh/edge_edge3.h"
23 : #include "libmesh/face_tri7.h"
24 : #include "libmesh/face_quad9.h"
25 : #include "libmesh/enum_io_package.h"
26 : #include "libmesh/enum_order.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 :
33 :
34 : // ------------------------------------------------------------
35 : // Pyramid18 class static member initializations
36 : const int Pyramid18::num_nodes;
37 : const int Pyramid18::nodes_per_side;
38 : const int Pyramid18::nodes_per_edge;
39 :
40 : const ReferenceElementTable<Pyramid18::num_sides, Pyramid18::nodes_per_side>
41 : Pyramid18::side_nodes_map = build_side_nodes<Pyramid18::num_sides, Pyramid18::nodes_per_side>(PYRAMID18);
42 :
43 : const ReferenceElementTable<Pyramid18::num_edges, Pyramid18::nodes_per_edge>
44 : Pyramid18::edge_nodes_map = pyramid_edge_nodes();
45 :
46 : // ------------------------------------------------------------
47 : // Pyramid18 class member functions
48 :
49 106704 : bool Pyramid18::is_vertex(const unsigned int i) const
50 : {
51 106704 : if (i < 5)
52 46680 : return true;
53 16002 : return false;
54 : }
55 :
56 :
57 :
58 14976 : bool Pyramid18::is_edge(const unsigned int i) const
59 : {
60 14976 : if (i < 5)
61 0 : return false;
62 14976 : if (i > 12)
63 5760 : return false;
64 2304 : return true;
65 : }
66 :
67 :
68 :
69 39130 : bool Pyramid18::is_face(const unsigned int i) const
70 : {
71 39130 : if (i > 12)
72 8710 : return true;
73 7800 : return false;
74 : }
75 :
76 :
77 :
78 30694 : bool Pyramid18::is_node_on_side(const unsigned int n,
79 : const unsigned int s) const
80 : {
81 7736 : libmesh_assert_less (s, n_sides());
82 7736 : return std::find(std::begin(side_nodes_map[s]),
83 7736 : std::end(side_nodes_map[s]),
84 30694 : n) != std::end(side_nodes_map[s]);
85 : }
86 :
87 : std::vector<unsigned>
88 16156 : Pyramid18::nodes_on_side(const unsigned int s) const
89 : {
90 4045 : libmesh_assert_less(s, n_sides());
91 16156 : auto trim = (s == 4) ? 0 : 2;
92 16156 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
93 : }
94 :
95 : std::vector<unsigned>
96 10752 : Pyramid18::nodes_on_edge(const unsigned int e) const
97 : {
98 2688 : libmesh_assert_less(e, n_edges());
99 10752 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
100 : }
101 :
102 55968 : bool Pyramid18::is_node_on_edge(const unsigned int n,
103 : const unsigned int e) const
104 : {
105 14856 : libmesh_assert_less (e, n_edges());
106 14856 : return std::find(std::begin(edge_nodes_map[e]),
107 14856 : std::end(edge_nodes_map[e]),
108 55968 : n) != std::end(edge_nodes_map[e]);
109 : }
110 :
111 :
112 :
113 18705 : bool Pyramid18::has_affine_map() const
114 : {
115 : // TODO: If the base is a parallelogram and all the triangular faces are planar,
116 : // the map should be linear, but I need to test this theory...
117 18705 : return false;
118 : }
119 :
120 :
121 :
122 2110665 : Order Pyramid18::default_order() const
123 : {
124 2110665 : return THIRD;
125 : }
126 :
127 :
128 :
129 0 : dof_id_type Pyramid18::key (const unsigned int s) const
130 : {
131 0 : libmesh_assert_less (s, this->n_sides());
132 :
133 0 : switch (s)
134 : {
135 0 : case 0: // triangular face 1
136 : case 1: // triangular face 2
137 : case 2: // triangular face 3
138 : case 3: // triangular face 4
139 0 : return this->compute_key (this->node_id(s+14));
140 :
141 0 : case 4: // the quad face at z=0
142 0 : return this->compute_key (this->node_id(13));
143 :
144 0 : default:
145 0 : libmesh_error_msg("Invalid side s = " << s);
146 : }
147 : }
148 :
149 :
150 :
151 7398 : unsigned int Pyramid18::local_side_node(unsigned int side,
152 : unsigned int side_node) const
153 : {
154 1860 : libmesh_assert_less (side, this->n_sides());
155 :
156 : // Never more than 9 nodes per side.
157 1860 : libmesh_assert_less(side_node, Pyramid18::nodes_per_side);
158 :
159 : // Some sides have 7 nodes.
160 1860 : libmesh_assert(side == 4 || side_node < 7);
161 :
162 7398 : return Pyramid18::side_nodes_map[side][side_node];
163 : }
164 :
165 :
166 :
167 44572 : unsigned int Pyramid18::local_edge_node(unsigned int edge,
168 : unsigned int edge_node) const
169 : {
170 11144 : libmesh_assert_less(edge, this->n_edges());
171 11144 : libmesh_assert_less(edge_node, Pyramid18::nodes_per_edge);
172 :
173 44572 : return Pyramid18::edge_nodes_map[edge][edge_node];
174 : }
175 :
176 :
177 :
178 43620 : std::unique_ptr<Elem> Pyramid18::build_side_ptr (const unsigned int i)
179 : {
180 13848 : libmesh_assert_less (i, this->n_sides());
181 :
182 43620 : std::unique_ptr<Elem> face;
183 :
184 43620 : switch (i)
185 : {
186 22695 : case 0: // triangular face 1
187 : case 1: // triangular face 2
188 : case 2: // triangular face 3
189 : case 3: // triangular face 4
190 : {
191 22695 : face = std::make_unique<Tri7>();
192 22695 : break;
193 : }
194 20925 : case 4: // the quad face at z=0
195 : {
196 20925 : face = std::make_unique<Quad9>();
197 20925 : break;
198 : }
199 0 : default:
200 0 : libmesh_error_msg("Invalid side i = " << i);
201 : }
202 :
203 : // Set the nodes
204 390810 : for (auto n : face->node_index_range())
205 457530 : face->set_node(n, this->node_ptr(Pyramid18::side_nodes_map[i][n]));
206 :
207 43620 : face->set_interior_parent(this);
208 29772 : face->inherit_data_from(*this);
209 :
210 43620 : return face;
211 0 : }
212 :
213 :
214 :
215 83424 : void Pyramid18::build_side_ptr (std::unique_ptr<Elem> & side,
216 : const unsigned int i)
217 : {
218 26736 : libmesh_assert_less (i, this->n_sides());
219 :
220 83424 : switch (i)
221 : {
222 20088 : case 0: // triangular face 1
223 : case 1: // triangular face 2
224 : case 2: // triangular face 3
225 : case 3: // triangular face 4
226 : {
227 62712 : if (!side.get() || side->type() != TRI7)
228 : {
229 298 : side = this->build_side_ptr(i);
230 199 : return;
231 : }
232 20038 : break;
233 : }
234 6648 : case 4: // the quad face at z=0
235 : {
236 20712 : if (!side.get() || side->type() != QUAD9)
237 : {
238 27562 : side = this->build_side_ptr(i);
239 20335 : return;
240 : }
241 94 : break;
242 : }
243 0 : default:
244 0 : libmesh_error_msg("Invalid side i = " << i);
245 : }
246 :
247 42758 : side->inherit_data_from(*this);
248 :
249 : // Set the nodes
250 503874 : for (auto n : side->node_index_range())
251 582096 : side->set_node(n, this->node_ptr(Pyramid18::side_nodes_map[i][n]));
252 : }
253 :
254 :
255 :
256 0 : std::unique_ptr<Elem> Pyramid18::build_edge_ptr (const unsigned int i)
257 : {
258 0 : return this->simple_build_edge_ptr<Edge3,Pyramid18>(i);
259 : }
260 :
261 :
262 :
263 0 : void Pyramid18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
264 : {
265 0 : this->simple_build_edge_ptr<Pyramid18>(edge, i, EDGE3);
266 0 : }
267 :
268 :
269 :
270 0 : void Pyramid18::connectivity(const unsigned int libmesh_dbg_var(sc),
271 : const IOPackage iop,
272 : std::vector<dof_id_type> & /*conn*/) const
273 : {
274 0 : libmesh_assert(_nodes);
275 0 : libmesh_assert_less (sc, this->n_sub_elem());
276 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
277 :
278 0 : switch (iop)
279 : {
280 0 : case TECPLOT:
281 : {
282 : // TODO
283 0 : libmesh_not_implemented();
284 : }
285 :
286 0 : case VTK:
287 : {
288 : // TODO
289 0 : libmesh_not_implemented();
290 : }
291 :
292 0 : default:
293 0 : libmesh_error_msg("Unsupported IO package " << iop);
294 : }
295 : }
296 :
297 :
298 :
299 116298 : unsigned int Pyramid18::n_second_order_adjacent_vertices (const unsigned int n) const
300 : {
301 33228 : switch (n)
302 : {
303 20448 : case 5:
304 : case 6:
305 : case 7:
306 : case 8:
307 : case 9:
308 : case 10:
309 : case 11:
310 : case 12:
311 20448 : return 2;
312 :
313 2556 : case 13:
314 2556 : return 4;
315 :
316 10224 : case 14:
317 : case 15:
318 : case 16:
319 : case 17:
320 10224 : return 3;
321 :
322 0 : default:
323 0 : libmesh_error_msg("Invalid node n = " << n);
324 : }
325 : }
326 :
327 :
328 286272 : unsigned short int Pyramid18::second_order_adjacent_vertex (const unsigned int n,
329 : const unsigned int v) const
330 : {
331 81792 : libmesh_assert_greater_equal (n, this->n_vertices());
332 81792 : libmesh_assert_less (n, this->n_nodes());
333 :
334 286272 : switch (n)
335 : {
336 143136 : case 5:
337 : case 6:
338 : case 7:
339 : case 8:
340 : case 9:
341 : case 10:
342 : case 11:
343 : case 12:
344 : {
345 40896 : libmesh_assert_less (v, 2);
346 :
347 : // This is the analog of the static, const arrays
348 : // {Hex,Prism,Tet10}::_second_order_adjacent_vertices
349 : // defined in the respective source files... possibly treat
350 : // this similarly once the Pyramid13 has been added?
351 143136 : constexpr unsigned short node_list[8][2] =
352 : {
353 : {0,1},
354 : {1,2},
355 : {2,3},
356 : {0,3},
357 : {0,4},
358 : {1,4},
359 : {2,4},
360 : {3,4}
361 : };
362 :
363 143136 : return node_list[n-5][v];
364 : }
365 :
366 : // mid-face node on bottom
367 10224 : case 13:
368 : {
369 10224 : libmesh_assert_less (v, 4);
370 :
371 : // The vertex nodes surrounding node 13 are 0, 1, 2, and 3.
372 : // Thus, the v'th node is simply = v.
373 35784 : return cast_int<unsigned short>(v);
374 : }
375 :
376 : // mid-face nodes on triangles
377 107352 : case 14:
378 : case 15:
379 : case 16:
380 : case 17:
381 : {
382 30672 : libmesh_assert_less (v, 3);
383 :
384 107352 : constexpr unsigned short node_list[4][3] =
385 : {
386 : {0,1,4},
387 : {1,2,4},
388 : {2,3,4},
389 : {0,3,4}
390 : };
391 :
392 107352 : return node_list[n-14][v];
393 : }
394 :
395 0 : default:
396 0 : libmesh_error_msg("Invalid n = " << n);
397 :
398 : }
399 : }
400 :
401 :
402 :
403 1188 : void Pyramid18::permute(unsigned int perm_num)
404 : {
405 312 : libmesh_assert_less (perm_num, 4);
406 :
407 2900 : for (unsigned int i = 0; i != perm_num; ++i)
408 : {
409 1712 : swap4nodes(0,1,2,3);
410 1712 : swap4nodes(5,6,7,8);
411 1712 : swap4nodes(9,10,11,12);
412 1712 : swap4nodes(14,15,16,17);
413 1264 : swap4neighbors(0,1,2,3);
414 : }
415 1188 : }
416 :
417 :
418 576 : void Pyramid18::flip(BoundaryInfo * boundary_info)
419 : {
420 144 : libmesh_assert(boundary_info);
421 :
422 576 : swap2nodes(0,1);
423 576 : swap2nodes(2,3);
424 576 : swap2nodes(6,8);
425 576 : swap2nodes(9,10);
426 576 : swap2nodes(11,12);
427 576 : swap2nodes(15,17);
428 144 : swap2neighbors(1,3);
429 576 : swap2boundarysides(1,3,boundary_info);
430 576 : swap2boundaryedges(1,3,boundary_info);
431 576 : swap2boundaryedges(4,5,boundary_info);
432 576 : swap2boundaryedges(6,7,boundary_info);
433 576 : }
434 :
435 :
436 960 : unsigned int Pyramid18::center_node_on_side(const unsigned short side) const
437 : {
438 240 : libmesh_assert_less (side, Pyramid18::num_sides);
439 960 : return side == 4 ? 13 : side+14;
440 : }
441 :
442 :
443 2880 : ElemType Pyramid18::side_type (const unsigned int s) const
444 : {
445 720 : libmesh_assert_less (s, 5);
446 2880 : if (s < 4)
447 2304 : return TRI7;
448 144 : return QUAD9;
449 : }
450 :
451 :
452 : } // namespace libMesh
|
