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