libMesh
fe_bernstein_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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  static_cast<Order>(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  {
121  // Compute quad shape functions as a tensor-product
122  auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
123 
124  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, p(0))*
125  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, p(1)));
126  }
127  // handle serendipity QUAD8 element separately
128  case QUAD8:
129  case QUADSHELL8:
130  {
131  libmesh_assert_less (totalorder, 3);
132 
133  const Real xi = p(0);
134  const Real eta = p(1);
135 
136  libmesh_assert_less (i, 8);
137 
138  // 0 1 2 3 4 5 6 7 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};
142 
143  //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta
144  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
145  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)
146  +scal[i]*
147  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)*
148  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta));
149 
150  }
151 
152  case TRI3:
153  case TRISHELL3:
154  libmesh_assert_less (totalorder, 2);
155  libmesh_fallthrough();
156  case TRI6:
157  case TRI7:
158  switch (totalorder)
159  {
160  case FIRST:
161  {
162  const Real x=p(0);
163  const Real y=p(1);
164  const Real r=1.-x-y;
165 
166  libmesh_assert_less (i, 3);
167 
168  switch(i)
169  {
170  case 0: return r; //f0,0,1
171  case 1: return x; //f0,1,1
172  case 2: return y; //f1,0,1
173 
174  default:
175  libmesh_error_msg("Invalid shape function index i = " << i);
176  }
177  }
178  case SECOND:
179  {
180  const Real x=p(0);
181  const Real y=p(1);
182  const Real r=1.-x-y;
183 
184  libmesh_assert_less (i, 6);
185 
186  switch(i)
187  {
188  case 0: return r*r;
189  case 1: return x*x;
190  case 2: return y*y;
191 
192  case 3: return 2.*x*r;
193  case 4: return 2.*x*y;
194  case 5: return 2.*r*y;
195 
196  default:
197  libmesh_error_msg("Invalid shape function index i = " << i);
198  }
199  }
200  case THIRD:
201  {
202  const Real x=p(0);
203  const Real y=p(1);
204  const Real r=1.-x-y;
205  libmesh_assert_less (i, 10);
206 
207  unsigned int shape=i;
208 
209 
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;
213 
214  switch(shape)
215  {
216  case 0: return r*r*r;
217  case 1: return x*x*x;
218  case 2: return y*y*y;
219 
220  case 3: return 3.*x*r*r;
221  case 4: return 3.*x*x*r;
222 
223  case 5: return 3.*y*x*x;
224  case 6: return 3.*y*y*x;
225 
226  case 7: return 3.*y*r*r;
227  case 8: return 3.*y*y*r;
228 
229  case 9: return 6.*x*y*r;
230 
231  default:
232  libmesh_error_msg("Invalid shape function index shape = " << shape);
233  }
234  }
235  case FOURTH:
236  {
237  const Real x=p(0);
238  const Real y=p(1);
239  const Real r=1-x-y;
240  unsigned int shape=i;
241 
242  libmesh_assert_less (i, 15);
243 
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;
247 
248 
249  switch(shape)
250  {
251  // point functions
252  case 0: return r*r*r*r;
253  case 1: return x*x*x*x;
254  case 2: return y*y*y*y;
255 
256  // edge functions
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;
260 
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;
264 
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;
268 
269  // inner functions
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;
273 
274  default:
275  libmesh_error_msg("Invalid shape function index shape = " << shape);
276  }
277  }
278  case FIFTH:
279  {
280  const Real x=p(0);
281  const Real y=p(1);
282  const Real r=1-x-y;
283  unsigned int shape=i;
284 
285  libmesh_assert_less (i, 21);
286 
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;
290 
291  switch(shape)
292  {
293  //point functions
294  case 0: return pow<5>(r);
295  case 1: return pow<5>(x);
296  case 2: return pow<5>(y);
297 
298  //edge functions
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;
303 
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;
308 
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;
313 
314  //inner functions
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;
321 
322  default:
323  libmesh_error_msg("Invalid shape function index shape = " << shape);
324  }
325  }
326  case SIXTH:
327  {
328  const Real x=p(0);
329  const Real y=p(1);
330  const Real r=1-x-y;
331  unsigned int shape=i;
332 
333  libmesh_assert_less (i, 28);
334 
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;
338 
339  switch(shape)
340  {
341  //point functions
342  case 0: return pow<6>(r);
343  case 1: return pow<6>(x);
344  case 2: return pow<6>(y);
345 
346  //edge functions
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;
352 
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;
358 
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;
364 
365  //inner functions
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;
376 
377  default:
378  libmesh_error_msg("Invalid shape function index shape = " << shape);
379  } // switch shape
380  }
381  default:
382  libmesh_error_msg("Invalid totalorder = " << totalorder);
383  } // switch order
384 
385  default:
386  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
387  } // switch type
388 }
389 
390 
391 template <>
393  const Order,
394  const unsigned int,
395  const Point &)
396 {
397  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
398  return 0.;
399 }
400 
401 
402 template <>
404  const Elem * elem,
405  const unsigned int i,
406  const Point & p,
407  const bool add_p_level)
408 {
409  return FE<2,BERNSTEIN>::shape(elem, fet.order, i, p, add_p_level);
410 }
411 
412 
413 
414 template <>
416  const Order order,
417  const unsigned int i,
418  const unsigned int j,
419  const Point & p,
420  const bool add_p_level)
421 {
422  libmesh_assert(elem);
423 
424  const ElemType type = elem->type();
425 
426  const Order totalorder =
427  static_cast<Order>(order + add_p_level * elem->p_level());
428 
429  switch (type)
430  {
431  // Hierarchic shape functions on the quadrilateral.
432  case QUAD4:
433  case QUAD9:
434  {
435  // Compute quad shape functions as a tensor-product
436  auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
437 
438  switch (j)
439  {
440  // d()/dxi
441  case 0:
442  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
443  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, p(1)));
444 
445  // d()/deta
446  case 1:
447  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, p(0))*
448  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
449 
450  default:
451  libmesh_error_msg("Invalid shape function derivative j = " << j);
452  }
453  }
454 
455  // Bernstein shape functions on the 8-noded quadrilateral
456  // is handled separately.
457  case QUAD8:
458  case QUADSHELL8:
459  {
460  libmesh_assert_less (totalorder, 3);
461 
462  const Real xi = p(0);
463  const Real eta = p(1);
464 
465  libmesh_assert_less (i, 8);
466 
467  // 0 1 2 3 4 5 6 7 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};
471  switch (j)
472  {
473  // d()/dxi
474  case 0:
475  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
476  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
477  +scal[i]*
478  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
479  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
480 
481  // d()/deta
482  case 1:
483  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
484  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
485  +scal[i]*
486  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
487  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
488 
489  default:
490  libmesh_error_msg("Invalid shape function derivative j = " << j);
491  }
492  }
493 
494  case TRI3:
495  case TRISHELL3:
496  libmesh_assert_less (totalorder, 2);
497  libmesh_fallthrough();
498  case TRI6:
499  case TRI7:
500  {
501  return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,BERNSTEIN>::shape);
502  }
503 
504  default:
505  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
506  }
507 }
508 
509 
510 
511 
512 template <>
514  const Order,
515  const unsigned int,
516  const unsigned int,
517  const Point &)
518 {
519  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
520  return 0.;
521 }
522 
523 template <>
525  const Elem * elem,
526  const unsigned int i,
527  const unsigned int j,
528  const Point & p,
529  const bool add_p_level)
530 {
531  return FE<2,BERNSTEIN>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
532 }
533 
534 
535 
536 
537 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
538 
539 
540 
541 template <>
543  const Order order,
544  const unsigned int i,
545  const unsigned int j,
546  const Point & p,
547  const bool add_p_level)
548 {
549  libmesh_assert(elem);
550 
551  const ElemType type = elem->type();
552 
553  const Order totalorder =
554  static_cast<Order>(order + add_p_level * elem->p_level());
555 
556  switch (type)
557  {
558  // Hierarchic shape functions on the quadrilateral.
559  case QUAD4:
560  case QUAD9:
561  {
562  // Compute quad shape functions as a tensor-product
563  auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
564 
565  switch (j)
566  {
567  // d^2() / dxi^2
568  case 0:
569  return (FE<1,BERNSTEIN>::shape_second_deriv(EDGE3, totalorder, i0, 0, p(0))*
570  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, p(1)));
571 
572  // d^2() / dxi deta
573  case 1:
574  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
575  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
576 
577  // d^2() / deta^2
578  case 2:
579  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, p(0))*
580  FE<1,BERNSTEIN>::shape_second_deriv(EDGE3, totalorder, i1, 0, p(1)));
581 
582  default:
583  libmesh_error_msg("Invalid shape function derivative j = " << j);
584  }
585  }
586 
587  // Going to be lazy again about the hard cases.
588  case TRI3:
589  case TRISHELL3:
590  libmesh_assert_less (totalorder, 2);
591  libmesh_fallthrough();
592  case QUAD8:
593  case QUADSHELL8:
594  case TRI6:
595  case TRI7:
596  {
597  return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
599  }
600 
601  default:
602  libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
603  }
604 }
605 
606 
607 template <>
609  const Order,
610  const unsigned int,
611  const unsigned int,
612  const Point &)
613 {
614  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
615  return 0.;
616 }
617 
618 
619 template <>
621  const Elem * elem,
622  const unsigned int i,
623  const unsigned int j,
624  const Point & p,
625  const bool add_p_level)
626 {
627  libmesh_assert(elem);
628  return FE<2,BERNSTEIN>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
629 }
630 
631 
632 
633 #endif
634 
635 } // namespace libMesh
636 
637 
638 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
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:2945
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
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:744
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)
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:851
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:2277
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)