20 #include "libmesh/libmesh_config.h"
21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
23 #include "libmesh/fe.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/number_lookups.h"
26 #include "libmesh/utility.h"
39 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
50 const bool add_p_level)
56 const Order totalorder =
57 static_cast<Order>(order + add_p_level * elem->
p_level());
72 const Real eta = p(1);
74 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
95 else if (i < totalorder + 3u)
96 { i0 = i - 2; i1 = 0; }
97 else if (i < 2u*totalorder + 2)
98 { i0 = 1; i1 = i - totalorder - 1; }
99 else if (i < 3u*totalorder + 1)
100 { i0 = i - 2u*totalorder; i1 = 1; }
101 else if (i < 4u*totalorder)
102 { i0 = 0; i1 = i - 3u*totalorder + 1; }
106 unsigned int basisnum = i - 4*totalorder;
114 if ((i>= 4 && i<= 4+ totalorder-2u) && elem->
point(0) > elem->
point(1)) i0=totalorder+2-i0;
115 else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->
point(1) > elem->
point(2)) i1=totalorder+2-i1;
116 else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->
point(3) > elem->
point(2)) i0=totalorder+2-i0;
117 else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->
point(0) > elem->
point(3)) i1=totalorder+2-i1;
127 libmesh_assert_less (totalorder, 3);
129 const Real xi = p(0);
130 const Real eta = p(1);
132 libmesh_assert_less (i, 8);
135 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
136 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
137 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
150 libmesh_assert_less (totalorder, 2);
151 libmesh_fallthrough();
161 libmesh_assert_less (i, 3);
170 libmesh_error_msg(
"Invalid shape function index i = " << i);
179 libmesh_assert_less (i, 6);
187 case 3:
return 2.*x*r;
188 case 4:
return 2.*x*y;
189 case 5:
return 2.*r*y;
192 libmesh_error_msg(
"Invalid shape function index i = " << i);
200 libmesh_assert_less (i, 10);
202 unsigned int shape=i;
205 if ((i==3||i==4) && elem->
point(0) > elem->
point(1)) shape=7-i;
206 if ((i==5||i==6) && elem->
point(1) > elem->
point(2)) shape=11-i;
207 if ((i==7||i==8) && elem->
point(0) > elem->
point(2)) shape=15-i;
211 case 0:
return r*r*r;
212 case 1:
return x*x*x;
213 case 2:
return y*y*y;
215 case 3:
return 3.*x*r*r;
216 case 4:
return 3.*x*x*r;
218 case 5:
return 3.*y*x*x;
219 case 6:
return 3.*y*y*x;
221 case 7:
return 3.*y*r*r;
222 case 8:
return 3.*y*y*r;
224 case 9:
return 6.*x*y*r;
227 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
235 unsigned int shape=i;
237 libmesh_assert_less (i, 15);
239 if ((i==3||i== 5) && elem->
point(0) > elem->
point(1)) shape=8-i;
240 if ((i==6||i== 8) && elem->
point(1) > elem->
point(2)) shape=14-i;
241 if ((i==9||i==11) && elem->
point(0) > elem->
point(2)) shape=20-i;
247 case 0:
return r*r*r*r;
248 case 1:
return x*x*x*x;
249 case 2:
return y*y*y*y;
252 case 3:
return 4.*x*r*r*r;
253 case 4:
return 6.*x*x*r*r;
254 case 5:
return 4.*x*x*x*r;
256 case 6:
return 4.*y*x*x*x;
257 case 7:
return 6.*y*y*x*x;
258 case 8:
return 4.*y*y*y*x;
260 case 9:
return 4.*y*r*r*r;
261 case 10:
return 6.*y*y*r*r;
262 case 11:
return 4.*y*y*y*r;
265 case 12:
return 12.*x*y*r*r;
266 case 13:
return 12.*x*x*y*r;
267 case 14:
return 12.*x*y*y*r;
270 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
278 unsigned int shape=i;
280 libmesh_assert_less (i, 21);
282 if ((i>= 3&&i<= 6) && elem->
point(0) > elem->
point(1)) shape=9-i;
283 if ((i>= 7&&i<=10) && elem->
point(1) > elem->
point(2)) shape=17-i;
284 if ((i>=11&&i<=14) && elem->
point(0) > elem->
point(2)) shape=25-i;
289 case 0:
return pow<5>(r);
290 case 1:
return pow<5>(x);
291 case 2:
return pow<5>(y);
294 case 3:
return 5.*x *pow<4>(r);
295 case 4:
return 10.*pow<2>(x)*pow<3>(r);
296 case 5:
return 10.*pow<3>(x)*pow<2>(r);
297 case 6:
return 5.*pow<4>(x)*r;
299 case 7:
return 5.*y *pow<4>(x);
300 case 8:
return 10.*pow<2>(y)*pow<3>(x);
301 case 9:
return 10.*pow<3>(y)*pow<2>(x);
302 case 10:
return 5.*pow<4>(y)*x;
304 case 11:
return 5.*y *pow<4>(r);
305 case 12:
return 10.*pow<2>(y)*pow<3>(r);
306 case 13:
return 10.*pow<3>(y)*pow<2>(r);
307 case 14:
return 5.*pow<4>(y)*r;
310 case 15:
return 20.*x*y*pow<3>(r);
311 case 16:
return 30.*x*pow<2>(y)*pow<2>(r);
312 case 17:
return 30.*pow<2>(x)*y*pow<2>(r);
313 case 18:
return 20.*x*pow<3>(y)*r;
314 case 19:
return 20.*pow<3>(x)*y*r;
315 case 20:
return 30.*pow<2>(x)*pow<2>(y)*r;
318 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
326 unsigned int shape=i;
328 libmesh_assert_less (i, 28);
330 if ((i>= 3&&i<= 7) && elem->
point(0) > elem->
point(1)) shape=10-i;
331 if ((i>= 8&&i<=12) && elem->
point(1) > elem->
point(2)) shape=20-i;
332 if ((i>=13&&i<=17) && elem->
point(0) > elem->
point(2)) shape=30-i;
337 case 0:
return pow<6>(r);
338 case 1:
return pow<6>(x);
339 case 2:
return pow<6>(y);
342 case 3:
return 6.*x *pow<5>(r);
343 case 4:
return 15.*pow<2>(x)*pow<4>(r);
344 case 5:
return 20.*pow<3>(x)*pow<3>(r);
345 case 6:
return 15.*pow<4>(x)*pow<2>(r);
346 case 7:
return 6.*pow<5>(x)*r;
348 case 8:
return 6.*y *pow<5>(x);
349 case 9:
return 15.*pow<2>(y)*pow<4>(x);
350 case 10:
return 20.*pow<3>(y)*pow<3>(x);
351 case 11:
return 15.*pow<4>(y)*pow<2>(x);
352 case 12:
return 6.*pow<5>(y)*x;
354 case 13:
return 6.*y *pow<5>(r);
355 case 14:
return 15.*pow<2>(y)*pow<4>(r);
356 case 15:
return 20.*pow<3>(y)*pow<3>(r);
357 case 16:
return 15.*pow<4>(y)*pow<2>(r);
358 case 17:
return 6.*pow<5>(y)*r;
361 case 18:
return 30.*x*y*pow<4>(r);
362 case 19:
return 60.*x*pow<2>(y)*pow<3>(r);
363 case 20:
return 60.* pow<2>(x)*y*pow<3>(r);
364 case 21:
return 60.*x*pow<3>(y)*pow<2>(r);
365 case 22:
return 60.*pow<3>(x)*y*pow<2>(r);
366 case 23:
return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
367 case 24:
return 30.*x*pow<4>(y)*r;
368 case 25:
return 60.*pow<2>(x)*pow<3>(y)*r;
369 case 26:
return 60.*pow<3>(x)*pow<2>(y)*r;
370 case 27:
return 30.*pow<4>(x)*y*r;
373 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
377 libmesh_error_msg(
"Invalid totalorder = " << totalorder);
381 libmesh_error_msg(
"ERROR: Unsupported element type = " << type);
394 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
403 const unsigned int i,
404 const unsigned int j,
406 const bool add_p_level)
412 const Order totalorder =
413 static_cast<Order>(order + add_p_level * elem->
p_level());
422 const Real xi = p(0);
423 const Real eta = p(1);
425 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
441 else if (i < totalorder + 3u)
442 { i0 = i - 2; i1 = 0; }
443 else if (i < 2u*totalorder + 2)
444 { i0 = 1; i1 = i - totalorder - 1; }
445 else if (i < 3u*totalorder + 1)
446 { i0 = i - 2u*totalorder; i1 = 1; }
447 else if (i < 4u*totalorder)
448 { i0 = 0; i1 = i - 3u*totalorder + 1; }
452 unsigned int basisnum = i - 4*totalorder;
460 if ((i>= 4 && i<= 4+ totalorder-2u) && elem->
point(0) > elem->
point(1)) i0=totalorder+2-i0;
461 else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->
point(1) > elem->
point(2)) i1=totalorder+2-i1;
462 else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->
point(3) > elem->
point(2)) i0=totalorder+2-i0;
463 else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->
point(0) > elem->
point(3)) i1=totalorder+2-i1;
478 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
487 libmesh_assert_less (totalorder, 3);
489 const Real xi = p(0);
490 const Real eta = p(1);
492 libmesh_assert_less (i, 8);
495 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
496 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
497 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
517 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
523 libmesh_assert_less (totalorder, 2);
524 libmesh_fallthrough();
529 const Real eps = 1.e-6;
536 const Point pp(p(0)+eps, p(1));
537 const Point pm(p(0)-eps, p(1));
546 const Point pp(p(0), p(1)+eps);
547 const Point pm(p(0), p(1)-eps);
555 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
560 libmesh_error_msg(
"ERROR: Unsupported element type = " << type);
566 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
574 static bool warning_given =
false;
577 libMesh::err <<
"Second derivatives for Bernstein elements "
578 <<
"are not yet implemented!"
581 warning_given =
true;
595 static bool warning_given =
false;
598 libMesh::err <<
"Second derivatives for Bernstein elements "
599 <<
"are not yet implemented!"
602 warning_given =
true;
611 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES