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/number_lookups.h" 27 #include "libmesh/utility.h" 28 #include "libmesh/enum_to_string.h" 40 std::pair<unsigned int, unsigned int> quad_i0_i1 (
const unsigned int i,
41 const Order totalorder,
44 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
60 else if (i < totalorder + 3u)
61 { i0 = i - 2; i1 = 0; }
62 else if (i < 2u*totalorder + 2)
63 { i0 = 1; i1 = i - totalorder - 1; }
64 else if (i < 3u*totalorder + 1)
65 { i0 = i - 2u*totalorder; i1 = 1; }
66 else if (i < 4u*totalorder)
67 { i0 = 0; i1 = i - 3u*totalorder + 1; }
71 unsigned int basisnum = i - 4*totalorder;
79 if ((i>= 4 && i<= 4+ totalorder-2u) && elem.
point(0) > elem.
point(1)) i0=totalorder+2-i0;
80 else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem.
point(1) > elem.
point(2)) i1=totalorder+2-i1;
81 else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem.
point(3) > elem.
point(2)) i0=totalorder+2-i0;
82 else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem.
point(0) > elem.
point(3)) i1=totalorder+2-i1;
84 return std::make_pair(i0, i1);
101 const bool add_p_level)
107 const Order totalorder =
108 static_cast<Order>(order + add_p_level * elem->
p_level());
122 auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
131 libmesh_assert_less (totalorder, 3);
133 const Real xi = p(0);
134 const Real eta = p(1);
136 libmesh_assert_less (i, 8);
139 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
140 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
141 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
154 libmesh_assert_less (totalorder, 2);
155 libmesh_fallthrough();
166 libmesh_assert_less (i, 3);
175 libmesh_error_msg(
"Invalid shape function index i = " << i);
184 libmesh_assert_less (i, 6);
192 case 3:
return 2.*x*r;
193 case 4:
return 2.*x*y;
194 case 5:
return 2.*r*y;
197 libmesh_error_msg(
"Invalid shape function index i = " << i);
205 libmesh_assert_less (i, 10);
207 unsigned int shape=i;
210 if ((i==3||i==4) && elem->
point(0) > elem->
point(1)) shape=7-i;
211 if ((i==5||i==6) && elem->
point(1) > elem->
point(2)) shape=11-i;
212 if ((i==7||i==8) && elem->
point(0) > elem->
point(2)) shape=15-i;
216 case 0:
return r*r*r;
217 case 1:
return x*x*x;
218 case 2:
return y*y*y;
220 case 3:
return 3.*x*r*r;
221 case 4:
return 3.*x*x*r;
223 case 5:
return 3.*y*x*x;
224 case 6:
return 3.*y*y*x;
226 case 7:
return 3.*y*r*r;
227 case 8:
return 3.*y*y*r;
229 case 9:
return 6.*x*y*r;
232 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
240 unsigned int shape=i;
242 libmesh_assert_less (i, 15);
244 if ((i==3||i== 5) && elem->
point(0) > elem->
point(1)) shape=8-i;
245 if ((i==6||i== 8) && elem->
point(1) > elem->
point(2)) shape=14-i;
246 if ((i==9||i==11) && elem->
point(0) > elem->
point(2)) shape=20-i;
252 case 0:
return r*r*r*r;
253 case 1:
return x*x*x*x;
254 case 2:
return y*y*y*y;
257 case 3:
return 4.*x*r*r*r;
258 case 4:
return 6.*x*x*r*r;
259 case 5:
return 4.*x*x*x*r;
261 case 6:
return 4.*y*x*x*x;
262 case 7:
return 6.*y*y*x*x;
263 case 8:
return 4.*y*y*y*x;
265 case 9:
return 4.*y*r*r*r;
266 case 10:
return 6.*y*y*r*r;
267 case 11:
return 4.*y*y*y*r;
270 case 12:
return 12.*x*y*r*r;
271 case 13:
return 12.*x*x*y*r;
272 case 14:
return 12.*x*y*y*r;
275 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
283 unsigned int shape=i;
285 libmesh_assert_less (i, 21);
287 if ((i>= 3&&i<= 6) && elem->
point(0) > elem->
point(1)) shape=9-i;
288 if ((i>= 7&&i<=10) && elem->
point(1) > elem->
point(2)) shape=17-i;
289 if ((i>=11&&i<=14) && elem->
point(0) > elem->
point(2)) shape=25-i;
294 case 0:
return pow<5>(r);
295 case 1:
return pow<5>(x);
296 case 2:
return pow<5>(y);
299 case 3:
return 5.*x *pow<4>(r);
300 case 4:
return 10.*pow<2>(x)*pow<3>(r);
301 case 5:
return 10.*pow<3>(x)*pow<2>(r);
302 case 6:
return 5.*pow<4>(x)*r;
304 case 7:
return 5.*y *pow<4>(x);
305 case 8:
return 10.*pow<2>(y)*pow<3>(x);
306 case 9:
return 10.*pow<3>(y)*pow<2>(x);
307 case 10:
return 5.*pow<4>(y)*x;
309 case 11:
return 5.*y *pow<4>(r);
310 case 12:
return 10.*pow<2>(y)*pow<3>(r);
311 case 13:
return 10.*pow<3>(y)*pow<2>(r);
312 case 14:
return 5.*pow<4>(y)*r;
315 case 15:
return 20.*x*y*pow<3>(r);
316 case 16:
return 30.*x*pow<2>(y)*pow<2>(r);
317 case 17:
return 30.*pow<2>(x)*y*pow<2>(r);
318 case 18:
return 20.*x*pow<3>(y)*r;
319 case 19:
return 20.*pow<3>(x)*y*r;
320 case 20:
return 30.*pow<2>(x)*pow<2>(y)*r;
323 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
331 unsigned int shape=i;
333 libmesh_assert_less (i, 28);
335 if ((i>= 3&&i<= 7) && elem->
point(0) > elem->
point(1)) shape=10-i;
336 if ((i>= 8&&i<=12) && elem->
point(1) > elem->
point(2)) shape=20-i;
337 if ((i>=13&&i<=17) && elem->
point(0) > elem->
point(2)) shape=30-i;
342 case 0:
return pow<6>(r);
343 case 1:
return pow<6>(x);
344 case 2:
return pow<6>(y);
347 case 3:
return 6.*x *pow<5>(r);
348 case 4:
return 15.*pow<2>(x)*pow<4>(r);
349 case 5:
return 20.*pow<3>(x)*pow<3>(r);
350 case 6:
return 15.*pow<4>(x)*pow<2>(r);
351 case 7:
return 6.*pow<5>(x)*r;
353 case 8:
return 6.*y *pow<5>(x);
354 case 9:
return 15.*pow<2>(y)*pow<4>(x);
355 case 10:
return 20.*pow<3>(y)*pow<3>(x);
356 case 11:
return 15.*pow<4>(y)*pow<2>(x);
357 case 12:
return 6.*pow<5>(y)*x;
359 case 13:
return 6.*y *pow<5>(r);
360 case 14:
return 15.*pow<2>(y)*pow<4>(r);
361 case 15:
return 20.*pow<3>(y)*pow<3>(r);
362 case 16:
return 15.*pow<4>(y)*pow<2>(r);
363 case 17:
return 6.*pow<5>(y)*r;
366 case 18:
return 30.*x*y*pow<4>(r);
367 case 19:
return 60.*x*pow<2>(y)*pow<3>(r);
368 case 20:
return 60.* pow<2>(x)*y*pow<3>(r);
369 case 21:
return 60.*x*pow<3>(y)*pow<2>(r);
370 case 22:
return 60.*pow<3>(x)*y*pow<2>(r);
371 case 23:
return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
372 case 24:
return 30.*x*pow<4>(y)*r;
373 case 25:
return 60.*pow<2>(x)*pow<3>(y)*r;
374 case 26:
return 60.*pow<3>(x)*pow<2>(y)*r;
375 case 27:
return 30.*pow<4>(x)*y*r;
378 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
382 libmesh_error_msg(
"Invalid totalorder = " << totalorder);
397 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
405 const unsigned int i,
407 const bool add_p_level)
417 const unsigned int i,
418 const unsigned int j,
420 const bool add_p_level)
426 const Order totalorder =
427 static_cast<Order>(order + add_p_level * elem->
p_level());
436 auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
451 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
460 libmesh_assert_less (totalorder, 3);
462 const Real xi = p(0);
463 const Real eta = p(1);
465 libmesh_assert_less (i, 8);
468 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
469 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
470 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
490 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
496 libmesh_assert_less (totalorder, 2);
497 libmesh_fallthrough();
519 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
526 const unsigned int i,
527 const unsigned int j,
529 const bool add_p_level)
537 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 544 const unsigned int i,
545 const unsigned int j,
547 const bool add_p_level)
553 const Order totalorder =
554 static_cast<Order>(order + add_p_level * elem->
p_level());
563 auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
583 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
590 libmesh_assert_less (totalorder, 2);
591 libmesh_fallthrough();
614 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
622 const unsigned int i,
623 const unsigned int j,
625 const bool add_p_level)
638 #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.
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...
const unsigned char square_number_column[]
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
A specific instantiation of the FEBase class.
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))
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const unsigned char square_number_row[]
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)