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