fe_bernstein_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
25 
26 #include "libmesh/fe.h"
27 #include "libmesh/elem.h"
28 
29 
30 namespace libMesh
31 {
32 
33 
34 
35 template <>
36 Real FE<3,BERNSTEIN>::shape(const ElemType,
37  const Order,
38  const unsigned int,
39  const Point &)
40 {
41  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
42  return 0.;
43 }
44 
45 
46 
47 template <>
48 Real FE<3,BERNSTEIN>::shape(const Elem * elem,
49  const Order order,
50  const unsigned int i,
51  const Point & p)
52 {
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 
1378 
1379 
1380 
1381 template <>
1382 Real FE<3,BERNSTEIN>::shape_deriv(const ElemType,
1383  const Order,
1384  const unsigned int,
1385  const unsigned int,
1386  const Point & )
1387 {
1388  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
1389  return 0.;
1390 }
1391 
1392 
1393 
1394 template <>
1395 Real FE<3,BERNSTEIN>::shape_deriv(const Elem * elem,
1396  const Order order,
1397  const unsigned int i,
1398  const unsigned int j,
1399  const Point & p)
1400 {
1401 
1402 #if LIBMESH_DIM == 3
1403  libmesh_assert(elem);
1404  const ElemType type = elem->type();
1405 
1406  const Order totalorder = static_cast<Order>(order + elem->p_level());
1407 
1408  libmesh_assert_less (j, 3);
1409 
1410  switch (totalorder)
1411  {
1412  // 1st order Bernstein.
1413  case FIRST:
1414  {
1415  switch (type)
1416  {
1417  // Bernstein shape functions on the tetrahedron.
1418  case TET4:
1419  case TET10:
1420  {
1421  // I have been lazy here and am using finite differences
1422  // to compute the derivatives!
1423  const Real eps = 1.e-6;
1424 
1425  libmesh_assert_less (i, 4);
1426  libmesh_assert_less (j, 3);
1427 
1428 
1429  switch (j)
1430  {
1431  // d()/dxi
1432  case 0:
1433  {
1434  const Point pp(p(0)+eps, p(1), p(2));
1435  const Point pm(p(0)-eps, p(1), p(2));
1436 
1437  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1438  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1439  }
1440 
1441  // d()/deta
1442  case 1:
1443  {
1444  const Point pp(p(0), p(1)+eps, p(2));
1445  const Point pm(p(0), p(1)-eps, p(2));
1446 
1447  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1448  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1449  }
1450  // d()/dzeta
1451  case 2:
1452  {
1453  const Point pp(p(0), p(1), p(2)+eps);
1454  const Point pm(p(0), p(1), p(2)-eps);
1455 
1456  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1457  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1458  }
1459  default:
1460  libmesh_error_msg("Invalid derivative index j = " << j);
1461  }
1462  }
1463 
1464 
1465 
1466 
1467  // Bernstein shape functions on the hexahedral.
1468  case HEX8:
1469  case HEX20:
1470  case HEX27:
1471  {
1472  libmesh_assert_less (i, 8);
1473 
1474  // Compute hex shape functions as a tensor-product
1475  const Real xi = p(0);
1476  const Real eta = p(1);
1477  const Real zeta = p(2);
1478 
1479  // The only way to make any sense of this
1480  // is to look at the mgflo/mg2/mgf documentation
1481  // and make the cut-out cube!
1482  // 0 1 2 3 4 5 6 7
1483  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
1484  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
1485  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
1486 
1487  switch (j)
1488  {
1489  // d()/dxi
1490  case 0:
1491  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1492  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1493  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1494 
1495  // d()/deta
1496  case 1:
1497  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1498  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1499  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1500 
1501  // d()/dzeta
1502  case 2:
1503  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1504  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1505  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1506 
1507  default:
1508  libmesh_error_msg("Invalid derivative index j = " << j);
1509  }
1510  }
1511 
1512  default:
1513  libmesh_error_msg("Invalid element type = " << type);
1514  }
1515  }
1516 
1517 
1518 
1519 
1520  case SECOND:
1521  {
1522  switch (type)
1523  {
1524  // Bernstein shape functions on the tetrahedron.
1525  case TET10:
1526  {
1527  // I have been lazy here and am using finite differences
1528  // to compute the derivatives!
1529  const Real eps = 1.e-6;
1530 
1531  libmesh_assert_less (i, 10);
1532  libmesh_assert_less (j, 3);
1533 
1534 
1535  switch (j)
1536  {
1537  // d()/dxi
1538  case 0:
1539  {
1540  const Point pp(p(0)+eps, p(1), p(2));
1541  const Point pm(p(0)-eps, p(1), p(2));
1542 
1543  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1544  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1545  }
1546 
1547  // d()/deta
1548  case 1:
1549  {
1550  const Point pp(p(0), p(1)+eps, p(2));
1551  const Point pm(p(0), p(1)-eps, p(2));
1552 
1553  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1554  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1555  }
1556  // d()/dzeta
1557  case 2:
1558  {
1559  const Point pp(p(0), p(1), p(2)+eps);
1560  const Point pm(p(0), p(1), p(2)-eps);
1561 
1562  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1563  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1564  }
1565  default:
1566  libmesh_error_msg("Invalid derivative index j = " << j);
1567  }
1568  }
1569 
1570  // Bernstein shape functions on the hexahedral.
1571  case HEX20:
1572  {
1573  libmesh_assert_less (i, 20);
1574 
1575  // Compute hex shape functions as a tensor-product
1576  const Real xi = p(0);
1577  const Real eta = p(1);
1578  const Real zeta = p(2);
1579 
1580  // The only way to make any sense of this
1581  // is to look at the mgflo/mg2/mgf documentation
1582  // and make the cut-out cube!
1583  // 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
1584  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};
1585  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};
1586  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};
1587  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};
1588  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};
1589  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};
1590  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};
1591  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};
1592  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};
1593  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};
1594 
1595  switch (j)
1596  {
1597  // d()/dxi
1598  case 0:
1599  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1600  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1601  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)
1602  +scal20[i]*
1603  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[20], 0, xi)*
1604  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[20], eta)*
1605  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[20], zeta)
1606  +scal21[i]*
1607  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[21], 0, xi)*
1608  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[21], eta)*
1609  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[21], zeta)
1610  +scal22[i]*
1611  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[22], 0, xi)*
1612  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[22], eta)*
1613  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[22], zeta)
1614  +scal23[i]*
1615  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[23], 0, xi)*
1616  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[23], eta)*
1617  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[23], zeta)
1618  +scal24[i]*
1619  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[24], 0, xi)*
1620  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[24], eta)*
1621  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[24], zeta)
1622  +scal25[i]*
1623  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[25], 0, xi)*
1624  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[25], eta)*
1625  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[25], zeta)
1626  +scal26[i]*
1627  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[26], 0, xi)*
1628  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[26], eta)*
1629  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[26], zeta));
1630 
1631  // d()/deta
1632  case 1:
1633  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1634  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1635  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta)
1636  +scal20[i]*
1637  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[20], xi)*
1638  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[20], 0, eta)*
1639  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[20], zeta)
1640  +scal21[i]*
1641  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[21], xi)*
1642  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[21], 0, eta)*
1643  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[21], zeta)
1644  +scal22[i]*
1645  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[22], xi)*
1646  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[22], 0, eta)*
1647  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[22], zeta)
1648  +scal23[i]*
1649  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[23], xi)*
1650  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[23], 0, eta)*
1651  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[23], zeta)
1652  +scal24[i]*
1653  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[24], xi)*
1654  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[24], 0, eta)*
1655  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[24], zeta)
1656  +scal25[i]*
1657  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[25], xi)*
1658  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[25], 0, eta)*
1659  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[25], zeta)
1660  +scal26[i]*
1661  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[26], xi)*
1662  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[26], 0, eta)*
1663  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[26], zeta));
1664 
1665  // d()/dzeta
1666  case 2:
1667  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1668  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1669  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta)
1670  +scal20[i]*
1671  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[20], xi)*
1672  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[20], eta)*
1673  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[20], 0, zeta)
1674  +scal21[i]*
1675  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[21], xi)*
1676  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[21], eta)*
1677  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[21], 0, zeta)
1678  +scal22[i]*
1679  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[22], xi)*
1680  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[22], eta)*
1681  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[22], 0, zeta)
1682  +scal23[i]*
1683  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[23], xi)*
1684  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[23], eta)*
1685  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[23], 0, zeta)
1686  +scal24[i]*
1687  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[24], xi)*
1688  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[24], eta)*
1689  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[24], 0, zeta)
1690  +scal25[i]*
1691  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[25], xi)*
1692  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[25], eta)*
1693  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[25], 0, zeta)
1694  +scal26[i]*
1695  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[26], xi)*
1696  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[26], eta)*
1697  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[26], 0, zeta));
1698 
1699  default:
1700  libmesh_error_msg("Invalid derivative index j = " << j);
1701  }
1702  }
1703 
1704  // Bernstein shape functions on the hexahedral.
1705  case HEX27:
1706  {
1707  libmesh_assert_less (i, 27);
1708 
1709  // Compute hex shape functions as a tensor-product
1710  const Real xi = p(0);
1711  const Real eta = p(1);
1712  const Real zeta = p(2);
1713 
1714  // The only way to make any sense of this
1715  // is to look at the mgflo/mg2/mgf documentation
1716  // and make the cut-out cube!
1717  // 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
1718  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};
1719  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};
1720  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};
1721 
1722  switch (j)
1723  {
1724  // d()/dxi
1725  case 0:
1726  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1727  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1728  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1729 
1730  // d()/deta
1731  case 1:
1732  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1733  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1734  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1735 
1736  // d()/dzeta
1737  case 2:
1738  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1739  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1740  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1741 
1742  default:
1743  libmesh_error_msg("Invalid derivative index j = " << j);
1744  }
1745  }
1746 
1747 
1748  default:
1749  libmesh_error_msg("Invalid element type = " << type);
1750  }
1751  }
1752 
1753 
1754 
1755  // 3rd-order Bernstein.
1756  case THIRD:
1757  {
1758  switch (type)
1759  {
1760 
1761  // // Bernstein shape functions derivatives.
1762  // case TET10:
1763  // {
1764  // // I have been lazy here and am using finite differences
1765  // // to compute the derivatives!
1766  // const Real eps = 1.e-6;
1767 
1768  // libmesh_assert_less (i, 20);
1769  // libmesh_assert_less (j, 3);
1770 
1771  // switch (j)
1772  // {
1773  // // d()/dxi
1774  // case 0:
1775  // {
1776  // const Point pp(p(0)+eps, p(1), p(2));
1777  // const Point pm(p(0)-eps, p(1), p(2));
1778 
1779  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1780  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1781  // }
1782 
1783  // // d()/deta
1784  // case 1:
1785  // {
1786  // const Point pp(p(0), p(1)+eps, p(2));
1787  // const Point pm(p(0), p(1)-eps, p(2));
1788 
1789  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1790  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1791  // }
1792  // // d()/dzeta
1793  // case 2:
1794  // {
1795  // const Point pp(p(0), p(1), p(2)+eps);
1796  // const Point pm(p(0), p(1), p(2)-eps);
1797 
1798  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1799  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1800  // }
1801  // default:
1802  // libmesh_error_msg("Invalid derivative index j = " << j);
1803  // }
1804 
1805 
1806  // }
1807 
1808 
1809  // Bernstein shape functions on the hexahedral.
1810  case HEX27:
1811  {
1812  // I have been lazy here and am using finite differences
1813  // to compute the derivatives!
1814  const Real eps = 1.e-6;
1815 
1816  libmesh_assert_less (i, 64);
1817  libmesh_assert_less (j, 3);
1818 
1819  switch (j)
1820  {
1821  // d()/dxi
1822  case 0:
1823  {
1824  const Point pp(p(0)+eps, p(1), p(2));
1825  const Point pm(p(0)-eps, p(1), p(2));
1826 
1827  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1828  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1829  }
1830 
1831  // d()/deta
1832  case 1:
1833  {
1834  const Point pp(p(0), p(1)+eps, p(2));
1835  const Point pm(p(0), p(1)-eps, p(2));
1836 
1837  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1838  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1839  }
1840  // d()/dzeta
1841  case 2:
1842  {
1843  const Point pp(p(0), p(1), p(2)+eps);
1844  const Point pm(p(0), p(1), p(2)-eps);
1845 
1846  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1847  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1848  }
1849  default:
1850  libmesh_error_msg("Invalid derivative index j = " << j);
1851  }
1852 
1853  }
1854 
1855  // // Compute hex shape functions as a tensor-product
1856  // const Real xi = p(0);
1857  // const Real eta = p(1);
1858  // const Real zeta = p(2);
1859  // Real xi_mapped = p(0);
1860  // Real eta_mapped = p(1);
1861  // Real zeta_mapped = p(2);
1862 
1863  // // The only way to make any sense of this
1864  // // is to look at the mgflo/mg2/mgf documentation
1865  // // and make the cut-out cube!
1866  // // 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
1867  // // 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
1868  // 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};
1869  // 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};
1870  // 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};
1871 
1872 
1873 
1874  // // handle the edge orientation
1875  // {
1876  // // Edge 0
1877  // if ((i1[i] == 0) && (i2[i] == 0))
1878  // {
1879  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
1880  // xi_mapped = -xi;
1881  // }
1882  // // Edge 1
1883  // else if ((i0[i] == 1) && (i2[i] == 0))
1884  // {
1885  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
1886  // eta_mapped = -eta;
1887  // }
1888  // // Edge 2
1889  // else if ((i1[i] == 1) && (i2[i] == 0))
1890  // {
1891  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
1892  // xi_mapped = -xi;
1893  // }
1894  // // Edge 3
1895  // else if ((i0[i] == 0) && (i2[i] == 0))
1896  // {
1897  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
1898  // eta_mapped = -eta;
1899  // }
1900  // // Edge 4
1901  // else if ((i0[i] == 0) && (i1[i] == 0))
1902  // {
1903  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
1904  // zeta_mapped = -zeta;
1905  // }
1906  // // Edge 5
1907  // else if ((i0[i] == 1) && (i1[i] == 0))
1908  // {
1909  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
1910  // zeta_mapped = -zeta;
1911  // }
1912  // // Edge 6
1913  // else if ((i0[i] == 1) && (i1[i] == 1))
1914  // {
1915  // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
1916  // zeta_mapped = -zeta;
1917  // }
1918  // // Edge 7
1919  // else if ((i0[i] == 0) && (i1[i] == 1))
1920  // {
1921  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
1922  // zeta_mapped = -zeta;
1923  // }
1924  // // Edge 8
1925  // else if ((i1[i] == 0) && (i2[i] == 1))
1926  // {
1927  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
1928  // xi_mapped = -xi;
1929  // }
1930  // // Edge 9
1931  // else if ((i0[i] == 1) && (i2[i] == 1))
1932  // {
1933  // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
1934  // eta_mapped = -eta;
1935  // }
1936  // // Edge 10
1937  // else if ((i1[i] == 1) && (i2[i] == 1))
1938  // {
1939  // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
1940  // xi_mapped = -xi;
1941  // }
1942  // // Edge 11
1943  // else if ((i0[i] == 0) && (i2[i] == 1))
1944  // {
1945  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
1946  // eta_mapped = -eta;
1947  // }
1948  // }
1949 
1950 
1951  // // handle the face orientation
1952  // {
1953  // // Face 0
1954  // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
1955  // {
1956  // const unsigned int min_node = std::min(elem->node_id(1),
1957  // std::min(elem->node_id(2),
1958  // std::min(elem->node_id(0),
1959  // elem->node_id(3))));
1960  // if (elem->node_id(0) == min_node)
1961  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
1962  // {
1963  // // Case 1
1964  // xi_mapped = xi;
1965  // eta_mapped = eta;
1966  // }
1967  // else
1968  // {
1969  // // Case 2
1970  // xi_mapped = eta;
1971  // eta_mapped = xi;
1972  // }
1973 
1974  // else if (elem->node_id(3) == min_node)
1975  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
1976  // {
1977  // // Case 3
1978  // xi_mapped = -eta;
1979  // eta_mapped = xi;
1980  // }
1981  // else
1982  // {
1983  // // Case 4
1984  // xi_mapped = xi;
1985  // eta_mapped = -eta;
1986  // }
1987 
1988  // else if (elem->node_id(2) == min_node)
1989  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
1990  // {
1991  // // Case 5
1992  // xi_mapped = -xi;
1993  // eta_mapped = -eta;
1994  // }
1995  // else
1996  // {
1997  // // Case 6
1998  // xi_mapped = -eta;
1999  // eta_mapped = -xi;
2000  // }
2001 
2002  // else if (elem->node_id(1) == min_node)
2003  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
2004  // {
2005  // // Case 7
2006  // xi_mapped = eta;
2007  // eta_mapped = -xi;
2008  // }
2009  // else
2010  // {
2011  // // Case 8
2012  // xi_mapped = -xi;
2013  // eta_mapped = eta;
2014  // }
2015  // }
2016 
2017 
2018  // // Face 1
2019  // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
2020  // {
2021  // const unsigned int min_node = std::min(elem->node_id(0),
2022  // std::min(elem->node_id(1),
2023  // std::min(elem->node_id(5),
2024  // elem->node_id(4))));
2025  // if (elem->node_id(0) == min_node)
2026  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
2027  // {
2028  // // Case 1
2029  // xi_mapped = xi;
2030  // zeta_mapped = zeta;
2031  // }
2032  // else
2033  // {
2034  // // Case 2
2035  // xi_mapped = zeta;
2036  // zeta_mapped = xi;
2037  // }
2038 
2039  // else if (elem->node_id(1) == min_node)
2040  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
2041  // {
2042  // // Case 3
2043  // xi_mapped = zeta;
2044  // zeta_mapped = -xi;
2045  // }
2046  // else
2047  // {
2048  // // Case 4
2049  // xi_mapped = -xi;
2050  // zeta_mapped = zeta;
2051  // }
2052 
2053  // else if (elem->node_id(5) == min_node)
2054  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
2055  // {
2056  // // Case 5
2057  // xi_mapped = -xi;
2058  // zeta_mapped = -zeta;
2059  // }
2060  // else
2061  // {
2062  // // Case 6
2063  // xi_mapped = -zeta;
2064  // zeta_mapped = -xi;
2065  // }
2066 
2067  // else if (elem->node_id(4) == min_node)
2068  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
2069  // {
2070  // // Case 7
2071  // xi_mapped = -xi;
2072  // zeta_mapped = zeta;
2073  // }
2074  // else
2075  // {
2076  // // Case 8
2077  // xi_mapped = xi;
2078  // zeta_mapped = -zeta;
2079  // }
2080  // }
2081 
2082 
2083  // // Face 2
2084  // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
2085  // {
2086  // const unsigned int min_node = std::min(elem->node_id(1),
2087  // std::min(elem->node_id(2),
2088  // std::min(elem->node_id(6),
2089  // elem->node_id(5))));
2090  // if (elem->node_id(1) == min_node)
2091  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
2092  // {
2093  // // Case 1
2094  // eta_mapped = eta;
2095  // zeta_mapped = zeta;
2096  // }
2097  // else
2098  // {
2099  // // Case 2
2100  // eta_mapped = zeta;
2101  // zeta_mapped = eta;
2102  // }
2103 
2104  // else if (elem->node_id(2) == min_node)
2105  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
2106  // {
2107  // // Case 3
2108  // eta_mapped = zeta;
2109  // zeta_mapped = -eta;
2110  // }
2111  // else
2112  // {
2113  // // Case 4
2114  // eta_mapped = -eta;
2115  // zeta_mapped = zeta;
2116  // }
2117 
2118  // else if (elem->node_id(6) == min_node)
2119  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
2120  // {
2121  // // Case 5
2122  // eta_mapped = -eta;
2123  // zeta_mapped = -zeta;
2124  // }
2125  // else
2126  // {
2127  // // Case 6
2128  // eta_mapped = -zeta;
2129  // zeta_mapped = -eta;
2130  // }
2131 
2132  // else if (elem->node_id(5) == min_node)
2133  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
2134  // {
2135  // // Case 7
2136  // eta_mapped = -zeta;
2137  // zeta_mapped = eta;
2138  // }
2139  // else
2140  // {
2141  // // Case 8
2142  // eta_mapped = eta;
2143  // zeta_mapped = -zeta;
2144  // }
2145  // }
2146 
2147 
2148  // // Face 3
2149  // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
2150  // {
2151  // const unsigned int min_node = std::min(elem->node_id(2),
2152  // std::min(elem->node_id(3),
2153  // std::min(elem->node_id(7),
2154  // elem->node_id(6))));
2155  // if (elem->node_id(3) == min_node)
2156  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
2157  // {
2158  // // Case 1
2159  // xi_mapped = xi;
2160  // zeta_mapped = zeta;
2161  // }
2162  // else
2163  // {
2164  // // Case 2
2165  // xi_mapped = zeta;
2166  // zeta_mapped = xi;
2167  // }
2168 
2169  // else if (elem->node_id(7) == min_node)
2170  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
2171  // {
2172  // // Case 3
2173  // xi_mapped = -zeta;
2174  // zeta_mapped = xi;
2175  // }
2176  // else
2177  // {
2178  // // Case 4
2179  // xi_mapped = xi;
2180  // zeta_mapped = -zeta;
2181  // }
2182 
2183  // else if (elem->node_id(6) == min_node)
2184  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
2185  // {
2186  // // Case 5
2187  // xi_mapped = -xi;
2188  // zeta_mapped = -zeta;
2189  // }
2190  // else
2191  // {
2192  // // Case 6
2193  // xi_mapped = -zeta;
2194  // zeta_mapped = -xi;
2195  // }
2196 
2197  // else if (elem->node_id(2) == min_node)
2198  // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
2199  // {
2200  // // Case 7
2201  // xi_mapped = zeta;
2202  // zeta_mapped = -xi;
2203  // }
2204  // else
2205  // {
2206  // // Case 8
2207  // xi_mapped = -xi;
2208  // zeta_mapped = zeta;
2209  // }
2210  // }
2211 
2212 
2213  // // Face 4
2214  // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
2215  // {
2216  // const unsigned int min_node = std::min(elem->node_id(3),
2217  // std::min(elem->node_id(0),
2218  // std::min(elem->node_id(4),
2219  // elem->node_id(7))));
2220  // if (elem->node_id(0) == min_node)
2221  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
2222  // {
2223  // // Case 1
2224  // eta_mapped = eta;
2225  // zeta_mapped = zeta;
2226  // }
2227  // else
2228  // {
2229  // // Case 2
2230  // eta_mapped = zeta;
2231  // zeta_mapped = eta;
2232  // }
2233 
2234  // else if (elem->node_id(4) == min_node)
2235  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
2236  // {
2237  // // Case 3
2238  // eta_mapped = -zeta;
2239  // zeta_mapped = eta;
2240  // }
2241  // else
2242  // {
2243  // // Case 4
2244  // eta_mapped = eta;
2245  // zeta_mapped = -zeta;
2246  // }
2247 
2248  // else if (elem->node_id(7) == min_node)
2249  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
2250  // {
2251  // // Case 5
2252  // eta_mapped = -eta;
2253  // zeta_mapped = -zeta;
2254  // }
2255  // else
2256  // {
2257  // // Case 6
2258  // eta_mapped = -zeta;
2259  // zeta_mapped = -eta;
2260  // }
2261 
2262  // else if (elem->node_id(3) == min_node)
2263  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
2264  // {
2265  // // Case 7
2266  // eta_mapped = zeta;
2267  // zeta_mapped = -eta;
2268  // }
2269  // else
2270  // {
2271  // // Case 8
2272  // eta_mapped = -eta;
2273  // zeta_mapped = zeta;
2274  // }
2275  // }
2276 
2277 
2278  // // Face 5
2279  // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
2280  // {
2281  // const unsigned int min_node = std::min(elem->node_id(4),
2282  // std::min(elem->node_id(5),
2283  // std::min(elem->node_id(6),
2284  // elem->node_id(7))));
2285  // if (elem->node_id(4) == min_node)
2286  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
2287  // {
2288  // // Case 1
2289  // xi_mapped = xi;
2290  // eta_mapped = eta;
2291  // }
2292  // else
2293  // {
2294  // // Case 2
2295  // xi_mapped = eta;
2296  // eta_mapped = xi;
2297  // }
2298 
2299  // else if (elem->node_id(5) == min_node)
2300  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
2301  // {
2302  // // Case 3
2303  // xi_mapped = eta;
2304  // eta_mapped = -xi;
2305  // }
2306  // else
2307  // {
2308  // // Case 4
2309  // xi_mapped = -xi;
2310  // eta_mapped = eta;
2311  // }
2312 
2313  // else if (elem->node_id(6) == min_node)
2314  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
2315  // {
2316  // // Case 5
2317  // xi_mapped = -xi;
2318  // eta_mapped = -eta;
2319  // }
2320  // else
2321  // {
2322  // // Case 6
2323  // xi_mapped = -eta;
2324  // eta_mapped = -xi;
2325  // }
2326 
2327  // else if (elem->node_id(7) == min_node)
2328  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
2329  // {
2330  // // Case 7
2331  // xi_mapped = -eta;
2332  // eta_mapped = xi;
2333  // }
2334  // else
2335  // {
2336  // // Case 8
2337  // xi_mapped = xi;
2338  // eta_mapped = eta;
2339  // }
2340  // }
2341 
2342 
2343  // }
2344 
2345 
2346 
2347  // libmesh_assert_less (j, 3);
2348 
2349  // switch (j)
2350  // {
2351  // // d()/dxi
2352  // case 0:
2353  // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
2354  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2355  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2356 
2357  // // d()/deta
2358  // case 1:
2359  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2360  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
2361  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2362 
2363  // // d()/dzeta
2364  // case 2:
2365  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2366  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2367  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
2368 
2369  // default:
2370  // libmesh_error_msg("Invalid derivative index j = " << j);
2371  // }
2372  // }
2373 
2374 
2375  default:
2376  libmesh_error_msg("Invalid element type = " << type);
2377  }
2378  }
2379 
2380  // 4th-order Bernstein.
2381  case FOURTH:
2382  {
2383  switch (type)
2384  {
2385 
2386  // Bernstein shape functions derivatives on the hexahedral.
2387  case HEX27:
2388  {
2389  const Real eps = 1.e-6;
2390 
2391  libmesh_assert_less (i, 125);
2392  libmesh_assert_less (j, 3);
2393 
2394  switch (j)
2395  {
2396  // d()/dxi
2397  case 0:
2398  {
2399  const Point pp(p(0)+eps, p(1), p(2));
2400  const Point pm(p(0)-eps, p(1), p(2));
2401 
2402  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
2403  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
2404  }
2405 
2406  // d()/deta
2407  case 1:
2408  {
2409  const Point pp(p(0), p(1)+eps, p(2));
2410  const Point pm(p(0), p(1)-eps, p(2));
2411 
2412  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
2413  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
2414  }
2415  // d()/dzeta
2416  case 2:
2417  {
2418  const Point pp(p(0), p(1), p(2)+eps);
2419  const Point pm(p(0), p(1), p(2)-eps);
2420 
2421  return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
2422  FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
2423  }
2424  default:
2425  libmesh_error_msg("Invalid derivative index j = " << j);
2426  }
2427  }
2428 
2429  // // Compute hex shape functions as a tensor-product
2430  // const Real xi = p(0);
2431  // const Real eta = p(1);
2432  // const Real zeta = p(2);
2433  // Real xi_mapped = p(0);
2434  // Real eta_mapped = p(1);
2435  // Real zeta_mapped = p(2);
2436 
2437  // // The only way to make any sense of this
2438  // // is to look at the mgflo/mg2/mgf documentation
2439  // // and make the cut-out cube!
2440  // // 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
2441  // // 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
2442  // 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};
2443  // 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};
2444  // 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};
2445 
2446 
2447 
2448  // // handle the edge orientation
2449  // {
2450  // // Edge 0
2451  // if ((i1[i] == 0) && (i2[i] == 0))
2452  // {
2453  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
2454  // xi_mapped = -xi;
2455  // }
2456  // // Edge 1
2457  // else if ((i0[i] == 1) && (i2[i] == 0))
2458  // {
2459  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
2460  // eta_mapped = -eta;
2461  // }
2462  // // Edge 2
2463  // else if ((i1[i] == 1) && (i2[i] == 0))
2464  // {
2465  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
2466  // xi_mapped = -xi;
2467  // }
2468  // // Edge 3
2469  // else if ((i0[i] == 0) && (i2[i] == 0))
2470  // {
2471  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
2472  // eta_mapped = -eta;
2473  // }
2474  // // Edge 4
2475  // else if ((i0[i] == 0) && (i1[i] == 0))
2476  // {
2477  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
2478  // zeta_mapped = -zeta;
2479  // }
2480  // // Edge 5
2481  // else if ((i0[i] == 1) && (i1[i] == 0))
2482  // {
2483  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
2484  // zeta_mapped = -zeta;
2485  // }
2486  // // Edge 6
2487  // else if ((i0[i] == 1) && (i1[i] == 1))
2488  // {
2489  // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
2490  // zeta_mapped = -zeta;
2491  // }
2492  // // Edge 7
2493  // else if ((i0[i] == 0) && (i1[i] == 1))
2494  // {
2495  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
2496  // zeta_mapped = -zeta;
2497  // }
2498  // // Edge 8
2499  // else if ((i1[i] == 0) && (i2[i] == 1))
2500  // {
2501  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
2502  // xi_mapped = -xi;
2503  // }
2504  // // Edge 9
2505  // else if ((i0[i] == 1) && (i2[i] == 1))
2506  // {
2507  // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
2508  // eta_mapped = -eta;
2509  // }
2510  // // Edge 10
2511  // else if ((i1[i] == 1) && (i2[i] == 1))
2512  // {
2513  // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
2514  // xi_mapped = -xi;
2515  // }
2516  // // Edge 11
2517  // else if ((i0[i] == 0) && (i2[i] == 1))
2518  // {
2519  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
2520  // eta_mapped = -eta;
2521  // }
2522  // }
2523 
2524 
2525  // // handle the face orientation
2526  // {
2527  // // Face 0
2528  // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
2529  // {
2530  // const unsigned int min_node = std::min(elem->node_id(1),
2531  // std::min(elem->node_id(2),
2532  // std::min(elem->node_id(0),
2533  // elem->node_id(3))));
2534  // if (elem->node_id(0) == min_node)
2535  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
2536  // {
2537  // // Case 1
2538  // xi_mapped = xi;
2539  // eta_mapped = eta;
2540  // }
2541  // else
2542  // {
2543  // // Case 2
2544  // xi_mapped = eta;
2545  // eta_mapped = xi;
2546  // }
2547 
2548  // else if (elem->node_id(3) == min_node)
2549  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
2550  // {
2551  // // Case 3
2552  // xi_mapped = -eta;
2553  // eta_mapped = xi;
2554  // }
2555  // else
2556  // {
2557  // // Case 4
2558  // xi_mapped = xi;
2559  // eta_mapped = -eta;
2560  // }
2561 
2562  // else if (elem->node_id(2) == min_node)
2563  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
2564  // {
2565  // // Case 5
2566  // xi_mapped = -xi;
2567  // eta_mapped = -eta;
2568  // }
2569  // else
2570  // {
2571  // // Case 6
2572  // xi_mapped = -eta;
2573  // eta_mapped = -xi;
2574  // }
2575 
2576  // else if (elem->node_id(1) == min_node)
2577  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
2578  // {
2579  // // Case 7
2580  // xi_mapped = eta;
2581  // eta_mapped = -xi;
2582  // }
2583  // else
2584  // {
2585  // // Case 8
2586  // xi_mapped = -xi;
2587  // eta_mapped = eta;
2588  // }
2589  // }
2590 
2591 
2592  // // Face 1
2593  // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
2594  // {
2595  // const unsigned int min_node = std::min(elem->node_id(0),
2596  // std::min(elem->node_id(1),
2597  // std::min(elem->node_id(5),
2598  // elem->node_id(4))));
2599  // if (elem->node_id(0) == min_node)
2600  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
2601  // {
2602  // // Case 1
2603  // xi_mapped = xi;
2604  // zeta_mapped = zeta;
2605  // }
2606  // else
2607  // {
2608  // // Case 2
2609  // xi_mapped = zeta;
2610  // zeta_mapped = xi;
2611  // }
2612 
2613  // else if (elem->node_id(1) == min_node)
2614  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
2615  // {
2616  // // Case 3
2617  // xi_mapped = zeta;
2618  // zeta_mapped = -xi;
2619  // }
2620  // else
2621  // {
2622  // // Case 4
2623  // xi_mapped = -xi;
2624  // zeta_mapped = zeta;
2625  // }
2626 
2627  // else if (elem->node_id(5) == min_node)
2628  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
2629  // {
2630  // // Case 5
2631  // xi_mapped = -xi;
2632  // zeta_mapped = -zeta;
2633  // }
2634  // else
2635  // {
2636  // // Case 6
2637  // xi_mapped = -zeta;
2638  // zeta_mapped = -xi;
2639  // }
2640 
2641  // else if (elem->node_id(4) == min_node)
2642  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
2643  // {
2644  // // Case 7
2645  // xi_mapped = -xi;
2646  // zeta_mapped = zeta;
2647  // }
2648  // else
2649  // {
2650  // // Case 8
2651  // xi_mapped = xi;
2652  // zeta_mapped = -zeta;
2653  // }
2654  // }
2655 
2656 
2657  // // Face 2
2658  // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
2659  // {
2660  // const unsigned int min_node = std::min(elem->node_id(1),
2661  // std::min(elem->node_id(2),
2662  // std::min(elem->node_id(6),
2663  // elem->node_id(5))));
2664  // if (elem->node_id(1) == min_node)
2665  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
2666  // {
2667  // // Case 1
2668  // eta_mapped = eta;
2669  // zeta_mapped = zeta;
2670  // }
2671  // else
2672  // {
2673  // // Case 2
2674  // eta_mapped = zeta;
2675  // zeta_mapped = eta;
2676  // }
2677 
2678  // else if (elem->node_id(2) == min_node)
2679  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
2680  // {
2681  // // Case 3
2682  // eta_mapped = zeta;
2683  // zeta_mapped = -eta;
2684  // }
2685  // else
2686  // {
2687  // // Case 4
2688  // eta_mapped = -eta;
2689  // zeta_mapped = zeta;
2690  // }
2691 
2692  // else if (elem->node_id(6) == min_node)
2693  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
2694  // {
2695  // // Case 5
2696  // eta_mapped = -eta;
2697  // zeta_mapped = -zeta;
2698  // }
2699  // else
2700  // {
2701  // // Case 6
2702  // eta_mapped = -zeta;
2703  // zeta_mapped = -eta;
2704  // }
2705 
2706  // else if (elem->node_id(5) == min_node)
2707  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
2708  // {
2709  // // Case 7
2710  // eta_mapped = -zeta;
2711  // zeta_mapped = eta;
2712  // }
2713  // else
2714  // {
2715  // // Case 8
2716  // eta_mapped = eta;
2717  // zeta_mapped = -zeta;
2718  // }
2719  // }
2720 
2721 
2722  // // Face 3
2723  // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
2724  // {
2725  // const unsigned int min_node = std::min(elem->node_id(2),
2726  // std::min(elem->node_id(3),
2727  // std::min(elem->node_id(7),
2728  // elem->node_id(6))));
2729  // if (elem->node_id(3) == min_node)
2730  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
2731  // {
2732  // // Case 1
2733  // xi_mapped = xi;
2734  // zeta_mapped = zeta;
2735  // }
2736  // else
2737  // {
2738  // // Case 2
2739  // xi_mapped = zeta;
2740  // zeta_mapped = xi;
2741  // }
2742 
2743  // else if (elem->node_id(7) == min_node)
2744  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
2745  // {
2746  // // Case 3
2747  // xi_mapped = -zeta;
2748  // zeta_mapped = xi;
2749  // }
2750  // else
2751  // {
2752  // // Case 4
2753  // xi_mapped = xi;
2754  // zeta_mapped = -zeta;
2755  // }
2756 
2757  // else if (elem->node_id(6) == min_node)
2758  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
2759  // {
2760  // // Case 5
2761  // xi_mapped = -xi;
2762  // zeta_mapped = -zeta;
2763  // }
2764  // else
2765  // {
2766  // // Case 6
2767  // xi_mapped = -zeta;
2768  // zeta_mapped = -xi;
2769  // }
2770 
2771  // else if (elem->node_id(2) == min_node)
2772  // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
2773  // {
2774  // // Case 7
2775  // xi_mapped = zeta;
2776  // zeta_mapped = -xi;
2777  // }
2778  // else
2779  // {
2780  // // Case 8
2781  // xi_mapped = -xi;
2782  // zeta_mapped = zeta;
2783  // }
2784  // }
2785 
2786 
2787  // // Face 4
2788  // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
2789  // {
2790  // const unsigned int min_node = std::min(elem->node_id(3),
2791  // std::min(elem->node_id(0),
2792  // std::min(elem->node_id(4),
2793  // elem->node_id(7))));
2794  // if (elem->node_id(0) == min_node)
2795  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
2796  // {
2797  // // Case 1
2798  // eta_mapped = eta;
2799  // zeta_mapped = zeta;
2800  // }
2801  // else
2802  // {
2803  // // Case 2
2804  // eta_mapped = zeta;
2805  // zeta_mapped = eta;
2806  // }
2807 
2808  // else if (elem->node_id(4) == min_node)
2809  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
2810  // {
2811  // // Case 3
2812  // eta_mapped = -zeta;
2813  // zeta_mapped = eta;
2814  // }
2815  // else
2816  // {
2817  // // Case 4
2818  // eta_mapped = eta;
2819  // zeta_mapped = -zeta;
2820  // }
2821 
2822  // else if (elem->node_id(7) == min_node)
2823  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
2824  // {
2825  // // Case 5
2826  // eta_mapped = -eta;
2827  // zeta_mapped = -zeta;
2828  // }
2829  // else
2830  // {
2831  // // Case 6
2832  // eta_mapped = -zeta;
2833  // zeta_mapped = -eta;
2834  // }
2835 
2836  // else if (elem->node_id(3) == min_node)
2837  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
2838  // {
2839  // // Case 7
2840  // eta_mapped = zeta;
2841  // zeta_mapped = -eta;
2842  // }
2843  // else
2844  // {
2845  // // Case 8
2846  // eta_mapped = -eta;
2847  // zeta_mapped = zeta;
2848  // }
2849  // }
2850 
2851 
2852  // // Face 5
2853  // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
2854  // {
2855  // const unsigned int min_node = std::min(elem->node_id(4),
2856  // std::min(elem->node_id(5),
2857  // std::min(elem->node_id(6),
2858  // elem->node_id(7))));
2859  // if (elem->node_id(4) == min_node)
2860  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
2861  // {
2862  // // Case 1
2863  // xi_mapped = xi;
2864  // eta_mapped = eta;
2865  // }
2866  // else
2867  // {
2868  // // Case 2
2869  // xi_mapped = eta;
2870  // eta_mapped = xi;
2871  // }
2872 
2873  // else if (elem->node_id(5) == min_node)
2874  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
2875  // {
2876  // // Case 3
2877  // xi_mapped = eta;
2878  // eta_mapped = -xi;
2879  // }
2880  // else
2881  // {
2882  // // Case 4
2883  // xi_mapped = -xi;
2884  // eta_mapped = eta;
2885  // }
2886 
2887  // else if (elem->node_id(6) == min_node)
2888  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
2889  // {
2890  // // Case 5
2891  // xi_mapped = -xi;
2892  // eta_mapped = -eta;
2893  // }
2894  // else
2895  // {
2896  // // Case 6
2897  // xi_mapped = -eta;
2898  // eta_mapped = -xi;
2899  // }
2900 
2901  // else if (elem->node_id(7) == min_node)
2902  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
2903  // {
2904  // // Case 7
2905  // xi_mapped = -eta;
2906  // eta_mapped = xi;
2907  // }
2908  // else
2909  // {
2910  // // Case 8
2911  // xi_mapped = xi;
2912  // eta_mapped = eta;
2913  // }
2914  // }
2915 
2916 
2917  // }
2918 
2919 
2920 
2921  // libmesh_assert_less (j, 3);
2922 
2923  // switch (j)
2924  // {
2925  // // d()/dxi
2926  // case 0:
2927  // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
2928  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2929  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2930 
2931  // // d()/deta
2932  // case 1:
2933  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2934  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
2935  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2936 
2937  // // d()/dzeta
2938  // case 2:
2939  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2940  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2941  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
2942 
2943  // default:
2944  // libmesh_error_msg("Invalid derivative index j = " << j);
2945  // }
2946  // }
2947 
2948 
2949  default:
2950  libmesh_error_msg("Invalid element type = " << type);
2951  }
2952  }
2953 
2954 
2955  default:
2956  libmesh_error_msg("Invalid totalorder = " << totalorder);
2957  }
2958 
2959 #endif
2960 }
2961 
2962 
2963 
2964 template <>
2965 Real FE<3,BERNSTEIN>::shape_second_deriv(const ElemType,
2966  const Order,
2967  const unsigned int,
2968  const unsigned int,
2969  const Point &)
2970 {
2971  static bool warning_given = false;
2972 
2973  if (!warning_given)
2974  libMesh::err << "Second derivatives for Bernstein elements "
2975  << "are not yet implemented!"
2976  << std::endl;
2977 
2978  warning_given = true;
2979  return 0.;
2980 }
2981 
2982 
2983 
2984 template <>
2985 Real FE<3,BERNSTEIN>::shape_second_deriv(const Elem *,
2986  const Order,
2987  const unsigned int,
2988  const unsigned int,
2989  const Point &)
2990 {
2991  static bool warning_given = false;
2992 
2993  if (!warning_given)
2994  libMesh::err << "Second derivatives for Bernstein elements "
2995  << "are not yet implemented!"
2996  << std::endl;
2997 
2998  warning_given = true;
2999  return 0.;
3000 }
3001 
3002 } // namespace libMesh
3003 
3004 
3005 
3006 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2555
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::err
OStreamProxy err(std::cerr)
Definition: libmesh_common.h:229
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libmesh_config.h
libMesh::Elem::type
virtual ElemType type() const =0
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
generated by doxygen