libMesh
fe_bernstein_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22 
23 // libmesh includes
24 #include "libmesh/fe.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/enum_to_string.h"
27 
28 
29 namespace {
30  using namespace libMesh;
31 
32  // Indices and coefficients for the HEX20
33  //
34  // The only way to make any sense of this
35  // is to look at the mgflo/mg2/mgf documentation
36  // and make the cut-out cube!
37  // 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
38  static const unsigned int hex20_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};
39  static const unsigned int hex20_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};
40  static const unsigned int hex20_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};
41  //To compute the hex20 shape functions the original shape functions for hex27 are used.
42  //hex20_scalx[i] tells how often the original x-th shape function has to be added to the original i-th shape function
43  //to compute the new i-th shape function for hex20
44  //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
45  // B_0^HEX20 = B_0^HEX27 + hex20_scal20[0]*B_20^HEX27 + hex20_scal21[0]*B_21^HEX27 + ...
46 
47  static const Real hex20_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};
48  static const Real hex20_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};
49  static const Real hex20_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};
50  static const Real hex20_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};
51  static const Real hex20_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};
52  static const Real hex20_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};
53  static const Real hex20_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};
54 
55  Point get_min_point(const Elem * elem,
56  unsigned int a, unsigned int b,
57  unsigned int c, unsigned int d)
58  {
59  return std::min(std::min(elem->point(a),elem->point(b)),
60  std::min(elem->point(c),elem->point(d)));
61  }
62 } // anonymous namespace
63 
64 
65 
66 namespace libMesh
67 {
68 
69 
70 LIBMESH_DEFAULT_VECTORIZED_FE(3,BERNSTEIN)
71 
72 
73 template <>
74 Real FE<3,BERNSTEIN>::shape(const Elem * elem,
75  const Order order,
76  const unsigned int i,
77  const Point & p,
78  const bool add_p_level)
79 {
80 
81 #if LIBMESH_DIM == 3
82 
83  libmesh_assert(elem);
84  const ElemType type = elem->type();
85 
86  const Order totalorder =
87  order + add_p_level*elem->p_level();
88 
89  auto hex_remap = [i, elem] (const Point & p_in,
90  const unsigned int * hex_i0,
91  const unsigned int * hex_i1,
92  const unsigned int * hex_i2) {
93  // Compute hex shape functions as a tensor-product
94  const Real xi = p_in(0);
95  const Real eta = p_in(1);
96  const Real zeta = p_in(2);
97  Point p_to_remap = p_in;
98  Real & xi_mapped = p_to_remap(0);
99  Real & eta_mapped = p_to_remap(1);
100  Real & zeta_mapped = p_to_remap(2);
101 
102  // handle the edge orientation
103  {
104  // Edge 0
105  if ((hex_i0[i] >= 2) && (hex_i1[i] == 0) && (hex_i2[i] == 0))
106  {
107  if (elem->positive_edge_orientation(0))
108  xi_mapped = -xi;
109  }
110  // Edge 1
111  else if ((hex_i0[i] == 1) && (hex_i1[i] >= 2) && (hex_i2[i] == 0))
112  {
113  if (elem->positive_edge_orientation(1))
114  eta_mapped = -eta;
115  }
116  // edge 2
117  else if ((hex_i0[i] >= 2) && (hex_i1[i] == 1) && (hex_i2[i] == 0))
118  {
119  if (!elem->positive_edge_orientation(2))
120  xi_mapped = -xi;
121  }
122  // edge 3
123  else if ((hex_i0[i] == 0) && (hex_i1[i] >= 2) && (hex_i2[i] == 0))
124  {
125  if (elem->positive_edge_orientation(3))
126  eta_mapped = -eta;
127  }
128  // edge 4
129  else if ((hex_i0[i] == 0) && (hex_i1[i] == 0) && (hex_i2[i] >=2 ))
130  {
131  if (elem->positive_edge_orientation(4))
132  zeta_mapped = -zeta;
133  }
134  // edge 5
135  else if ((hex_i0[i] == 1) && (hex_i1[i] == 0) && (hex_i2[i] >=2 ))
136  {
137  if (elem->positive_edge_orientation(5))
138  zeta_mapped = -zeta;
139  }
140  // edge 6
141  else if ((hex_i0[i] == 1) && (hex_i1[i] == 1) && (hex_i2[i] >=2 ))
142  {
143  if (elem->positive_edge_orientation(6))
144  zeta_mapped = -zeta;
145  }
146  // edge 7
147  else if ((hex_i0[i] == 0) && (hex_i1[i] == 1) && (hex_i2[i] >=2 ))
148  {
149  if (elem->positive_edge_orientation(7))
150  zeta_mapped = -zeta;
151  }
152  // edge 8
153  else if ((hex_i0[i] >=2 ) && (hex_i1[i] == 0) && (hex_i2[i] == 1))
154  {
155  if (elem->positive_edge_orientation(8))
156  xi_mapped = -xi;
157  }
158  // edge 9
159  else if ((hex_i0[i] == 1) && (hex_i1[i] >=2 ) && (hex_i2[i] == 1))
160  {
161  if (elem->positive_edge_orientation(9))
162  eta_mapped = -eta;
163  }
164  // edge 10
165  else if ((hex_i0[i] >=2 ) && (hex_i1[i] == 1) && (hex_i2[i] == 1))
166  {
167  if (!elem->positive_edge_orientation(10))
168  xi_mapped = -xi;
169  }
170  // edge 11
171  else if ((hex_i0[i] == 0) && (hex_i1[i] >=2 ) && (hex_i2[i] == 1))
172  {
173  if (elem->positive_edge_orientation(11))
174  eta_mapped = -eta;
175  }
176  }
177 
178 
179  // handle the face orientation
180  {
181  // face 0
182  if ((hex_i2[i] == 0) && (hex_i0[i] >= 2) && (hex_i1[i] >= 2))
183  {
184  const Point min_point = get_min_point(elem, 1, 2, 0, 3);
185 
186  if (elem->point(0) == min_point)
187  if (elem->positive_face_orientation(0))
188  {
189  // case 1
190  xi_mapped = xi;
191  eta_mapped = eta;
192  }
193  else
194  {
195  // case 2
196  xi_mapped = eta;
197  eta_mapped = xi;
198  }
199 
200  else if (elem->point(3) == min_point)
201  if (elem->positive_face_orientation(0))
202  {
203  // case 3
204  xi_mapped = -eta;
205  eta_mapped = xi;
206  }
207  else
208  {
209  // case 4
210  xi_mapped = xi;
211  eta_mapped = -eta;
212  }
213 
214  else if (elem->point(2) == min_point)
215  if (elem->positive_face_orientation(0))
216  {
217  // case 5
218  xi_mapped = -xi;
219  eta_mapped = -eta;
220  }
221  else
222  {
223  // case 6
224  xi_mapped = -eta;
225  eta_mapped = -xi;
226  }
227 
228  else if (elem->point(1) == min_point)
229  {
230  if (elem->positive_face_orientation(0))
231  {
232  // case 7
233  xi_mapped = eta;
234  eta_mapped = -xi;
235  }
236  else
237  {
238  // Case 8
239  xi_mapped = -xi;
240  eta_mapped = eta;
241  }
242  }
243  }
244 
245 
246  // Face 1
247  else if ((hex_i1[i] == 0) && (hex_i0[i] >= 2) && (hex_i2[i] >= 2))
248  {
249  const Point min_point = get_min_point(elem, 0, 1, 5, 4);
250 
251  if (elem->point(0) == min_point)
252  if (!elem->positive_face_orientation(1))
253  {
254  // Case 1
255  xi_mapped = xi;
256  zeta_mapped = zeta;
257  }
258  else
259  {
260  // Case 2
261  xi_mapped = zeta;
262  zeta_mapped = xi;
263  }
264 
265  else if (elem->point(1) == min_point)
266  if (!elem->positive_face_orientation(1))
267  {
268  // Case 3
269  xi_mapped = zeta;
270  zeta_mapped = -xi;
271  }
272  else
273  {
274  // Case 4
275  xi_mapped = -xi;
276  zeta_mapped = zeta;
277  }
278 
279  else if (elem->point(5) == min_point)
280  if (!elem->positive_face_orientation(1))
281  {
282  // Case 5
283  xi_mapped = -xi;
284  zeta_mapped = -zeta;
285  }
286  else
287  {
288  // Case 6
289  xi_mapped = -zeta;
290  zeta_mapped = -xi;
291  }
292 
293  else if (elem->point(4) == min_point)
294  {
295  if (!elem->positive_face_orientation(1))
296  {
297  // Case 7
298  xi_mapped = -xi;
299  zeta_mapped = zeta;
300  }
301  else
302  {
303  // Case 8
304  xi_mapped = xi;
305  zeta_mapped = -zeta;
306  }
307  }
308  }
309 
310 
311  // Face 2
312  else if ((hex_i0[i] == 1) && (hex_i1[i] >= 2) && (hex_i2[i] >= 2))
313  {
314  const Point min_point = get_min_point(elem, 1, 2, 6, 5);
315 
316  if (elem->point(1) == min_point)
317  if (!elem->positive_face_orientation(2))
318  {
319  // Case 1
320  eta_mapped = eta;
321  zeta_mapped = zeta;
322  }
323  else
324  {
325  // Case 2
326  eta_mapped = zeta;
327  zeta_mapped = eta;
328  }
329 
330  else if (elem->point(2) == min_point)
331  if (!elem->positive_face_orientation(2))
332  {
333  // Case 3
334  eta_mapped = zeta;
335  zeta_mapped = -eta;
336  }
337  else
338  {
339  // Case 4
340  eta_mapped = -eta;
341  zeta_mapped = zeta;
342  }
343 
344  else if (elem->point(6) == min_point)
345  if (!elem->positive_face_orientation(2))
346  {
347  // Case 5
348  eta_mapped = -eta;
349  zeta_mapped = -zeta;
350  }
351  else
352  {
353  // Case 6
354  eta_mapped = -zeta;
355  zeta_mapped = -eta;
356  }
357 
358  else if (elem->point(5) == min_point)
359  {
360  if (!elem->positive_face_orientation(2))
361  {
362  // Case 7
363  eta_mapped = -zeta;
364  zeta_mapped = eta;
365  }
366  else
367  {
368  // Case 8
369  eta_mapped = eta;
370  zeta_mapped = -zeta;
371  }
372  }
373  }
374 
375 
376  // Face 3
377  else if ((hex_i1[i] == 1) && (hex_i0[i] >= 2) && (hex_i2[i] >= 2))
378  {
379  const Point min_point = get_min_point(elem, 2, 3, 7, 6);
380 
381  if (elem->point(3) == min_point)
382  if (elem->positive_face_orientation(3))
383  {
384  // Case 1
385  xi_mapped = xi;
386  zeta_mapped = zeta;
387  }
388  else
389  {
390  // Case 2
391  xi_mapped = zeta;
392  zeta_mapped = xi;
393  }
394 
395  else if (elem->point(7) == min_point)
396  if (elem->positive_face_orientation(3))
397  {
398  // Case 3
399  xi_mapped = -zeta;
400  zeta_mapped = xi;
401  }
402  else
403  {
404  // Case 4
405  xi_mapped = xi;
406  zeta_mapped = -zeta;
407  }
408 
409  else if (elem->point(6) == min_point)
410  if (elem->positive_face_orientation(3))
411  {
412  // Case 5
413  xi_mapped = -xi;
414  zeta_mapped = -zeta;
415  }
416  else
417  {
418  // Case 6
419  xi_mapped = -zeta;
420  zeta_mapped = -xi;
421  }
422 
423  else if (elem->point(2) == min_point)
424  {
425  if (elem->positive_face_orientation(3))
426  {
427  // Case 7
428  xi_mapped = zeta;
429  zeta_mapped = -xi;
430  }
431  else
432  {
433  // Case 8
434  xi_mapped = -xi;
435  zeta_mapped = zeta;
436  }
437  }
438  }
439 
440 
441  // Face 4
442  else if ((hex_i0[i] == 0) && (hex_i1[i] >= 2) && (hex_i2[i] >= 2))
443  {
444  const Point min_point = get_min_point(elem, 3, 0, 4, 7);
445 
446  if (elem->point(0) == min_point)
447  if (elem->positive_face_orientation(4))
448  {
449  // Case 1
450  eta_mapped = eta;
451  zeta_mapped = zeta;
452  }
453  else
454  {
455  // Case 2
456  eta_mapped = zeta;
457  zeta_mapped = eta;
458  }
459 
460  else if (elem->point(4) == min_point)
461  if (elem->positive_face_orientation(4))
462  {
463  // Case 3
464  eta_mapped = -zeta;
465  zeta_mapped = eta;
466  }
467  else
468  {
469  // Case 4
470  eta_mapped = eta;
471  zeta_mapped = -zeta;
472  }
473 
474  else if (elem->point(7) == min_point)
475  if (elem->positive_face_orientation(4))
476  {
477  // Case 5
478  eta_mapped = -eta;
479  zeta_mapped = -zeta;
480  }
481  else
482  {
483  // Case 6
484  eta_mapped = -zeta;
485  zeta_mapped = -eta;
486  }
487 
488  else if (elem->point(3) == min_point)
489  {
490  if (elem->positive_face_orientation(4))
491  {
492  // Case 7
493  eta_mapped = zeta;
494  zeta_mapped = -eta;
495  }
496  else
497  {
498  // Case 8
499  eta_mapped = -eta;
500  zeta_mapped = zeta;
501  }
502  }
503  }
504 
505 
506  // Face 5
507  else if ((hex_i2[i] == 1) && (hex_i0[i] >= 2) && (hex_i1[i] >= 2))
508  {
509  const Point min_point = get_min_point(elem, 4, 5, 6, 7);
510 
511  if (elem->point(4) == min_point)
512  if (!elem->positive_face_orientation(5))
513  {
514  // Case 1
515  xi_mapped = xi;
516  eta_mapped = eta;
517  }
518  else
519  {
520  // Case 2
521  xi_mapped = eta;
522  eta_mapped = xi;
523  }
524 
525  else if (elem->point(5) == min_point)
526  if (!elem->positive_face_orientation(5))
527  {
528  // Case 3
529  xi_mapped = eta;
530  eta_mapped = -xi;
531  }
532  else
533  {
534  // Case 4
535  xi_mapped = -xi;
536  eta_mapped = eta;
537  }
538 
539  else if (elem->point(6) == min_point)
540  if (!elem->positive_face_orientation(5))
541  {
542  // Case 5
543  xi_mapped = -xi;
544  eta_mapped = -eta;
545  }
546  else
547  {
548  // Case 6
549  xi_mapped = -eta;
550  eta_mapped = -xi;
551  }
552 
553  else if (elem->point(7) == min_point)
554  {
555  if (!elem->positive_face_orientation(5))
556  {
557  // Case 7
558  xi_mapped = -eta;
559  eta_mapped = xi;
560  }
561  else
562  {
563  // Case 8
564  xi_mapped = xi;
565  eta_mapped = eta;
566  }
567  }
568  }
569  }
570 
571  return p_to_remap;
572  };
573 
574  switch (totalorder)
575  {
576  // 1st order Bernstein.
577  case FIRST:
578  {
579  switch (type)
580  {
581 
582  // Bernstein shape functions on the tetrahedron.
583  case TET4:
584  case TET10:
585  case TET14:
586  {
587  libmesh_assert_less (i, 4);
588 
589  // Area coordinates
590  const Real zeta1 = p(0);
591  const Real zeta2 = p(1);
592  const Real zeta3 = p(2);
593  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
594 
595  switch(i)
596  {
597  case 0: return zeta0;
598  case 1: return zeta1;
599  case 2: return zeta2;
600  case 3: return zeta3;
601 
602  default:
603  libmesh_error_msg("Invalid shape function index i = " << i);
604  }
605  }
606 
607  // Bernstein shape functions on the hexahedral.
608  case HEX8:
609  case HEX20:
610  case HEX27:
611  {
612  libmesh_assert_less (i, 8);
613 
614  // Compute hex shape functions as a tensor-product
615  const Real xi = p(0);
616  const Real eta = p(1);
617  const Real zeta = p(2);
618 
619  // The only way to make any sense of this
620  // is to look at the mgflo/mg2/mgf documentation
621  // and make the cut-out cube!
622  // 0 1 2 3 4 5 6 7
623  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
624  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
625  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
626 
627  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
628  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
629  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
630  }
631 
632 
633  default:
634  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
635  }
636  }
637 
638 
639 
640 
641  case SECOND:
642  {
643  switch (type)
644  {
645 
646  // Bernstein shape functions on the tetrahedron.
647  case TET10:
648  case TET14:
649  {
650  libmesh_assert_less (i, 10);
651 
652  // Area coordinates
653  const Real zeta1 = p(0);
654  const Real zeta2 = p(1);
655  const Real zeta3 = p(2);
656  const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
657 
658  switch(i)
659  {
660  case 0: return zeta0*zeta0;
661  case 1: return zeta1*zeta1;
662  case 2: return zeta2*zeta2;
663  case 3: return zeta3*zeta3;
664  case 4: return 2.*zeta0*zeta1;
665  case 5: return 2.*zeta1*zeta2;
666  case 6: return 2.*zeta0*zeta2;
667  case 7: return 2.*zeta3*zeta0;
668  case 8: return 2.*zeta1*zeta3;
669  case 9: return 2.*zeta2*zeta3;
670 
671  default:
672  libmesh_error_msg("Invalid shape function index i = " << i);
673  }
674  }
675 
676  // Bernstein shape functions on the 20-noded hexahedral.
677  case HEX20:
678  {
679  libmesh_assert_less (i, 20);
680 
681  // Compute hex shape functions as a tensor-product
682  const Real xi = p(0);
683  const Real eta = p(1);
684  const Real zeta = p(2);
685 
686  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[i], xi)*
687  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[i], eta)*
688  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[i], zeta)
689  +hex20_scal20[i]*
690  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[20], xi)*
691  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[20], eta)*
692  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[20], zeta)
693  +hex20_scal21[i]*
694  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[21], xi)*
695  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[21], eta)*
696  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[21], zeta)
697  +hex20_scal22[i]*
698  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[22], xi)*
699  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[22], eta)*
700  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[22], zeta)
701  +hex20_scal23[i]*
702  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[23], xi)*
703  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[23], eta)*
704  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[23], zeta)
705  +hex20_scal24[i]*
706  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[24], xi)*
707  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[24], eta)*
708  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[24], zeta)
709  +hex20_scal25[i]*
710  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[25], xi)*
711  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[25], eta)*
712  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[25], zeta)
713  +hex20_scal26[i]*
714  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i0[26], xi)*
715  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i1[26], eta)*
716  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, hex20_i2[26], zeta));
717  }
718 
719  // Bernstein shape functions on the hexahedral.
720  case HEX27:
721  {
722  libmesh_assert_less (i, 27);
723 
724  // Compute hex shape functions as a tensor-product
725  const Real xi = p(0);
726  const Real eta = p(1);
727  const Real zeta = p(2);
728 
729  // The only way to make any sense of this
730  // is to look at the mgflo/mg2/mgf documentation
731  // and make the cut-out cube!
732  // 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
733  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};
734  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};
735  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};
736 
737  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
738  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
739  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
740  }
741 
742 
743  default:
744  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
745  }
746 
747  }
748 
749 
750 
751  // 3rd-order Bernstein.
752  case THIRD:
753  {
754  switch (type)
755  {
756 
757  // // Bernstein shape functions on the tetrahedron.
758  // case TET10:
759  // {
760  // libmesh_assert_less (i, 20);
761 
762  // // Area coordinates
763  // const Real zeta1 = p(0);
764  // const Real zeta2 = p(1);
765  // const Real zeta3 = p(2);
766  // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
767 
768 
769  // unsigned int shape=i;
770 
771  // // handle the edge orientation
772 
773  // if ((i== 4||i== 5) && elem->node_id(0) > elem->node_id(1))shape= 9-i; //Edge 0
774  // if ((i== 6||i== 7) && elem->node_id(1) > elem->node_id(2))shape=13-i; //Edge 1
775  // if ((i== 8||i== 9) && elem->node_id(0) > elem->node_id(2))shape=17-i; //Edge 2
776  // if ((i==10||i==11) && elem->node_id(0) > elem->node_id(3))shape=21-i; //Edge 3
777  // if ((i==12||i==13) && elem->node_id(1) > elem->node_id(3))shape=25-i; //Edge 4
778  // if ((i==14||i==15) && elem->node_id(2) > elem->node_id(3))shape=29-i; //Edge 5
779 
780  // // No need to handle face orientation in 3rd order.
781 
782 
783  // switch(shape)
784  // {
785  // //point function
786  // case 0: return zeta0*zeta0*zeta0;
787  // case 1: return zeta1*zeta1*zeta1;
788  // case 2: return zeta2*zeta2*zeta2;
789  // case 3: return zeta3*zeta3*zeta3;
790 
791  // //edge functions
792  // case 4: return 3.*zeta0*zeta0*zeta1;
793  // case 5: return 3.*zeta1*zeta1*zeta0;
794 
795  // case 6: return 3.*zeta1*zeta1*zeta2;
796  // case 7: return 3.*zeta2*zeta2*zeta1;
797 
798  // case 8: return 3.*zeta0*zeta0*zeta2;
799  // case 9: return 3.*zeta2*zeta2*zeta0;
800 
801  // case 10: return 3.*zeta0*zeta0*zeta3;
802  // case 11: return 3.*zeta3*zeta3*zeta0;
803 
804  // case 12: return 3.*zeta1*zeta1*zeta3;
805  // case 13: return 3.*zeta3*zeta3*zeta1;
806 
807  // case 14: return 3.*zeta2*zeta2*zeta3;
808  // case 15: return 3.*zeta3*zeta3*zeta2;
809 
810  // //face functions
811  // case 16: return 6.*zeta0*zeta1*zeta2;
812  // case 17: return 6.*zeta0*zeta1*zeta3;
813  // case 18: return 6.*zeta1*zeta2*zeta3;
814  // case 19: return 6.*zeta2*zeta0*zeta3;
815 
816  // default:
817  // libmesh_error_msg("Invalid shape function index i = " << i);
818  // }
819  // }
820 
821 
822  // Bernstein shape functions on the hexahedral.
823  case HEX27:
824  {
825  libmesh_assert_less (i, 64);
826 
827  // The only way to make any sense of this
828  // is to look at the mgflo/mg2/mgf documentation
829  // and make the cut-out cube!
830  // 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
831  // 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
832  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};
833  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};
834  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};
835 
836  const Point p_mapped = hex_remap(p, i0, i1, i2);
837  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], p_mapped(0))*
838  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], p_mapped(1))*
839  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], p_mapped(2)));
840  }
841 
842  default:
843  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
844  }
845 
846  }//case THIRD
847 
848 
849  // 4th-order Bernstein.
850  case FOURTH:
851  {
852  switch (type)
853  {
854 
855  // Bernstein shape functions on the hexahedral.
856  case HEX27:
857  {
858  libmesh_assert_less (i, 125);
859 
860  // The only way to make any sense of this
861  // is to look at the mgflo/mg2/mgf documentation
862  // and make the cut-out cube!
863  // 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
864  // 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
865  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};
866  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};
867  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};
868 
869  const Point p_mapped = hex_remap(p, i0, i1, i2);
870  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], p_mapped(0))*
871  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], p_mapped(1))*
872  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], p_mapped(2)));
873  }
874 
875 
876  default:
877  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
878  }
879  }
880 
881 
882  default:
883  libmesh_error_msg("Invalid totalorder = " << totalorder);
884  }
885 #else // LIBMESH_DIM != 3
886  libmesh_ignore(elem, order, i, p, add_p_level);
887  libmesh_not_implemented();
888 #endif
889 }
890 
891 
892 
893 template <>
894 Real FE<3,BERNSTEIN>::shape(const ElemType,
895  const Order,
896  const unsigned int,
897  const Point &)
898 {
899  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
900  return 0.;
901 }
902 
903 
904 template <>
905 Real FE<3,BERNSTEIN>::shape(const FEType fet,
906  const Elem * elem,
907  const unsigned int i,
908  const Point & p,
909  const bool add_p_level)
910 {
911  return FE<3,BERNSTEIN>::shape(elem, fet.order, i, p, add_p_level);
912 }
913 
914 template <>
915 Real FE<3,BERNSTEIN>::shape_deriv(const Elem * elem,
916  const Order order,
917  const unsigned int i,
918  const unsigned int j,
919  const Point & p,
920  const bool add_p_level)
921 {
922 
923 #if LIBMESH_DIM == 3
924  libmesh_assert(elem);
925  const ElemType type = elem->type();
926 
927  const Order totalorder =
928  order + add_p_level*elem->p_level();
929 
930  libmesh_assert_less (j, 3);
931 
932  switch (totalorder)
933  {
934  // 1st order Bernstein.
935  case FIRST:
936  {
937  switch (type)
938  {
939  // Bernstein shape functions on the tetrahedron.
940  case TET4:
941  case TET10:
942  case TET14:
943  {
944  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
945  }
946 
947 
948  // Bernstein shape functions on the hexahedral.
949  case HEX8:
950  case HEX20:
951  case HEX27:
952  {
953  libmesh_assert_less (i, 8);
954 
955  // Compute hex shape functions as a tensor-product
956  const Real xi = p(0);
957  const Real eta = p(1);
958  const Real zeta = p(2);
959 
960  // The only way to make any sense of this
961  // is to look at the mgflo/mg2/mgf documentation
962  // and make the cut-out cube!
963  // 0 1 2 3 4 5 6 7
964  static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
965  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
966  static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
967 
968  switch (j)
969  {
970  // d()/dxi
971  case 0:
972  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
973  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
974  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
975 
976  // d()/deta
977  case 1:
978  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
979  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
980  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
981 
982  // d()/dzeta
983  case 2:
984  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
985  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
986  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
987 
988  default:
989  libmesh_error_msg("Invalid derivative index j = " << j);
990  }
991  }
992 
993  default:
994  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
995  }
996  }
997 
998 
999 
1000 
1001  case SECOND:
1002  {
1003  switch (type)
1004  {
1005  // Bernstein shape functions on the tetrahedron.
1006  case TET10:
1007  case TET14:
1008  {
1009  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
1010  }
1011 
1012  // Bernstein shape functions on the hexahedral.
1013  case HEX20:
1014  {
1015  libmesh_assert_less (i, 20);
1016 
1017  // Compute hex shape functions as a tensor-product
1018  const Real xi = p(0);
1019  const Real eta = p(1);
1020  const Real zeta = p(2);
1021 
1022  switch (j)
1023  {
1024  // d()/dxi
1025  case 0:
1026  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[i], 0, xi)*
1027  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[i], eta)*
1028  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[i], zeta)
1029  +hex20_scal20[i]*
1030  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[20], 0, xi)*
1031  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[20], eta)*
1032  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[20], zeta)
1033  +hex20_scal21[i]*
1034  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[21], 0, xi)*
1035  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[21], eta)*
1036  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[21], zeta)
1037  +hex20_scal22[i]*
1038  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[22], 0, xi)*
1039  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[22], eta)*
1040  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[22], zeta)
1041  +hex20_scal23[i]*
1042  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[23], 0, xi)*
1043  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[23], eta)*
1044  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[23], zeta)
1045  +hex20_scal24[i]*
1046  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[24], 0, xi)*
1047  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[24], eta)*
1048  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[24], zeta)
1049  +hex20_scal25[i]*
1050  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[25], 0, xi)*
1051  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[25], eta)*
1052  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[25], zeta)
1053  +hex20_scal26[i]*
1054  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i0[26], 0, xi)*
1055  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[26], eta)*
1056  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[26], zeta));
1057 
1058  // d()/deta
1059  case 1:
1060  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[i], xi)*
1061  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[i], 0, eta)*
1062  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[i], zeta)
1063  +hex20_scal20[i]*
1064  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[20], xi)*
1065  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[20], 0, eta)*
1066  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[20], zeta)
1067  +hex20_scal21[i]*
1068  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[21], xi)*
1069  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[21], 0, eta)*
1070  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[21], zeta)
1071  +hex20_scal22[i]*
1072  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[22], xi)*
1073  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[22], 0, eta)*
1074  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[22], zeta)
1075  +hex20_scal23[i]*
1076  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[23], xi)*
1077  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[23], 0, eta)*
1078  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[23], zeta)
1079  +hex20_scal24[i]*
1080  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[24], xi)*
1081  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[24], 0, eta)*
1082  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[24], zeta)
1083  +hex20_scal25[i]*
1084  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[25], xi)*
1085  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[25], 0, eta)*
1086  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[25], zeta)
1087  +hex20_scal26[i]*
1088  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[26], xi)*
1089  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i1[26], 0, eta)*
1090  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i2[26], zeta));
1091 
1092  // d()/dzeta
1093  case 2:
1094  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[i], xi)*
1095  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[i], eta)*
1096  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[i], 0, zeta)
1097  +hex20_scal20[i]*
1098  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[20], xi)*
1099  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[20], eta)*
1100  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[20], 0, zeta)
1101  +hex20_scal21[i]*
1102  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[21], xi)*
1103  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[21], eta)*
1104  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[21], 0, zeta)
1105  +hex20_scal22[i]*
1106  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[22], xi)*
1107  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[22], eta)*
1108  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[22], 0, zeta)
1109  +hex20_scal23[i]*
1110  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[23], xi)*
1111  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[23], eta)*
1112  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[23], 0, zeta)
1113  +hex20_scal24[i]*
1114  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[24], xi)*
1115  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[24], eta)*
1116  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[24], 0, zeta)
1117  +hex20_scal25[i]*
1118  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[25], xi)*
1119  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[25], eta)*
1120  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[25], 0, zeta)
1121  +hex20_scal26[i]*
1122  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i0[26], xi)*
1123  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, hex20_i1[26], eta)*
1124  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, hex20_i2[26], 0, zeta));
1125 
1126  default:
1127  libmesh_error_msg("Invalid derivative index j = " << j);
1128  }
1129  }
1130 
1131  // Bernstein shape functions on the hexahedral.
1132  case HEX27:
1133  {
1134  libmesh_assert_less (i, 27);
1135 
1136  // Compute hex shape functions as a tensor-product
1137  const Real xi = p(0);
1138  const Real eta = p(1);
1139  const Real zeta = p(2);
1140 
1141  // The only way to make any sense of this
1142  // is to look at the mgflo/mg2/mgf documentation
1143  // and make the cut-out cube!
1144  // 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
1145  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};
1146  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};
1147  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};
1148 
1149  switch (j)
1150  {
1151  // d()/dxi
1152  case 0:
1153  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1154  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1155  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1156 
1157  // d()/deta
1158  case 1:
1159  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1160  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)*
1161  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta));
1162 
1163  // d()/dzeta
1164  case 2:
1165  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1166  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)*
1167  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));
1168 
1169  default:
1170  libmesh_error_msg("Invalid derivative index j = " << j);
1171  }
1172  }
1173 
1174 
1175  default:
1176  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1177  }
1178  }
1179 
1180 
1181 
1182  // 3rd-order Bernstein.
1183  case THIRD:
1184  {
1185  switch (type)
1186  {
1187 
1188  // // Bernstein shape functions derivatives.
1189  // case TET10:
1190  // {
1191  // // I have been lazy here and am using finite differences
1192  // // to compute the derivatives!
1193  // const Real eps = 1.e-4;
1194 
1195  // libmesh_assert_less (i, 20);
1196  // libmesh_assert_less (j, 3);
1197 
1198  // switch (j)
1199  // {
1200  // // d()/dxi
1201  // case 0:
1202  // {
1203  // const Point pp(p(0)+eps, p(1), p(2));
1204  // const Point pm(p(0)-eps, p(1), p(2));
1205 
1206  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1207  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1208  // }
1209 
1210  // // d()/deta
1211  // case 1:
1212  // {
1213  // const Point pp(p(0), p(1)+eps, p(2));
1214  // const Point pm(p(0), p(1)-eps, p(2));
1215 
1216  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1217  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1218  // }
1219  // // d()/dzeta
1220  // case 2:
1221  // {
1222  // const Point pp(p(0), p(1), p(2)+eps);
1223  // const Point pm(p(0), p(1), p(2)-eps);
1224 
1225  // return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
1226  // FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1227  // }
1228  // default:
1229  // libmesh_error_msg("Invalid derivative index j = " << j);
1230  // }
1231 
1232 
1233  // }
1234 
1235 
1236  // Bernstein shape functions on the hexahedral.
1237  case HEX27:
1238  {
1239  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
1240  }
1241 
1242  // // Compute hex shape functions as a tensor-product
1243  // const Real xi = p(0);
1244  // const Real eta = p(1);
1245  // const Real zeta = p(2);
1246  // Real xi_mapped = p(0);
1247  // Real eta_mapped = p(1);
1248  // Real zeta_mapped = p(2);
1249 
1250  // // The only way to make any sense of this
1251  // // is to look at the mgflo/mg2/mgf documentation
1252  // // and make the cut-out cube!
1253  // // 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
1254  // // 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
1255  // 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};
1256  // 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};
1257  // 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};
1258 
1259 
1260 
1261  // // handle the edge orientation
1262  // {
1263  // // Edge 0
1264  // if ((i1[i] == 0) && (i2[i] == 0))
1265  // {
1266  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
1267  // xi_mapped = -xi;
1268  // }
1269  // // Edge 1
1270  // else if ((i0[i] == 1) && (i2[i] == 0))
1271  // {
1272  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
1273  // eta_mapped = -eta;
1274  // }
1275  // // Edge 2
1276  // else if ((i1[i] == 1) && (i2[i] == 0))
1277  // {
1278  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
1279  // xi_mapped = -xi;
1280  // }
1281  // // Edge 3
1282  // else if ((i0[i] == 0) && (i2[i] == 0))
1283  // {
1284  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
1285  // eta_mapped = -eta;
1286  // }
1287  // // Edge 4
1288  // else if ((i0[i] == 0) && (i1[i] == 0))
1289  // {
1290  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
1291  // zeta_mapped = -zeta;
1292  // }
1293  // // Edge 5
1294  // else if ((i0[i] == 1) && (i1[i] == 0))
1295  // {
1296  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
1297  // zeta_mapped = -zeta;
1298  // }
1299  // // Edge 6
1300  // else if ((i0[i] == 1) && (i1[i] == 1))
1301  // {
1302  // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
1303  // zeta_mapped = -zeta;
1304  // }
1305  // // Edge 7
1306  // else if ((i0[i] == 0) && (i1[i] == 1))
1307  // {
1308  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
1309  // zeta_mapped = -zeta;
1310  // }
1311  // // Edge 8
1312  // else if ((i1[i] == 0) && (i2[i] == 1))
1313  // {
1314  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
1315  // xi_mapped = -xi;
1316  // }
1317  // // Edge 9
1318  // else if ((i0[i] == 1) && (i2[i] == 1))
1319  // {
1320  // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
1321  // eta_mapped = -eta;
1322  // }
1323  // // Edge 10
1324  // else if ((i1[i] == 1) && (i2[i] == 1))
1325  // {
1326  // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
1327  // xi_mapped = -xi;
1328  // }
1329  // // Edge 11
1330  // else if ((i0[i] == 0) && (i2[i] == 1))
1331  // {
1332  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
1333  // eta_mapped = -eta;
1334  // }
1335  // }
1336 
1337 
1338  // // handle the face orientation
1339  // {
1340  // // Face 0
1341  // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
1342  // {
1343  // const unsigned int min_node = std::min(elem->node_id(1),
1344  // std::min(elem->node_id(2),
1345  // std::min(elem->node_id(0),
1346  // elem->node_id(3))));
1347  // if (elem->node_id(0) == min_node)
1348  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
1349  // {
1350  // // Case 1
1351  // xi_mapped = xi;
1352  // eta_mapped = eta;
1353  // }
1354  // else
1355  // {
1356  // // Case 2
1357  // xi_mapped = eta;
1358  // eta_mapped = xi;
1359  // }
1360 
1361  // else if (elem->node_id(3) == min_node)
1362  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
1363  // {
1364  // // Case 3
1365  // xi_mapped = -eta;
1366  // eta_mapped = xi;
1367  // }
1368  // else
1369  // {
1370  // // Case 4
1371  // xi_mapped = xi;
1372  // eta_mapped = -eta;
1373  // }
1374 
1375  // else if (elem->node_id(2) == min_node)
1376  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
1377  // {
1378  // // Case 5
1379  // xi_mapped = -xi;
1380  // eta_mapped = -eta;
1381  // }
1382  // else
1383  // {
1384  // // Case 6
1385  // xi_mapped = -eta;
1386  // eta_mapped = -xi;
1387  // }
1388 
1389  // else if (elem->node_id(1) == min_node)
1390  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
1391  // {
1392  // // Case 7
1393  // xi_mapped = eta;
1394  // eta_mapped = -xi;
1395  // }
1396  // else
1397  // {
1398  // // Case 8
1399  // xi_mapped = -xi;
1400  // eta_mapped = eta;
1401  // }
1402  // }
1403 
1404 
1405  // // Face 1
1406  // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
1407  // {
1408  // const unsigned int min_node = std::min(elem->node_id(0),
1409  // std::min(elem->node_id(1),
1410  // std::min(elem->node_id(5),
1411  // elem->node_id(4))));
1412  // if (elem->node_id(0) == min_node)
1413  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
1414  // {
1415  // // Case 1
1416  // xi_mapped = xi;
1417  // zeta_mapped = zeta;
1418  // }
1419  // else
1420  // {
1421  // // Case 2
1422  // xi_mapped = zeta;
1423  // zeta_mapped = xi;
1424  // }
1425 
1426  // else if (elem->node_id(1) == min_node)
1427  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
1428  // {
1429  // // Case 3
1430  // xi_mapped = zeta;
1431  // zeta_mapped = -xi;
1432  // }
1433  // else
1434  // {
1435  // // Case 4
1436  // xi_mapped = -xi;
1437  // zeta_mapped = zeta;
1438  // }
1439 
1440  // else if (elem->node_id(5) == min_node)
1441  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
1442  // {
1443  // // Case 5
1444  // xi_mapped = -xi;
1445  // zeta_mapped = -zeta;
1446  // }
1447  // else
1448  // {
1449  // // Case 6
1450  // xi_mapped = -zeta;
1451  // zeta_mapped = -xi;
1452  // }
1453 
1454  // else if (elem->node_id(4) == min_node)
1455  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
1456  // {
1457  // // Case 7
1458  // xi_mapped = -xi;
1459  // zeta_mapped = zeta;
1460  // }
1461  // else
1462  // {
1463  // // Case 8
1464  // xi_mapped = xi;
1465  // zeta_mapped = -zeta;
1466  // }
1467  // }
1468 
1469 
1470  // // Face 2
1471  // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
1472  // {
1473  // const unsigned int min_node = std::min(elem->node_id(1),
1474  // std::min(elem->node_id(2),
1475  // std::min(elem->node_id(6),
1476  // elem->node_id(5))));
1477  // if (elem->node_id(1) == min_node)
1478  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
1479  // {
1480  // // Case 1
1481  // eta_mapped = eta;
1482  // zeta_mapped = zeta;
1483  // }
1484  // else
1485  // {
1486  // // Case 2
1487  // eta_mapped = zeta;
1488  // zeta_mapped = eta;
1489  // }
1490 
1491  // else if (elem->node_id(2) == min_node)
1492  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
1493  // {
1494  // // Case 3
1495  // eta_mapped = zeta;
1496  // zeta_mapped = -eta;
1497  // }
1498  // else
1499  // {
1500  // // Case 4
1501  // eta_mapped = -eta;
1502  // zeta_mapped = zeta;
1503  // }
1504 
1505  // else if (elem->node_id(6) == min_node)
1506  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
1507  // {
1508  // // Case 5
1509  // eta_mapped = -eta;
1510  // zeta_mapped = -zeta;
1511  // }
1512  // else
1513  // {
1514  // // Case 6
1515  // eta_mapped = -zeta;
1516  // zeta_mapped = -eta;
1517  // }
1518 
1519  // else if (elem->node_id(5) == min_node)
1520  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
1521  // {
1522  // // Case 7
1523  // eta_mapped = -zeta;
1524  // zeta_mapped = eta;
1525  // }
1526  // else
1527  // {
1528  // // Case 8
1529  // eta_mapped = eta;
1530  // zeta_mapped = -zeta;
1531  // }
1532  // }
1533 
1534 
1535  // // Face 3
1536  // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
1537  // {
1538  // const unsigned int min_node = std::min(elem->node_id(2),
1539  // std::min(elem->node_id(3),
1540  // std::min(elem->node_id(7),
1541  // elem->node_id(6))));
1542  // if (elem->node_id(3) == min_node)
1543  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
1544  // {
1545  // // Case 1
1546  // xi_mapped = xi;
1547  // zeta_mapped = zeta;
1548  // }
1549  // else
1550  // {
1551  // // Case 2
1552  // xi_mapped = zeta;
1553  // zeta_mapped = xi;
1554  // }
1555 
1556  // else if (elem->node_id(7) == min_node)
1557  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
1558  // {
1559  // // Case 3
1560  // xi_mapped = -zeta;
1561  // zeta_mapped = xi;
1562  // }
1563  // else
1564  // {
1565  // // Case 4
1566  // xi_mapped = xi;
1567  // zeta_mapped = -zeta;
1568  // }
1569 
1570  // else if (elem->node_id(6) == min_node)
1571  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
1572  // {
1573  // // Case 5
1574  // xi_mapped = -xi;
1575  // zeta_mapped = -zeta;
1576  // }
1577  // else
1578  // {
1579  // // Case 6
1580  // xi_mapped = -zeta;
1581  // zeta_mapped = -xi;
1582  // }
1583 
1584  // else if (elem->node_id(2) == min_node)
1585  // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
1586  // {
1587  // // Case 7
1588  // xi_mapped = zeta;
1589  // zeta_mapped = -xi;
1590  // }
1591  // else
1592  // {
1593  // // Case 8
1594  // xi_mapped = -xi;
1595  // zeta_mapped = zeta;
1596  // }
1597  // }
1598 
1599 
1600  // // Face 4
1601  // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
1602  // {
1603  // const unsigned int min_node = std::min(elem->node_id(3),
1604  // std::min(elem->node_id(0),
1605  // std::min(elem->node_id(4),
1606  // elem->node_id(7))));
1607  // if (elem->node_id(0) == min_node)
1608  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
1609  // {
1610  // // Case 1
1611  // eta_mapped = eta;
1612  // zeta_mapped = zeta;
1613  // }
1614  // else
1615  // {
1616  // // Case 2
1617  // eta_mapped = zeta;
1618  // zeta_mapped = eta;
1619  // }
1620 
1621  // else if (elem->node_id(4) == min_node)
1622  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
1623  // {
1624  // // Case 3
1625  // eta_mapped = -zeta;
1626  // zeta_mapped = eta;
1627  // }
1628  // else
1629  // {
1630  // // Case 4
1631  // eta_mapped = eta;
1632  // zeta_mapped = -zeta;
1633  // }
1634 
1635  // else if (elem->node_id(7) == min_node)
1636  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
1637  // {
1638  // // Case 5
1639  // eta_mapped = -eta;
1640  // zeta_mapped = -zeta;
1641  // }
1642  // else
1643  // {
1644  // // Case 6
1645  // eta_mapped = -zeta;
1646  // zeta_mapped = -eta;
1647  // }
1648 
1649  // else if (elem->node_id(3) == min_node)
1650  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
1651  // {
1652  // // Case 7
1653  // eta_mapped = zeta;
1654  // zeta_mapped = -eta;
1655  // }
1656  // else
1657  // {
1658  // // Case 8
1659  // eta_mapped = -eta;
1660  // zeta_mapped = zeta;
1661  // }
1662  // }
1663 
1664 
1665  // // Face 5
1666  // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
1667  // {
1668  // const unsigned int min_node = std::min(elem->node_id(4),
1669  // std::min(elem->node_id(5),
1670  // std::min(elem->node_id(6),
1671  // elem->node_id(7))));
1672  // if (elem->node_id(4) == min_node)
1673  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
1674  // {
1675  // // Case 1
1676  // xi_mapped = xi;
1677  // eta_mapped = eta;
1678  // }
1679  // else
1680  // {
1681  // // Case 2
1682  // xi_mapped = eta;
1683  // eta_mapped = xi;
1684  // }
1685 
1686  // else if (elem->node_id(5) == min_node)
1687  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
1688  // {
1689  // // Case 3
1690  // xi_mapped = eta;
1691  // eta_mapped = -xi;
1692  // }
1693  // else
1694  // {
1695  // // Case 4
1696  // xi_mapped = -xi;
1697  // eta_mapped = eta;
1698  // }
1699 
1700  // else if (elem->node_id(6) == min_node)
1701  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
1702  // {
1703  // // Case 5
1704  // xi_mapped = -xi;
1705  // eta_mapped = -eta;
1706  // }
1707  // else
1708  // {
1709  // // Case 6
1710  // xi_mapped = -eta;
1711  // eta_mapped = -xi;
1712  // }
1713 
1714  // else if (elem->node_id(7) == min_node)
1715  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
1716  // {
1717  // // Case 7
1718  // xi_mapped = -eta;
1719  // eta_mapped = xi;
1720  // }
1721  // else
1722  // {
1723  // // Case 8
1724  // xi_mapped = xi;
1725  // eta_mapped = eta;
1726  // }
1727  // }
1728 
1729 
1730  // }
1731 
1732 
1733 
1734  // libmesh_assert_less (j, 3);
1735 
1736  // switch (j)
1737  // {
1738  // // d()/dxi
1739  // case 0:
1740  // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
1741  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
1742  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
1743 
1744  // // d()/deta
1745  // case 1:
1746  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
1747  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
1748  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
1749 
1750  // // d()/dzeta
1751  // case 2:
1752  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
1753  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
1754  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
1755 
1756  // default:
1757  // libmesh_error_msg("Invalid derivative index j = " << j);
1758  // }
1759 
1760 
1761  default:
1762  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1763  }
1764  }
1765 
1766  // 4th-order Bernstein.
1767  case FOURTH:
1768  {
1769  switch (type)
1770  {
1771 
1772  // Bernstein shape functions derivatives on the hexahedral.
1773  case HEX27:
1774  {
1775  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,BERNSTEIN>::shape);
1776  }
1777 
1778  // // Compute hex shape functions as a tensor-product
1779  // const Real xi = p(0);
1780  // const Real eta = p(1);
1781  // const Real zeta = p(2);
1782  // Real xi_mapped = p(0);
1783  // Real eta_mapped = p(1);
1784  // Real zeta_mapped = p(2);
1785 
1786  // // The only way to make any sense of this
1787  // // is to look at the mgflo/mg2/mgf documentation
1788  // // and make the cut-out cube!
1789  // // 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
1790  // // 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
1791  // 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};
1792  // 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};
1793  // 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};
1794 
1795 
1796 
1797  // // handle the edge orientation
1798  // {
1799  // // Edge 0
1800  // if ((i1[i] == 0) && (i2[i] == 0))
1801  // {
1802  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(1)))
1803  // xi_mapped = -xi;
1804  // }
1805  // // Edge 1
1806  // else if ((i0[i] == 1) && (i2[i] == 0))
1807  // {
1808  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(2)))
1809  // eta_mapped = -eta;
1810  // }
1811  // // Edge 2
1812  // else if ((i1[i] == 1) && (i2[i] == 0))
1813  // {
1814  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(2)))
1815  // xi_mapped = -xi;
1816  // }
1817  // // Edge 3
1818  // else if ((i0[i] == 0) && (i2[i] == 0))
1819  // {
1820  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(3)))
1821  // eta_mapped = -eta;
1822  // }
1823  // // Edge 4
1824  // else if ((i0[i] == 0) && (i1[i] == 0))
1825  // {
1826  // if (elem->node_id(0) != std::min(elem->node_id(0), elem->node_id(4)))
1827  // zeta_mapped = -zeta;
1828  // }
1829  // // Edge 5
1830  // else if ((i0[i] == 1) && (i1[i] == 0))
1831  // {
1832  // if (elem->node_id(1) != std::min(elem->node_id(1), elem->node_id(5)))
1833  // zeta_mapped = -zeta;
1834  // }
1835  // // Edge 6
1836  // else if ((i0[i] == 1) && (i1[i] == 1))
1837  // {
1838  // if (elem->node_id(2) != std::min(elem->node_id(2), elem->node_id(6)))
1839  // zeta_mapped = -zeta;
1840  // }
1841  // // Edge 7
1842  // else if ((i0[i] == 0) && (i1[i] == 1))
1843  // {
1844  // if (elem->node_id(3) != std::min(elem->node_id(3), elem->node_id(7)))
1845  // zeta_mapped = -zeta;
1846  // }
1847  // // Edge 8
1848  // else if ((i1[i] == 0) && (i2[i] == 1))
1849  // {
1850  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(5)))
1851  // xi_mapped = -xi;
1852  // }
1853  // // Edge 9
1854  // else if ((i0[i] == 1) && (i2[i] == 1))
1855  // {
1856  // if (elem->node_id(5) != std::min(elem->node_id(5), elem->node_id(6)))
1857  // eta_mapped = -eta;
1858  // }
1859  // // Edge 10
1860  // else if ((i1[i] == 1) && (i2[i] == 1))
1861  // {
1862  // if (elem->node_id(7) != std::min(elem->node_id(7), elem->node_id(6)))
1863  // xi_mapped = -xi;
1864  // }
1865  // // Edge 11
1866  // else if ((i0[i] == 0) && (i2[i] == 1))
1867  // {
1868  // if (elem->node_id(4) != std::min(elem->node_id(4), elem->node_id(7)))
1869  // eta_mapped = -eta;
1870  // }
1871  // }
1872 
1873 
1874  // // handle the face orientation
1875  // {
1876  // // Face 0
1877  // if ((i2[i] == 0) && (i0[i] >= 2) && (i1[i] >= 2))
1878  // {
1879  // const unsigned int min_node = std::min(elem->node_id(1),
1880  // std::min(elem->node_id(2),
1881  // std::min(elem->node_id(0),
1882  // elem->node_id(3))));
1883  // if (elem->node_id(0) == min_node)
1884  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(3)))
1885  // {
1886  // // Case 1
1887  // xi_mapped = xi;
1888  // eta_mapped = eta;
1889  // }
1890  // else
1891  // {
1892  // // Case 2
1893  // xi_mapped = eta;
1894  // eta_mapped = xi;
1895  // }
1896 
1897  // else if (elem->node_id(3) == min_node)
1898  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(2)))
1899  // {
1900  // // Case 3
1901  // xi_mapped = -eta;
1902  // eta_mapped = xi;
1903  // }
1904  // else
1905  // {
1906  // // Case 4
1907  // xi_mapped = xi;
1908  // eta_mapped = -eta;
1909  // }
1910 
1911  // else if (elem->node_id(2) == min_node)
1912  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(1)))
1913  // {
1914  // // Case 5
1915  // xi_mapped = -xi;
1916  // eta_mapped = -eta;
1917  // }
1918  // else
1919  // {
1920  // // Case 6
1921  // xi_mapped = -eta;
1922  // eta_mapped = -xi;
1923  // }
1924 
1925  // else if (elem->node_id(1) == min_node)
1926  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(0)))
1927  // {
1928  // // Case 7
1929  // xi_mapped = eta;
1930  // eta_mapped = -xi;
1931  // }
1932  // else
1933  // {
1934  // // Case 8
1935  // xi_mapped = -xi;
1936  // eta_mapped = eta;
1937  // }
1938  // }
1939 
1940 
1941  // // Face 1
1942  // else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
1943  // {
1944  // const unsigned int min_node = std::min(elem->node_id(0),
1945  // std::min(elem->node_id(1),
1946  // std::min(elem->node_id(5),
1947  // elem->node_id(4))));
1948  // if (elem->node_id(0) == min_node)
1949  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(4)))
1950  // {
1951  // // Case 1
1952  // xi_mapped = xi;
1953  // zeta_mapped = zeta;
1954  // }
1955  // else
1956  // {
1957  // // Case 2
1958  // xi_mapped = zeta;
1959  // zeta_mapped = xi;
1960  // }
1961 
1962  // else if (elem->node_id(1) == min_node)
1963  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(0)))
1964  // {
1965  // // Case 3
1966  // xi_mapped = zeta;
1967  // zeta_mapped = -xi;
1968  // }
1969  // else
1970  // {
1971  // // Case 4
1972  // xi_mapped = -xi;
1973  // zeta_mapped = zeta;
1974  // }
1975 
1976  // else if (elem->node_id(5) == min_node)
1977  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(1)))
1978  // {
1979  // // Case 5
1980  // xi_mapped = -xi;
1981  // zeta_mapped = -zeta;
1982  // }
1983  // else
1984  // {
1985  // // Case 6
1986  // xi_mapped = -zeta;
1987  // zeta_mapped = -xi;
1988  // }
1989 
1990  // else if (elem->node_id(4) == min_node)
1991  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(5)))
1992  // {
1993  // // Case 7
1994  // xi_mapped = -xi;
1995  // zeta_mapped = zeta;
1996  // }
1997  // else
1998  // {
1999  // // Case 8
2000  // xi_mapped = xi;
2001  // zeta_mapped = -zeta;
2002  // }
2003  // }
2004 
2005 
2006  // // Face 2
2007  // else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
2008  // {
2009  // const unsigned int min_node = std::min(elem->node_id(1),
2010  // std::min(elem->node_id(2),
2011  // std::min(elem->node_id(6),
2012  // elem->node_id(5))));
2013  // if (elem->node_id(1) == min_node)
2014  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(5)))
2015  // {
2016  // // Case 1
2017  // eta_mapped = eta;
2018  // zeta_mapped = zeta;
2019  // }
2020  // else
2021  // {
2022  // // Case 2
2023  // eta_mapped = zeta;
2024  // zeta_mapped = eta;
2025  // }
2026 
2027  // else if (elem->node_id(2) == min_node)
2028  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(1)))
2029  // {
2030  // // Case 3
2031  // eta_mapped = zeta;
2032  // zeta_mapped = -eta;
2033  // }
2034  // else
2035  // {
2036  // // Case 4
2037  // eta_mapped = -eta;
2038  // zeta_mapped = zeta;
2039  // }
2040 
2041  // else if (elem->node_id(6) == min_node)
2042  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(2)))
2043  // {
2044  // // Case 5
2045  // eta_mapped = -eta;
2046  // zeta_mapped = -zeta;
2047  // }
2048  // else
2049  // {
2050  // // Case 6
2051  // eta_mapped = -zeta;
2052  // zeta_mapped = -eta;
2053  // }
2054 
2055  // else if (elem->node_id(5) == min_node)
2056  // if (elem->node_id(1) == std::min(elem->node_id(1), elem->node_id(6)))
2057  // {
2058  // // Case 7
2059  // eta_mapped = -zeta;
2060  // zeta_mapped = eta;
2061  // }
2062  // else
2063  // {
2064  // // Case 8
2065  // eta_mapped = eta;
2066  // zeta_mapped = -zeta;
2067  // }
2068  // }
2069 
2070 
2071  // // Face 3
2072  // else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
2073  // {
2074  // const unsigned int min_node = std::min(elem->node_id(2),
2075  // std::min(elem->node_id(3),
2076  // std::min(elem->node_id(7),
2077  // elem->node_id(6))));
2078  // if (elem->node_id(3) == min_node)
2079  // if (elem->node_id(2) == std::min(elem->node_id(2), elem->node_id(7)))
2080  // {
2081  // // Case 1
2082  // xi_mapped = xi;
2083  // zeta_mapped = zeta;
2084  // }
2085  // else
2086  // {
2087  // // Case 2
2088  // xi_mapped = zeta;
2089  // zeta_mapped = xi;
2090  // }
2091 
2092  // else if (elem->node_id(7) == min_node)
2093  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(6)))
2094  // {
2095  // // Case 3
2096  // xi_mapped = -zeta;
2097  // zeta_mapped = xi;
2098  // }
2099  // else
2100  // {
2101  // // Case 4
2102  // xi_mapped = xi;
2103  // zeta_mapped = -zeta;
2104  // }
2105 
2106  // else if (elem->node_id(6) == min_node)
2107  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(2)))
2108  // {
2109  // // Case 5
2110  // xi_mapped = -xi;
2111  // zeta_mapped = -zeta;
2112  // }
2113  // else
2114  // {
2115  // // Case 6
2116  // xi_mapped = -zeta;
2117  // zeta_mapped = -xi;
2118  // }
2119 
2120  // else if (elem->node_id(2) == min_node)
2121  // if (elem->node_id(6) == std::min(elem->node_id(3), elem->node_id(6)))
2122  // {
2123  // // Case 7
2124  // xi_mapped = zeta;
2125  // zeta_mapped = -xi;
2126  // }
2127  // else
2128  // {
2129  // // Case 8
2130  // xi_mapped = -xi;
2131  // zeta_mapped = zeta;
2132  // }
2133  // }
2134 
2135 
2136  // // Face 4
2137  // else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
2138  // {
2139  // const unsigned int min_node = std::min(elem->node_id(3),
2140  // std::min(elem->node_id(0),
2141  // std::min(elem->node_id(4),
2142  // elem->node_id(7))));
2143  // if (elem->node_id(0) == min_node)
2144  // if (elem->node_id(3) == std::min(elem->node_id(3), elem->node_id(4)))
2145  // {
2146  // // Case 1
2147  // eta_mapped = eta;
2148  // zeta_mapped = zeta;
2149  // }
2150  // else
2151  // {
2152  // // Case 2
2153  // eta_mapped = zeta;
2154  // zeta_mapped = eta;
2155  // }
2156 
2157  // else if (elem->node_id(4) == min_node)
2158  // if (elem->node_id(0) == std::min(elem->node_id(0), elem->node_id(7)))
2159  // {
2160  // // Case 3
2161  // eta_mapped = -zeta;
2162  // zeta_mapped = eta;
2163  // }
2164  // else
2165  // {
2166  // // Case 4
2167  // eta_mapped = eta;
2168  // zeta_mapped = -zeta;
2169  // }
2170 
2171  // else if (elem->node_id(7) == min_node)
2172  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(3)))
2173  // {
2174  // // Case 5
2175  // eta_mapped = -eta;
2176  // zeta_mapped = -zeta;
2177  // }
2178  // else
2179  // {
2180  // // Case 6
2181  // eta_mapped = -zeta;
2182  // zeta_mapped = -eta;
2183  // }
2184 
2185  // else if (elem->node_id(3) == min_node)
2186  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(0)))
2187  // {
2188  // // Case 7
2189  // eta_mapped = zeta;
2190  // zeta_mapped = -eta;
2191  // }
2192  // else
2193  // {
2194  // // Case 8
2195  // eta_mapped = -eta;
2196  // zeta_mapped = zeta;
2197  // }
2198  // }
2199 
2200 
2201  // // Face 5
2202  // else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
2203  // {
2204  // const unsigned int min_node = std::min(elem->node_id(4),
2205  // std::min(elem->node_id(5),
2206  // std::min(elem->node_id(6),
2207  // elem->node_id(7))));
2208  // if (elem->node_id(4) == min_node)
2209  // if (elem->node_id(5) == std::min(elem->node_id(5), elem->node_id(7)))
2210  // {
2211  // // Case 1
2212  // xi_mapped = xi;
2213  // eta_mapped = eta;
2214  // }
2215  // else
2216  // {
2217  // // Case 2
2218  // xi_mapped = eta;
2219  // eta_mapped = xi;
2220  // }
2221 
2222  // else if (elem->node_id(5) == min_node)
2223  // if (elem->node_id(6) == std::min(elem->node_id(6), elem->node_id(4)))
2224  // {
2225  // // Case 3
2226  // xi_mapped = eta;
2227  // eta_mapped = -xi;
2228  // }
2229  // else
2230  // {
2231  // // Case 4
2232  // xi_mapped = -xi;
2233  // eta_mapped = eta;
2234  // }
2235 
2236  // else if (elem->node_id(6) == min_node)
2237  // if (elem->node_id(7) == std::min(elem->node_id(7), elem->node_id(5)))
2238  // {
2239  // // Case 5
2240  // xi_mapped = -xi;
2241  // eta_mapped = -eta;
2242  // }
2243  // else
2244  // {
2245  // // Case 6
2246  // xi_mapped = -eta;
2247  // eta_mapped = -xi;
2248  // }
2249 
2250  // else if (elem->node_id(7) == min_node)
2251  // if (elem->node_id(4) == std::min(elem->node_id(4), elem->node_id(6)))
2252  // {
2253  // // Case 7
2254  // xi_mapped = -eta;
2255  // eta_mapped = xi;
2256  // }
2257  // else
2258  // {
2259  // // Case 8
2260  // xi_mapped = xi;
2261  // eta_mapped = eta;
2262  // }
2263  // }
2264 
2265 
2266  // }
2267 
2268 
2269 
2270  // libmesh_assert_less (j, 3);
2271 
2272  // switch (j)
2273  // {
2274  // // d()/dxi
2275  // case 0:
2276  // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi_mapped)*
2277  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2278  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2279 
2280  // // d()/deta
2281  // case 1:
2282  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2283  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta_mapped)*
2284  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i2[i], zeta_mapped));
2285 
2286  // // d()/dzeta
2287  // case 2:
2288  // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi_mapped)*
2289  // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta_mapped)*
2290  // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta_mapped));
2291 
2292  // default:
2293  // libmesh_error_msg("Invalid derivative index j = " << j);
2294  // }
2295 
2296 
2297  default:
2298  libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
2299  }
2300  }
2301 
2302 
2303  default:
2304  libmesh_error_msg("Invalid totalorder = " << totalorder);
2305  }
2306 
2307 #else // LIBMESH_DIM != 3
2308  libmesh_ignore(elem, order, i, j, p, add_p_level);
2309  libmesh_not_implemented();
2310 #endif
2311 }
2312 
2313 
2314 template <>
2315 Real FE<3,BERNSTEIN>::shape_deriv(const ElemType,
2316  const Order,
2317  const unsigned int,
2318  const unsigned int,
2319  const Point & )
2320 {
2321  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
2322  return 0.;
2323 }
2324 
2325 template <>
2326 Real FE<3,BERNSTEIN>::shape_deriv(const FEType fet,
2327  const Elem * elem,
2328  const unsigned int i,
2329  const unsigned int j,
2330  const Point & p,
2331  const bool add_p_level)
2332 {
2333  return FE<3,BERNSTEIN>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
2334 }
2335 
2336 
2337 
2338 
2339 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2340 
2341 template <>
2342 Real FE<3,BERNSTEIN>::shape_second_deriv(const Elem * elem,
2343  const Order order,
2344  const unsigned int i,
2345  const unsigned int j,
2346  const Point & p,
2347  const bool add_p_level)
2348 {
2349  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
2350  FE<3,BERNSTEIN>::shape_deriv);
2351 }
2352 
2353 
2354 template <>
2355 Real FE<3,BERNSTEIN>::shape_second_deriv(const ElemType,
2356  const Order,
2357  const unsigned int,
2358  const unsigned int,
2359  const Point &)
2360 {
2361  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge and face orientation is needed.");
2362  return 0.;
2363 }
2364 
2365 
2366 template <>
2367 Real FE<3,BERNSTEIN>::shape_second_deriv(const FEType fet,
2368  const Elem * elem,
2369  const unsigned int i,
2370  const unsigned int j,
2371  const Point & p,
2372  const bool add_p_level)
2373 {
2374  return FE<3,BERNSTEIN>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
2375 }
2376 
2377 
2378 
2379 
2380 #endif
2381 
2382 } // namespace libMesh
2383 
2384 
2385 
2386 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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:3109
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::fe_fdm_deriv
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven&#39;t be...
Definition: fe.C:911
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition: fe_bernstein_shape_0D.C:30
libMesh::FE
A specific instantiation of the FEBase class.
Definition: fe.h:127
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition: elem.C:3604
libMesh::fe_fdm_second_deriv
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
Definition: fe.C:969
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::positive_face_orientation
bool positive_face_orientation(const unsigned int i) const
Definition: elem.C:3613
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2454
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 on Thu Oct 30 2025 20:01:50 for libMesh by   doxygen 1.8.14