libMesh
fe_bernstein_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22 
23 // libmesh includes
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"
29 
30 
31 namespace {
32 
33 using namespace libMesh;
34 
35 /*
36  * Helper function for finding indices, when computing quad shape
37  * functions and derivatives via a tensor-product of 1D functions and
38  * derivatives.
39  */
40 std::pair<unsigned int, unsigned int> quad_i0_i1 (const unsigned int i,
41  const Order totalorder,
42  const Elem & elem)
43 {
44  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
45 
46  unsigned int i0, i1;
47 
48  // Vertex DoFs
49  if (i == 0)
50  { i0 = 0; i1 = 0; }
51  else if (i == 1)
52  { i0 = 1; i1 = 0; }
53  else if (i == 2)
54  { i0 = 1; i1 = 1; }
55  else if (i == 3)
56  { i0 = 0; i1 = 1; }
57 
58 
59  // Edge DoFs
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; }
68  // Interior DoFs
69  else
70  {
71  unsigned int basisnum = i - 4*totalorder;
72  i0 = square_number_column[basisnum] + 2;
73  i1 = square_number_row[basisnum] + 2;
74  }
75 
76 
77  // Flip odd degree of freedom values if necessary
78  // to keep continuity on sides
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;
83 
84  return std::make_pair(i0, i1);
85 }
86 
87 }
88 
89 namespace libMesh
90 {
91 
92 
94 
95 
96 template <>
97 Real FE<2,BERNSTEIN>::shape(const Elem * elem,
98  const Order order,
99  const unsigned int i,
100  const Point & p,
101  const bool add_p_level)
102 {
103  libmesh_assert(elem);
104 
105  const ElemType type = elem->type();
106 
107  const Order totalorder =
108  order + add_p_level*elem->p_level();
109 
110  // Declare that we are using our own special power function
111  // from the Utility namespace. This saves typing later.
112  using Utility::pow;
113 
114  switch (type)
115  {
116  // Hierarchic shape functions on the quadrilateral.
117  case QUAD4:
118  case QUADSHELL4:
119  case QUAD9:
120  case QUADSHELL9:
121  {
122  // Compute quad shape functions as a tensor-product
123  auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
124 
125  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, p(0))*
126  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, p(1)));
127  }
128  // handle serendipity QUAD8 element separately
129  case QUAD8:
130  case QUADSHELL8:
131  {
132  libmesh_assert_less (totalorder, 3);
133 
134  const Real xi = p(0);
135  const Real eta = p(1);
136 
137  libmesh_assert_less (i, 8);
138 
139  // 0 1 2 3 4 5 6 7 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};
143 
144  //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta
145  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
146  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)
147  +scal[i]*
148  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)*
149  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta));
150 
151  }
152 
153  case TRI3:
154  case TRISHELL3:
155  libmesh_assert_less (totalorder, 2);
156  libmesh_fallthrough();
157  case TRI6:
158  case TRI7:
159  switch (totalorder)
160  {
161  case FIRST:
162  {
163  const Real x=p(0);
164  const Real y=p(1);
165  const Real r=1.-x-y;
166 
167  libmesh_assert_less (i, 3);
168 
169  switch(i)
170  {
171  case 0: return r; //f0,0,1
172  case 1: return x; //f0,1,1
173  case 2: return y; //f1,0,1
174 
175  default:
176  libmesh_error_msg("Invalid shape function index i = " << i);
177  }
178  }
179  case SECOND:
180  {
181  const Real x=p(0);
182  const Real y=p(1);
183  const Real r=1.-x-y;
184 
185  libmesh_assert_less (i, 6);
186 
187  switch(i)
188  {
189  case 0: return r*r;
190  case 1: return x*x;
191  case 2: return y*y;
192 
193  case 3: return 2.*x*r;
194  case 4: return 2.*x*y;
195  case 5: return 2.*r*y;
196 
197  default:
198  libmesh_error_msg("Invalid shape function index i = " << i);
199  }
200  }
201  case THIRD:
202  {
203  const Real x=p(0);
204  const Real y=p(1);
205  const Real r=1.-x-y;
206  libmesh_assert_less (i, 10);
207 
208  unsigned int shape=i;
209 
210 
211  if ((i==3||i==4) && elem->positive_edge_orientation(0)) shape=7-i;
212  if ((i==5||i==6) && elem->positive_edge_orientation(1)) shape=11-i;
213  if ((i==7||i==8) && !elem->positive_edge_orientation(2)) shape=15-i;
214 
215  switch(shape)
216  {
217  case 0: return r*r*r;
218  case 1: return x*x*x;
219  case 2: return y*y*y;
220 
221  case 3: return 3.*x*r*r;
222  case 4: return 3.*x*x*r;
223 
224  case 5: return 3.*y*x*x;
225  case 6: return 3.*y*y*x;
226 
227  case 7: return 3.*y*r*r;
228  case 8: return 3.*y*y*r;
229 
230  case 9: return 6.*x*y*r;
231 
232  default:
233  libmesh_error_msg("Invalid shape function index shape = " << shape);
234  }
235  }
236  case FOURTH:
237  {
238  const Real x=p(0);
239  const Real y=p(1);
240  const Real r=1-x-y;
241  unsigned int shape=i;
242 
243  libmesh_assert_less (i, 15);
244 
245  if ((i==3||i== 5) && elem->positive_edge_orientation(0)) shape=8-i;
246  if ((i==6||i== 8) && elem->positive_edge_orientation(1)) shape=14-i;
247  if ((i==9||i==11) && !elem->positive_edge_orientation(2)) shape=20-i;
248 
249 
250  switch(shape)
251  {
252  // point functions
253  case 0: return r*r*r*r;
254  case 1: return x*x*x*x;
255  case 2: return y*y*y*y;
256 
257  // edge functions
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;
261 
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;
265 
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;
269 
270  // inner functions
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;
274 
275  default:
276  libmesh_error_msg("Invalid shape function index shape = " << shape);
277  }
278  }
279  case FIFTH:
280  {
281  const Real x=p(0);
282  const Real y=p(1);
283  const Real r=1-x-y;
284  unsigned int shape=i;
285 
286  libmesh_assert_less (i, 21);
287 
288  if ((i>= 3&&i<= 6) && elem->positive_edge_orientation(0)) shape=9-i;
289  if ((i>= 7&&i<=10) && elem->positive_edge_orientation(1)) shape=17-i;
290  if ((i>=11&&i<=14) && !elem->positive_edge_orientation(2)) shape=25-i;
291 
292  switch(shape)
293  {
294  //point functions
295  case 0: return pow<5>(r);
296  case 1: return pow<5>(x);
297  case 2: return pow<5>(y);
298 
299  //edge functions
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;
304 
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;
309 
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;
314 
315  //inner functions
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;
322 
323  default:
324  libmesh_error_msg("Invalid shape function index shape = " << shape);
325  }
326  }
327  case SIXTH:
328  {
329  const Real x=p(0);
330  const Real y=p(1);
331  const Real r=1-x-y;
332  unsigned int shape=i;
333 
334  libmesh_assert_less (i, 28);
335 
336  if ((i>= 3&&i<= 7) && elem->positive_edge_orientation(0)) shape=10-i;
337  if ((i>= 8&&i<=12) && elem->positive_edge_orientation(1)) shape=20-i;
338  if ((i>=13&&i<=17) && !elem->positive_edge_orientation(2)) shape=30-i;
339 
340  switch(shape)
341  {
342  //point functions
343  case 0: return pow<6>(r);
344  case 1: return pow<6>(x);
345  case 2: return pow<6>(y);
346 
347  //edge functions
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;
353 
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;
359 
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;
365 
366  //inner functions
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;
377 
378  default:
379  libmesh_error_msg("Invalid shape function index shape = " << shape);
380  } // switch shape
381  }
382  default:
383  libmesh_error_msg("Invalid totalorder = " << totalorder);
384  } // switch order
385 
386  default:
387  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
388  } // switch type
389 }
390 
391 
392 template <>
394  const Order,
395  const unsigned int,
396  const Point &)
397 {
398  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
399  return 0.;
400 }
401 
402 
403 template <>
405  const Elem * elem,
406  const unsigned int i,
407  const Point & p,
408  const bool add_p_level)
409 {
410  return FE<2,BERNSTEIN>::shape(elem, fet.order, i, p, add_p_level);
411 }
412 
413 
414 
415 template <>
417  const Order order,
418  const unsigned int i,
419  const unsigned int j,
420  const Point & p,
421  const bool add_p_level)
422 {
423  libmesh_assert(elem);
424 
425  const ElemType type = elem->type();
426 
427  const Order totalorder =
428  order + add_p_level*elem->p_level();
429 
430  switch (type)
431  {
432  // Hierarchic shape functions on the quadrilateral.
433  case QUAD4:
434  case QUAD9:
435  case QUADSHELL9:
436  {
437  // Compute quad shape functions as a tensor-product
438  auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
439 
440  switch (j)
441  {
442  // d()/dxi
443  case 0:
444  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
445  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, p(1)));
446 
447  // d()/deta
448  case 1:
449  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, p(0))*
450  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
451 
452  default:
453  libmesh_error_msg("Invalid shape function derivative j = " << j);
454  }
455  }
456 
457  // Bernstein shape functions on the 8-noded quadrilateral
458  // is handled separately.
459  case QUAD8:
460  case QUADSHELL8:
461  {
462  libmesh_assert_less (totalorder, 3);
463 
464  const Real xi = p(0);
465  const Real eta = p(1);
466 
467  libmesh_assert_less (i, 8);
468 
469  // 0 1 2 3 4 5 6 7 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};
473  switch (j)
474  {
475  // d()/dxi
476  case 0:
477  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
478  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
479  +scal[i]*
480  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
481  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
482 
483  // d()/deta
484  case 1:
485  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
486  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
487  +scal[i]*
488  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
489  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
490 
491  default:
492  libmesh_error_msg("Invalid shape function derivative j = " << j);
493  }
494  }
495 
496  case TRI3:
497  case TRISHELL3:
498  libmesh_assert_less (totalorder, 2);
499  libmesh_fallthrough();
500  case TRI6:
501  case TRI7:
502  {
503  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,BERNSTEIN>::shape);
504  }
505 
506  default:
507  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
508  }
509 }
510 
511 
512 
513 
514 template <>
516  const Order,
517  const unsigned int,
518  const unsigned int,
519  const Point &)
520 {
521  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
522  return 0.;
523 }
524 
525 template <>
527  const Elem * elem,
528  const unsigned int i,
529  const unsigned int j,
530  const Point & p,
531  const bool add_p_level)
532 {
533  return FE<2,BERNSTEIN>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
534 }
535 
536 
537 
538 
539 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
540 
541 
542 
543 template <>
545  const Order order,
546  const unsigned int i,
547  const unsigned int j,
548  const Point & p,
549  const bool add_p_level)
550 {
551  libmesh_assert(elem);
552 
553  const ElemType type = elem->type();
554 
555  const Order totalorder =
556  order + add_p_level*elem->p_level();
557 
558  switch (type)
559  {
560  // Hierarchic shape functions on the quadrilateral.
561  case QUAD4:
562  case QUAD9:
563  case QUADSHELL9:
564  {
565  // Compute quad shape functions as a tensor-product
566  auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
567 
568  switch (j)
569  {
570  // d^2() / dxi^2
571  case 0:
572  return (FE<1,BERNSTEIN>::shape_second_deriv(EDGE3, totalorder, i0, 0, p(0))*
573  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, p(1)));
574 
575  // d^2() / dxi deta
576  case 1:
577  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
578  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
579 
580  // d^2() / deta^2
581  case 2:
582  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, p(0))*
583  FE<1,BERNSTEIN>::shape_second_deriv(EDGE3, totalorder, i1, 0, p(1)));
584 
585  default:
586  libmesh_error_msg("Invalid shape function derivative j = " << j);
587  }
588  }
589 
590  // Going to be lazy again about the hard cases.
591  case TRI3:
592  case TRISHELL3:
593  libmesh_assert_less (totalorder, 2);
594  libmesh_fallthrough();
595  case QUAD8:
596  case QUADSHELL8:
597  case TRI6:
598  case TRI7:
599  {
600  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
602  }
603 
604  default:
605  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
606  }
607 }
608 
609 
610 template <>
612  const Order,
613  const unsigned int,
614  const unsigned int,
615  const Point &)
616 {
617  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
618  return 0.;
619 }
620 
621 
622 template <>
624  const Elem * elem,
625  const unsigned int i,
626  const unsigned int j,
627  const Point & p,
628  const bool add_p_level)
629 {
630  libmesh_assert(elem);
631  return FE<2,BERNSTEIN>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
632 }
633 
634 
635 
636 #endif
637 
638 } // namespace libMesh
639 
640 
641 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
ElemType
Defines an enum for geometric element types.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
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.
Definition: elem.h:94
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
Definition: elem.h:3108
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
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&#39;t be...
Definition: fe.C:911
const unsigned char square_number_column[]
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
A specific instantiation of the FEBase class.
Definition: fe.h:127
T pow(const T &x)
Definition: utility.h:328
libmesh_assert(ctx)
bool positive_edge_orientation(const unsigned int i) const
Definition: elem.C:3589
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))
Definition: fe.C:969
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.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2453
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)