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/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 : // ------------------------------------------------------------
40 : // Pyramid18 class member functions
41 :
42 756144 : bool Pyramid18::is_vertex(const unsigned int i) const
43 : {
44 756144 : if (i < 5)
45 314856 : return true;
46 24498 : return false;
47 : }
48 :
49 :
50 :
51 119808 : bool Pyramid18::is_edge(const unsigned int i) const
52 : {
53 119808 : if (i < 5)
54 0 : return false;
55 119808 : if (i > 12)
56 46080 : return false;
57 6144 : return true;
58 : }
59 :
60 :
61 :
62 147948 : bool Pyramid18::is_face(const unsigned int i) const
63 : {
64 147948 : if (i > 12)
65 31572 : return true;
66 9360 : return false;
67 : }
68 :
69 :
70 :
71 97380 : bool Pyramid18::is_node_on_side(const unsigned int n,
72 : const unsigned int s) const
73 : {
74 8064 : libmesh_assert_less (s, n_sides());
75 8064 : return std::find(std::begin(side_nodes_map[s]),
76 97380 : std::end(side_nodes_map[s]),
77 97380 : n) != std::end(side_nodes_map[s]);
78 : }
79 :
80 : std::vector<unsigned>
81 195455 : Pyramid18::nodes_on_side(const unsigned int s) const
82 : {
83 16237 : libmesh_assert_less(s, n_sides());
84 195455 : auto trim = (s == 4) ? 0 : 2;
85 195455 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
86 : }
87 :
88 : std::vector<unsigned>
89 92160 : Pyramid18::nodes_on_edge(const unsigned int e) const
90 : {
91 7680 : libmesh_assert_less(e, n_edges());
92 99840 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
93 : }
94 :
95 340704 : bool Pyramid18::is_node_on_edge(const unsigned int n,
96 : const unsigned int e) const
97 : {
98 14856 : libmesh_assert_less (e, n_edges());
99 14856 : return std::find(std::begin(edge_nodes_map[e]),
100 340704 : std::end(edge_nodes_map[e]),
101 340704 : n) != std::end(edge_nodes_map[e]);
102 : }
103 :
104 :
105 :
106 77505 : bool Pyramid18::has_affine_map() const
107 : {
108 : // TODO: If the base is a parallelogram and all the triangular faces are planar,
109 : // the map should be linear, but I need to test this theory...
110 77505 : return false;
111 : }
112 :
113 :
114 :
115 24856859 : Order Pyramid18::default_order() const
116 : {
117 24856859 : return THIRD;
118 : }
119 :
120 :
121 :
122 0 : dof_id_type Pyramid18::key (const unsigned int s) const
123 : {
124 0 : libmesh_assert_less (s, this->n_sides());
125 :
126 0 : switch (s)
127 : {
128 0 : case 0: // triangular face 1
129 : case 1: // triangular face 2
130 : case 2: // triangular face 3
131 : case 3: // triangular face 4
132 0 : return this->compute_key (this->node_id(s+14));
133 :
134 0 : case 4: // the quad face at z=0
135 0 : return this->compute_key (this->node_id(13));
136 :
137 0 : default:
138 0 : libmesh_error_msg("Invalid side s = " << s);
139 : }
140 : }
141 :
142 :
143 :
144 40422 : unsigned int Pyramid18::local_side_node(unsigned int side,
145 : unsigned int side_node) const
146 : {
147 3204 : libmesh_assert_less (side, this->n_sides());
148 :
149 : // Never more than 9 nodes per side.
150 3204 : libmesh_assert_less(side_node, Pyramid18::nodes_per_side);
151 :
152 : // Some sides have 7 nodes.
153 3204 : libmesh_assert(side == 4 || side_node < 7);
154 :
155 40422 : return Pyramid18::side_nodes_map[side][side_node];
156 : }
157 :
158 :
159 :
160 669916 : unsigned int Pyramid18::local_edge_node(unsigned int edge,
161 : unsigned int edge_node) const
162 : {
163 55752 : libmesh_assert_less(edge, this->n_edges());
164 55752 : libmesh_assert_less(edge_node, Pyramid18::nodes_per_edge);
165 :
166 669916 : return Pyramid18::edge_nodes_map[edge][edge_node];
167 : }
168 :
169 :
170 :
171 441036 : std::unique_ptr<Elem> Pyramid18::build_side_ptr (const unsigned int i)
172 : {
173 9360 : libmesh_assert_less (i, this->n_sides());
174 :
175 441036 : std::unique_ptr<Elem> face;
176 :
177 441036 : switch (i)
178 : {
179 223323 : case 0: // triangular face 1
180 : case 1: // triangular face 2
181 : case 2: // triangular face 3
182 : case 3: // triangular face 4
183 : {
184 223323 : face = std::make_unique<Tri7>();
185 223323 : break;
186 : }
187 217713 : case 4: // the quad face at z=0
188 : {
189 217713 : face = std::make_unique<Quad9>();
190 217713 : break;
191 : }
192 0 : default:
193 0 : libmesh_error_msg("Invalid side i = " << i);
194 : }
195 :
196 : // Set the nodes
197 3963714 : for (auto n : face->node_index_range())
198 3597114 : face->set_node(n, this->node_ptr(Pyramid18::side_nodes_map[i][n]));
199 :
200 441036 : face->set_interior_parent(this);
201 431676 : face->inherit_data_from(*this);
202 :
203 441036 : return face;
204 0 : }
205 :
206 :
207 :
208 869424 : void Pyramid18::build_side_ptr (std::unique_ptr<Elem> & side,
209 : const unsigned int i)
210 : {
211 17760 : libmesh_assert_less (i, this->n_sides());
212 :
213 869424 : switch (i)
214 : {
215 13356 : case 0: // triangular face 1
216 : case 1: // triangular face 2
217 : case 2: // triangular face 3
218 : case 3: // triangular face 4
219 : {
220 652500 : if (!side.get() || side->type() != TRI7)
221 : {
222 1194 : side = this->build_side_ptr(i);
223 647 : return;
224 : }
225 13306 : break;
226 : }
227 4404 : case 4: // the quad face at z=0
228 : {
229 216924 : if (!side.get() || side->type() != QUAD9)
230 : {
231 423066 : side = this->build_side_ptr(i);
232 215843 : return;
233 : }
234 94 : break;
235 : }
236 0 : default:
237 0 : libmesh_error_msg("Invalid side i = " << i);
238 : }
239 :
240 639534 : side->inherit_data_from(*this);
241 :
242 : // Set the nodes
243 5225634 : for (auto n : side->node_index_range())
244 4666688 : side->set_node(n, this->node_ptr(Pyramid18::side_nodes_map[i][n]));
245 : }
246 :
247 :
248 :
249 4608 : std::unique_ptr<Elem> Pyramid18::build_edge_ptr (const unsigned int i)
250 : {
251 4608 : return this->simple_build_edge_ptr<Edge3,Pyramid18>(i);
252 : }
253 :
254 :
255 :
256 0 : void Pyramid18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
257 : {
258 0 : this->simple_build_edge_ptr<Pyramid18>(edge, i, EDGE3);
259 0 : }
260 :
261 :
262 :
263 0 : void Pyramid18::connectivity(const unsigned int libmesh_dbg_var(sc),
264 : const IOPackage iop,
265 : std::vector<dof_id_type> & /*conn*/) const
266 : {
267 0 : libmesh_assert(_nodes);
268 0 : libmesh_assert_less (sc, this->n_sub_elem());
269 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
270 :
271 0 : switch (iop)
272 : {
273 0 : case TECPLOT:
274 : {
275 : // TODO
276 0 : libmesh_not_implemented();
277 : }
278 :
279 0 : case VTK:
280 : {
281 : // TODO
282 0 : libmesh_not_implemented();
283 : }
284 :
285 0 : default:
286 0 : libmesh_error_msg("Unsupported IO package " << iop);
287 : }
288 : }
289 :
290 :
291 :
292 1412190 : unsigned int Pyramid18::n_second_order_adjacent_vertices (const unsigned int n) const
293 : {
294 39780 : switch (n)
295 : {
296 24480 : case 5:
297 : case 6:
298 : case 7:
299 : case 8:
300 : case 9:
301 : case 10:
302 : case 11:
303 : case 12:
304 24480 : return 2;
305 :
306 3060 : case 13:
307 3060 : return 4;
308 :
309 12240 : case 14:
310 : case 15:
311 : case 16:
312 : case 17:
313 12240 : return 3;
314 :
315 0 : default:
316 0 : libmesh_error_msg("Invalid node n = " << n);
317 : }
318 : }
319 :
320 :
321 3476160 : unsigned short int Pyramid18::second_order_adjacent_vertex (const unsigned int n,
322 : const unsigned int v) const
323 : {
324 97920 : libmesh_assert_greater_equal (n, this->n_vertices());
325 97920 : libmesh_assert_less (n, this->n_nodes());
326 :
327 3476160 : switch (n)
328 : {
329 1738080 : case 5:
330 : case 6:
331 : case 7:
332 : case 8:
333 : case 9:
334 : case 10:
335 : case 11:
336 : case 12:
337 : {
338 48960 : libmesh_assert_less (v, 2);
339 :
340 : // This is the analog of the static, const arrays
341 : // {Hex,Prism,Tet10}::_second_order_adjacent_vertices
342 : // defined in the respective source files... possibly treat
343 : // this similarly once the Pyramid13 has been added?
344 1738080 : constexpr unsigned short node_list[8][2] =
345 : {
346 : {0,1},
347 : {1,2},
348 : {2,3},
349 : {0,3},
350 : {0,4},
351 : {1,4},
352 : {2,4},
353 : {3,4}
354 : };
355 :
356 1738080 : return node_list[n-5][v];
357 : }
358 :
359 : // mid-face node on bottom
360 12240 : case 13:
361 : {
362 12240 : libmesh_assert_less (v, 4);
363 :
364 : // The vertex nodes surrounding node 13 are 0, 1, 2, and 3.
365 : // Thus, the v'th node is simply = v.
366 434520 : return cast_int<unsigned short>(v);
367 : }
368 :
369 : // mid-face nodes on triangles
370 1303560 : case 14:
371 : case 15:
372 : case 16:
373 : case 17:
374 : {
375 36720 : libmesh_assert_less (v, 3);
376 :
377 1303560 : constexpr unsigned short node_list[4][3] =
378 : {
379 : {0,1,4},
380 : {1,2,4},
381 : {2,3,4},
382 : {0,3,4}
383 : };
384 :
385 1303560 : return node_list[n-14][v];
386 : }
387 :
388 0 : default:
389 0 : libmesh_error_msg("Invalid n = " << n);
390 :
391 : }
392 : }
393 :
394 :
395 :
396 10224 : void Pyramid18::permute(unsigned int perm_num)
397 : {
398 720 : libmesh_assert_less (perm_num, 4);
399 :
400 25008 : for (unsigned int i = 0; i != perm_num; ++i)
401 : {
402 14784 : swap4nodes(0,1,2,3);
403 14784 : swap4nodes(5,6,7,8);
404 14784 : swap4nodes(9,10,11,12);
405 14784 : swap4nodes(14,15,16,17);
406 13728 : swap4neighbors(0,1,2,3);
407 : }
408 10224 : }
409 :
410 :
411 1728 : void Pyramid18::flip(BoundaryInfo * boundary_info)
412 : {
413 144 : libmesh_assert(boundary_info);
414 :
415 1728 : swap2nodes(0,1);
416 1728 : swap2nodes(2,3);
417 1728 : swap2nodes(6,8);
418 1728 : swap2nodes(9,10);
419 1728 : swap2nodes(11,12);
420 1728 : swap2nodes(15,17);
421 144 : swap2neighbors(1,3);
422 1728 : swap2boundarysides(1,3,boundary_info);
423 1728 : swap2boundaryedges(1,3,boundary_info);
424 1728 : swap2boundaryedges(4,5,boundary_info);
425 1728 : swap2boundaryedges(6,7,boundary_info);
426 1728 : }
427 :
428 :
429 2880 : unsigned int Pyramid18::center_node_on_side(const unsigned short side) const
430 : {
431 240 : libmesh_assert_less (side, Pyramid18::num_sides);
432 2880 : return side == 4 ? 13 : side+14;
433 : }
434 :
435 :
436 56448 : ElemType Pyramid18::side_type (const unsigned int s) const
437 : {
438 4704 : libmesh_assert_less (s, 5);
439 56448 : if (s < 4)
440 46080 : return TRI7;
441 864 : return QUAD9;
442 : }
443 :
444 :
445 : } // namespace libMesh
|
