libMesh
fe_bernstein_shape_3D.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 // C++ includes
20 
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
25 
26 #include "libmesh/fe.h"
27 #include "libmesh/elem.h"
28 
29 
30 namespace libMesh
31 {
32 
33 
34 
35 template <>
36 Real FE<3,BERNSTEIN>::shape(const ElemType,
37  const Order,
38  const unsigned int,
39  const Point &)
40 {
41  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
42  return 0.;
43 }
44 
45 
46 
47 template <>
48 Real FE<3,BERNSTEIN>::shape(const Elem * elem,
49  const Order order,
50  const unsigned int i,
51  const Point & p,
52  const bool add_p_level)
53 {
54 
55 #if LIBMESH_DIM == 3
56 
57  libmesh_assert(elem);
58  const ElemType type = elem->type();
59 
60  const Order totalorder =
61  static_cast<Order>(order + add_p_level * elem->p_level());
62 
63  switch (totalorder)
64  {
65  // 1st order Bernstein.
66  case FIRST:
67  {
68  switch (type)
69  {
70 
71  // Bernstein shape functions on the tetrahedron.
72  case TET4:
73  case TET10:
74  {
75  libmesh_assert_less (i, 4);
76 
77  // Area coordinates
78  const Real zeta1 = p(0);
79  const Real zeta2 = p(1);
80  const Real zeta3 = p(2);
81  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
82 
83  switch(i)
84  {
85  case 0: return zeta0;
86  case 1: return zeta1;
87  case 2: return zeta2;
88  case 3: return zeta3;
89 
90  default:
91  libmesh_error_msg("Invalid shape function index i = " << i);
92  }
93  }
94 
95  // Bernstein shape functions on the hexahedral.
96  case HEX8:
97  case HEX20:
98  case HEX27:
99  {
100  libmesh_assert_less (i, 8);
101 
102  // Compute hex shape functions as a tensor-product
103  const Real xi = p(0);
104  const Real eta = p(1);
105  const Real zeta = p(2);
106 
107  // The only way to make any sense of this
108  // is to look at the mgflo/mg2/mgf documentation
109  // and make the cut-out cube!
110  // 0 1 2 3 4 5 6 7
111  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
112  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
113  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
114 
115  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
116  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
117  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
118  }
119 
120 
121  default:
122  libmesh_error_msg("Invalid element type = " << type);
123  }
124  }
125 
126 
127 
128 
129  case SECOND:
130  {
131  switch (type)
132  {
133 
134  // Bernstein shape functions on the tetrahedron.
135  case TET10:
136  {
137  libmesh_assert_less (i, 10);
138 
139  // Area coordinates
140  const Real zeta1 = p(0);
141  const Real zeta2 = p(1);
142  const Real zeta3 = p(2);
143  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
144 
145  switch(i)
146  {
147  case 0: return zeta0*zeta0;
148  case 1: return zeta1*zeta1;
149  case 2: return zeta2*zeta2;
150  case 3: return zeta3*zeta3;
151  case 4: return 2.*zeta0*zeta1;
152  case 5: return 2.*zeta1*zeta2;
153  case 6: return 2.*zeta0*zeta2;
154  case 7: return 2.*zeta3*zeta0;
155  case 8: return 2.*zeta1*zeta3;
156  case 9: return 2.*zeta2*zeta3;
157 
158  default:
159  libmesh_error_msg("Invalid shape function index i = " << i);
160  }
161  }
162 
163  // Bernstein shape functions on the 20-noded hexahedral.
164  case HEX20:
165  {
166  libmesh_assert_less (i, 20);
167 
168  // Compute hex shape functions as a tensor-product
169  const Real xi = p(0);
170  const Real eta = p(1);
171  const Real zeta = p(2);
172 
173  // The only way to make any sense of this
174  // is to look at the mgflo/mg2/mgf documentation
175  // and make the cut-out cube!
176  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
177  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
178  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
179  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
180  //To compute the hex20 shape functions the original shape functions for hex27 are used.
181  //scalx[i] tells how often the original x-th shape function has to be added to the original i-th shape function
182  //to compute the new i-th shape function for hex20
183  //example: B_0^HEX20 = B_0^HEX27 - 0.25*B_20^HEX27 - 0.25*B_21^HEX27 + 0*B_22^HEX27 + 0*B_23^HEX27 - 0.25*B_24^HEX27 + 0*B_25^HEX27 - 0.25*B_26^HEX27
184  // B_0^HEX20 = B_0^HEX27 + scal20[0]*B_20^HEX27 + scal21[0]*B_21^HEX27 + ...
185  static const Real scal20[] = {-0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0};
186  static const Real scal21[] = {-0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0};
187  static const Real scal22[] = {0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0};
188  static const Real scal23[] = {0, 0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0};
189  static const Real scal24[] = {-0.25, 0, 0, -0.25, -0.25, 0, 0, -0.25, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0, 0.5};
190  static const Real scal25[] = {0, 0, 0, 0, -0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5};
191  static const Real scal26[] = {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
192 
193  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
194  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
195  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta)
196  +scal20[i]*
197  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[20], xi)*
198  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[20], eta)*
199  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[20], zeta)
200  +scal21[i]*
201  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[21], xi)*
202  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[21], eta)*
203  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[21], zeta)
204  +scal22[i]*
205  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[22], xi)*
206  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[22], eta)*
207  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[22], zeta)
208  +scal23[i]*
209  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[23], xi)*
210  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[23], eta)*
211  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[23], zeta)
212  +scal24[i]*
213  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[24], xi)*
214  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[24], eta)*
215  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[24], zeta)
216  +scal25[i]*
217  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[25], xi)*
218  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[25], eta)*
219  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[25], zeta)
220  +scal26[i]*
221  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[26], xi)*
222  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[26], eta)*
223  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[26], zeta));
224  }
225 
226  // Bernstein shape functions on the hexahedral.
227  case HEX27:
228  {
229  libmesh_assert_less (i, 27);
230 
231  // Compute hex shape functions as a tensor-product
232  const Real xi = p(0);
233  const Real eta = p(1);
234  const Real zeta = p(2);
235 
236  // The only way to make any sense of this
237  // is to look at the mgflo/mg2/mgf documentation
238  // and make the cut-out cube!
239  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
240  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
241  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
242  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
243 
244  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
245  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
246  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
247  }
248 
249 
250  default:
251  libmesh_error_msg("Invalid element type = " << type);
252  }
253 
254  }
255 
256 
257 
258  // 3rd-order Bernstein.
259  case THIRD:
260  {
261  switch (type)
262  {
263 
264  // // Bernstein shape functions on the tetrahedron.
265  // case TET10:
266  // {
267  // libmesh_assert_less (i, 20);
268 
269  // // Area coordinates
270  // const Real zeta1 = p(0);
271  // const Real zeta2 = p(1);
272  // const Real zeta3 = p(2);
273  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
274 
275 
276  // unsigned int shape=i;
277 
278  // // handle the edge orientation
279 
280  // if ((i== 4||i== 5) && elem->node_id(0) > elem->node_id(1))shape= 9-i; //Edge 0
281  // if ((i== 6||i== 7) && elem->node_id(1) > elem->node_id(2))shape=13-i; //Edge 1
282  // if ((i== 8||i== 9) && elem->node_id(0) > elem->node_id(2))shape=17-i; //Edge 2
283  // if ((i==10||i==11) && elem->node_id(0) > elem->node_id(3))shape=21-i; //Edge 3
284  // if ((i==12||i==13) && elem->node_id(1) > elem->node_id(3))shape=25-i; //Edge 4
285  // if ((i==14||i==15) && elem->node_id(2) > elem->node_id(3))shape=29-i; //Edge 5
286 
287  // // No need to handle face orientation in 3rd order.
288 
289 
290  // switch(shape)
291  // {
292  // //point function
293  // case 0: return zeta0*zeta0*zeta0;
294  // case 1: return zeta1*zeta1*zeta1;
295  // case 2: return zeta2*zeta2*zeta2;
296  // case 3: return zeta3*zeta3*zeta3;
297 
298  // //edge functions
299  // case 4: return 3.*zeta0*zeta0*zeta1;
300  // case 5: return 3.*zeta1*zeta1*zeta0;
301 
302  // case 6: return 3.*zeta1*zeta1*zeta2;
303  // case 7: return 3.*zeta2*zeta2*zeta1;
304 
305  // case 8: return 3.*zeta0*zeta0*zeta2;
306  // case 9: return 3.*zeta2*zeta2*zeta0;
307 
308  // case 10: return 3.*zeta0*zeta0*zeta3;
309  // case 11: return 3.*zeta3*zeta3*zeta0;
310 
311  // case 12: return 3.*zeta1*zeta1*zeta3;
312  // case 13: return 3.*zeta3*zeta3*zeta1;
313 
314  // case 14: return 3.*zeta2*zeta2*zeta3;
315  // case 15: return 3.*zeta3*zeta3*zeta2;
316 
317  // //face functions
318  // case 16: return 6.*zeta0*zeta1*zeta2;
319  // case 17: return 6.*zeta0*zeta1*zeta3;
320  // case 18: return 6.*zeta1*zeta2*zeta3;
321  // case 19: return 6.*zeta2*zeta0*zeta3;
322 
323  // default:
324  // libmesh_error_msg("Invalid shape function index i = " << i);
325  // }
326  // }
327 
328 
329  // Bernstein shape functions on the hexahedral.
330  case HEX27:
331  {
332  libmesh_assert_less (i, 64);
333 
334  // Compute hex shape functions as a tensor-product
335  const Real xi = p(0);
336  const Real eta = p(1);
337  const Real zeta = p(2);
338  Real xi_mapped = p(0);
339  Real eta_mapped = p(1);
340  Real zeta_mapped = p(2);
341 
342  // The only way to make any sense of this
343  // is to look at the mgflo/mg2/mgf documentation
344  // and make the cut-out cube!
345  // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
346  // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63
347  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
348  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
349  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
350 
351 
352 
353  // handle the edge orientation
354  {
355  // Edge 0
356  if ((i0[i] >= 2) && (i1[i] == 0) && (i2[i] == 0))
357  {
358  if (elem->point(0) != std::min(elem->point(0), elem->point(1)))
359  xi_mapped = -xi;
360  }
361  // Edge 1
362  else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] == 0))
363  {
364  if (elem->point(1) != std::min(elem->point(1), elem->point(2)))
365  eta_mapped = -eta;
366  }
367  // Edge 2
368  else if ((i0[i] >= 2) && (i1[i] == 1) && (i2[i] == 0))
369  {
370  if (elem->point(3) != std::min(elem->point(3), elem->point(2)))
371  xi_mapped = -xi;
372  }
373  // Edge 3
374  else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] == 0))
375  {
376  if (elem->point(0) != std::min(elem->point(0), elem->point(3)))
377  eta_mapped = -eta;
378  }
379  // Edge 4
380  else if ((i0[i] == 0) && (i1[i] == 0) && (i2[i] >=2 ))
381  {
382  if (elem->point(0) != std::min(elem->point(0), elem->point(4)))
383  zeta_mapped = -zeta;
384  }
385  // Edge 5
386  else if ((i0[i] == 1) && (i1[i] == 0) && (i2[i] >=2 ))
387  {
388  if (elem->point(1) != std::min(elem->point(1), elem->point(5)))
389  zeta_mapped = -zeta;
390  }
391  // Edge 6
392  else if ((i0[i] == 1) && (i1[i] == 1) && (i2[i] >=2 ))
393  {
394  if (elem->point(2) != std::min(elem->point(2), elem->point(6)))
395  zeta_mapped = -zeta;
396  }
397  // Edge 7
398  else if ((i0[i] == 0) && (i1[i] == 1) && (i2[i] >=2 ))
399  {
400  if (elem->point(3) != std::min(elem->point(3), elem->point(7)))
401  zeta_mapped = -zeta;
402  }
403  // Edge 8
404  else if ((i0[i] >=2 ) && (i1[i] == 0) && (i2[i] == 1))
405  {
406  if (elem->point(4) != std::min(elem->point(4), elem->point(5)))
407  xi_mapped = -xi;
408  }
409  // Edge 9
410  else if ((i0[i] == 1) && (i1[i] >=2 ) && (i2[i] == 1))
411  {
412  if (elem->point(5) != std::min(elem->point(5), elem->point(6)))
413  eta_mapped = -eta;
414  }
415  // Edge 10
416  else if ((i0[i] >=2 ) && (i1[i] == 1) && (i2[i] == 1))
417  {
418  if (elem->point(7) != std::min(elem->point(7), elem->point(6)))
419  xi_mapped = -xi;
420  }
421  // Edge 11
422  else if ((i0[i] == 0) && (i1[i] >=2 ) && (i2[i] == 1))
423  {
424  if (elem->point(4) != std::min(elem->point(4), elem->point(7)))
425  eta_mapped = -eta;
426  }
427  }
428 
429 
430  // handle the face orientation
431  {
432  // Face 0
433  if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
434  {
435  const Point min_point = std::min(elem->point(1),
436  std::min(elem->point(2),
437  std::min(elem->point(0),
438  elem->point(3))));
439  if (elem->point(0) == min_point)
440  if (elem->point(1) == std::min(elem->point(1), elem->point(3)))
441  {
442  // Case 1
443  xi_mapped = xi;
444  eta_mapped = eta;
445  }
446  else
447  {
448  // Case 2
449  xi_mapped = eta;
450  eta_mapped = xi;
451  }
452 
453  else if (elem->point(3) == min_point)
454  if (elem->point(0) == std::min(elem->point(0), elem->point(2)))
455  {
456  // Case 3
457  xi_mapped = -eta;
458  eta_mapped = xi;
459  }
460  else
461  {
462  // Case 4
463  xi_mapped = xi;
464  eta_mapped = -eta;
465  }
466 
467  else if (elem->point(2) == min_point)
468  if (elem->point(3) == std::min(elem->point(3), elem->point(1)))
469  {
470  // Case 5
471  xi_mapped = -xi;
472  eta_mapped = -eta;
473  }
474  else
475  {
476  // Case 6
477  xi_mapped = -eta;
478  eta_mapped = -xi;
479  }
480 
481  else if (elem->point(1) == min_point)
482  {
483  if (elem->point(2) == std::min(elem->point(2), elem->point(0)))
484  {
485  // Case 7
486  xi_mapped = eta;
487  eta_mapped = -xi;
488  }
489  else
490  {
491  // Case 8
492  xi_mapped = -xi;
493  eta_mapped = eta;
494  }
495  }
496  }
497 
498 
499  // Face 1
500  else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
501  {
502  const Point min_point = std::min(elem->point(0),
503  std::min(elem->point(1),
504  std::min(elem->point(5),
505  elem->point(4))));
506  if (elem->point(0) == min_point)
507  if (elem->point(1) == std::min(elem->point(1), elem->point(4)))
508  {
509  // Case 1
510  xi_mapped = xi;
511  zeta_mapped = zeta;
512  }
513  else
514  {
515  // Case 2
516  xi_mapped = zeta;
517  zeta_mapped = xi;
518  }
519 
520  else if (elem->point(1) == min_point)
521  if (elem->point(5) == std::min(elem->point(5), elem->point(0)))
522  {
523  // Case 3
524  xi_mapped = zeta;
525  zeta_mapped = -xi;
526  }
527  else
528  {
529  // Case 4
530  xi_mapped = -xi;
531  zeta_mapped = zeta;
532  }
533 
534  else if (elem->point(5) == min_point)
535  if (elem->point(4) == std::min(elem->point(4), elem->point(1)))
536  {
537  // Case 5
538  xi_mapped = -xi;
539  zeta_mapped = -zeta;
540  }
541  else
542  {
543  // Case 6
544  xi_mapped = -zeta;
545  zeta_mapped = -xi;
546  }
547 
548  else if (elem->point(4) == min_point)
549  {
550  if (elem->point(0) == std::min(elem->point(0), elem->point(5)))
551  {
552  // Case 7
553  xi_mapped = -xi;
554  zeta_mapped = zeta;
555  }
556  else
557  {
558  // Case 8
559  xi_mapped = xi;
560  zeta_mapped = -zeta;
561  }
562  }
563  }
564 
565 
566  // Face 2
567  else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
568  {
569  const Point min_point = std::min(elem->point(1),
570  std::min(elem->point(2),
571  std::min(elem->point(6),
572  elem->point(5))));
573  if (elem->point(1) == min_point)
574  if (elem->point(2) == std::min(elem->point(2), elem->point(5)))
575  {
576  // Case 1
577  eta_mapped = eta;
578  zeta_mapped = zeta;
579  }
580  else
581  {
582  // Case 2
583  eta_mapped = zeta;
584  zeta_mapped = eta;
585  }
586 
587  else if (elem->point(2) == min_point)
588  if (elem->point(6) == std::min(elem->point(6), elem->point(1)))
589  {
590  // Case 3
591  eta_mapped = zeta;
592  zeta_mapped = -eta;
593  }
594  else
595  {
596  // Case 4
597  eta_mapped = -eta;
598  zeta_mapped = zeta;
599  }
600 
601  else if (elem->point(6) == min_point)
602  if (elem->point(5) == std::min(elem->point(5), elem->point(2)))
603  {
604  // Case 5
605  eta_mapped = -eta;
606  zeta_mapped = -zeta;
607  }
608  else
609  {
610  // Case 6
611  eta_mapped = -zeta;
612  zeta_mapped = -eta;
613  }
614 
615  else if (elem->point(5) == min_point)
616  {
617  if (elem->point(1) == std::min(elem->point(1), elem->point(6)))
618  {
619  // Case 7
620  eta_mapped = -zeta;
621  zeta_mapped = eta;
622  }
623  else
624  {
625  // Case 8
626  eta_mapped = eta;
627  zeta_mapped = -zeta;
628  }
629  }
630  }
631 
632 
633  // Face 3
634  else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
635  {
636  const Point min_point = std::min(elem->point(2),
637  std::min(elem->point(3),
638  std::min(elem->point(7),
639  elem->point(6))));
640  if (elem->point(3) == min_point)
641  if (elem->point(2) == std::min(elem->point(2), elem->point(7)))
642  {
643  // Case 1
644  xi_mapped = xi;
645  zeta_mapped = zeta;
646  }
647  else
648  {
649  // Case 2
650  xi_mapped = zeta;
651  zeta_mapped = xi;
652  }
653 
654  else if (elem->point(7) == min_point)
655  if (elem->point(3) == std::min(elem->point(3), elem->point(6)))
656  {
657  // Case 3
658  xi_mapped = -zeta;
659  zeta_mapped = xi;
660  }
661  else
662  {
663  // Case 4
664  xi_mapped = xi;
665  zeta_mapped = -zeta;
666  }
667 
668  else if (elem->point(6) == min_point)
669  if (elem->point(7) == std::min(elem->point(7), elem->point(2)))
670  {
671  // Case 5
672  xi_mapped = -xi;
673  zeta_mapped = -zeta;
674  }
675  else
676  {
677  // Case 6
678  xi_mapped = -zeta;
679  zeta_mapped = -xi;
680  }
681 
682  else if (elem->point(2) == min_point)
683  {
684  if (elem->point(6) == std::min(elem->point(3), elem->point(6)))
685  {
686  // Case 7
687  xi_mapped = zeta;
688  zeta_mapped = -xi;
689  }
690  else
691  {
692  // Case 8
693  xi_mapped = -xi;
694  zeta_mapped = zeta;
695  }
696  }
697  }
698 
699 
700  // Face 4
701  else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
702  {
703  const Point min_point = std::min(elem->point(3),
704  std::min(elem->point(0),
705  std::min(elem->point(4),
706  elem->point(7))));
707  if (elem->point(0) == min_point)
708  if (elem->point(3) == std::min(elem->point(3), elem->point(4)))
709  {
710  // Case 1
711  eta_mapped = eta;
712  zeta_mapped = zeta;
713  }
714  else
715  {
716  // Case 2
717  eta_mapped = zeta;
718  zeta_mapped = eta;
719  }
720 
721  else if (elem->point(4) == min_point)
722  if (elem->point(0) == std::min(elem->point(0), elem->point(7)))
723  {
724  // Case 3
725  eta_mapped = -zeta;
726  zeta_mapped = eta;
727  }
728  else
729  {
730  // Case 4
731  eta_mapped = eta;
732  zeta_mapped = -zeta;
733  }
734 
735  else if (elem->point(7) == min_point)
736  if (elem->point(4) == std::min(elem->point(4), elem->point(3)))
737  {
738  // Case 5
739  eta_mapped = -eta;
740  zeta_mapped = -zeta;
741  }
742  else
743  {
744  // Case 6
745  eta_mapped = -zeta;
746  zeta_mapped = -eta;
747  }
748 
749  else if (elem->point(3) == min_point)
750  {
751  if (elem->point(7) == std::min(elem->point(7), elem->point(0)))
752  {
753  // Case 7
754  eta_mapped = zeta;
755  zeta_mapped = -eta;
756  }
757  else
758  {
759  // Case 8
760  eta_mapped = -eta;
761  zeta_mapped = zeta;
762  }
763  }
764  }
765 
766 
767  // Face 5
768  else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
769  {
770  const Point min_point = std::min(elem->point(4),
771  std::min(elem->point(5),
772  std::min(elem->point(6),
773  elem->point(7))));
774  if (elem->point(4) == min_point)
775  if (elem->point(5) == std::min(elem->point(5), elem->point(7)))
776  {
777  // Case 1
778  xi_mapped = xi;
779  eta_mapped = eta;
780  }
781  else
782  {
783  // Case 2
784  xi_mapped = eta;
785  eta_mapped = xi;
786  }
787 
788  else if (elem->point(5) == min_point)
789  if (elem->point(6) == std::min(elem->point(6), elem->point(4)))
790  {
791  // Case 3
792  xi_mapped = eta;
793  eta_mapped = -xi;
794  }
795  else
796  {
797  // Case 4
798  xi_mapped = -xi;
799  eta_mapped = eta;
800  }
801 
802  else if (elem->point(6) == min_point)
803  if (elem->point(7) == std::min(elem->point(7), elem->point(5)))
804  {
805  // Case 5
806  xi_mapped = -xi;
807  eta_mapped = -eta;
808  }
809  else
810  {
811  // Case 6
812  xi_mapped = -eta;
813  eta_mapped = -xi;
814  }
815 
816  else if (elem->point(7) == min_point)
817  {
818  if (elem->point(4) == std::min(elem->point(4), elem->point(6)))
819  {
820  // Case 7
821  xi_mapped = -eta;
822  eta_mapped = xi;
823  }
824  else
825  {
826  // Case 8
827  xi_mapped = xi;
828  eta_mapped = eta;
829  }
830  }
831  }
832  }
833 
834  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi_mapped)*
835  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta_mapped)*
836  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta_mapped));
837  }
838 
839 
840  default:
841  libmesh_error_msg("Invalid element type = " << type);
842  } //case HEX27
843 
844  }//case THIRD
845 
846 
847  // 4th-order Bernstein.
848  case FOURTH:
849  {
850  switch (type)
851  {
852 
853  // Bernstein shape functions on the hexahedral.
854  case HEX27:
855  {
856  libmesh_assert_less (i, 125);
857 
858  // Compute hex shape functions as a tensor-product
859  const Real xi = p(0);
860  const Real eta = p(1);
861  const Real zeta = p(2);
862  Real xi_mapped = p(0);
863  Real eta_mapped = p(1);
864  Real zeta_mapped = p(2);
865 
866  // The only way to make any sense of this
867  // is to look at the mgflo/mg2/mgf documentation
868  // and make the cut-out cube!
869  // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
870  // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
871  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4};
872  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
873  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4};
874 
875 
876 
877  // handle the edge orientation
878  {
879  // Edge 0
880  if ((i0[i] >= 2) && (i1[i] == 0) && (i2[i] == 0))
881  {
882  if (elem->point(0) != std::min(elem->point(0), elem->point(1)))
883  xi_mapped = -xi;
884  }
885  // Edge 1
886  else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] == 0))
887  {
888  if (elem->point(1) != std::min(elem->point(1), elem->point(2)))
889  eta_mapped = -eta;
890  }
891  // Edge 2
892  else if ((i0[i] >= 2) && (i1[i] == 1) && (i2[i] == 0))
893  {
894  if (elem->point(3) != std::min(elem->point(3), elem->point(2)))
895  xi_mapped = -xi;
896  }
897  // Edge 3
898  else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] == 0))
899  {
900  if (elem->point(0) != std::min(elem->point(0), elem->point(3)))
901  eta_mapped = -eta;
902  }
903  // Edge 4
904  else if ((i0[i] == 0) && (i1[i] == 0) && (i2[i] >=2 ))
905  {
906  if (elem->point(0) != std::min(elem->point(0), elem->point(4)))
907  zeta_mapped = -zeta;
908  }
909  // Edge 5
910  else if ((i0[i] == 1) && (i1[i] == 0) && (i2[i] >=2 ))
911  {
912  if (elem->point(1) != std::min(elem->point(1), elem->point(5)))
913  zeta_mapped = -zeta;
914  }
915  // Edge 6
916  else if ((i0[i] == 1) && (i1[i] == 1) && (i2[i] >=2 ))
917  {
918  if (elem->point(2) != std::min(elem->point(2), elem->point(6)))
919  zeta_mapped = -zeta;
920  }
921  // Edge 7
922  else if ((i0[i] == 0) && (i1[i] == 1) && (i2[i] >=2 ))
923  {
924  if (elem->point(3) != std::min(elem->point(3), elem->point(7)))
925  zeta_mapped = -zeta;
926  }
927  // Edge 8
928  else if ((i0[i] >=2 ) && (i1[i] == 0) && (i2[i] == 1))
929  {
930  if (elem->point(4) != std::min(elem->point(4), elem->point(5)))
931  xi_mapped = -xi;
932  }
933  // Edge 9
934  else if ((i0[i] == 1) && (i1[i] >=2 ) && (i2[i] == 1))
935  {
936  if (elem->point(5) != std::min(elem->point(5), elem->point(6)))
937  eta_mapped = -eta;
938  }
939  // Edge 10
940  else if ((i0[i] >=2 ) && (i1[i] == 1) && (i2[i] == 1))
941  {
942  if (elem->point(7) != std::min(elem->point(7), elem->point(6)))
943  xi_mapped = -xi;
944  }
945  // Edge 11
946  else if ((i0[i] == 0) && (i1[i] >=2 ) && (i2[i] == 1))
947  {
948  if (elem->point(4) != std::min(elem->point(4), elem->point(7)))
949  eta_mapped = -eta;
950  }
951  }
952 
953 
954  // handle the face orientation
955  {
956  // Face 0
957  if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
958  {
959  const Point min_point = std::min(elem->point(1),
960  std::min(elem->point(2),
961  std::min(elem->point(0),
962  elem->point(3))));
963  if (elem->point(0) == min_point)
964  if (elem->point(1) == std::min(elem->point(1), elem->point(3)))
965  {
966  // Case 1
967  xi_mapped = xi;
968  eta_mapped = eta;
969  }
970  else
971  {
972  // Case 2
973  xi_mapped = eta;
974  eta_mapped = xi;
975  }
976 
977  else if (elem->point(3) == min_point)
978  if (elem->point(0) == std::min(elem->point(0), elem->point(2)))
979  {
980  // Case 3
981  xi_mapped = -eta;
982  eta_mapped = xi;
983  }
984  else
985  {
986  // Case 4
987  xi_mapped = xi;
988  eta_mapped = -eta;
989  }
990 
991  else if (elem->point(2) == min_point)
992  if (elem->point(3) == std::min(elem->point(3), elem->point(1)))
993  {
994  // Case 5
995  xi_mapped = -xi;
996  eta_mapped = -eta;
997  }
998  else
999  {
1000  // Case 6
1001  xi_mapped = -eta;
1002  eta_mapped = -xi;
1003  }
1004 
1005  else if (elem->point(1) == min_point)
1006  {
1007  if (elem->point(2) == std::min(elem->point(2), elem->point(0)))
1008  {
1009  // Case 7
1010  xi_mapped = eta;
1011  eta_mapped = -xi;
1012  }
1013  else
1014  {
1015  // Case 8
1016  xi_mapped = -xi;
1017  eta_mapped = eta;
1018  }
1019  }
1020  }
1021 
1022 
1023  // Face 1
1024  else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
1025  {
1026  const Point min_point = std::min(elem->point(0),
1027  std::min(elem->point(1),
1028  std::min(elem->point(5),
1029  elem->point(4))));
1030  if (elem->point(0) == min_point)
1031  if (elem->point(1) == std::min(elem->point(1), elem->point(4)))
1032  {
1033  // Case 1
1034  xi_mapped = xi;
1035  zeta_mapped = zeta;
1036  }
1037  else
1038  {
1039  // Case 2
1040  xi_mapped = zeta;
1041  zeta_mapped = xi;
1042  }
1043 
1044  else if (elem->point(1) == min_point)
1045  if (elem->point(5) == std::min(elem->point(5), elem->point(0)))
1046  {
1047  // Case 3
1048  xi_mapped = zeta;
1049  zeta_mapped = -xi;
1050  }
1051  else
1052  {
1053  // Case 4
1054  xi_mapped = -xi;
1055  zeta_mapped = zeta;
1056  }
1057 
1058  else if (elem->point(5) == min_point)
1059  if (elem->point(4) == std::min(elem->point(4), elem->point(1)))
1060  {
1061  // Case 5
1062  xi_mapped = -xi;
1063  zeta_mapped = -zeta;
1064  }
1065  else
1066  {
1067  // Case 6
1068  xi_mapped = -zeta;
1069  zeta_mapped = -xi;
1070  }
1071 
1072  else if (elem->point(4) == min_point)
1073  {
1074  if (elem->point(0) == std::min(elem->point(0), elem->point(5)))
1075  {
1076  // Case 7
1077  xi_mapped = -xi;
1078  zeta_mapped = zeta;
1079  }
1080  else
1081  {
1082  // Case 8
1083  xi_mapped = xi;
1084  zeta_mapped = -zeta;
1085  }
1086  }
1087  }
1088 
1089 
1090  // Face 2
1091  else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
1092  {
1093  const Point min_point = std::min(elem->point(1),
1094  std::min(elem->point(2),
1095  std::min(elem->point(6),
1096  elem->point(5))));
1097  if (elem->point(1) == min_point)
1098  if (elem->point(2) == std::min(elem->point(2), elem->point(5)))
1099  {
1100  // Case 1
1101  eta_mapped = eta;
1102  zeta_mapped = zeta;
1103  }
1104  else
1105  {
1106  // Case 2
1107  eta_mapped = zeta;
1108  zeta_mapped = eta;
1109  }
1110 
1111  else if (elem->point(2) == min_point)
1112  if (elem->point(6) == std::min(elem->point(6), elem->point(1)))
1113  {
1114  // Case 3
1115  eta_mapped = zeta;
1116  zeta_mapped = -eta;
1117  }
1118  else
1119  {
1120  // Case 4
1121  eta_mapped = -eta;
1122  zeta_mapped = zeta;
1123  }
1124 
1125  else if (elem->point(6) == min_point)
1126  if (elem->point(5) == std::min(elem->point(5), elem->point(2)))
1127  {
1128  // Case 5
1129  eta_mapped = -eta;
1130  zeta_mapped = -zeta;
1131  }
1132  else
1133  {
1134  // Case 6
1135  eta_mapped = -zeta;
1136  zeta_mapped = -eta;
1137  }
1138 
1139  else if (elem->point(5) == min_point)
1140  {
1141  if (elem->point(1) == std::min(elem->point(1), elem->point(6)))
1142  {
1143  // Case 7
1144  eta_mapped = -zeta;
1145  zeta_mapped = eta;
1146  }
1147  else
1148  {
1149  // Case 8
1150  eta_mapped = eta;
1151  zeta_mapped = -zeta;
1152  }
1153  }
1154  }
1155 
1156 
1157  // Face 3
1158  else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
1159  {
1160  const Point min_point = std::min(elem->point(2),
1161  std::min(elem->point(3),
1162  std::min(elem->point(7),
1163  elem->point(6))));
1164  if (elem->point(3) == min_point)
1165  if (elem->point(2) == std::min(elem->point(2), elem->point(7)))
1166  {
1167  // Case 1
1168  xi_mapped = xi;
1169  zeta_mapped = zeta;
1170  }
1171  else
1172  {
1173  // Case 2
1174  xi_mapped = zeta;
1175  zeta_mapped = xi;
1176  }
1177 
1178  else if (elem->point(7) == min_point)
1179  if (elem->point(3) == std::min(elem->point(3), elem->point(6)))
1180  {
1181  // Case 3
1182  xi_mapped = -zeta;
1183  zeta_mapped = xi;
1184  }
1185  else
1186  {
1187  // Case 4
1188  xi_mapped = xi;
1189  zeta_mapped = -zeta;
1190  }
1191 
1192  else if (elem->point(6) == min_point)
1193  if (elem->point(7) == std::min(elem->point(7), elem->point(2)))
1194  {
1195  // Case 5
1196  xi_mapped = -xi;
1197  zeta_mapped = -zeta;
1198  }
1199  else
1200  {
1201  // Case 6
1202  xi_mapped = -zeta;
1203  zeta_mapped = -xi;
1204  }
1205 
1206  else if (elem->point(2) == min_point)
1207  {
1208  if (elem->point(6) == std::min(elem->point(3), elem->point(6)))
1209  {
1210  // Case 7
1211  xi_mapped = zeta;
1212  zeta_mapped = -xi;
1213  }
1214  else
1215  {
1216  // Case 8
1217  xi_mapped = -xi;
1218  zeta_mapped = zeta;
1219  }
1220  }
1221  }
1222 
1223 
1224  // Face 4
1225  else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
1226  {
1227  const Point min_point = std::min(elem->point(3),
1228  std::min(elem->point(0),
1229  std::min(elem->point(4),
1230  elem->point(7))));
1231  if (elem->point(0) == min_point)
1232  if (elem->point(3) == std::min(elem->point(3), elem->point(4)))
1233  {
1234  // Case 1
1235  eta_mapped = eta;
1236  zeta_mapped = zeta;
1237  }
1238  else
1239  {
1240  // Case 2
1241  eta_mapped = zeta;
1242  zeta_mapped = eta;
1243  }
1244 
1245  else if (elem->point(4) == min_point)
1246  if (elem->point(0) == std::min(elem->point(0), elem->point(7)))
1247  {
1248  // Case 3
1249  eta_mapped = -zeta;
1250  zeta_mapped = eta;
1251  }
1252  else
1253  {
1254  // Case 4
1255  eta_mapped = eta;
1256  zeta_mapped = -zeta;
1257  }
1258 
1259  else if (elem->point(7) == min_point)
1260  if (elem->point(4) == std::min(elem->point(4), elem->point(3)))
1261  {
1262  // Case 5
1263  eta_mapped = -eta;
1264  zeta_mapped = -zeta;
1265  }
1266  else
1267  {
1268  // Case 6
1269  eta_mapped = -zeta;
1270  zeta_mapped = -eta;
1271  }
1272 
1273  else if (elem->point(3) == min_point)
1274  {
1275  if (elem->point(7) == std::min(elem->point(7), elem->point(0)))
1276  {
1277  // Case 7
1278  eta_mapped = zeta;
1279  zeta_mapped = -eta;
1280  }
1281  else
1282  {
1283  // Case 8
1284  eta_mapped = -eta;
1285  zeta_mapped = zeta;
1286  }
1287  }
1288  }
1289 
1290 
1291  // Face 5
1292  else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
1293  {
1294  const Point min_point = std::min(elem->point(4),
1295  std::min(elem->point(5),
1296  std::min(elem->point(6),
1297  elem->point(7))));
1298  if (elem->point(4) == min_point)
1299  if (elem->point(5) == std::min(elem->point(5), elem->point(7)))
1300  {
1301  // Case 1
1302  xi_mapped = xi;
1303  eta_mapped = eta;
1304  }
1305  else
1306  {
1307  // Case 2
1308  xi_mapped = eta;
1309  eta_mapped = xi;
1310  }
1311 
1312  else if (elem->point(5) == min_point)
1313  if (elem->point(6) == std::min(elem->point(6), elem->point(4)))
1314  {
1315  // Case 3
1316  xi_mapped = eta;
1317  eta_mapped = -xi;
1318  }
1319  else
1320  {
1321  // Case 4
1322  xi_mapped = -xi;
1323  eta_mapped = eta;
1324  }
1325 
1326  else if (elem->point(6) == min_point)
1327  if (elem->point(7) == std::min(elem->point(7), elem->point(5)))
1328  {
1329  // Case 5
1330  xi_mapped = -xi;
1331  eta_mapped = -eta;
1332  }
1333  else
1334  {
1335  // Case 6
1336  xi_mapped = -eta;
1337  eta_mapped = -xi;
1338  }
1339 
1340  else if (elem->point(7) == min_point)
1341  {
1342  if (elem->point(4) == std::min(elem->point(4), elem->point(6)))
1343  {
1344  // Case 7
1345  xi_mapped = -eta;
1346  eta_mapped = xi;
1347  }
1348  else
1349  {
1350  // Case 8
1351  xi_mapped = xi;
1352  eta_mapped = eta;
1353  }
1354  }
1355  }
1356 
1357 
1358  }
1359 
1360 
1361  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi_mapped)*
1362  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta_mapped)*
1363  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta_mapped));
1364  }
1365 
1366 
1367  default:
1368  libmesh_error_msg("Invalid element type = " << type);
1369  }
1370  }
1371 
1372 
1373  default:
1374  libmesh_error_msg("Invalid totalorder = " << totalorder);
1375  }
1376 #else // LIBMESH_DIM != 3
1377  libmesh_ignore(elem, order, i, p, add_p_level);
1378  libmesh_not_implemented();
1379 #endif
1380 }
1381 
1382 
1383 
1384 
1385 template <>
1386 Real FE<3,BERNSTEIN>::shape_deriv(const ElemType,
1387  const Order,
1388  const unsigned int,
1389  const unsigned int,
1390  const Point & )
1391 {
1392  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
1393  return 0.;
1394 }
1395 
1396 
1397 
1398 template <>
1399 Real FE<3,BERNSTEIN>::shape_deriv(const Elem * elem,
1400  const Order order,
1401  const unsigned int i,
1402  const unsigned int j,
1403  const Point & p,
1404  const bool add_p_level)
1405 {
1406 
1407 #if LIBMESH_DIM == 3
1408  libmesh_assert(elem);
1409  const ElemType type = elem->type();
1410 
1411  const Order totalorder =
1412  static_cast<Order>(order + add_p_level * elem->p_level());
1413 
1414  libmesh_assert_less (j, 3);
1415 
1416  switch (totalorder)
1417  {
1418  // 1st order Bernstein.
1419  case FIRST:
1420  {
1421  switch (type)
1422  {
1423  // Bernstein shape functions on the tetrahedron.
1424  case TET4:
1425  case TET10:
1426  {
1427  // I have been lazy here and am using finite differences
1428  // to compute the derivatives!
1429  const Real eps = 1.e-6;
1430 
1431  libmesh_assert_less (i, 4);
1432  libmesh_assert_less (j, 3);
1433 
1434 
1435  switch (j)
1436  {
1437  // d()/dxi
1438  case 0:
1439  {
1440  const Point pp(p(0)+eps, p(1), p(2));
1441  const Point pm(p(0)-eps, p(1), p(2));
1442 
1443  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1444  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1445  }
1446 
1447  // d()/deta
1448  case 1:
1449  {
1450  const Point pp(p(0), p(1)+eps, p(2));
1451  const Point pm(p(0), p(1)-eps, p(2));
1452 
1453  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1454  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1455  }
1456  // d()/dzeta
1457  case 2:
1458  {
1459  const Point pp(p(0), p(1), p(2)+eps);
1460  const Point pm(p(0), p(1), p(2)-eps);
1461 
1462  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1463  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1464  }
1465  default:
1466  libmesh_error_msg("Invalid derivative index j = " << j);
1467  }
1468  }
1469 
1470 
1471 
1472 
1473  // Bernstein shape functions on the hexahedral.
1474  case HEX8:
1475  case HEX20:
1476  case HEX27:
1477  {
1478  libmesh_assert_less (i, 8);
1479 
1480  // Compute hex shape functions as a tensor-product
1481  const Real xi = p(0);
1482  const Real eta = p(1);
1483  const Real zeta = p(2);
1484 
1485  // The only way to make any sense of this
1486  // is to look at the mgflo/mg2/mgf documentation
1487  // and make the cut-out cube!
1488  // 0 1 2 3 4 5 6 7
1489  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
1490  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
1491  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
1492 
1493  switch (j)
1494  {
1495  // d()/dxi
1496  case 0:
1497  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1498  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1499  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1500 
1501  // d()/deta
1502  case 1:
1503  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1504  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1505  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1506 
1507  // d()/dzeta
1508  case 2:
1509  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1510  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1511  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1512 
1513  default:
1514  libmesh_error_msg("Invalid derivative index j = " << j);
1515  }
1516  }
1517 
1518  default:
1519  libmesh_error_msg("Invalid element type = " << type);
1520  }
1521  }
1522 
1523 
1524 
1525 
1526  case SECOND:
1527  {
1528  switch (type)
1529  {
1530  // Bernstein shape functions on the tetrahedron.
1531  case TET10:
1532  {
1533  // I have been lazy here and am using finite differences
1534  // to compute the derivatives!
1535  const Real eps = 1.e-6;
1536 
1537  libmesh_assert_less (i, 10);
1538  libmesh_assert_less (j, 3);
1539 
1540 
1541  switch (j)
1542  {
1543  // d()/dxi
1544  case 0:
1545  {
1546  const Point pp(p(0)+eps, p(1), p(2));
1547  const Point pm(p(0)-eps, p(1), p(2));
1548 
1549  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1550  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1551  }
1552 
1553  // d()/deta
1554  case 1:
1555  {
1556  const Point pp(p(0), p(1)+eps, p(2));
1557  const Point pm(p(0), p(1)-eps, p(2));
1558 
1559  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1560  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1561  }
1562  // d()/dzeta
1563  case 2:
1564  {
1565  const Point pp(p(0), p(1), p(2)+eps);
1566  const Point pm(p(0), p(1), p(2)-eps);
1567 
1568  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1569  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1570  }
1571  default:
1572  libmesh_error_msg("Invalid derivative index j = " << j);
1573  }
1574  }
1575 
1576  // Bernstein shape functions on the hexahedral.
1577  case HEX20:
1578  {
1579  libmesh_assert_less (i, 20);
1580 
1581  // Compute hex shape functions as a tensor-product
1582  const Real xi = p(0);
1583  const Real eta = p(1);
1584  const Real zeta = p(2);
1585 
1586  // The only way to make any sense of this
1587  // is to look at the mgflo/mg2/mgf documentation
1588  // and make the cut-out cube!
1589  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
1590  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
1591  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
1592  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
1593  static const Real scal20[] = {-0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0};
1594  static const Real scal21[] = {-0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0};
1595  static const Real scal22[] = {0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0};
1596  static const Real scal23[] = {0, 0, -0.25, -0.25, 0, 0, -0.25, -0.25, 0, 0, 0.5, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0};
1597  static const Real scal24[] = {-0.25, 0, 0, -0.25, -0.25, 0, 0, -0.25, 0, 0, 0, 0.5, 0.5, 0, 0, 0.5, 0, 0, 0, 0.5};
1598  static const Real scal25[] = {0, 0, 0, 0, -0.25, -0.25, -0.25, -0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5};
1599  static const Real scal26[] = {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
1600 
1601  switch (j)
1602  {
1603  // d()/dxi
1604  case 0:
1605  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1606  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1607  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)
1608  +scal20[i]*
1609  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[20], 0, xi)*
1610  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[20], eta)*
1611  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[20], zeta)
1612  +scal21[i]*
1613  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[21], 0, xi)*
1614  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[21], eta)*
1615  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[21], zeta)
1616  +scal22[i]*
1617  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[22], 0, xi)*
1618  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[22], eta)*
1619  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[22], zeta)
1620  +scal23[i]*
1621  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[23], 0, xi)*
1622  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[23], eta)*
1623  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[23], zeta)
1624  +scal24[i]*
1625  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[24], 0, xi)*
1626  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[24], eta)*
1627  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[24], zeta)
1628  +scal25[i]*
1629  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[25], 0, xi)*
1630  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[25], eta)*
1631  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[25], zeta)
1632  +scal26[i]*
1633  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[26], 0, xi)*
1634  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[26], eta)*
1635  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[26], zeta));
1636 
1637  // d()/deta
1638  case 1:
1639  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1640  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1641  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)
1642  +scal20[i]*
1643  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[20], xi)*
1644  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[20], 0, eta)*
1645  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[20], zeta)
1646  +scal21[i]*
1647  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[21], xi)*
1648  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[21], 0, eta)*
1649  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[21], zeta)
1650  +scal22[i]*
1651  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[22], xi)*
1652  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[22], 0, eta)*
1653  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[22], zeta)
1654  +scal23[i]*
1655  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[23], xi)*
1656  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[23], 0, eta)*
1657  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[23], zeta)
1658  +scal24[i]*
1659  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[24], xi)*
1660  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[24], 0, eta)*
1661  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[24], zeta)
1662  +scal25[i]*
1663  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[25], xi)*
1664  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[25], 0, eta)*
1665  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[25], zeta)
1666  +scal26[i]*
1667  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[26], xi)*
1668  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[26], 0, eta)*
1669  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[26], zeta));
1670 
1671  // d()/dzeta
1672  case 2:
1673  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1674  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1675  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta)
1676  +scal20[i]*
1677  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[20], xi)*
1678  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[20], eta)*
1679  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[20], 0, zeta)
1680  +scal21[i]*
1681  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[21], xi)*
1682  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[21], eta)*
1683  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[21], 0, zeta)
1684  +scal22[i]*
1685  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[22], xi)*
1686  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[22], eta)*
1687  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[22], 0, zeta)
1688  +scal23[i]*
1689  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[23], xi)*
1690  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[23], eta)*
1691  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[23], 0, zeta)
1692  +scal24[i]*
1693  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[24], xi)*
1694  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[24], eta)*
1695  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[24], 0, zeta)
1696  +scal25[i]*
1697  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[25], xi)*
1698  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[25], eta)*
1699  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[25], 0, zeta)
1700  +scal26[i]*
1701  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[26], xi)*
1702  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[26], eta)*
1703  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[26], 0, zeta));
1704 
1705  default:
1706  libmesh_error_msg("Invalid derivative index j = " << j);
1707  }
1708  }
1709 
1710  // Bernstein shape functions on the hexahedral.
1711  case HEX27:
1712  {
1713  libmesh_assert_less (i, 27);
1714 
1715  // Compute hex shape functions as a tensor-product
1716  const Real xi = p(0);
1717  const Real eta = p(1);
1718  const Real zeta = p(2);
1719 
1720  // The only way to make any sense of this
1721  // is to look at the mgflo/mg2/mgf documentation
1722  // and make the cut-out cube!
1723  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
1724  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
1725  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
1726  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
1727 
1728  switch (j)
1729  {
1730  // d()/dxi
1731  case 0:
1732  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1733  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1734  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1735 
1736  // d()/deta
1737  case 1:
1738  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1739  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1740  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1741 
1742  // d()/dzeta
1743  case 2:
1744  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1745  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1746  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1747 
1748  default:
1749  libmesh_error_msg("Invalid derivative index j = " << j);
1750  }
1751  }
1752 
1753 
1754  default:
1755  libmesh_error_msg("Invalid element type = " << type);
1756  }
1757  }
1758 
1759 
1760 
1761  // 3rd-order Bernstein.
1762  case THIRD:
1763  {
1764  switch (type)
1765  {
1766 
1767  // // Bernstein shape functions derivatives.
1768  // case TET10:
1769  // {
1770  // // I have been lazy here and am using finite differences
1771  // // to compute the derivatives!
1772  // const Real eps = 1.e-6;
1773 
1774  // libmesh_assert_less (i, 20);
1775  // libmesh_assert_less (j, 3);
1776 
1777  // switch (j)
1778  // {
1779  // // d()/dxi
1780  // case 0:
1781  // {
1782  // const Point pp(p(0)+eps, p(1), p(2));
1783  // const Point pm(p(0)-eps, p(1), p(2));
1784 
1785  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1786  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1787  // }
1788 
1789  // // d()/deta
1790  // case 1:
1791  // {
1792  // const Point pp(p(0), p(1)+eps, p(2));
1793  // const Point pm(p(0), p(1)-eps, p(2));
1794 
1795  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1796  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1797  // }
1798  // // d()/dzeta
1799  // case 2:
1800  // {
1801  // const Point pp(p(0), p(1), p(2)+eps);
1802  // const Point pm(p(0), p(1), p(2)-eps);
1803 
1804  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1805  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1806  // }
1807  // default:
1808  // libmesh_error_msg("Invalid derivative index j = " << j);
1809  // }
1810 
1811 
1812  // }
1813 
1814 
1815  // Bernstein shape functions on the hexahedral.
1816  case HEX27:
1817  {
1818  // I have been lazy here and am using finite differences
1819  // to compute the derivatives!
1820  const Real eps = 1.e-6;
1821 
1822  libmesh_assert_less (i, 64);
1823  libmesh_assert_less (j, 3);
1824 
1825  switch (j)
1826  {
1827  // d()/dxi
1828  case 0:
1829  {
1830  const Point pp(p(0)+eps, p(1), p(2));
1831  const Point pm(p(0)-eps, p(1), p(2));
1832 
1833  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1834  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1835  }
1836 
1837  // d()/deta
1838  case 1:
1839  {
1840  const Point pp(p(0), p(1)+eps, p(2));
1841  const Point pm(p(0), p(1)-eps, p(2));
1842 
1843  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1844  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1845  }
1846  // d()/dzeta
1847  case 2:
1848  {
1849  const Point pp(p(0), p(1), p(2)+eps);
1850  const Point pm(p(0), p(1), p(2)-eps);
1851 
1852  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1853  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1854  }
1855  default:
1856  libmesh_error_msg("Invalid derivative index j = " << j);
1857  }
1858 
1859  }
1860 
1861  // // Compute hex shape functions as a tensor-product
1862  // const Real xi = p(0);
1863  // const Real eta = p(1);
1864  // const Real zeta = p(2);
1865  // Real xi_mapped = p(0);
1866  // Real eta_mapped = p(1);
1867  // Real zeta_mapped = p(2);
1868 
1869  // // The only way to make any sense of this
1870  // // is to look at the mgflo/mg2/mgf documentation
1871  // // and make the cut-out cube!
1872  // // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
1873  // // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63
1874  // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
1875  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
1876  // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
1877 
1878 
1879 
1880  // // handle the edge orientation
1881  // {
1882  // // Edge 0
1883  // if ((i1[i] == 0) && (i2[i] == 0))
1884  // {
1885  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
1886  // xi_mapped = -xi;
1887  // }
1888  // // Edge 1
1889  // else if ((i0[i] == 1) && (i2[i] == 0))
1890  // {
1891  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
1892  // eta_mapped = -eta;
1893  // }
1894  // // Edge 2
1895  // else if ((i1[i] == 1) && (i2[i] == 0))
1896  // {
1897  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
1898  // xi_mapped = -xi;
1899  // }
1900  // // Edge 3
1901  // else if ((i0[i] == 0) && (i2[i] == 0))
1902  // {
1903  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
1904  // eta_mapped = -eta;
1905  // }
1906  // // Edge 4
1907  // else if ((i0[i] == 0) && (i1[i] == 0))
1908  // {
1909  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
1910  // zeta_mapped = -zeta;
1911  // }
1912  // // Edge 5
1913  // else if ((i0[i] == 1) && (i1[i] == 0))
1914  // {
1915  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
1916  // zeta_mapped = -zeta;
1917  // }
1918  // // Edge 6
1919  // else if ((i0[i] == 1) && (i1[i] == 1))
1920  // {
1921  // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
1922  // zeta_mapped = -zeta;
1923  // }
1924  // // Edge 7
1925  // else if ((i0[i] == 0) && (i1[i] == 1))
1926  // {
1927  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
1928  // zeta_mapped = -zeta;
1929  // }
1930  // // Edge 8
1931  // else if ((i1[i] == 0) && (i2[i] == 1))
1932  // {
1933  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
1934  // xi_mapped = -xi;
1935  // }
1936  // // Edge 9
1937  // else if ((i0[i] == 1) && (i2[i] == 1))
1938  // {
1939  // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
1940  // eta_mapped = -eta;
1941  // }
1942  // // Edge 10
1943  // else if ((i1[i] == 1) && (i2[i] == 1))
1944  // {
1945  // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
1946  // xi_mapped = -xi;
1947  // }
1948  // // Edge 11
1949  // else if ((i0[i] == 0) && (i2[i] == 1))
1950  // {
1951  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
1952  // eta_mapped = -eta;
1953  // }
1954  // }
1955 
1956 
1957  // // handle the face orientation
1958  // {
1959  // // Face 0
1960  // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
1961  // {
1962  // const unsigned int min_node = std::min(elem->node_id(1),
1963  // std::min(elem->node_id(2),
1964  // std::min(elem->node_id(0),
1965  // elem->node_id(3))));
1966  // if (elem->node_id(0) == min_node)
1967  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
1968  // {
1969  // // Case 1
1970  // xi_mapped = xi;
1971  // eta_mapped = eta;
1972  // }
1973  // else
1974  // {
1975  // // Case 2
1976  // xi_mapped = eta;
1977  // eta_mapped = xi;
1978  // }
1979 
1980  // else if (elem->node_id(3) == min_node)
1981  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
1982  // {
1983  // // Case 3
1984  // xi_mapped = -eta;
1985  // eta_mapped = xi;
1986  // }
1987  // else
1988  // {
1989  // // Case 4
1990  // xi_mapped = xi;
1991  // eta_mapped = -eta;
1992  // }
1993 
1994  // else if (elem->node_id(2) == min_node)
1995  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
1996  // {
1997  // // Case 5
1998  // xi_mapped = -xi;
1999  // eta_mapped = -eta;
2000  // }
2001  // else
2002  // {
2003  // // Case 6
2004  // xi_mapped = -eta;
2005  // eta_mapped = -xi;
2006  // }
2007 
2008  // else if (elem->node_id(1) == min_node)
2009  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
2010  // {
2011  // // Case 7
2012  // xi_mapped = eta;
2013  // eta_mapped = -xi;
2014  // }
2015  // else
2016  // {
2017  // // Case 8
2018  // xi_mapped = -xi;
2019  // eta_mapped = eta;
2020  // }
2021  // }
2022 
2023 
2024  // // Face 1
2025  // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
2026  // {
2027  // const unsigned int min_node = std::min(elem->node_id(0),
2028  // std::min(elem->node_id(1),
2029  // std::min(elem->node_id(5),
2030  // elem->node_id(4))));
2031  // if (elem->node_id(0) == min_node)
2032  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
2033  // {
2034  // // Case 1
2035  // xi_mapped = xi;
2036  // zeta_mapped = zeta;
2037  // }
2038  // else
2039  // {
2040  // // Case 2
2041  // xi_mapped = zeta;
2042  // zeta_mapped = xi;
2043  // }
2044 
2045  // else if (elem->node_id(1) == min_node)
2046  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
2047  // {
2048  // // Case 3
2049  // xi_mapped = zeta;
2050  // zeta_mapped = -xi;
2051  // }
2052  // else
2053  // {
2054  // // Case 4
2055  // xi_mapped = -xi;
2056  // zeta_mapped = zeta;
2057  // }
2058 
2059  // else if (elem->node_id(5) == min_node)
2060  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
2061  // {
2062  // // Case 5
2063  // xi_mapped = -xi;
2064  // zeta_mapped = -zeta;
2065  // }
2066  // else
2067  // {
2068  // // Case 6
2069  // xi_mapped = -zeta;
2070  // zeta_mapped = -xi;
2071  // }
2072 
2073  // else if (elem->node_id(4) == min_node)
2074  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
2075  // {
2076  // // Case 7
2077  // xi_mapped = -xi;
2078  // zeta_mapped = zeta;
2079  // }
2080  // else
2081  // {
2082  // // Case 8
2083  // xi_mapped = xi;
2084  // zeta_mapped = -zeta;
2085  // }
2086  // }
2087 
2088 
2089  // // Face 2
2090  // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
2091  // {
2092  // const unsigned int min_node = std::min(elem->node_id(1),
2093  // std::min(elem->node_id(2),
2094  // std::min(elem->node_id(6),
2095  // elem->node_id(5))));
2096  // if (elem->node_id(1) == min_node)
2097  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
2098  // {
2099  // // Case 1
2100  // eta_mapped = eta;
2101  // zeta_mapped = zeta;
2102  // }
2103  // else
2104  // {
2105  // // Case 2
2106  // eta_mapped = zeta;
2107  // zeta_mapped = eta;
2108  // }
2109 
2110  // else if (elem->node_id(2) == min_node)
2111  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
2112  // {
2113  // // Case 3
2114  // eta_mapped = zeta;
2115  // zeta_mapped = -eta;
2116  // }
2117  // else
2118  // {
2119  // // Case 4
2120  // eta_mapped = -eta;
2121  // zeta_mapped = zeta;
2122  // }
2123 
2124  // else if (elem->node_id(6) == min_node)
2125  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
2126  // {
2127  // // Case 5
2128  // eta_mapped = -eta;
2129  // zeta_mapped = -zeta;
2130  // }
2131  // else
2132  // {
2133  // // Case 6
2134  // eta_mapped = -zeta;
2135  // zeta_mapped = -eta;
2136  // }
2137 
2138  // else if (elem->node_id(5) == min_node)
2139  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
2140  // {
2141  // // Case 7
2142  // eta_mapped = -zeta;
2143  // zeta_mapped = eta;
2144  // }
2145  // else
2146  // {
2147  // // Case 8
2148  // eta_mapped = eta;
2149  // zeta_mapped = -zeta;
2150  // }
2151  // }
2152 
2153 
2154  // // Face 3
2155  // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
2156  // {
2157  // const unsigned int min_node = std::min(elem->node_id(2),
2158  // std::min(elem->node_id(3),
2159  // std::min(elem->node_id(7),
2160  // elem->node_id(6))));
2161  // if (elem->node_id(3) == min_node)
2162  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
2163  // {
2164  // // Case 1
2165  // xi_mapped = xi;
2166  // zeta_mapped = zeta;
2167  // }
2168  // else
2169  // {
2170  // // Case 2
2171  // xi_mapped = zeta;
2172  // zeta_mapped = xi;
2173  // }
2174 
2175  // else if (elem->node_id(7) == min_node)
2176  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
2177  // {
2178  // // Case 3
2179  // xi_mapped = -zeta;
2180  // zeta_mapped = xi;
2181  // }
2182  // else
2183  // {
2184  // // Case 4
2185  // xi_mapped = xi;
2186  // zeta_mapped = -zeta;
2187  // }
2188 
2189  // else if (elem->node_id(6) == min_node)
2190  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
2191  // {
2192  // // Case 5
2193  // xi_mapped = -xi;
2194  // zeta_mapped = -zeta;
2195  // }
2196  // else
2197  // {
2198  // // Case 6
2199  // xi_mapped = -zeta;
2200  // zeta_mapped = -xi;
2201  // }
2202 
2203  // else if (elem->node_id(2) == min_node)
2204  // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
2205  // {
2206  // // Case 7
2207  // xi_mapped = zeta;
2208  // zeta_mapped = -xi;
2209  // }
2210  // else
2211  // {
2212  // // Case 8
2213  // xi_mapped = -xi;
2214  // zeta_mapped = zeta;
2215  // }
2216  // }
2217 
2218 
2219  // // Face 4
2220  // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
2221  // {
2222  // const unsigned int min_node = std::min(elem->node_id(3),
2223  // std::min(elem->node_id(0),
2224  // std::min(elem->node_id(4),
2225  // elem->node_id(7))));
2226  // if (elem->node_id(0) == min_node)
2227  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
2228  // {
2229  // // Case 1
2230  // eta_mapped = eta;
2231  // zeta_mapped = zeta;
2232  // }
2233  // else
2234  // {
2235  // // Case 2
2236  // eta_mapped = zeta;
2237  // zeta_mapped = eta;
2238  // }
2239 
2240  // else if (elem->node_id(4) == min_node)
2241  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
2242  // {
2243  // // Case 3
2244  // eta_mapped = -zeta;
2245  // zeta_mapped = eta;
2246  // }
2247  // else
2248  // {
2249  // // Case 4
2250  // eta_mapped = eta;
2251  // zeta_mapped = -zeta;
2252  // }
2253 
2254  // else if (elem->node_id(7) == min_node)
2255  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
2256  // {
2257  // // Case 5
2258  // eta_mapped = -eta;
2259  // zeta_mapped = -zeta;
2260  // }
2261  // else
2262  // {
2263  // // Case 6
2264  // eta_mapped = -zeta;
2265  // zeta_mapped = -eta;
2266  // }
2267 
2268  // else if (elem->node_id(3) == min_node)
2269  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
2270  // {
2271  // // Case 7
2272  // eta_mapped = zeta;
2273  // zeta_mapped = -eta;
2274  // }
2275  // else
2276  // {
2277  // // Case 8
2278  // eta_mapped = -eta;
2279  // zeta_mapped = zeta;
2280  // }
2281  // }
2282 
2283 
2284  // // Face 5
2285  // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
2286  // {
2287  // const unsigned int min_node = std::min(elem->node_id(4),
2288  // std::min(elem->node_id(5),
2289  // std::min(elem->node_id(6),
2290  // elem->node_id(7))));
2291  // if (elem->node_id(4) == min_node)
2292  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
2293  // {
2294  // // Case 1
2295  // xi_mapped = xi;
2296  // eta_mapped = eta;
2297  // }
2298  // else
2299  // {
2300  // // Case 2
2301  // xi_mapped = eta;
2302  // eta_mapped = xi;
2303  // }
2304 
2305  // else if (elem->node_id(5) == min_node)
2306  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
2307  // {
2308  // // Case 3
2309  // xi_mapped = eta;
2310  // eta_mapped = -xi;
2311  // }
2312  // else
2313  // {
2314  // // Case 4
2315  // xi_mapped = -xi;
2316  // eta_mapped = eta;
2317  // }
2318 
2319  // else if (elem->node_id(6) == min_node)
2320  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
2321  // {
2322  // // Case 5
2323  // xi_mapped = -xi;
2324  // eta_mapped = -eta;
2325  // }
2326  // else
2327  // {
2328  // // Case 6
2329  // xi_mapped = -eta;
2330  // eta_mapped = -xi;
2331  // }
2332 
2333  // else if (elem->node_id(7) == min_node)
2334  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
2335  // {
2336  // // Case 7
2337  // xi_mapped = -eta;
2338  // eta_mapped = xi;
2339  // }
2340  // else
2341  // {
2342  // // Case 8
2343  // xi_mapped = xi;
2344  // eta_mapped = eta;
2345  // }
2346  // }
2347 
2348 
2349  // }
2350 
2351 
2352 
2353  // libmesh_assert_less (j, 3);
2354 
2355  // switch (j)
2356  // {
2357  // // d()/dxi
2358  // case 0:
2359  // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
2360  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2361  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2362 
2363  // // d()/deta
2364  // case 1:
2365  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2366  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
2367  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2368 
2369  // // d()/dzeta
2370  // case 2:
2371  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2372  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2373  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
2374 
2375  // default:
2376  // libmesh_error_msg("Invalid derivative index j = " << j);
2377  // }
2378  // }
2379 
2380 
2381  default:
2382  libmesh_error_msg("Invalid element type = " << type);
2383  }
2384  }
2385 
2386  // 4th-order Bernstein.
2387  case FOURTH:
2388  {
2389  switch (type)
2390  {
2391 
2392  // Bernstein shape functions derivatives on the hexahedral.
2393  case HEX27:
2394  {
2395  const Real eps = 1.e-6;
2396 
2397  libmesh_assert_less (i, 125);
2398  libmesh_assert_less (j, 3);
2399 
2400  switch (j)
2401  {
2402  // d()/dxi
2403  case 0:
2404  {
2405  const Point pp(p(0)+eps, p(1), p(2));
2406  const Point pm(p(0)-eps, p(1), p(2));
2407 
2408  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
2409  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
2410  }
2411 
2412  // d()/deta
2413  case 1:
2414  {
2415  const Point pp(p(0), p(1)+eps, p(2));
2416  const Point pm(p(0), p(1)-eps, p(2));
2417 
2418  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
2419  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
2420  }
2421  // d()/dzeta
2422  case 2:
2423  {
2424  const Point pp(p(0), p(1), p(2)+eps);
2425  const Point pm(p(0), p(1), p(2)-eps);
2426 
2427  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
2428  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
2429  }
2430  default:
2431  libmesh_error_msg("Invalid derivative index j = " << j);
2432  }
2433  }
2434 
2435  // // Compute hex shape functions as a tensor-product
2436  // const Real xi = p(0);
2437  // const Real eta = p(1);
2438  // const Real zeta = p(2);
2439  // Real xi_mapped = p(0);
2440  // Real eta_mapped = p(1);
2441  // Real zeta_mapped = p(2);
2442 
2443  // // The only way to make any sense of this
2444  // // is to look at the mgflo/mg2/mgf documentation
2445  // // and make the cut-out cube!
2446  // // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
2447  // // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
2448  // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4};
2449  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
2450  // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4};
2451 
2452 
2453 
2454  // // handle the edge orientation
2455  // {
2456  // // Edge 0
2457  // if ((i1[i] == 0) && (i2[i] == 0))
2458  // {
2459  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
2460  // xi_mapped = -xi;
2461  // }
2462  // // Edge 1
2463  // else if ((i0[i] == 1) && (i2[i] == 0))
2464  // {
2465  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
2466  // eta_mapped = -eta;
2467  // }
2468  // // Edge 2
2469  // else if ((i1[i] == 1) && (i2[i] == 0))
2470  // {
2471  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
2472  // xi_mapped = -xi;
2473  // }
2474  // // Edge 3
2475  // else if ((i0[i] == 0) && (i2[i] == 0))
2476  // {
2477  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
2478  // eta_mapped = -eta;
2479  // }
2480  // // Edge 4
2481  // else if ((i0[i] == 0) && (i1[i] == 0))
2482  // {
2483  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
2484  // zeta_mapped = -zeta;
2485  // }
2486  // // Edge 5
2487  // else if ((i0[i] == 1) && (i1[i] == 0))
2488  // {
2489  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
2490  // zeta_mapped = -zeta;
2491  // }
2492  // // Edge 6
2493  // else if ((i0[i] == 1) && (i1[i] == 1))
2494  // {
2495  // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
2496  // zeta_mapped = -zeta;
2497  // }
2498  // // Edge 7
2499  // else if ((i0[i] == 0) && (i1[i] == 1))
2500  // {
2501  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
2502  // zeta_mapped = -zeta;
2503  // }
2504  // // Edge 8
2505  // else if ((i1[i] == 0) && (i2[i] == 1))
2506  // {
2507  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
2508  // xi_mapped = -xi;
2509  // }
2510  // // Edge 9
2511  // else if ((i0[i] == 1) && (i2[i] == 1))
2512  // {
2513  // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
2514  // eta_mapped = -eta;
2515  // }
2516  // // Edge 10
2517  // else if ((i1[i] == 1) && (i2[i] == 1))
2518  // {
2519  // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
2520  // xi_mapped = -xi;
2521  // }
2522  // // Edge 11
2523  // else if ((i0[i] == 0) && (i2[i] == 1))
2524  // {
2525  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
2526  // eta_mapped = -eta;
2527  // }
2528  // }
2529 
2530 
2531  // // handle the face orientation
2532  // {
2533  // // Face 0
2534  // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
2535  // {
2536  // const unsigned int min_node = std::min(elem->node_id(1),
2537  // std::min(elem->node_id(2),
2538  // std::min(elem->node_id(0),
2539  // elem->node_id(3))));
2540  // if (elem->node_id(0) == min_node)
2541  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
2542  // {
2543  // // Case 1
2544  // xi_mapped = xi;
2545  // eta_mapped = eta;
2546  // }
2547  // else
2548  // {
2549  // // Case 2
2550  // xi_mapped = eta;
2551  // eta_mapped = xi;
2552  // }
2553 
2554  // else if (elem->node_id(3) == min_node)
2555  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
2556  // {
2557  // // Case 3
2558  // xi_mapped = -eta;
2559  // eta_mapped = xi;
2560  // }
2561  // else
2562  // {
2563  // // Case 4
2564  // xi_mapped = xi;
2565  // eta_mapped = -eta;
2566  // }
2567 
2568  // else if (elem->node_id(2) == min_node)
2569  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
2570  // {
2571  // // Case 5
2572  // xi_mapped = -xi;
2573  // eta_mapped = -eta;
2574  // }
2575  // else
2576  // {
2577  // // Case 6
2578  // xi_mapped = -eta;
2579  // eta_mapped = -xi;
2580  // }
2581 
2582  // else if (elem->node_id(1) == min_node)
2583  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
2584  // {
2585  // // Case 7
2586  // xi_mapped = eta;
2587  // eta_mapped = -xi;
2588  // }
2589  // else
2590  // {
2591  // // Case 8
2592  // xi_mapped = -xi;
2593  // eta_mapped = eta;
2594  // }
2595  // }
2596 
2597 
2598  // // Face 1
2599  // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
2600  // {
2601  // const unsigned int min_node = std::min(elem->node_id(0),
2602  // std::min(elem->node_id(1),
2603  // std::min(elem->node_id(5),
2604  // elem->node_id(4))));
2605  // if (elem->node_id(0) == min_node)
2606  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
2607  // {
2608  // // Case 1
2609  // xi_mapped = xi;
2610  // zeta_mapped = zeta;
2611  // }
2612  // else
2613  // {
2614  // // Case 2
2615  // xi_mapped = zeta;
2616  // zeta_mapped = xi;
2617  // }
2618 
2619  // else if (elem->node_id(1) == min_node)
2620  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
2621  // {
2622  // // Case 3
2623  // xi_mapped = zeta;
2624  // zeta_mapped = -xi;
2625  // }
2626  // else
2627  // {
2628  // // Case 4
2629  // xi_mapped = -xi;
2630  // zeta_mapped = zeta;
2631  // }
2632 
2633  // else if (elem->node_id(5) == min_node)
2634  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
2635  // {
2636  // // Case 5
2637  // xi_mapped = -xi;
2638  // zeta_mapped = -zeta;
2639  // }
2640  // else
2641  // {
2642  // // Case 6
2643  // xi_mapped = -zeta;
2644  // zeta_mapped = -xi;
2645  // }
2646 
2647  // else if (elem->node_id(4) == min_node)
2648  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
2649  // {
2650  // // Case 7
2651  // xi_mapped = -xi;
2652  // zeta_mapped = zeta;
2653  // }
2654  // else
2655  // {
2656  // // Case 8
2657  // xi_mapped = xi;
2658  // zeta_mapped = -zeta;
2659  // }
2660  // }
2661 
2662 
2663  // // Face 2
2664  // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
2665  // {
2666  // const unsigned int min_node = std::min(elem->node_id(1),
2667  // std::min(elem->node_id(2),
2668  // std::min(elem->node_id(6),
2669  // elem->node_id(5))));
2670  // if (elem->node_id(1) == min_node)
2671  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
2672  // {
2673  // // Case 1
2674  // eta_mapped = eta;
2675  // zeta_mapped = zeta;
2676  // }
2677  // else
2678  // {
2679  // // Case 2
2680  // eta_mapped = zeta;
2681  // zeta_mapped = eta;
2682  // }
2683 
2684  // else if (elem->node_id(2) == min_node)
2685  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
2686  // {
2687  // // Case 3
2688  // eta_mapped = zeta;
2689  // zeta_mapped = -eta;
2690  // }
2691  // else
2692  // {
2693  // // Case 4
2694  // eta_mapped = -eta;
2695  // zeta_mapped = zeta;
2696  // }
2697 
2698  // else if (elem->node_id(6) == min_node)
2699  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
2700  // {
2701  // // Case 5
2702  // eta_mapped = -eta;
2703  // zeta_mapped = -zeta;
2704  // }
2705  // else
2706  // {
2707  // // Case 6
2708  // eta_mapped = -zeta;
2709  // zeta_mapped = -eta;
2710  // }
2711 
2712  // else if (elem->node_id(5) == min_node)
2713  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
2714  // {
2715  // // Case 7
2716  // eta_mapped = -zeta;
2717  // zeta_mapped = eta;
2718  // }
2719  // else
2720  // {
2721  // // Case 8
2722  // eta_mapped = eta;
2723  // zeta_mapped = -zeta;
2724  // }
2725  // }
2726 
2727 
2728  // // Face 3
2729  // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
2730  // {
2731  // const unsigned int min_node = std::min(elem->node_id(2),
2732  // std::min(elem->node_id(3),
2733  // std::min(elem->node_id(7),
2734  // elem->node_id(6))));
2735  // if (elem->node_id(3) == min_node)
2736  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
2737  // {
2738  // // Case 1
2739  // xi_mapped = xi;
2740  // zeta_mapped = zeta;
2741  // }
2742  // else
2743  // {
2744  // // Case 2
2745  // xi_mapped = zeta;
2746  // zeta_mapped = xi;
2747  // }
2748 
2749  // else if (elem->node_id(7) == min_node)
2750  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
2751  // {
2752  // // Case 3
2753  // xi_mapped = -zeta;
2754  // zeta_mapped = xi;
2755  // }
2756  // else
2757  // {
2758  // // Case 4
2759  // xi_mapped = xi;
2760  // zeta_mapped = -zeta;
2761  // }
2762 
2763  // else if (elem->node_id(6) == min_node)
2764  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
2765  // {
2766  // // Case 5
2767  // xi_mapped = -xi;
2768  // zeta_mapped = -zeta;
2769  // }
2770  // else
2771  // {
2772  // // Case 6
2773  // xi_mapped = -zeta;
2774  // zeta_mapped = -xi;
2775  // }
2776 
2777  // else if (elem->node_id(2) == min_node)
2778  // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
2779  // {
2780  // // Case 7
2781  // xi_mapped = zeta;
2782  // zeta_mapped = -xi;
2783  // }
2784  // else
2785  // {
2786  // // Case 8
2787  // xi_mapped = -xi;
2788  // zeta_mapped = zeta;
2789  // }
2790  // }
2791 
2792 
2793  // // Face 4
2794  // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
2795  // {
2796  // const unsigned int min_node = std::min(elem->node_id(3),
2797  // std::min(elem->node_id(0),
2798  // std::min(elem->node_id(4),
2799  // elem->node_id(7))));
2800  // if (elem->node_id(0) == min_node)
2801  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
2802  // {
2803  // // Case 1
2804  // eta_mapped = eta;
2805  // zeta_mapped = zeta;
2806  // }
2807  // else
2808  // {
2809  // // Case 2
2810  // eta_mapped = zeta;
2811  // zeta_mapped = eta;
2812  // }
2813 
2814  // else if (elem->node_id(4) == min_node)
2815  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
2816  // {
2817  // // Case 3
2818  // eta_mapped = -zeta;
2819  // zeta_mapped = eta;
2820  // }
2821  // else
2822  // {
2823  // // Case 4
2824  // eta_mapped = eta;
2825  // zeta_mapped = -zeta;
2826  // }
2827 
2828  // else if (elem->node_id(7) == min_node)
2829  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
2830  // {
2831  // // Case 5
2832  // eta_mapped = -eta;
2833  // zeta_mapped = -zeta;
2834  // }
2835  // else
2836  // {
2837  // // Case 6
2838  // eta_mapped = -zeta;
2839  // zeta_mapped = -eta;
2840  // }
2841 
2842  // else if (elem->node_id(3) == min_node)
2843  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
2844  // {
2845  // // Case 7
2846  // eta_mapped = zeta;
2847  // zeta_mapped = -eta;
2848  // }
2849  // else
2850  // {
2851  // // Case 8
2852  // eta_mapped = -eta;
2853  // zeta_mapped = zeta;
2854  // }
2855  // }
2856 
2857 
2858  // // Face 5
2859  // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
2860  // {
2861  // const unsigned int min_node = std::min(elem->node_id(4),
2862  // std::min(elem->node_id(5),
2863  // std::min(elem->node_id(6),
2864  // elem->node_id(7))));
2865  // if (elem->node_id(4) == min_node)
2866  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
2867  // {
2868  // // Case 1
2869  // xi_mapped = xi;
2870  // eta_mapped = eta;
2871  // }
2872  // else
2873  // {
2874  // // Case 2
2875  // xi_mapped = eta;
2876  // eta_mapped = xi;
2877  // }
2878 
2879  // else if (elem->node_id(5) == min_node)
2880  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
2881  // {
2882  // // Case 3
2883  // xi_mapped = eta;
2884  // eta_mapped = -xi;
2885  // }
2886  // else
2887  // {
2888  // // Case 4
2889  // xi_mapped = -xi;
2890  // eta_mapped = eta;
2891  // }
2892 
2893  // else if (elem->node_id(6) == min_node)
2894  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
2895  // {
2896  // // Case 5
2897  // xi_mapped = -xi;
2898  // eta_mapped = -eta;
2899  // }
2900  // else
2901  // {
2902  // // Case 6
2903  // xi_mapped = -eta;
2904  // eta_mapped = -xi;
2905  // }
2906 
2907  // else if (elem->node_id(7) == min_node)
2908  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
2909  // {
2910  // // Case 7
2911  // xi_mapped = -eta;
2912  // eta_mapped = xi;
2913  // }
2914  // else
2915  // {
2916  // // Case 8
2917  // xi_mapped = xi;
2918  // eta_mapped = eta;
2919  // }
2920  // }
2921 
2922 
2923  // }
2924 
2925 
2926 
2927  // libmesh_assert_less (j, 3);
2928 
2929  // switch (j)
2930  // {
2931  // // d()/dxi
2932  // case 0:
2933  // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
2934  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2935  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2936 
2937  // // d()/deta
2938  // case 1:
2939  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2940  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
2941  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2942 
2943  // // d()/dzeta
2944  // case 2:
2945  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2946  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2947  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
2948 
2949  // default:
2950  // libmesh_error_msg("Invalid derivative index j = " << j);
2951  // }
2952  // }
2953 
2954 
2955  default:
2956  libmesh_error_msg("Invalid element type = " << type);
2957  }
2958  }
2959 
2960 
2961  default:
2962  libmesh_error_msg("Invalid totalorder = " << totalorder);
2963  }
2964 
2965 #else // LIBMESH_DIM != 3
2966  libmesh_ignore(elem, order, i, j, p, add_p_level);
2967  libmesh_not_implemented();
2968 #endif
2969 }
2970 
2971 
2972 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2973 
2974 template <>
2975 Real FE<3,BERNSTEIN>::shape_second_deriv(const ElemType,
2976  const Order,
2977  const unsigned int,
2978  const unsigned int,
2979  const Point &)
2980 {
2981  static bool warning_given = false;
2982 
2983  if (!warning_given)
2984  libMesh::err << "Second derivatives for Bernstein elements "
2985  << "are not yet implemented!"
2986  << std::endl;
2987 
2988  warning_given = true;
2989  return 0.;
2990 }
2991 
2992 
2993 
2994 template <>
2995 Real FE<3,BERNSTEIN>::shape_second_deriv(const Elem *,
2996  const Order,
2997  const unsigned int,
2998  const unsigned int,
2999  const Point &,
3000  const bool)
3001 {
3002  static bool warning_given = false;
3003 
3004  if (!warning_given)
3005  libMesh::err << "Second derivatives for Bernstein elements "
3006  << "are not yet implemented!"
3007  << std::endl;
3008 
3009  warning_given = true;
3010  return 0.;
3011 }
3012 
3013 #endif
3014 
3015 } // namespace libMesh
3016 
3017 
3018 
3019 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::HEX20
Definition: enum_elem_type.h:48
libMesh::HEX8
Definition: enum_elem_type.h:47
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::TET10
Definition: enum_elem_type.h:46
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::SECOND
Definition: enum_order.h:43
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1955
libMesh::TET4
Definition: enum_elem_type.h:45
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::HEX27
Definition: enum_elem_type.h:49
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FE
A specific instantiation of the FEBase class.
Definition: fe.h:95
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:526
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::FOURTH
Definition: enum_order.h:45
libMesh::EDGE3
Definition: enum_elem_type.h:36
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::THIRD
Definition: enum_order.h:44
libMesh::err
OStreamProxy err
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::FIRST
Definition: enum_order.h:42
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
Generated on Sat Jan 25 2020 12:06:57 for libMesh by   doxygen 1.8.16