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"
40std::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)
This is the base class from which all geometric element types are derived.
const Point & point(const unsigned int i) const
bool positive_edge_orientation(const unsigned int i) const
virtual ElemType type() const =0
unsigned int p_level() const
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
A specific instantiation of the FEBase class.
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
A Point defines a location in LIBMESH_DIM dimensional Real space.
std::string enum_to_string(const T e)
The libMesh namespace provides an interface to certain functionality in the library.
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
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...
ElemType
Defines an enum for geometric element types.
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))
const unsigned char square_number_row[]
const unsigned char square_number_column[]
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real