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