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_pyramid14.h"
21 : #include "libmesh/fe_reference_element_traits.h"
22 : #include "libmesh/edge_edge3.h"
23 : #include "libmesh/face_tri6.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 : // Pyramid14 class static member initializations
36 : const int Pyramid14::num_nodes;
37 : const int Pyramid14::nodes_per_side;
38 : const int Pyramid14::nodes_per_edge;
39 :
40 : const ReferenceElementTable<Pyramid14::num_sides, Pyramid14::nodes_per_side>
41 : Pyramid14::side_nodes_map = build_side_nodes<Pyramid14::num_sides, Pyramid14::nodes_per_side>(PYRAMID14);
42 :
43 : const ReferenceElementTable<Pyramid14::num_edges, Pyramid14::nodes_per_edge>
44 : Pyramid14::edge_nodes_map = pyramid_edge_nodes();
45 :
46 : // ------------------------------------------------------------
47 : // Pyramid14 class member functions
48 :
49 289968 : bool Pyramid14::is_vertex(const unsigned int i) const
50 : {
51 289968 : if (i < 5)
52 153768 : return true;
53 37608 : return false;
54 : }
55 :
56 :
57 :
58 10368 : bool Pyramid14::is_edge(const unsigned int i) const
59 : {
60 10368 : if (i < 5)
61 0 : return false;
62 10368 : if (i == 13)
63 1152 : return false;
64 2304 : return true;
65 : }
66 :
67 :
68 :
69 23952 : bool Pyramid14::is_face(const unsigned int i) const
70 : {
71 23952 : if (i == 13)
72 2112 : return true;
73 5460 : return false;
74 : }
75 :
76 :
77 :
78 39370 : bool Pyramid14::is_node_on_side(const unsigned int n,
79 : const unsigned int s) const
80 : {
81 9860 : libmesh_assert_less (s, n_sides());
82 9860 : return std::find(std::begin(side_nodes_map[s]),
83 9860 : std::end(side_nodes_map[s]),
84 39370 : n) != std::end(side_nodes_map[s]);
85 : }
86 :
87 : std::vector<unsigned>
88 14891 : Pyramid14::nodes_on_side(const unsigned int s) const
89 : {
90 3730 : libmesh_assert_less(s, n_sides());
91 14891 : auto trim = (s == 4) ? 0 : 3;
92 14891 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
93 : }
94 :
95 : std::vector<unsigned>
96 10808 : Pyramid14::nodes_on_edge(const unsigned int e) const
97 : {
98 2704 : libmesh_assert_less(e, n_edges());
99 10808 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
100 : }
101 :
102 202960 : bool Pyramid14::is_node_on_edge(const unsigned int n,
103 : const unsigned int e) const
104 : {
105 56144 : libmesh_assert_less (e, n_edges());
106 56144 : return std::find(std::begin(edge_nodes_map[e]),
107 56144 : std::end(edge_nodes_map[e]),
108 202960 : n) != std::end(edge_nodes_map[e]);
109 : }
110 :
111 18680 : bool Pyramid14::has_affine_map() const
112 : {
113 : // TODO: If the base is a parallelogram and all the triangular faces are planar,
114 : // the map should be linear, but I need to test this theory...
115 18680 : return false;
116 : }
117 :
118 :
119 :
120 1379536 : Order Pyramid14::default_order() const
121 : {
122 1379536 : return SECOND;
123 : }
124 :
125 :
126 :
127 0 : dof_id_type Pyramid14::key (const unsigned int s) const
128 : {
129 0 : libmesh_assert_less (s, this->n_sides());
130 :
131 0 : switch (s)
132 : {
133 0 : case 0: // triangular face 1
134 : case 1: // triangular face 2
135 : case 2: // triangular face 3
136 : case 3: // triangular face 4
137 0 : return Pyramid::key(s);
138 :
139 0 : case 4: // the quad face at z=0
140 0 : return this->compute_key (this->node_id(13));
141 :
142 0 : default:
143 0 : libmesh_error_msg("Invalid side s = " << s);
144 : }
145 : }
146 :
147 :
148 :
149 6574 : unsigned int Pyramid14::local_side_node(unsigned int side,
150 : unsigned int side_node) const
151 : {
152 1652 : libmesh_assert_less (side, this->n_sides());
153 :
154 : // Never more than 9 nodes per side.
155 1652 : libmesh_assert_less(side_node, Pyramid14::nodes_per_side);
156 :
157 : // Some sides have 6 nodes.
158 1652 : libmesh_assert(side == 4 || side_node < 6);
159 :
160 6574 : return Pyramid14::side_nodes_map[side][side_node];
161 : }
162 :
163 :
164 :
165 44572 : unsigned int Pyramid14::local_edge_node(unsigned int edge,
166 : unsigned int edge_node) const
167 : {
168 11144 : libmesh_assert_less(edge, this->n_edges());
169 11144 : libmesh_assert_less(edge_node, Pyramid14::nodes_per_edge);
170 :
171 44572 : return Pyramid14::edge_nodes_map[edge][edge_node];
172 : }
173 :
174 :
175 :
176 54245 : std::unique_ptr<Elem> Pyramid14::build_side_ptr (const unsigned int i)
177 : {
178 17254 : libmesh_assert_less (i, this->n_sides());
179 :
180 54245 : std::unique_ptr<Elem> face;
181 :
182 54245 : switch (i)
183 : {
184 28046 : case 0: // triangular face 1
185 : case 1: // triangular face 2
186 : case 2: // triangular face 3
187 : case 3: // triangular face 4
188 : {
189 28046 : face = std::make_unique<Tri6>();
190 28046 : break;
191 : }
192 26199 : case 4: // the quad face at z=0
193 : {
194 26199 : face = std::make_unique<Quad9>();
195 26199 : break;
196 : }
197 0 : default:
198 0 : libmesh_error_msg("Invalid side i = " << i);
199 : }
200 :
201 : // Set the nodes
202 458312 : for (auto n : face->node_index_range())
203 532773 : face->set_node(n, this->node_ptr(Pyramid14::side_nodes_map[i][n]));
204 :
205 54245 : face->set_interior_parent(this);
206 36991 : face->inherit_data_from(*this);
207 :
208 54245 : return face;
209 0 : }
210 :
211 :
212 :
213 104429 : void Pyramid14::build_side_ptr (std::unique_ptr<Elem> & side,
214 : const unsigned int i)
215 : {
216 33478 : libmesh_assert_less (i, this->n_sides());
217 :
218 104429 : switch (i)
219 : {
220 25136 : case 0: // triangular face 1
221 : case 1: // triangular face 2
222 : case 2: // triangular face 3
223 : case 3: // triangular face 4
224 : {
225 78436 : if (!side.get() || side->type() != TRI6)
226 : {
227 308 : side = this->build_side_ptr(i);
228 206 : return;
229 : }
230 25084 : break;
231 : }
232 8342 : case 4: // the quad face at z=0
233 : {
234 25993 : if (!side.get() || side->type() != QUAD9)
235 : {
236 34686 : side = this->build_side_ptr(i);
237 25581 : return;
238 : }
239 104 : break;
240 : }
241 0 : default:
242 0 : libmesh_error_msg("Invalid side i = " << i);
243 : }
244 :
245 53454 : side->inherit_data_from(*this);
246 :
247 : // Set the nodes
248 551730 : for (auto n : side->node_index_range())
249 624528 : side->set_node(n, this->node_ptr(Pyramid14::side_nodes_map[i][n]));
250 : }
251 :
252 :
253 :
254 168 : std::unique_ptr<Elem> Pyramid14::build_edge_ptr (const unsigned int i)
255 : {
256 168 : return this->simple_build_edge_ptr<Edge3,Pyramid14>(i);
257 : }
258 :
259 :
260 :
261 0 : void Pyramid14::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
262 : {
263 0 : this->simple_build_edge_ptr<Pyramid14>(edge, i, EDGE3);
264 0 : }
265 :
266 :
267 :
268 0 : void Pyramid14::connectivity(const unsigned int libmesh_dbg_var(sc),
269 : const IOPackage iop,
270 : std::vector<dof_id_type> & /*conn*/) const
271 : {
272 0 : libmesh_assert(_nodes);
273 0 : libmesh_assert_less (sc, this->n_sub_elem());
274 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
275 :
276 0 : switch (iop)
277 : {
278 0 : case TECPLOT:
279 : {
280 : // TODO
281 0 : libmesh_not_implemented();
282 : }
283 :
284 0 : case VTK:
285 : {
286 : // TODO
287 0 : libmesh_not_implemented();
288 : }
289 :
290 0 : default:
291 0 : libmesh_error_msg("Unsupported IO package " << iop);
292 : }
293 : }
294 :
295 :
296 :
297 97902 : unsigned int Pyramid14::n_second_order_adjacent_vertices (const unsigned int n) const
298 : {
299 97902 : switch (n)
300 : {
301 24864 : case 5:
302 : case 6:
303 : case 7:
304 : case 8:
305 : case 9:
306 : case 10:
307 : case 11:
308 : case 12:
309 24864 : return 2;
310 :
311 10878 : case 13:
312 10878 : return 4;
313 :
314 0 : default:
315 0 : libmesh_error_msg("Invalid node n = " << n);
316 : }
317 : }
318 :
319 :
320 217560 : unsigned short int Pyramid14::second_order_adjacent_vertex (const unsigned int n,
321 : const unsigned int v) const
322 : {
323 62160 : libmesh_assert_greater_equal (n, this->n_vertices());
324 62160 : libmesh_assert_less (n, this->n_nodes());
325 :
326 217560 : switch (n)
327 : {
328 174048 : case 5:
329 : case 6:
330 : case 7:
331 : case 8:
332 : case 9:
333 : case 10:
334 : case 11:
335 : case 12:
336 : {
337 49728 : libmesh_assert_less (v, 2);
338 :
339 : // This is the analog of the static, const arrays
340 : // {Hex,Prism,Tet10}::_second_order_adjacent_vertices
341 : // defined in the respective source files... possibly treat
342 : // this similarly once the Pyramid13 has been added?
343 174048 : unsigned short node_list[8][2] =
344 : {
345 : {0,1},
346 : {1,2},
347 : {2,3},
348 : {0,3},
349 : {0,4},
350 : {1,4},
351 : {2,4},
352 : {3,4}
353 : };
354 :
355 174048 : return node_list[n-5][v];
356 : }
357 :
358 : // mid-face node on bottom
359 12432 : case 13:
360 : {
361 12432 : libmesh_assert_less (v, 4);
362 :
363 : // The vertex nodes surrounding node 13 are 0, 1, 2, and 3.
364 : // Thus, the v'th node is simply = v.
365 43512 : return cast_int<unsigned short>(v);
366 : }
367 :
368 0 : default:
369 0 : libmesh_error_msg("Invalid n = " << n);
370 :
371 : }
372 : }
373 :
374 :
375 :
376 201 : Real Pyramid14::volume () const
377 : {
378 : // This specialization is good for Lagrange mappings only in general
379 201 : if (this->mapping_type() != LAGRANGE_MAP)
380 0 : return this->Elem::volume();
381 :
382 : // Make copies of our points. It makes the subsequent calculations a bit
383 : // shorter and avoids dereferencing the same pointer multiple times.
384 : Point
385 551 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5), x6 = point(6),
386 501 : x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13);
387 :
388 : // dx/dxi and dx/deta have 15 components while dx/dzeta has 20.
389 : // These are copied directly from the output of a Python script.
390 : Point dx_dxi[15] =
391 : {
392 52 : x6/2 - x8/2,
393 52 : x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - 3*x6/2 + 3*x8/2 - x9,
394 52 : -x0/2 + x1/2 - 2*x10 - 2*x11 + 2*x12 + x2/2 - x3/2 + 3*x6/2 - 3*x8/2 + 2*x9,
395 52 : x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
396 52 : x0/4 - x1/4 + x2/4 - x3/4,
397 52 : -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
398 52 : x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
399 52 : -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
400 52 : x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
401 52 : -2*x13 + x6 + x8,
402 52 : 4*x13 - 2*x6 - 2*x8,
403 52 : -2*x13 + x6 + x8,
404 52 : -x0/2 - x1/2 + x2/2 + x3/2 + x5 - x7,
405 52 : x0/2 + x1/2 - x2/2 - x3/2 - x5 + x7,
406 52 : x0/2 + x1/2 + 2*x13 + x2/2 + x3/2 - x5 - x6 - x7 - x8
407 832 : };
408 :
409 : // dx/dxi and dx/deta have 15 components while dx/dzeta has 20.
410 : // These are copied directly from the output of a Python script.
411 : Point dx_deta[15] =
412 : {
413 52 : -x5/2 + x7/2,
414 52 : x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + 3*x5/2 - 3*x7/2 - x9,
415 52 : -x0/2 - x1/2 + 2*x10 - 2*x11 - 2*x12 + x2/2 + x3/2 - 3*x5/2 + 3*x7/2 + 2*x9,
416 52 : x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
417 52 : -2*x13 + x5 + x7,
418 52 : 4*x13 - 2*x5 - 2*x7,
419 52 : -2*x13 + x5 + x7,
420 52 : x0/4 - x1/4 + x2/4 - x3/4,
421 52 : -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
422 52 : x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
423 52 : -x0/2 + x1/2 + x2/2 - x3/2 - x6 + x8,
424 52 : x0/2 - x1/2 - x2/2 + x3/2 + x6 - x8,
425 52 : -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
426 52 : x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2,
427 52 : x0/2 + x1/2 + 2*x13 + x2/2 + x3/2 - x5 - x6 - x7 - x8
428 832 : };
429 :
430 : // dx/dxi and dx/deta have 15 components while dx/dzeta has 20.
431 : // These are copied directly from the output of a Python script.
432 : Point dx_dzeta[20] =
433 : {
434 52 : -x0/4 - x1/4 + x10 + x11 + x12 - 2*x13 - x2/4 - x3/4 - x4 + x9,
435 52 : 5*x0/4 + 5*x1/4 - 5*x10 - 5*x11 - 5*x12 + 8*x13 + 5*x2/4 + 5*x3/4 + 7*x4 - 5*x9,
436 52 : -9*x0/4 - 9*x1/4 + 9*x10 + 9*x11 + 9*x12 - 12*x13 - 9*x2/4 - 9*x3/4 - 15*x4 + 9*x9,
437 52 : 7*x0/4 + 7*x1/4 - 7*x10 - 7*x11 - 7*x12 + 8*x13 + 7*x2/4 + 7*x3/4 + 13*x4 - 7*x9,
438 52 : -x0/2 - x1/2 + 2*x10 + 2*x11 + 2*x12 - 2*x13 - x2/2 - x3/2 - 4*x4 + 2*x9,
439 52 : x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
440 52 : -3*x0/4 - 3*x1/4 + 3*x10 - 3*x11 - 3*x12 + 3*x2/4 + 3*x3/4 - 3*x5/2 + 3*x7/2 + 3*x9,
441 52 : 3*x0/4 + 3*x1/4 - 3*x10 + 3*x11 + 3*x12 - 3*x2/4 - 3*x3/4 + 3*x5/2 - 3*x7/2 - 3*x9,
442 52 : -x0/4 - x1/4 + x10 - x11 - x12 + x2/4 + x3/4 - x5/2 + x7/2 + x9,
443 52 : x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
444 52 : -3*x0/4 + 3*x1/4 - 3*x10 - 3*x11 + 3*x12 + 3*x2/4 - 3*x3/4 + 3*x6/2 - 3*x8/2 + 3*x9,
445 52 : 3*x0/4 - 3*x1/4 + 3*x10 + 3*x11 - 3*x12 - 3*x2/4 + 3*x3/4 - 3*x6/2 + 3*x8/2 - 3*x9,
446 52 : -x0/4 + x1/4 - x10 - x11 + x12 + x2/4 - x3/4 + x6/2 - x8/2 + x9,
447 52 : -x0/4 + x1/4 - x10 + x11 - x12 - x2/4 + x3/4 + x9,
448 52 : x0/4 - x1/4 + x10 - x11 + x12 + x2/4 - x3/4 - x9,
449 52 : -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
450 52 : x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
451 52 : -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
452 52 : x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2,
453 52 : x0/2 + x1/2 + 2*x13 + x2/2 + x3/2 - x5 - x6 - x7 - x8
454 1092 : };
455 :
456 : // The (xi, eta, zeta) exponents for each of the dx_dxi terms
457 : static const int dx_dxi_exponents[15][3] =
458 : {
459 : {0, 0, 0},
460 : {0, 0, 1},
461 : {0, 0, 2},
462 : {0, 0, 3},
463 : {0, 1, 0},
464 : {0, 1, 1},
465 : {0, 1, 2},
466 : {0, 2, 0},
467 : {0, 2, 1},
468 : {1, 0, 0},
469 : {1, 0, 1},
470 : {1, 0, 2},
471 : {1, 1, 0},
472 : {1, 1, 1},
473 : {1, 2, 0}
474 : };
475 :
476 : // The (xi, eta, zeta) exponents for each of the dx_deta terms
477 : static const int dx_deta_exponents[15][3] =
478 : {
479 : {0, 0, 0},
480 : {0, 0, 1},
481 : {0, 0, 2},
482 : {0, 0, 3},
483 : {0, 1, 0},
484 : {0, 1, 1},
485 : {0, 1, 2},
486 : {1, 0, 0},
487 : {1, 0, 1},
488 : {1, 0, 2},
489 : {1, 1, 0},
490 : {1, 1, 1},
491 : {2, 0, 0},
492 : {2, 0, 1},
493 : {2, 1, 0}
494 : };
495 :
496 : // The (xi, eta, zeta) exponents for each of the dx_dzeta terms
497 : static const int dx_dzeta_exponents[20][3] =
498 : {
499 : {0, 0, 0},
500 : {0, 0, 1},
501 : {0, 0, 2},
502 : {0, 0, 3},
503 : {0, 0, 4},
504 : {0, 1, 0},
505 : {0, 1, 1},
506 : {0, 1, 2},
507 : {0, 1, 3},
508 : {1, 0, 0},
509 : {1, 0, 1},
510 : {1, 0, 2},
511 : {1, 0, 3},
512 : {1, 1, 0},
513 : {1, 1, 1},
514 : {1, 2, 0},
515 : {1, 2, 1},
516 : {2, 1, 0},
517 : {2, 1, 1},
518 : {2, 2, 0},
519 : };
520 :
521 : // Number of points in the quadrature rule
522 52 : const int N = 27;
523 :
524 : // Parameters of the quadrature rule
525 : static const Real
526 : // Parameters used for (xi, eta) quadrature points.
527 : a1 = -7.1805574131988893873307823958101e-01_R,
528 : a2 = -5.0580870785392503961340276902425e-01_R,
529 : a3 = -2.2850430565396735359574351631314e-01_R,
530 : // Parameters used for zeta quadrature points.
531 : b1 = 7.2994024073149732155837979012003e-02_R,
532 : b2 = 3.4700376603835188472176354340395e-01_R,
533 : b3 = 7.0500220988849838312239847758405e-01_R,
534 : // There are 9 unique weight values since there are three
535 : // for each of the three unique zeta values.
536 : w1 = 4.8498876871878584357513834016440e-02_R,
537 : w2 = 4.5137737425884574692441981593901e-02_R,
538 : w3 = 9.2440441384508327195915094925393e-03_R,
539 : w4 = 7.7598202995005734972022134426305e-02_R,
540 : w5 = 7.2220379881415319507907170550242e-02_R,
541 : w6 = 1.4790470621521332351346415188063e-02_R,
542 : w7 = 1.2415712479200917595523541508209e-01_R,
543 : w8 = 1.1555260781026451121265147288039e-01_R,
544 : w9 = 2.3664752994434131762154264300901e-02_R;
545 :
546 : // The points and weights of the 3x3x3 quadrature rule
547 : static const Real xi[N][3] =
548 : {// ^0 ^1 ^2
549 : { 1., a1, a1*a1},
550 : { 1., a2, a2*a2},
551 : { 1., a3, a3*a3},
552 : { 1., a1, a1*a1},
553 : { 1., a2, a2*a2},
554 : { 1., a3, a3*a3},
555 : { 1., a1, a1*a1},
556 : { 1., a2, a2*a2},
557 : { 1., a3, a3*a3},
558 : { 1., 0., 0. },
559 : { 1., 0., 0. },
560 : { 1., 0., 0. },
561 : { 1., 0., 0. },
562 : { 1., 0., 0. },
563 : { 1., 0., 0. },
564 : { 1., 0., 0. },
565 : { 1., 0., 0. },
566 : { 1., 0., 0. },
567 : { 1., -a1, a1*a1},
568 : { 1., -a2, a2*a2},
569 : { 1., -a3, a3*a3},
570 : { 1., -a1, a1*a1},
571 : { 1., -a2, a2*a2},
572 : { 1., -a3, a3*a3},
573 : { 1., -a1, a1*a1},
574 : { 1., -a2, a2*a2},
575 : { 1., -a3, a3*a3}
576 : };
577 :
578 : static const Real eta[N][3] =
579 : {// ^0 ^1 ^2
580 : { 1., a1, a1*a1},
581 : { 1., a2, a2*a2},
582 : { 1., a3, a3*a3},
583 : { 1., 0., 0. },
584 : { 1., 0., 0. },
585 : { 1., 0., 0. },
586 : { 1., -a1, a1*a1},
587 : { 1., -a2, a2*a2},
588 : { 1., -a3, a3*a3},
589 : { 1., a1, a1*a1},
590 : { 1., a2, a2*a2},
591 : { 1., a3, a3*a3},
592 : { 1., 0., 0. },
593 : { 1., 0., 0. },
594 : { 1., 0., 0. },
595 : { 1., -a1, a1*a1},
596 : { 1., -a2, a2*a2},
597 : { 1., -a3, a3*a3},
598 : { 1., a1, a1*a1},
599 : { 1., a2, a2*a2},
600 : { 1., a3, a3*a3},
601 : { 1., 0., 0. },
602 : { 1., 0., 0. },
603 : { 1., 0., 0. },
604 : { 1., -a1, a1*a1},
605 : { 1., -a2, a2*a2},
606 : { 1., -a3, a3*a3}
607 : };
608 :
609 : static const Real zeta[N][5] =
610 : {// ^0 ^1 ^2 ^3 ^4
611 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
612 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
613 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
614 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
615 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
616 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
617 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
618 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
619 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
620 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
621 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
622 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
623 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
624 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
625 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
626 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
627 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
628 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
629 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
630 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
631 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
632 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
633 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
634 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
635 : { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
636 : { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
637 : { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3}
638 : };
639 :
640 : static const Real w[N] = {w1, w2, w3, w4, w5, w6, // 0-5
641 : w1, w2, w3, w4, w5, w6, // 6-11
642 : w7, w8, w9, w4, w5, w6, // 12-17
643 : w1, w2, w3, w4, w5, w6, // 18-23
644 : w1, w2, w3}; // 24-26
645 :
646 52 : Real vol = 0.;
647 5628 : for (int q=0; q<N; ++q)
648 : {
649 : // Compute denominators for the current q.
650 : Real
651 5427 : den2 = (1. - zeta[q][1])*(1. - zeta[q][1]),
652 5427 : den3 = den2*(1. - zeta[q][1]);
653 :
654 : // Compute dx/dxi and dx/deta at the current q.
655 1404 : Point dx_dxi_q, dx_deta_q;
656 86832 : for (int c=0; c<15; ++c)
657 : {
658 : dx_dxi_q +=
659 122715 : xi[q][dx_dxi_exponents[c][0]]*
660 80595 : eta[q][dx_dxi_exponents[c][1]]*
661 81405 : zeta[q][dx_dxi_exponents[c][2]]*dx_dxi[c];
662 :
663 : dx_deta_q +=
664 122715 : xi[q][dx_deta_exponents[c][0]]*
665 80595 : eta[q][dx_deta_exponents[c][1]]*
666 81405 : zeta[q][dx_deta_exponents[c][2]]*dx_deta[c];
667 : }
668 :
669 : // Compute dx/dzeta at the current q.
670 1404 : Point dx_dzeta_q;
671 113967 : for (int c=0; c<20; ++c)
672 : {
673 : dx_dzeta_q +=
674 163620 : xi[q][dx_dzeta_exponents[c][0]]*
675 107460 : eta[q][dx_dzeta_exponents[c][1]]*
676 108540 : zeta[q][dx_dzeta_exponents[c][2]]*dx_dzeta[c];
677 : }
678 :
679 : // Scale everything appropriately
680 1404 : dx_dxi_q /= den2;
681 1404 : dx_deta_q /= den2;
682 1404 : dx_dzeta_q /= den3;
683 :
684 : // Compute scalar triple product, multiply by weight, and accumulate volume.
685 6777 : vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
686 : }
687 :
688 52 : return vol;
689 : }
690 :
691 :
692 1188 : void Pyramid14::permute(unsigned int perm_num)
693 : {
694 312 : libmesh_assert_less (perm_num, 4);
695 :
696 2900 : for (unsigned int i = 0; i != perm_num; ++i)
697 : {
698 1712 : swap4nodes(0,1,2,3);
699 1712 : swap4nodes(5,6,7,8);
700 1712 : swap4nodes(9,10,11,12);
701 1264 : swap4neighbors(0,1,2,3);
702 : }
703 1188 : }
704 :
705 :
706 576 : void Pyramid14::flip(BoundaryInfo * boundary_info)
707 : {
708 144 : libmesh_assert(boundary_info);
709 :
710 576 : swap2nodes(0,1);
711 576 : swap2nodes(2,3);
712 576 : swap2nodes(6,8);
713 576 : swap2nodes(9,10);
714 576 : swap2nodes(11,12);
715 144 : swap2neighbors(1,3);
716 576 : swap2boundarysides(1,3,boundary_info);
717 576 : swap2boundaryedges(1,3,boundary_info);
718 576 : swap2boundaryedges(4,5,boundary_info);
719 576 : swap2boundaryedges(6,7,boundary_info);
720 576 : }
721 :
722 :
723 960 : unsigned int Pyramid14::center_node_on_side(const unsigned short side) const
724 : {
725 240 : libmesh_assert_less (side, Pyramid14::num_sides);
726 960 : return side == 4 ? 13 : invalid_uint;
727 : }
728 :
729 :
730 2880 : ElemType Pyramid14::side_type (const unsigned int s) const
731 : {
732 720 : libmesh_assert_less (s, 5);
733 2880 : if (s < 4)
734 2304 : return TRI6;
735 144 : return QUAD9;
736 : }
737 :
738 :
739 : } // namespace libMesh
|
