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 order + add_p_level*elem->
p_level();
123 auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
132 libmesh_assert_less (totalorder, 3);
134 const Real xi = p(0);
135 const Real eta = p(1);
137 libmesh_assert_less (i, 8);
140 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
141 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
142 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
155 libmesh_assert_less (totalorder, 2);
156 libmesh_fallthrough();
167 libmesh_assert_less (i, 3);
176 libmesh_error_msg(
"Invalid shape function index i = " << i);
185 libmesh_assert_less (i, 6);
193 case 3:
return 2.*x*r;
194 case 4:
return 2.*x*y;
195 case 5:
return 2.*r*y;
198 libmesh_error_msg(
"Invalid shape function index i = " << i);
206 libmesh_assert_less (i, 10);
208 unsigned int shape=i;
217 case 0:
return r*r*r;
218 case 1:
return x*x*x;
219 case 2:
return y*y*y;
221 case 3:
return 3.*x*r*r;
222 case 4:
return 3.*x*x*r;
224 case 5:
return 3.*y*x*x;
225 case 6:
return 3.*y*y*x;
227 case 7:
return 3.*y*r*r;
228 case 8:
return 3.*y*y*r;
230 case 9:
return 6.*x*y*r;
233 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
241 unsigned int shape=i;
243 libmesh_assert_less (i, 15);
253 case 0:
return r*r*r*r;
254 case 1:
return x*x*x*x;
255 case 2:
return y*y*y*y;
258 case 3:
return 4.*x*r*r*r;
259 case 4:
return 6.*x*x*r*r;
260 case 5:
return 4.*x*x*x*r;
262 case 6:
return 4.*y*x*x*x;
263 case 7:
return 6.*y*y*x*x;
264 case 8:
return 4.*y*y*y*x;
266 case 9:
return 4.*y*r*r*r;
267 case 10:
return 6.*y*y*r*r;
268 case 11:
return 4.*y*y*y*r;
271 case 12:
return 12.*x*y*r*r;
272 case 13:
return 12.*x*x*y*r;
273 case 14:
return 12.*x*y*y*r;
276 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
284 unsigned int shape=i;
286 libmesh_assert_less (i, 21);
295 case 0:
return pow<5>(r);
296 case 1:
return pow<5>(x);
297 case 2:
return pow<5>(y);
300 case 3:
return 5.*x *pow<4>(r);
301 case 4:
return 10.*pow<2>(x)*pow<3>(r);
302 case 5:
return 10.*pow<3>(x)*pow<2>(r);
303 case 6:
return 5.*pow<4>(x)*r;
305 case 7:
return 5.*y *pow<4>(x);
306 case 8:
return 10.*pow<2>(y)*pow<3>(x);
307 case 9:
return 10.*pow<3>(y)*pow<2>(x);
308 case 10:
return 5.*pow<4>(y)*x;
310 case 11:
return 5.*y *pow<4>(r);
311 case 12:
return 10.*pow<2>(y)*pow<3>(r);
312 case 13:
return 10.*pow<3>(y)*pow<2>(r);
313 case 14:
return 5.*pow<4>(y)*r;
316 case 15:
return 20.*x*y*pow<3>(r);
317 case 16:
return 30.*x*pow<2>(y)*pow<2>(r);
318 case 17:
return 30.*pow<2>(x)*y*pow<2>(r);
319 case 18:
return 20.*x*pow<3>(y)*r;
320 case 19:
return 20.*pow<3>(x)*y*r;
321 case 20:
return 30.*pow<2>(x)*pow<2>(y)*r;
324 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
332 unsigned int shape=i;
334 libmesh_assert_less (i, 28);
343 case 0:
return pow<6>(r);
344 case 1:
return pow<6>(x);
345 case 2:
return pow<6>(y);
348 case 3:
return 6.*x *pow<5>(r);
349 case 4:
return 15.*pow<2>(x)*pow<4>(r);
350 case 5:
return 20.*pow<3>(x)*pow<3>(r);
351 case 6:
return 15.*pow<4>(x)*pow<2>(r);
352 case 7:
return 6.*pow<5>(x)*r;
354 case 8:
return 6.*y *pow<5>(x);
355 case 9:
return 15.*pow<2>(y)*pow<4>(x);
356 case 10:
return 20.*pow<3>(y)*pow<3>(x);
357 case 11:
return 15.*pow<4>(y)*pow<2>(x);
358 case 12:
return 6.*pow<5>(y)*x;
360 case 13:
return 6.*y *pow<5>(r);
361 case 14:
return 15.*pow<2>(y)*pow<4>(r);
362 case 15:
return 20.*pow<3>(y)*pow<3>(r);
363 case 16:
return 15.*pow<4>(y)*pow<2>(r);
364 case 17:
return 6.*pow<5>(y)*r;
367 case 18:
return 30.*x*y*pow<4>(r);
368 case 19:
return 60.*x*pow<2>(y)*pow<3>(r);
369 case 20:
return 60.* pow<2>(x)*y*pow<3>(r);
370 case 21:
return 60.*x*pow<3>(y)*pow<2>(r);
371 case 22:
return 60.*pow<3>(x)*y*pow<2>(r);
372 case 23:
return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
373 case 24:
return 30.*x*pow<4>(y)*r;
374 case 25:
return 60.*pow<2>(x)*pow<3>(y)*r;
375 case 26:
return 60.*pow<3>(x)*pow<2>(y)*r;
376 case 27:
return 30.*pow<4>(x)*y*r;
379 libmesh_error_msg(
"Invalid shape function index shape = " << shape);
383 libmesh_error_msg(
"Invalid totalorder = " << totalorder);
398 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
406 const unsigned int i,
408 const bool add_p_level)
418 const unsigned int i,
419 const unsigned int j,
421 const bool add_p_level)
427 const Order totalorder =
428 order + add_p_level*elem->
p_level();
438 auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
453 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
462 libmesh_assert_less (totalorder, 3);
464 const Real xi = p(0);
465 const Real eta = p(1);
467 libmesh_assert_less (i, 8);
470 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
471 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
472 static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
492 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
498 libmesh_assert_less (totalorder, 2);
499 libmesh_fallthrough();
521 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
528 const unsigned int i,
529 const unsigned int j,
531 const bool add_p_level)
539 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 546 const unsigned int i,
547 const unsigned int j,
549 const bool add_p_level)
555 const Order totalorder =
556 order + add_p_level*elem->
p_level();
566 auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
586 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
593 libmesh_assert_less (totalorder, 2);
594 libmesh_fallthrough();
617 libmesh_error_msg(
"Bernstein polynomials require the element type \nbecause edge orientation is needed.");
625 const unsigned int i,
626 const unsigned int j,
628 const bool add_p_level)
641 #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.
bool positive_edge_orientation(const unsigned int i) const
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)