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