19 #include "libmesh/libmesh_config.h" 21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES 24 #include "libmesh/fe.h" 25 #include "libmesh/elem.h" 26 #include "libmesh/utility.h" 27 #include "libmesh/enum_to_string.h" 49 const unsigned int totalorder,
52 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
55 if (i < 4 || i >= 4*totalorder)
58 const int edge = (i-4) / (totalorder-1);
59 libmesh_assert_less(edge, 4);
60 libmesh_assert_greater_equal(edge, 0);
62 const int edge_end = (edge+1)%4;
64 const int edge_i = i - 4 - edge*(totalorder-1);
69 ((elem->
point(edge) > elem->
point(edge_end)) ==
91 const bool add_p_level)
97 const Order totalorder =
static_cast<Order>(order + add_p_level * elem->
p_level());
117 const Real l1 = 1-p(0)-p(1);
118 const Real l2 = p(0);
119 const Real l3 = p(1);
121 libmesh_assert_less (i, 6);
129 case 3:
return l1*l2*(-4.*sqrt6);
130 case 4:
return l2*l3*(-4.*sqrt6);
131 case 5:
return l3*l1*(-4.*sqrt6);
134 libmesh_error_msg(
"Invalid i = " << i);
145 const Real xi = p(0);
146 const Real eta = p(1);
148 libmesh_assert_less (i, 9);
151 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
152 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
175 Real l1 = 1-p(0)-p(1);
181 libmesh_assert_less (i, 10);
184 if (i==4 && (elem->
point(0) > elem->
point(1)))f=-1;
185 if (i==6 && (elem->
point(1) > elem->
point(2)))f=-1;
186 if (i==8 && (elem->
point(2) > elem->
point(0)))f=-1;
197 case 3:
return l1*l2*(-4.*sqrt6);
198 case 4:
return f*l1*l2*(-4.*sqrt10)*(l2-l1);
200 case 5:
return l2*l3*(-4.*sqrt6);
201 case 6:
return f*l2*l3*(-4.*sqrt10)*(l3-l2);
203 case 7:
return l3*l1*(-4.*sqrt6);
204 case 8:
return f*l3*l1*(-4.*sqrt10)*(l1-l3);
207 case 9:
return l1*l2*l3;
210 libmesh_error_msg(
"Invalid i = " << i);
220 const Real xi = p(0);
221 const Real eta = p(1);
223 libmesh_assert_less (i, 16);
226 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
227 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
229 const Real f = quad_flip(elem, totalorder, i);
252 Real l1 = 1-p(0)-p(1);
258 libmesh_assert_less (i, 15);
261 if (i== 4 && (elem->
point(0) > elem->
point(1)))f=-1;
262 if (i== 7 && (elem->
point(1) > elem->
point(2)))f=-1;
263 if (i==10 && (elem->
point(2) > elem->
point(0)))f=-1;
274 case 3:
return l1*l2*(-4.*sqrt6);
275 case 4:
return f*l1*l2*(-4.*sqrt10)*(l2-l1);
276 case 5:
return l1*l2*(-sqrt14)*(5.*pow<2>(l2-l1)-1);
278 case 6:
return l2*l3*(-4.*sqrt6);
279 case 7:
return f*l2*l3*(-4.*sqrt10)*(l3-l2);
280 case 8:
return l2*l3*(-sqrt14)*(5.*pow<2>(l3-l2)-1);
282 case 9:
return l3*l1*(-4.*sqrt6);
283 case 10:
return f*l3*l1*(-4.*sqrt10)*(l1-l3);
284 case 11:
return l3*l1*(-sqrt14)*(5.*pow<2>(l1-l3)-1);
287 case 12:
return l1*l2*l3;
289 case 13:
return l1*l2*l3*(l2-l1);
290 case 14:
return l1*l2*l3*(2*l3-1);
293 libmesh_error_msg(
"Invalid i = " << i);
303 const Real xi = p(0);
304 const Real eta = p(1);
306 libmesh_assert_less (i, 25);
309 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
310 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
312 const Real f = quad_flip(elem, totalorder, i);
335 Real l1 = 1-p(0)-p(1);
340 const Real y=2.*l3-1;
344 libmesh_assert_less (i, 21);
347 if ((i== 4||i== 6) && (elem->
point(0) > elem->
point(1)))f=-1;
348 if ((i== 8||i==10) && (elem->
point(1) > elem->
point(2)))f=-1;
349 if ((i==12||i==14) && (elem->
point(2) > elem->
point(0)))f=-1;
360 case 3:
return l1*l2*(-4.*sqrt6);
361 case 4:
return f*l1*l2*(-4.*sqrt10)*(l2-l1);
362 case 5:
return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
363 case 6:
return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
365 case 7:
return l2*l3*(-4.*sqrt6);
366 case 8:
return f*l2*l3*(-4.*sqrt10)*(l3-l2);
367 case 9:
return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
368 case 10:
return -f*sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
370 case 11:
return l3*l1*(-4.*sqrt6);
371 case 12:
return f*l3*l1*(-4.*sqrt10)*(l1-l3);
372 case 13:
return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
373 case 14:
return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
376 case 15:
return l1*l2*l3;
378 case 16:
return l1*l2*l3*x;
379 case 17:
return l1*l2*l3*y;
381 case 18:
return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
382 case 19:
return l1*l2*l3*x*y;
383 case 20:
return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
386 libmesh_error_msg(
"Invalid i = " << i);
395 const Real xi = p(0);
396 const Real eta = p(1);
398 libmesh_assert_less (i, 36);
401 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
402 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
404 const Real f = quad_flip(elem, totalorder, i);
427 Real l1 = 1-p(0)-p(1);
432 const Real y=2.*l3-1;
436 libmesh_assert_less (i, 28);
439 if ((i== 4||i== 6) && (elem->
point(0) > elem->
point(1)))f=-1;
440 if ((i== 9||i==11) && (elem->
point(1) > elem->
point(2)))f=-1;
441 if ((i==14||i==16) && (elem->
point(2) > elem->
point(0)))f=-1;
452 case 3:
return l1*l2*(-4.*sqrt6);
453 case 4:
return f*l1*l2*(-4.*sqrt10)*(l2-l1);
454 case 5:
return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
455 case 6:
return f*(-sqrt2)*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
456 case 7:
return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
458 case 8:
return l2*l3*(-4.*sqrt6);
459 case 9:
return f*l2*l3*(-4.*sqrt10)*(l3-l2);
460 case 10:
return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
461 case 11:
return f*(-sqrt2)*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
462 case 12:
return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
464 case 13:
return l3*l1*(-4.*sqrt6);
465 case 14:
return f*l3*l1*(-4.*sqrt10)*(l1-l3);
466 case 15:
return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
467 case 16:
return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
468 case 17:
return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
473 case 18:
return l1*l2*l3;
475 case 19:
return l1*l2*l3*x;
476 case 20:
return l1*l2*l3*y;
478 case 21:
return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2);
479 case 22:
return l1*l2*l3*(l2-l1)*(2.0*l3-1.0);
480 case 23:
return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3);
481 case 24:
return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
482 case 25:
return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
483 case 26:
return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
484 case 27:
return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
488 libmesh_error_msg(
"Invalid i = " << i);
497 const Real xi = p(0);
498 const Real eta = p(1);
500 libmesh_assert_less (i, 49);
503 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
504 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
506 const Real f = quad_flip(elem, totalorder, i);
531 Real l1 = 1-p(0)-p(1);
536 const Real y=2.*l3-1.;
540 libmesh_assert_less (i, 36);
543 if ((i>= 4&&i<= 8) && (elem->
point(0) > elem->
point(1)))f=-1;
544 if ((i>=10&&i<=14) && (elem->
point(1) > elem->
point(2)))f=-1;
545 if ((i>=16&&i<=20) && (elem->
point(2) > elem->
point(0)))f=-1;
556 case 3:
return l1*l2*(-4.*sqrt6);
557 case 4:
return f*l1*l2*(-4.*sqrt10)*(l2-l1);
559 case 5:
return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
560 case 6:
return f*-sqrt2*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
561 case 7:
return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
562 case 8:
return f*-sqrt26*l1*l2*((-0.25E1+(15.0-0.165E2*l1*l1)*l1*l1)*l1+(0.25E1+(-45.0+0.825E2*l1*l1)*l1*l1+((45.0-0.165E3*l1*l1)*l1+(-15.0+0.165E3*l1*l1+(-0.825E2*l1+0.165E2*l2)*l2)*l2)*l2)*l2);
564 case 9:
return l2*l3*(-4.*sqrt6);
565 case 10:
return f*l2*l3*(-4.*sqrt10)*(l3-l2);
567 case 11:
return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
568 case 12:
return f*-sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
569 case 13:
return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
570 case 14:
return f*-sqrt26*l2*l3*((0.25E1+(-15.0+0.165E2*l3*l3)*l3*l3)*l3+(-0.25E1+(45.0-0.825E2*l3*l3)*l3*l3+((-45.0+0.165E3*l3*l3)*l3+(15.0-0.165E3*l3*l3+(0.825E2*l3-0.165E2*l2)*l2)*l2)*l2)*l2);
572 case 15:
return l3*l1*(-4.*sqrt6);
573 case 16:
return f*l3*l1*(-4.*sqrt10)*(l1-l3);
575 case 17:
return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
576 case 18:
return -f*sqrt2*l3*l1*((9.-21.*l3*l3)*l3+(-9.+63.*l3*l3+(-63.*l3+21.*l1)*l1)*l1);
577 case 19:
return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
578 case 20:
return f*-sqrt26*l3*l1*((-0.25E1+(15.0-0.165E2*l3*l3)*l3*l3)*l3+(0.25E1+(-45.0+0.825E2*l3*l3)*l3*l3+((45.0-0.165E3*l3*l3)*l3+(-15.0+0.165E3*l3*l3+(-0.825E2*l3+0.165E2*l1)*l1)*l1)*l1)*l1);
582 case 21:
return l1*l2*l3;
584 case 22:
return l1*l2*l3*x;
585 case 23:
return l1*l2*l3*y;
587 case 24:
return l1*l2*l3*0.5*(3.*pow<2>(x)-1.);
588 case 25:
return l1*l2*l3*x*y;
589 case 26:
return l1*l2*l3*0.5*(3.*pow<2>(y)-1.);
591 case 27:
return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
592 case 28:
return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
593 case 29:
return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
594 case 30:
return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
595 case 31:
return 0.125*l1*l2*l3*((-15.0+(70.0-63.0*l1*l1)*l1*l1)*l1+(15.0+(-210.0+315.0*l1*l1)*l1*l1+((210.0-630.0*l1*l1)*l1+(-70.0+630.0*l1*l1+(-315.0*l1+63.0*l2)*l2)*l2)*l2)*l2);
596 case 32:
return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2)*(2.0*l3-1.0);
597 case 33:
return 0.25*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0+(-12.0+12.0*l3)*l3);
598 case 34:
return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3)*(l2-l1);
599 case 35:
return 0.125*l1*l2*l3*(-8.0+(240.0+(-1680.0+(4480.0+(-5040.0+2016.0*l3)*l3)*l3)*l3)*l3);
602 libmesh_error_msg(
"Invalid i = " << i);
611 const Real xi = p(0);
612 const Real eta = p(1);
614 libmesh_assert_less (i, 64);
617 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
618 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
620 const Real f = quad_flip(elem, totalorder, i);
637 libmesh_error_msg(
"ERROR: Unsupported polynomial order!");
651 libmesh_error_msg(
"Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
660 const unsigned int i,
662 const bool add_p_level)
674 const unsigned int i,
675 const unsigned int j,
677 const bool add_p_level)
683 const Order totalorder =
static_cast<Order>(order + add_p_level * elem->
p_level());
710 const Real xi = p(0);
711 const Real eta = p(1);
713 libmesh_assert_less (i, 9);
716 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
717 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
732 libmesh_error_msg(
"Invalid j = " << j);
761 const Real xi = p(0);
762 const Real eta = p(1);
764 libmesh_assert_less (i, 16);
767 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
768 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
770 const Real f = quad_flip(elem, totalorder, i);
785 libmesh_error_msg(
"Invalid j = " << j);
816 const Real xi = p(0);
817 const Real eta = p(1);
819 libmesh_assert_less (i, 25);
822 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
823 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
825 const Real f = quad_flip(elem, totalorder, i);
840 libmesh_error_msg(
"Invalid j = " << j);
871 const Real xi = p(0);
872 const Real eta = p(1);
874 libmesh_assert_less (i, 36);
877 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
878 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
880 const Real f = quad_flip(elem, totalorder, i);
895 libmesh_error_msg(
"Invalid j = " << j);
924 const Real xi = p(0);
925 const Real eta = p(1);
927 libmesh_assert_less (i, 49);
930 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
931 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
933 const Real f = quad_flip(elem, totalorder, i);
948 libmesh_error_msg(
"Invalid j = " << j);
977 const Real xi = p(0);
978 const Real eta = p(1);
980 libmesh_assert_less (i, 64);
983 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
984 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
986 const Real f = quad_flip(elem, totalorder, i);
1001 libmesh_error_msg(
"Invalid j = " << j);
1014 libmesh_error_msg(
"ERROR: Unsupported polynomial order!");
1029 libmesh_error_msg(
"Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
1037 const unsigned int i,
1038 const unsigned int j,
1040 const bool add_p_level)
1048 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 1057 static bool warning_given =
false;
1060 libMesh::err <<
"Second derivatives for Szabab elements " 1061 <<
" are not yet implemented!" 1064 warning_given =
true;
1078 static bool warning_given =
false;
1081 libMesh::err <<
"Second derivatives for Szabab elements " 1082 <<
" are not yet implemented!" 1085 warning_given =
true;
1098 static bool warning_given =
false;
1101 libMesh::err <<
"Second derivatives for Szabab elements " 1102 <<
" are not yet implemented!" 1105 warning_given =
true;
1113 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
ElemType
Defines an enum for geometric element types.
Order
defines an enum for polynomial orders.
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
This is the base class from which all geometric element types are derived.
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
OrderWrapper order
The approximation order of the element.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
The libMesh namespace provides an interface to certain functionality in the library.
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't be...
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
A specific instantiation of the FEBase class.
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
const Point & point(const unsigned int i) const
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)