libMesh
Loading...
Searching...
No Matches
fe_bernstein_shape_2D.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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
31namespace {
32
33using 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 */
40std::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
89namespace libMesh
90{
91
92
94
95
96template <>
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
392template <>
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
403template <>
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
415template <>
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
514template <>
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
525template <>
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
543template <>
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
610template <>
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
622template <>
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
This is the base class from which all geometric element types are derived.
Definition elem.h:96
const Point & point(const unsigned int i) const
Definition elem.h:2462
bool positive_edge_orientation(const unsigned int i) const
Definition elem.C:3634
virtual ElemType type() const =0
unsigned int p_level() const
Definition elem.h:3125
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
A specific instantiation of the FEBase class.
Definition fe.h:128
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.
Definition point.h:40
T pow(const T &x)
Definition utility.h:296
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...
Definition fe.C:911
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))
Definition fe.C:969
libmesh_assert(ctx)
const unsigned char square_number_row[]
const unsigned char square_number_column[]
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real