libMesh
cell_pyramid13.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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 (auto n : face->node_index_range())
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  // This specialization is good for Lagrange mappings only
359  if (this->mapping_type() != LAGRANGE_MAP)
360  return this->Elem::volume();
361 
362  // Make copies of our points. It makes the subsequent calculations a bit
363  // shorter and avoids dereferencing the same pointer multiple times.
364  Point
365  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4), x5 = point(5), x6 = point(6),
366  x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11), x12 = point(12);
367 
368  // dx/dxi and dx/deta have 14 components while dx/dzeta has 19.
369  // These are copied directly from the output of a Python script.
370  Point dx_dxi[14] =
371  {
372  x6/2 - x8/2,
373  x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - 3*x6/2 + 3*x8/2 - x9,
374  -x0/2 + x1/2 - 2*x10 - 2*x11 + 2*x12 + x2/2 - x3/2 + 3*x6/2 - 3*x8/2 + 2*x9,
375  x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
376  x0/4 - x1/4 + x2/4 - x3/4,
377  -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
378  x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
379  -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
380  x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
381  x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
382  -x0 - x1 - x2 - x3 + 2*x5 + 2*x7,
383  x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
384  -x0/2 - x1/2 + x2/2 + x3/2 + x5 - x7,
385  x0/2 + x1/2 - x2/2 - x3/2 - x5 + x7
386  };
387 
388  // dx/dxi and dx/deta have 14 components while dx/dzeta has 19.
389  // These are copied directly from the output of a Python script.
390  Point dx_deta[14] =
391  {
392  -x5/2 + x7/2,
393  x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + 3*x5/2 - 3*x7/2 - x9,
394  -x0/2 - x1/2 + 2*x10 - 2*x11 - 2*x12 + x2/2 + x3/2 - 3*x5/2 + 3*x7/2 + 2*x9,
395  x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
396  x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
397  -x0 - x1 - x2 - x3 + 2*x6 + 2*x8,
398  x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
399  x0/4 - x1/4 + x2/4 - x3/4,
400  -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
401  x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
402  -x0/2 + x1/2 + x2/2 - x3/2 - x6 + x8,
403  x0/2 - x1/2 - x2/2 + x3/2 + x6 - x8,
404  -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
405  x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
406  };
407 
408  // dx/dxi and dx/deta have 14 components while dx/dzeta has 19.
409  // These are copied directly from the output of a Python script.
410  Point dx_dzeta[19] =
411  {
412  x0/4 + x1/4 + x10 + x11 + x12 + x2/4 + x3/4 - x4 - x5 - x6 - x7 - x8 + x9,
413  -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,
414  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,
415  -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,
416  2*x10 + 2*x11 + 2*x12 - 4*x4 - x5 - x6 - x7 - x8 + 2*x9,
417  x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
418  -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,
419  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,
420  -x0/4 - x1/4 + x10 - x11 - x12 + x2/4 + x3/4 - x5/2 + x7/2 + x9,
421  x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
422  -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,
423  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,
424  -x0/4 + x1/4 - x10 - x11 + x12 + x2/4 - x3/4 + x6/2 - x8/2 + x9,
425  -x0/4 + x1/4 - x10 + x11 - x12 - x2/4 + x3/4 + x9,
426  x0/4 - x1/4 + x10 - x11 + x12 + x2/4 - x3/4 - x9,
427  -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
428  x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
429  -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
430  x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
431  };
432 
433  // The (xi, eta, zeta) exponents for each of the dx_dxi terms
434  static const int dx_dxi_exponents[14][3] =
435  {
436  {0, 0, 0},
437  {0, 0, 1},
438  {0, 0, 2},
439  {0, 0, 3},
440  {0, 1, 0},
441  {0, 1, 1},
442  {0, 1, 2},
443  {0, 2, 0},
444  {0, 2, 1},
445  {1, 0, 0},
446  {1, 0, 1},
447  {1, 0, 2},
448  {1, 1, 0},
449  {1, 1, 1}
450  };
451 
452  // The (xi, eta, zeta) exponents for each of the dx_deta terms
453  static const int dx_deta_exponents[14][3] =
454  {
455  {0, 0, 0},
456  {0, 0, 1},
457  {0, 0, 2},
458  {0, 0, 3},
459  {0, 1, 0},
460  {0, 1, 1},
461  {0, 1, 2},
462  {1, 0, 0},
463  {1, 0, 1},
464  {1, 0, 2},
465  {1, 1, 0},
466  {1, 1, 1},
467  {2, 0, 0},
468  {2, 0, 1}
469  };
470 
471  // The (xi, eta, zeta) exponents for each of the dx_dzeta terms
472  static const int dx_dzeta_exponents[19][3] =
473  {
474  {0, 0, 0},
475  {0, 0, 1},
476  {0, 0, 2},
477  {0, 0, 3},
478  {0, 0, 4},
479  {0, 1, 0},
480  {0, 1, 1},
481  {0, 1, 2},
482  {0, 1, 3},
483  {1, 0, 0},
484  {1, 0, 1},
485  {1, 0, 2},
486  {1, 0, 3},
487  {1, 1, 0},
488  {1, 1, 1},
489  {1, 2, 0},
490  {1, 2, 1},
491  {2, 1, 0},
492  {2, 1, 1},
493  };
494 
495  // Number of points in the quadrature rule
496  const int N = 27;
497 
498  // Parameters of the quadrature rule
499  static const Real
500  // Parameters used for (xi, eta) quadrature points.
501  a1 = Real(-7.1805574131988893873307823958101e-01L),
502  a2 = Real(-5.0580870785392503961340276902425e-01L),
503  a3 = Real(-2.2850430565396735359574351631314e-01L),
504  // Parameters used for zeta quadrature points.
505  b1 = Real( 7.2994024073149732155837979012003e-02L),
506  b2 = Real( 3.4700376603835188472176354340395e-01L),
507  b3 = Real( 7.0500220988849838312239847758405e-01L),
508  // There are 9 unique weight values since there are three
509  // for each of the three unique zeta values.
510  w1 = Real( 4.8498876871878584357513834016440e-02L),
511  w2 = Real( 4.5137737425884574692441981593901e-02L),
512  w3 = Real( 9.2440441384508327195915094925393e-03L),
513  w4 = Real( 7.7598202995005734972022134426305e-02L),
514  w5 = Real( 7.2220379881415319507907170550242e-02L),
515  w6 = Real( 1.4790470621521332351346415188063e-02L),
516  w7 = Real( 1.2415712479200917595523541508209e-01L),
517  w8 = Real( 1.1555260781026451121265147288039e-01L),
518  w9 = Real( 2.3664752994434131762154264300901e-02L);
519 
520  // The points and weights of the 3x3x3 quadrature rule
521  static const Real xi[N][3] =
522  {// ^0 ^1 ^2
523  { 1., a1, a1*a1},
524  { 1., a2, a2*a2},
525  { 1., a3, a3*a3},
526  { 1., a1, a1*a1},
527  { 1., a2, a2*a2},
528  { 1., a3, a3*a3},
529  { 1., a1, a1*a1},
530  { 1., a2, a2*a2},
531  { 1., a3, a3*a3},
532  { 1., 0., 0. },
533  { 1., 0., 0. },
534  { 1., 0., 0. },
535  { 1., 0., 0. },
536  { 1., 0., 0. },
537  { 1., 0., 0. },
538  { 1., 0., 0. },
539  { 1., 0., 0. },
540  { 1., 0., 0. },
541  { 1., -a1, a1*a1},
542  { 1., -a2, a2*a2},
543  { 1., -a3, a3*a3},
544  { 1., -a1, a1*a1},
545  { 1., -a2, a2*a2},
546  { 1., -a3, a3*a3},
547  { 1., -a1, a1*a1},
548  { 1., -a2, a2*a2},
549  { 1., -a3, a3*a3}
550  };
551 
552  static const Real eta[N][3] =
553  {// ^0 ^1 ^2
554  { 1., a1, a1*a1},
555  { 1., a2, a2*a2},
556  { 1., a3, a3*a3},
557  { 1., 0., 0. },
558  { 1., 0., 0. },
559  { 1., 0., 0. },
560  { 1., -a1, a1*a1},
561  { 1., -a2, a2*a2},
562  { 1., -a3, a3*a3},
563  { 1., a1, a1*a1},
564  { 1., a2, a2*a2},
565  { 1., a3, a3*a3},
566  { 1., 0., 0. },
567  { 1., 0., 0. },
568  { 1., 0., 0. },
569  { 1., -a1, a1*a1},
570  { 1., -a2, a2*a2},
571  { 1., -a3, a3*a3},
572  { 1., a1, a1*a1},
573  { 1., a2, a2*a2},
574  { 1., a3, a3*a3},
575  { 1., 0., 0. },
576  { 1., 0., 0. },
577  { 1., 0., 0. },
578  { 1., -a1, a1*a1},
579  { 1., -a2, a2*a2},
580  { 1., -a3, a3*a3}
581  };
582 
583  static const Real zeta[N][5] =
584  {// ^0 ^1 ^2 ^3 ^4
585  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
586  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
587  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
588  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
589  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
590  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
591  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
592  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
593  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
594  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
595  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
596  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
597  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
598  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
599  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
600  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
601  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
602  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
603  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
604  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
605  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
606  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
607  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
608  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
609  { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
610  { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
611  { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3}
612  };
613 
614  static const Real w[N] = {w1, w2, w3, w4, w5, w6, // 0-5
615  w1, w2, w3, w4, w5, w6, // 6-11
616  w7, w8, w9, w4, w5, w6, // 12-17
617  w1, w2, w3, w4, w5, w6, // 18-23
618  w1, w2, w3}; // 24-26
619 
620  Real vol = 0.;
621  for (int q=0; q<N; ++q)
622  {
623  // Compute denominators for the current q.
624  Real
625  den2 = (1. - zeta[q][1])*(1. - zeta[q][1]),
626  den3 = den2*(1. - zeta[q][1]);
627 
628  // Compute dx/dxi and dx/deta at the current q.
629  Point dx_dxi_q, dx_deta_q;
630  for (int c=0; c<14; ++c)
631  {
632  dx_dxi_q +=
633  xi[q][dx_dxi_exponents[c][0]]*
634  eta[q][dx_dxi_exponents[c][1]]*
635  zeta[q][dx_dxi_exponents[c][2]]*dx_dxi[c];
636 
637  dx_deta_q +=
638  xi[q][dx_deta_exponents[c][0]]*
639  eta[q][dx_deta_exponents[c][1]]*
640  zeta[q][dx_deta_exponents[c][2]]*dx_deta[c];
641  }
642 
643  // Compute dx/dzeta at the current q.
644  Point dx_dzeta_q;
645  for (int c=0; c<19; ++c)
646  {
647  dx_dzeta_q +=
648  xi[q][dx_dzeta_exponents[c][0]]*
649  eta[q][dx_dzeta_exponents[c][1]]*
650  zeta[q][dx_dzeta_exponents[c][2]]*dx_dzeta[c];
651  }
652 
653  // Scale everything appropriately
654  dx_dxi_q /= den2;
655  dx_deta_q /= den2;
656  dx_dzeta_q /= den3;
657 
658  // Compute scalar triple product, multiply by weight, and accumulate volume.
659  vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
660  }
661 
662  return vol;
663 }
664 
665 } // namespace libMesh
libMesh::Pyramid13::nodes_per_edge
static const int nodes_per_edge
Definition: cell_pyramid13.h:193
libMesh::Pyramid13::n_nodes
virtual unsigned int n_nodes() const override
Definition: cell_pyramid13.h:89
libMesh::Pyramid13::edge_nodes_map
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
Definition: cell_pyramid13.h:205
libMesh::IOPackage
IOPackage
libMesh interfaces with several different software packages for the purposes of creating,...
Definition: enum_io_package.h:37
libMesh::Pyramid13::side_nodes_map
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
Definition: cell_pyramid13.h:199
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
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::Order
Order
Definition: enum_order.h:40
libMesh::INVALID_IO_PACKAGE
Definition: enum_io_package.h:48
end
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end.
Definition: variant_filter_iterator.h:343
libMesh::Pyramid13::default_order
virtual Order default_order() const override
Definition: cell_pyramid13.C:131
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::connectivity
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: cell_pyramid13.C:262
libMesh::triple_product
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1106
libMesh::Pyramid::n_edges
virtual unsigned int n_edges() const override
Definition: cell_pyramid.h:93
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
libMesh::SECOND
Definition: enum_order.h:43
libMesh::Elem::_nodes
Node ** _nodes
Pointers to the nodes we are connected to.
Definition: elem.h:1743
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1955
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Pyramid::n_vertices
virtual unsigned int n_vertices() const override
Definition: cell_pyramid.h:88
libMesh::Elem::volume
virtual Real volume() const
Definition: elem.C:2617
libMesh::Pyramid13::num_nodes
static const int num_nodes
Geometric constants for Pyramid13.
Definition: cell_pyramid13.h:188
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::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition: elem.h:2524
libMesh::Pyramid13::volume
virtual Real volume() const override
Specialization for computing the volume of a Pyramid13.
Definition: cell_pyramid13.C:356
libMesh::Pyramid13::nodes_per_side
static const int nodes_per_side
Definition: cell_pyramid13.h:192
libMesh::Pyramid13::num_edges
static const int num_edges
Definition: cell_pyramid13.h:190
libMesh::Pyramid13::n_sub_elem
virtual unsigned int n_sub_elem() const override
FIXME: we don't yet have a refinement pattern for pyramids...
Definition: cell_pyramid13.h:100
libMesh::Pyramid::n_sides
virtual unsigned int n_sides() const override
Definition: cell_pyramid.h:83
libMesh::Pyramid13::is_edge
virtual bool is_edge(const unsigned int i) const override
Definition: cell_pyramid13.C:78
libMesh::Pyramid13::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true) override
Builds a QUAD8 or TRI6 coincident with face i.
Definition: cell_pyramid13.C:154
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::Pyramid13::is_face
virtual bool is_face(const unsigned int i) const override
Definition: cell_pyramid13.C:87
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::LAGRANGE_MAP
Definition: enum_elem_type.h:83
libMesh::Pyramid13::has_affine_map
virtual bool has_affine_map() const override
Definition: cell_pyramid13.C:122
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition: elem.h:2069
libMesh::Pyramid13::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 coincident with edge i.
Definition: cell_pyramid13.C:253
libMesh::Pyramid13::num_children
static const int num_children
Definition: cell_pyramid13.h:191
libMesh::VTK
Definition: enum_io_package.h:42
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::num_sides
static const int num_sides
Definition: cell_pyramid13.h:189
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2009
libMesh::Pyramid13::is_vertex
virtual bool is_vertex(const unsigned int i) const override
Definition: cell_pyramid13.C:69
libMesh::TECPLOT
Definition: enum_io_package.h:39
libMesh::QUAD8
Definition: enum_elem_type.h:42
Generated on Sat Jan 25 2020 12:06:59 for libMesh by   doxygen 1.8.16