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