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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo