cell_pyramid13.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2012 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/side.h"
21 #include "libmesh/cell_pyramid13.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::num_sides;
38 const int Pyramid13::num_edges;
39 const int Pyramid13::num_children;
40 const int Pyramid13::nodes_per_side;
41 const int Pyramid13::nodes_per_edge;
42 
43 const unsigned int Pyramid13::side_nodes_map[Pyramid13::num_sides][Pyramid13::nodes_per_side] =
44  {
45  {0, 1, 4, 5, 10, 9, 99, 99}, // Side 0 (front)
46  {1, 2, 4, 6, 11, 10, 99, 99}, // Side 1 (right)
47  {2, 3, 4, 7, 12, 11, 99, 99}, // Side 2 (back)
48  {3, 0, 4, 8, 9, 12, 99, 99}, // Side 3 (left)
49  {0, 3, 2, 1, 8, 7, 6, 5} // Side 4 (base)
50  };
51 
52 const unsigned int Pyramid13::edge_nodes_map[Pyramid13::num_edges][Pyramid13::nodes_per_edge] =
53  {
54  {0, 1, 5}, // Edge 0
55  {1, 2, 6}, // Edge 1
56  {2, 3, 7}, // Edge 2
57  {0, 3, 8}, // Edge 3
58  {0, 4, 9}, // Edge 4
59  {1, 4, 10}, // Edge 5
60  {2, 4, 11}, // Edge 6
61  {3, 4, 12} // Edge 7
62  };
63 
64 
65 
66 // ------------------------------------------------------------
67 // Pyramid13 class member functions
68 
69 bool Pyramid13::is_vertex(const unsigned int i) const
70 {
71  if (i < 5)
72  return true;
73  return false;
74 }
75 
76 
77 
78 bool Pyramid13::is_edge(const unsigned int i) const
79 {
80  if (i < 5)
81  return false;
82  return true;
83 }
84 
85 
86 
87 bool Pyramid13::is_face(const unsigned int) const
88 {
89  return false;
90 }
91 
92 
93 
94 bool Pyramid13::is_node_on_side(const unsigned int n,
95  const unsigned int s) const
96 {
97  libmesh_assert_less (s, n_sides());
98  return std::find(std::begin(side_nodes_map[s]),
99  std::end(side_nodes_map[s]),
100  n) != std::end(side_nodes_map[s]);
101 }
102 
103 std::vector<unsigned>
104 Pyramid13::nodes_on_side(const unsigned int s) const
105 {
106  libmesh_assert_less(s, n_sides());
107  auto trim = (s == 4) ? 0 : 2;
108  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
109 }
110 
111 bool Pyramid13::is_node_on_edge(const unsigned int n,
112  const unsigned int e) const
113 {
114  libmesh_assert_less (e, n_edges());
115  return std::find(std::begin(edge_nodes_map[e]),
116  std::end(edge_nodes_map[e]),
117  n) != std::end(edge_nodes_map[e]);
118 }
119 
120 
121 
122 bool Pyramid13::has_affine_map() const
123 {
124  // TODO: If the base is a parallelogram and all the triangular faces are planar,
125  // the map should be linear, but I need to test this theory...
126  return false;
127 }
128 
129 
130 
131 Order Pyramid13::default_order() const
132 {
133  return SECOND;
134 }
135 
136 
137 
138 unsigned int Pyramid13::which_node_am_i(unsigned int side,
139  unsigned int side_node) const
140 {
141  libmesh_assert_less (side, this->n_sides());
142 
143  // Never more than 8 nodes per side.
144  libmesh_assert_less(side_node, Pyramid13::nodes_per_side);
145 
146  // Some sides have 6 nodes.
147  libmesh_assert(side == 4 || side_node < 6);
148 
149  return Pyramid13::side_nodes_map[side][side_node];
150 }
151 
152 
153 
154 std::unique_ptr<Elem> Pyramid13::build_side_ptr (const unsigned int i, bool proxy)
155 {
156  libmesh_assert_less (i, this->n_sides());
157 
158  if (proxy)
159  {
160  switch (i)
161  {
162  case 0:
163  case 1:
164  case 2:
165  case 3:
166  return libmesh_make_unique<Side<Tri6,Pyramid13>>(this,i);
167 
168  case 4:
169  return libmesh_make_unique<Side<Quad8,Pyramid13>>(this,i);
170 
171  default:
172  libmesh_error_msg("Invalid side i = " << i);
173  }
174  }
175 
176  else
177  {
178  // Return value
179  std::unique_ptr<Elem> face;
180 
181  switch (i)
182  {
183  case 0: // triangular face 1
184  case 1: // triangular face 2
185  case 2: // triangular face 3
186  case 3: // triangular face 4
187  {
188  face = libmesh_make_unique<Tri6>();
189  break;
190  }
191  case 4: // the quad face at z=0
192  {
193  face = libmesh_make_unique<Quad8>();
194  break;
195  }
196  default:
197  libmesh_error_msg("Invalid side i = " << i);
198  }
199 
200  face->subdomain_id() = this->subdomain_id();
201 
202  // Set the nodes
203  for (unsigned n=0; n<face->n_nodes(); ++n)
204  face->set_node(n) = this->node_ptr(Pyramid13::side_nodes_map[i][n]);
205 
206  return face;
207  }
208 }
209 
210 
211 
212 void Pyramid13::build_side_ptr (std::unique_ptr<Elem> & side,
213  const unsigned int i)
214 {
215  libmesh_assert_less (i, this->n_sides());
216 
217  switch (i)
218  {
219  case 0: // triangular face 1
220  case 1: // triangular face 2
221  case 2: // triangular face 3
222  case 3: // triangular face 4
223  {
224  if (!side.get() || side->type() != TRI6)
225  {
226  side = this->build_side_ptr(i, false);
227  return;
228  }
229  break;
230  }
231  case 4: // the quad face at z=0
232  {
233  if (!side.get() || side->type() != QUAD8)
234  {
235  side = this->build_side_ptr(i, false);
236  return;
237  }
238  break;
239  }
240  default:
241  libmesh_error_msg("Invalid side i = " << i);
242  }
243 
244  side->subdomain_id() = this->subdomain_id();
245 
246  // Set the nodes
247  for (auto n : side->node_index_range())
248  side->set_node(n) = this->node_ptr(Pyramid13::side_nodes_map[i][n]);
249 }
250 
251 
252 
253 std::unique_ptr<Elem> Pyramid13::build_edge_ptr (const unsigned int i)
254 {
255  libmesh_assert_less (i, this->n_edges());
256 
257  return libmesh_make_unique<SideEdge<Edge3,Pyramid13>>(this,i);
258 }
259 
260 
261 
262 void Pyramid13::connectivity(const unsigned int libmesh_dbg_var(sc),
263  const IOPackage iop,
264  std::vector<dof_id_type> & /*conn*/) const
265 {
266  libmesh_assert(_nodes);
267  libmesh_assert_less (sc, this->n_sub_elem());
268  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
269 
270  switch (iop)
271  {
272  case TECPLOT:
273  {
274  // TODO
275  libmesh_not_implemented();
276  }
277 
278  case VTK:
279  {
280  // TODO
281  libmesh_not_implemented();
282  }
283 
284  default:
285  libmesh_error_msg("Unsupported IO package " << iop);
286  }
287 }
288 
289 
290 
291 unsigned int Pyramid13::n_second_order_adjacent_vertices (const unsigned int n) const
292 {
293  switch (n)
294  {
295  case 5:
296  case 6:
297  case 7:
298  case 8:
299  case 9:
300  case 10:
301  case 11:
302  case 12:
303  return 2;
304 
305  default:
306  libmesh_error_msg("Invalid node n = " << n);
307  }
308 }
309 
310 
311 unsigned short int Pyramid13::second_order_adjacent_vertex (const unsigned int n,
312  const unsigned int v) const
313 {
314  libmesh_assert_greater_equal (n, this->n_vertices());
315  libmesh_assert_less (n, this->n_nodes());
316 
317  switch (n)
318  {
319  case 5:
320  case 6:
321  case 7:
322  case 8:
323  case 9:
324  case 10:
325  case 11:
326  case 12:
327  {
328  libmesh_assert_less (v, 2);
329 
330  // This is the analog of the static, const arrays
331  // {Hex,Prism,Tet10}::_second_order_adjacent_vertices
332  // defined in the respective source files...
333  unsigned short node_list[8][2] =
334  {
335  {0,1},
336  {1,2},
337  {2,3},
338  {0,3},
339  {0,4},
340  {1,4},
341  {2,4},
342  {3,4}
343  };
344 
345  return node_list[n-5][v];
346  }
347 
348  default:
349  libmesh_error_msg("Invalid n = " << n);
350 
351  }
352 }
353 
354 
355 
356 Real Pyramid13::volume () const
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 = Real(-7.1805574131988893873307823958101e-01L),
498  a2 = Real(-5.0580870785392503961340276902425e-01L),
499  a3 = Real(-2.2850430565396735359574351631314e-01L),
500  // Parameters used for zeta quadrature points.
501  b1 = Real( 7.2994024073149732155837979012003e-02L),
502  b2 = Real( 3.4700376603835188472176354340395e-01L),
503  b3 = Real( 7.0500220988849838312239847758405e-01L),
504  // There are 9 unique weight values since there are three
505  // for each of the three unique zeta values.
506  w1 = Real( 4.8498876871878584357513834016440e-02L),
507  w2 = Real( 4.5137737425884574692441981593901e-02L),
508  w3 = Real( 9.2440441384508327195915094925393e-03L),
509  w4 = Real( 7.7598202995005734972022134426305e-02L),
510  w5 = Real( 7.2220379881415319507907170550242e-02L),
511  w6 = Real( 1.4790470621521332351346415188063e-02L),
512  w7 = Real( 1.2415712479200917595523541508209e-01L),
513  w8 = Real( 1.1555260781026451121265147288039e-01L),
514  w9 = Real( 2.3664752994434131762154264300901e-02L);
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 } // namespace libMesh
libMesh::Pyramid13::edge_nodes_map
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
Definition: cell_pyramid13.h:205
libMesh::Pyramid13::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
Definition: cell_pyramid13.h:199
libMesh::Pyramid13::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: cell_pyramid13.C:104
libMesh::Pyramid13::has_affine_map
virtual bool has_affine_map() const override
Definition: cell_pyramid13.C:122
libMesh::Elem::_nodes
Node ** _nodes
Definition: elem.h:1695
libMesh::Pyramid13::n_nodes
virtual unsigned int n_nodes() const override
Definition: cell_pyramid13.h:89
libMesh::Pyramid13::volume
virtual Real volume() const override
Definition: cell_pyramid13.C:356
libMesh::Pyramid13::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: cell_pyramid13.C:94
libMesh::Pyramid13::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: cell_pyramid13.C:262
cell_pyramid13.h
face_tri6.h
libMesh::Pyramid13::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const override
Definition: cell_pyramid13.C:291
libMesh::Pyramid13::nodes_per_edge
static const int nodes_per_edge
Definition: cell_pyramid13.h:193
side
unsigned short int side
Definition: xdr_io.C:50
libMesh::Pyramid13::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Definition: cell_pyramid13.C:253
libMesh::Pyramid::n_sides
virtual unsigned int n_sides() const override
Definition: cell_pyramid.h:83
libMesh::Pyramid13::nodes_per_side
static const int nodes_per_side
Definition: cell_pyramid13.h:192
end
IterBase * end
Definition: variant_filter_iterator.h:343
libMesh::Pyramid13::num_children
static const int num_children
Definition: cell_pyramid13.h:191
libMesh::Pyramid::n_vertices
virtual unsigned int n_vertices() const override
Definition: cell_pyramid.h:88
libMesh::Pyramid13::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition: cell_pyramid13.C:78
libMesh::Pyramid13::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition: cell_pyramid13.C:69
face_quad8.h
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1054
edge_edge3.h
libMesh::Pyramid13::n_sub_elem
virtual unsigned int n_sub_elem() const override
Definition: cell_pyramid13.h:100
libMesh::Pyramid13::which_node_am_i
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const override
Definition: cell_pyramid13.C:138
libMesh::Pyramid13::num_sides
static const int num_sides
Definition: cell_pyramid13.h:189
libMesh::Pyramid13::num_nodes
static const int num_nodes
Definition: cell_pyramid13.h:188
libMesh::Pyramid13::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy) override
Definition: cell_pyramid13.C:154
enum_io_package.h
libMesh::Pyramid13::is_face
virtual bool is_face(const unsigned int i) const override
Definition: cell_pyramid13.C:87
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::Pyramid13::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition: cell_pyramid13.C:311
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
libMesh::Pyramid::n_edges
virtual unsigned int n_edges() const override
Definition: cell_pyramid.h:93
libMesh::Pyramid13::num_edges
static const int num_edges
Definition: cell_pyramid13.h:190
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1957
libMesh::Pyramid13::default_order
virtual Order default_order() const override
Definition: cell_pyramid13.C:131
enum_order.h
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
libMesh::Elem::side
std::unique_ptr< Elem > side(const unsigned int i) const
Definition: elem.h:2202
libMesh::Pyramid13::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override
Definition: cell_pyramid13.C:111
generated by doxygen