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