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