Loading [MathJax]/extensions/tex2jax.js
libMesh
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
cell_pyramid13.C
Go to the documentation of this file.
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;
38 
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 
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 bool Pyramid13::is_vertex(const unsigned int i) const
64 {
65  if (i < 5)
66  return true;
67  return false;
68 }
69 
70 
71 
72 bool Pyramid13::is_edge(const unsigned int i) const
73 {
74  if (i < 5)
75  return false;
76  return true;
77 }
78 
79 
80 
81 bool Pyramid13::is_face(const unsigned int) const
82 {
83  return false;
84 }
85 
86 
87 
88 bool Pyramid13::is_node_on_side(const unsigned int n,
89  const unsigned int s) const
90 {
91  libmesh_assert_less (s, n_sides());
92  return std::find(std::begin(side_nodes_map[s]),
93  std::end(side_nodes_map[s]),
94  n) != std::end(side_nodes_map[s]);
95 }
96 
97 std::vector<unsigned>
98 Pyramid13::nodes_on_side(const unsigned int s) const
99 {
100  libmesh_assert_less(s, n_sides());
101  auto trim = (s == 4) ? 0 : 2;
102  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
103 }
104 
105 std::vector<unsigned>
106 Pyramid13::nodes_on_edge(const unsigned int e) const
107 {
108  libmesh_assert_less(e, n_edges());
109  return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
110 }
111 
112 bool Pyramid13::is_node_on_edge(const unsigned int n,
113  const unsigned int e) const
114 {
115  libmesh_assert_less (e, n_edges());
116  return std::find(std::begin(edge_nodes_map[e]),
117  std::end(edge_nodes_map[e]),
118  n) != std::end(edge_nodes_map[e]);
119 }
120 
121 
122 
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  return false;
128 }
129 
130 
131 
133 {
134  return SECOND;
135 }
136 
137 
138 
139 unsigned int Pyramid13::local_side_node(unsigned int side,
140  unsigned int side_node) const
141 {
142  libmesh_assert_less (side, this->n_sides());
143 
144  // Never more than 8 nodes per side.
145  libmesh_assert_less(side_node, Pyramid13::nodes_per_side);
146 
147  // Some sides have 6 nodes.
148  libmesh_assert(side == 4 || side_node < 6);
149 
150  return Pyramid13::side_nodes_map[side][side_node];
151 }
152 
153 
154 
155 unsigned int Pyramid13::local_edge_node(unsigned int edge,
156  unsigned int edge_node) const
157 {
158  libmesh_assert_less(edge, this->n_edges());
159  libmesh_assert_less(edge_node, Pyramid13::nodes_per_edge);
160 
161  return Pyramid13::edge_nodes_map[edge][edge_node];
162 }
163 
164 
165 
166 std::unique_ptr<Elem> Pyramid13::build_side_ptr (const unsigned int i)
167 {
168  libmesh_assert_less (i, this->n_sides());
169 
170  std::unique_ptr<Elem> face;
171 
172  switch (i)
173  {
174  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  face = std::make_unique<Tri6>();
180  break;
181  }
182  case 4: // the quad face at z=0
183  {
184  face = std::make_unique<Quad8>();
185  break;
186  }
187  default:
188  libmesh_error_msg("Invalid side i = " << i);
189  }
190 
191  // Set the nodes
192  for (auto n : face->node_index_range())
193  face->set_node(n, this->node_ptr(Pyramid13::side_nodes_map[i][n]));
194 
195  face->set_interior_parent(this);
196  face->inherit_data_from(*this);
197 
198  return face;
199 }
200 
201 
202 
203 void Pyramid13::build_side_ptr (std::unique_ptr<Elem> & side,
204  const unsigned int i)
205 {
206  libmesh_assert_less (i, this->n_sides());
207 
208  switch (i)
209  {
210  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  if (!side.get() || side->type() != TRI6)
216  {
217  side = this->build_side_ptr(i);
218  return;
219  }
220  break;
221  }
222  case 4: // the quad face at z=0
223  {
224  if (!side.get() || side->type() != QUAD8)
225  {
226  side = this->build_side_ptr(i);
227  return;
228  }
229  break;
230  }
231  default:
232  libmesh_error_msg("Invalid side i = " << i);
233  }
234 
235  side->inherit_data_from(*this);
236 
237  // Set the nodes
238  for (auto n : side->node_index_range())
239  side->set_node(n, this->node_ptr(Pyramid13::side_nodes_map[i][n]));
240 }
241 
242 
243 
244 std::unique_ptr<Elem> Pyramid13::build_edge_ptr (const unsigned int i)
245 {
246  return this->simple_build_edge_ptr<Edge3,Pyramid13>(i);
247 }
248 
249 
250 
251 void Pyramid13::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
252 {
253  this->simple_build_edge_ptr<Pyramid13>(edge, i, EDGE3);
254 }
255 
256 
257 
258 void Pyramid13::connectivity(const unsigned int libmesh_dbg_var(sc),
259  const IOPackage iop,
260  std::vector<dof_id_type> & /*conn*/) const
261 {
263  libmesh_assert_less (sc, this->n_sub_elem());
264  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
265 
266  switch (iop)
267  {
268  case TECPLOT:
269  {
270  // TODO
271  libmesh_not_implemented();
272  }
273 
274  case VTK:
275  {
276  // TODO
277  libmesh_not_implemented();
278  }
279 
280  default:
281  libmesh_error_msg("Unsupported IO package " << iop);
282  }
283 }
284 
285 
286 
287 unsigned int Pyramid13::n_second_order_adjacent_vertices (const unsigned int n) const
288 {
289  switch (n)
290  {
291  case 5:
292  case 6:
293  case 7:
294  case 8:
295  case 9:
296  case 10:
297  case 11:
298  case 12:
299  return 2;
300 
301  default:
302  libmesh_error_msg("Invalid node n = " << n);
303  }
304 }
305 
306 
307 unsigned short int Pyramid13::second_order_adjacent_vertex (const unsigned int n,
308  const unsigned int v) const
309 {
310  libmesh_assert_greater_equal (n, this->n_vertices());
311  libmesh_assert_less (n, this->n_nodes());
312 
313  switch (n)
314  {
315  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  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  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  return node_list[n-5][v];
342  }
343 
344  default:
345  libmesh_error_msg("Invalid n = " << n);
346 
347  }
348 }
349 
350 
351 
353 {
354  // This specialization is good for Lagrange mappings only
355  if (this->mapping_type() != LAGRANGE_MAP)
356  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  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5), x6 = point(6),
362  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  x6/2 - x8/2,
369  x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - 3*x6/2 + 3*x8/2 - x9,
370  -x0/2 + x1/2 - 2*x10 - 2*x11 + 2*x12 + x2/2 - x3/2 + 3*x6/2 - 3*x8/2 + 2*x9,
371  x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
372  x0/4 - x1/4 + x2/4 - x3/4,
373  -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
374  x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
375  -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
376  x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
377  x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
378  -x0 - x1 - x2 - x3 + 2*x5 + 2*x7,
379  x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
380  -x0/2 - x1/2 + x2/2 + x3/2 + x5 - x7,
381  x0/2 + x1/2 - x2/2 - x3/2 - x5 + x7
382  };
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  -x5/2 + x7/2,
389  x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + 3*x5/2 - 3*x7/2 - x9,
390  -x0/2 - x1/2 + 2*x10 - 2*x11 - 2*x12 + x2/2 + x3/2 - 3*x5/2 + 3*x7/2 + 2*x9,
391  x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
392  x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
393  -x0 - x1 - x2 - x3 + 2*x6 + 2*x8,
394  x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
395  x0/4 - x1/4 + x2/4 - x3/4,
396  -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
397  x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
398  -x0/2 + x1/2 + x2/2 - x3/2 - x6 + x8,
399  x0/2 - x1/2 - x2/2 + x3/2 + x6 - x8,
400  -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
401  x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
402  };
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  x0/4 + x1/4 + x10 + x11 + x12 + x2/4 + x3/4 - x4 - x5 - x6 - x7 - x8 + x9,
409  -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  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  -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  2*x10 + 2*x11 + 2*x12 - 4*x4 - x5 - x6 - x7 - x8 + 2*x9,
413  x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
414  -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  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  -x0/4 - x1/4 + x10 - x11 - x12 + x2/4 + x3/4 - x5/2 + x7/2 + x9,
417  x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
418  -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  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  -x0/4 + x1/4 - x10 - x11 + x12 + x2/4 - x3/4 + x6/2 - x8/2 + x9,
421  -x0/4 + x1/4 - x10 + x11 - x12 - x2/4 + x3/4 + x9,
422  x0/4 - x1/4 + x10 - x11 + x12 + x2/4 - x3/4 - x9,
423  -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
424  x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
425  -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
426  x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
427  };
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  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  Real vol = 0.;
617  for (int q=0; q<N; ++q)
618  {
619  // Compute denominators for the current q.
620  Real
621  den2 = (1. - zeta[q][1])*(1. - zeta[q][1]),
622  den3 = den2*(1. - zeta[q][1]);
623 
624  // Compute dx/dxi and dx/deta at the current q.
625  Point dx_dxi_q, dx_deta_q;
626  for (int c=0; c<14; ++c)
627  {
628  dx_dxi_q +=
629  xi[q][dx_dxi_exponents[c][0]]*
630  eta[q][dx_dxi_exponents[c][1]]*
631  zeta[q][dx_dxi_exponents[c][2]]*dx_dxi[c];
632 
633  dx_deta_q +=
634  xi[q][dx_deta_exponents[c][0]]*
635  eta[q][dx_deta_exponents[c][1]]*
636  zeta[q][dx_deta_exponents[c][2]]*dx_deta[c];
637  }
638 
639  // Compute dx/dzeta at the current q.
640  Point dx_dzeta_q;
641  for (int c=0; c<19; ++c)
642  {
643  dx_dzeta_q +=
644  xi[q][dx_dzeta_exponents[c][0]]*
645  eta[q][dx_dzeta_exponents[c][1]]*
646  zeta[q][dx_dzeta_exponents[c][2]]*dx_dzeta[c];
647  }
648 
649  // Scale everything appropriately
650  dx_dxi_q /= den2;
651  dx_deta_q /= den2;
652  dx_dzeta_q /= den3;
653 
654  // Compute scalar triple product, multiply by weight, and accumulate volume.
655  vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
656  }
657 
658  return vol;
659 }
660 
661 
662 void Pyramid13::permute(unsigned int perm_num)
663 {
664  libmesh_assert_less (perm_num, 4);
665 
666  for (unsigned int i = 0; i != perm_num; ++i)
667  {
668  swap4nodes(0,1,2,3);
669  swap4nodes(5,6,7,8);
670  swap4nodes(9,10,11,12);
671  swap4neighbors(0,1,2,3);
672  }
673 }
674 
675 
676 void Pyramid13::flip(BoundaryInfo * boundary_info)
677 {
678  libmesh_assert(boundary_info);
679 
680  swap2nodes(0,1);
681  swap2nodes(2,3);
682  swap2nodes(6,8);
683  swap2nodes(9,10);
684  swap2nodes(11,12);
685  swap2neighbors(1,3);
686  swap2boundarysides(1,3,boundary_info);
687  swap2boundaryedges(1,3,boundary_info);
688  swap2boundaryedges(4,5,boundary_info);
689  swap2boundaryedges(6,7,boundary_info);
690 }
691 
692 
693 ElemType Pyramid13::side_type (const unsigned int s) const
694 {
695  libmesh_assert_less (s, 5);
696  if (s < 4)
697  return TRI6;
698  return QUAD8;
699 }
700 
701 
702 } // namespace libMesh
ElemType
Defines an enum for geometric element types.
void swap2boundaryedges(unsigned short e1, unsigned short e2, BoundaryInfo *boundary_info) const
Swaps two edges in boundary_info, if it is non-null.
Definition: elem.C:3550
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
virtual bool has_affine_map() const override
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
Node ** _nodes
Pointers to the nodes we are connected to.
Definition: elem.h:2245
virtual unsigned int n_nodes() const override
virtual Real volume() const override
Specialization for computing the volume of a Pyramid13.
static const int num_edges
Definition: cell_pyramid.h:78
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const override
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
static const int nodes_per_edge
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 coincident with edge i.
void swap2boundarysides(unsigned short s1, unsigned short s2, BoundaryInfo *boundary_info) const
Swaps two sides in boundary_info, if it is non-null.
Definition: elem.C:3534
virtual unsigned int n_sides() const override
Definition: cell_pyramid.h:90
static const int nodes_per_side
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int n_vertices() const override
Definition: cell_pyramid.h:95
virtual bool is_edge(const unsigned int i) const override
virtual bool is_vertex(const unsigned int i) const override
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1029
void swap4nodes(unsigned int n1, unsigned int n2, unsigned int n3, unsigned int n4)
Swaps four node_ptrs, "rotating" them.
Definition: elem.h:2143
ElemMappingType mapping_type() const
Definition: elem.h:3120
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
Definition: elem.h:2092
virtual unsigned int n_sub_elem() const override
FIXME: we don&#39;t yet have a refinement pattern for pyramids...
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libmesh_assert(ctx)
static const int num_nodes
Geometric constants for Pyramid13.
virtual void flip(BoundaryInfo *) override final
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
static const int num_sides
Geometric constants for all Pyramids.
Definition: cell_pyramid.h:77
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD8 or TRI6 coincident with face i.
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
Definition: elem.h:2102
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
virtual bool is_face(const unsigned int i) const override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
virtual unsigned int n_edges() const override
Definition: cell_pyramid.h:100
virtual void permute(unsigned int perm_num) override final
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
ElemType side_type(const unsigned int s) const override final
virtual Real volume() const
Definition: elem.C:3429
void swap4neighbors(unsigned int n1, unsigned int n2, unsigned int n3, unsigned int n4)
Swaps four neighbor_ptrs, "rotating" them.
Definition: elem.h:2153
virtual Order default_order() const override
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2453
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override