libMesh
Loading...
Searching...
No Matches
fe_szabab_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// Local includes
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/utility.h"
27#include "libmesh/enum_to_string.h"
28
29// C++ includes
30#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
31#include <cmath> // for std::sqrt
32
33// Anonymous namespace to hold static std::sqrt values
34namespace
35{
36using namespace libMesh;
37
38static const Real sqrt2 = std::sqrt(2.);
39static const Real sqrt6 = std::sqrt(6.);
40static const Real sqrt10 = std::sqrt(10.);
41static const Real sqrt14 = std::sqrt(14.);
42static const Real sqrt22 = std::sqrt(22.);
43static const Real sqrt26 = std::sqrt(26.);
44
45
46// Take care of edge orientation, this is needed at edge shapes with
47// (y=0)-asymmetric 1D shapes
48Real quad_flip(const Elem * elem,
49 const unsigned int totalorder,
50 const unsigned int i)
51{
52 libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
53
54 // vertex and interior shape functions don't flip
55 if (i < 4 || i >= 4*totalorder)
56 return 1;
57
58 const int edge = (i-4) / (totalorder-1);
59 libmesh_assert_less(edge, 4);
60 libmesh_assert_greater_equal(edge, 0);
61
62 const int edge_i = i - 4 - edge*(totalorder-1);
63
64 // The "natural" orientation is p1>p0 on edge 0,
65 // p2>p1 on e1, p2>p3 on e2, p3>p0 on e3
66 if (edge_i%2 &&
67 (elem->positive_edge_orientation(edge) ==
68 (edge < 2)))
69 return -1;
70
71 return 1;
72}
73
74} // anonymous namespace
75
76
77namespace libMesh
78{
79
80
82
83
84template <>
85Real FE<2,SZABAB>::shape(const Elem * elem,
86 const Order order,
87 const unsigned int i,
88 const Point & p,
89 const bool add_p_level)
90{
91 libmesh_assert(elem);
92
93 const ElemType type = elem->type();
94
95 const Order totalorder = order + add_p_level*elem->p_level();
96
97 // Declare that we are using our own special power function
98 // from the Utility namespace. This saves typing later.
99 using Utility::pow;
100
101 switch (totalorder)
102 {
103 // 1st & 2nd-order Szabo-Babuska.
104 case FIRST:
105 case SECOND:
106 {
107 switch (type)
108 {
109
110 // Szabo-Babuska shape functions on the triangle.
111 case TRI3:
112 case TRI6:
113 case TRI7:
114 {
115 const Real l1 = 1-p(0)-p(1);
116 const Real l2 = p(0);
117 const Real l3 = p(1);
118
119 libmesh_assert_less (i, 6);
120
121 switch (i)
122 {
123 case 0: return l1;
124 case 1: return l2;
125 case 2: return l3;
126
127 case 3: return l1*l2*(-4.*sqrt6);
128 case 4: return l2*l3*(-4.*sqrt6);
129 case 5: return l3*l1*(-4.*sqrt6);
130
131 default:
132 libmesh_error_msg("Invalid i = " << i);
133 }
134 }
135
136
137 // Szabo-Babuska shape functions on the quadrilateral.
138 case QUAD4:
139 case QUAD8:
140 case QUAD9:
141 {
142 // Compute quad shape functions as a tensor-product
143 const Real xi = p(0);
144 const Real eta = p(1);
145
146 libmesh_assert_less (i, 9);
147
148 // 0 1 2 3 4 5 6 7 8
149 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
150 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
151
152 return (FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
153 FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
154
155 }
156
157 default:
158 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
159 }
160 }
161
162
163 // 3rd-order Szabo-Babuska.
164 case THIRD:
165 {
166 switch (type)
167 {
168
169 // Szabo-Babuska shape functions on the triangle.
170 case TRI6:
171 case TRI7:
172 {
173 Real l1 = 1-p(0)-p(1);
174 Real l2 = p(0);
175 Real l3 = p(1);
176
177 Real f=1;
178
179 libmesh_assert_less (i, 10);
180
181
182 if (i==4 && elem->positive_edge_orientation(0))f=-1;
183 if (i==6 && elem->positive_edge_orientation(1))f=-1;
184 if (i==8 && elem->positive_edge_orientation(2))f=-1;
185
186
187 switch (i)
188 {
189 //nodal modes
190 case 0: return l1;
191 case 1: return l2;
192 case 2: return l3;
193
194 //side modes
195 case 3: return l1*l2*(-4.*sqrt6);
196 case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
197
198 case 5: return l2*l3*(-4.*sqrt6);
199 case 6: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
200
201 case 7: return l3*l1*(-4.*sqrt6);
202 case 8: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
203
204 //internal modes
205 case 9: return l1*l2*l3;
206
207 default:
208 libmesh_error_msg("Invalid i = " << i);
209 }
210 }
211
212
213 // Szabo-Babuska shape functions on the quadrilateral.
214 case QUAD8:
215 case QUAD9:
216 {
217 // Compute quad shape functions as a tensor-product
218 const Real xi = p(0);
219 const Real eta = p(1);
220
221 libmesh_assert_less (i, 16);
222
223 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
224 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
225 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
226
227 const Real f = quad_flip(elem, totalorder, i);
228
229 return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
230 FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
231 }
232
233 default:
234 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
235 }
236 }
237
238
239
240
241 // 4th-order Szabo-Babuska.
242 case FOURTH:
243 {
244 switch (type)
245 {
246 // Szabo-Babuska shape functions on the triangle.
247 case TRI6:
248 case TRI7:
249 {
250 Real l1 = 1-p(0)-p(1);
251 Real l2 = p(0);
252 Real l3 = p(1);
253
254 Real f=1;
255
256 libmesh_assert_less (i, 15);
257
258
259 if (i== 4 && elem->positive_edge_orientation(0))f=-1;
260 if (i== 7 && elem->positive_edge_orientation(1))f=-1;
261 if (i==10 && elem->positive_edge_orientation(2))f=-1;
262
263
264 switch (i)
265 {
266 //nodal modes
267 case 0: return l1;
268 case 1: return l2;
269 case 2: return l3;
270
271 //side modes
272 case 3: return l1*l2*(-4.*sqrt6);
273 case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
274 case 5: return l1*l2*(-sqrt14)*(5.*pow<2>(l2-l1)-1);
275
276 case 6: return l2*l3*(-4.*sqrt6);
277 case 7: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
278 case 8: return l2*l3*(-sqrt14)*(5.*pow<2>(l3-l2)-1);
279
280 case 9: return l3*l1*(-4.*sqrt6);
281 case 10: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
282 case 11: return l3*l1*(-sqrt14)*(5.*pow<2>(l1-l3)-1);
283
284 //internal modes
285 case 12: return l1*l2*l3;
286
287 case 13: return l1*l2*l3*(l2-l1);
288 case 14: return l1*l2*l3*(2*l3-1);
289
290 default:
291 libmesh_error_msg("Invalid i = " << i);
292 }
293 }
294
295
296 // Szabo-Babuska shape functions on the quadrilateral.
297 case QUAD8:
298 case QUAD9:
299 {
300 // Compute quad shape functions as a tensor-product
301 const Real xi = p(0);
302 const Real eta = p(1);
303
304 libmesh_assert_less (i, 25);
305
306 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
307 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
308 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
309
310 const Real f = quad_flip(elem, totalorder, i);
311
312 return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
313 FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
314 }
315
316 default:
317 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
318 }
319 }
320
321
322
323
324 // 5th-order Szabo-Babuska.
325 case FIFTH:
326 {
327 switch (type)
328 {
329 // Szabo-Babuska shape functions on the triangle.
330 case TRI6:
331 case TRI7:
332 {
333 Real l1 = 1-p(0)-p(1);
334 Real l2 = p(0);
335 Real l3 = p(1);
336
337 const Real x=l2-l1;
338 const Real y=2.*l3-1;
339
340 Real f=1;
341
342 libmesh_assert_less (i, 21);
343
344
345 if ((i== 4||i== 6) && elem->positive_edge_orientation(0))f=-1;
346 if ((i== 8||i==10) && elem->positive_edge_orientation(1))f=-1;
347 if ((i==12||i==14) && elem->positive_edge_orientation(2))f=-1;
348
349
350 switch (i)
351 {
352 //nodal modes
353 case 0: return l1;
354 case 1: return l2;
355 case 2: return l3;
356
357 //side modes
358 case 3: return l1*l2*(-4.*sqrt6);
359 case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
360 case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
361 case 6: return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
362
363 case 7: return l2*l3*(-4.*sqrt6);
364 case 8: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
365 case 9: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
366 case 10: return -f*sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
367
368 case 11: return l3*l1*(-4.*sqrt6);
369 case 12: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
370 case 13: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
371 case 14: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
372
373 //internal modes
374 case 15: return l1*l2*l3;
375
376 case 16: return l1*l2*l3*x;
377 case 17: return l1*l2*l3*y;
378
379 case 18: return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
380 case 19: return l1*l2*l3*x*y;
381 case 20: return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
382
383 default:
384 libmesh_error_msg("Invalid i = " << i);
385 }
386 } // case TRI6/TRI7
387
388 // Szabo-Babuska shape functions on the quadrilateral.
389 case QUAD8:
390 case QUAD9:
391 {
392 // Compute quad shape functions as a tensor-product
393 const Real xi = p(0);
394 const Real eta = p(1);
395
396 libmesh_assert_less (i, 36);
397
398 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
399 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
400 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
401
402 const Real f = quad_flip(elem, totalorder, i);
403
404 return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
405 FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
406
407 } // case QUAD8/QUAD9
408
409 default:
410 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
411
412 } // switch type
413
414 } // case FIFTH
415
416 // 6th-order Szabo-Babuska.
417 case SIXTH:
418 {
419 switch (type)
420 {
421 // Szabo-Babuska shape functions on the triangle.
422 case TRI6:
423 case TRI7:
424 {
425 Real l1 = 1-p(0)-p(1);
426 Real l2 = p(0);
427 Real l3 = p(1);
428
429 const Real x=l2-l1;
430 const Real y=2.*l3-1;
431
432 Real f=1;
433
434 libmesh_assert_less (i, 28);
435
436
437 if ((i== 4||i== 6) && elem->positive_edge_orientation(0))f=-1;
438 if ((i== 9||i==11) && elem->positive_edge_orientation(1))f=-1;
439 if ((i==14||i==16) && elem->positive_edge_orientation(2))f=-1;
440
441
442 switch (i)
443 {
444 //nodal modes
445 case 0: return l1;
446 case 1: return l2;
447 case 2: return l3;
448
449 //side modes
450 case 3: return l1*l2*(-4.*sqrt6);
451 case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
452 case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
453 case 6: return f*(-sqrt2)*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
454 case 7: return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
455
456 case 8: return l2*l3*(-4.*sqrt6);
457 case 9: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
458 case 10: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
459 case 11: return f*(-sqrt2)*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
460 case 12: return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
461
462 case 13: return l3*l1*(-4.*sqrt6);
463 case 14: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
464 case 15: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
465 case 16: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
466 case 17: return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
467
468
469
470 //internal modes
471 case 18: return l1*l2*l3;
472
473 case 19: return l1*l2*l3*x;
474 case 20: return l1*l2*l3*y;
475
476 case 21: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2);
477 case 22: return l1*l2*l3*(l2-l1)*(2.0*l3-1.0);
478 case 23: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3);
479 case 24: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
480 case 25: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
481 case 26: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
482 case 27: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
483
484
485 default:
486 libmesh_error_msg("Invalid i = " << i);
487 }
488 } // case TRI6/TRI7
489
490 // Szabo-Babuska shape functions on the quadrilateral.
491 case QUAD8:
492 case QUAD9:
493 {
494 // Compute quad shape functions as a tensor-product
495 const Real xi = p(0);
496 const Real eta = p(1);
497
498 libmesh_assert_less (i, 49);
499
500 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
501 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
502 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
503
504 const Real f = quad_flip(elem, totalorder, i);
505
506 return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
507 FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
508
509 } // case QUAD8/QUAD9
510
511 default:
512 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
513
514 } // switch type
515
516 } // case SIXTH
517
518
519 // 7th-order Szabo-Babuska.
520 case SEVENTH:
521 {
522 switch (type)
523 {
524 // Szabo-Babuska shape functions on the triangle.
525 case TRI6:
526 case TRI7:
527 {
528
529 Real l1 = 1-p(0)-p(1);
530 Real l2 = p(0);
531 Real l3 = p(1);
532
533 const Real x=l2-l1;
534 const Real y=2.*l3-1.;
535
536 Real f=1;
537
538 libmesh_assert_less (i, 36);
539
540
541 if ((i>= 4&&i<= 8) && elem->positive_edge_orientation(0))f=-1;
542 if ((i>=10&&i<=14) && elem->positive_edge_orientation(1))f=-1;
543 if ((i>=16&&i<=20) && elem->positive_edge_orientation(2))f=-1;
544
545
546 switch (i)
547 {
548 //nodal modes
549 case 0: return l1;
550 case 1: return l2;
551 case 2: return l3;
552
553 //side modes
554 case 3: return l1*l2*(-4.*sqrt6);
555 case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
556
557 case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
558 case 6: return f*-sqrt2*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
559 case 7: return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
560 case 8: return f*-sqrt26*l1*l2*((-0.25E1+(15.0-0.165E2*l1*l1)*l1*l1)*l1+(0.25E1+(-45.0+0.825E2*l1*l1)*l1*l1+((45.0-0.165E3*l1*l1)*l1+(-15.0+0.165E3*l1*l1+(-0.825E2*l1+0.165E2*l2)*l2)*l2)*l2)*l2);
561
562 case 9: return l2*l3*(-4.*sqrt6);
563 case 10: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
564
565 case 11: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
566 case 12: return f*-sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
567 case 13: return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
568 case 14: return f*-sqrt26*l2*l3*((0.25E1+(-15.0+0.165E2*l3*l3)*l3*l3)*l3+(-0.25E1+(45.0-0.825E2*l3*l3)*l3*l3+((-45.0+0.165E3*l3*l3)*l3+(15.0-0.165E3*l3*l3+(0.825E2*l3-0.165E2*l2)*l2)*l2)*l2)*l2);
569
570 case 15: return l3*l1*(-4.*sqrt6);
571 case 16: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
572
573 case 17: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
574 case 18: return -f*sqrt2*l3*l1*((9.-21.*l3*l3)*l3+(-9.+63.*l3*l3+(-63.*l3+21.*l1)*l1)*l1);
575 case 19: return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
576 case 20: return f*-sqrt26*l3*l1*((-0.25E1+(15.0-0.165E2*l3*l3)*l3*l3)*l3+(0.25E1+(-45.0+0.825E2*l3*l3)*l3*l3+((45.0-0.165E3*l3*l3)*l3+(-15.0+0.165E3*l3*l3+(-0.825E2*l3+0.165E2*l1)*l1)*l1)*l1)*l1);
577
578
579 //internal modes
580 case 21: return l1*l2*l3;
581
582 case 22: return l1*l2*l3*x;
583 case 23: return l1*l2*l3*y;
584
585 case 24: return l1*l2*l3*0.5*(3.*pow<2>(x)-1.);
586 case 25: return l1*l2*l3*x*y;
587 case 26: return l1*l2*l3*0.5*(3.*pow<2>(y)-1.);
588
589 case 27: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
590 case 28: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
591 case 29: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
592 case 30: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
593 case 31: return 0.125*l1*l2*l3*((-15.0+(70.0-63.0*l1*l1)*l1*l1)*l1+(15.0+(-210.0+315.0*l1*l1)*l1*l1+((210.0-630.0*l1*l1)*l1+(-70.0+630.0*l1*l1+(-315.0*l1+63.0*l2)*l2)*l2)*l2)*l2);
594 case 32: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2)*(2.0*l3-1.0);
595 case 33: return 0.25*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0+(-12.0+12.0*l3)*l3);
596 case 34: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3)*(l2-l1);
597 case 35: return 0.125*l1*l2*l3*(-8.0+(240.0+(-1680.0+(4480.0+(-5040.0+2016.0*l3)*l3)*l3)*l3)*l3);
598
599 default:
600 libmesh_error_msg("Invalid i = " << i);
601 }
602 } // case TRI6/TRI7
603
604 // Szabo-Babuska shape functions on the quadrilateral.
605 case QUAD8:
606 case QUAD9:
607 {
608 // Compute quad shape functions as a tensor-product
609 const Real xi = p(0);
610 const Real eta = p(1);
611
612 libmesh_assert_less (i, 64);
613
614 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
615 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
616 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
617
618 const Real f = quad_flip(elem, totalorder, i);
619
620 return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
621 FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
622
623 } // case QUAD8/QUAD9
624
625 default:
626 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
627
628 } // switch type
629
630 } // case SEVENTH
631
632
633 // by default throw an error
634 default:
635 libmesh_error_msg("ERROR: Unsupported polynomial order!");
636 } // switch order
637}
638
639
640
641
642
643template <>
645 const Order,
646 const unsigned int,
647 const Point &)
648{
649 libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
650 return 0.;
651}
652
653
654
655template <>
657 const Elem * elem,
658 const unsigned int i,
659 const Point & p,
660 const bool add_p_level)
661{
662 return FE<2,SZABAB>::shape(elem, fet.order, i, p, add_p_level);
663}
664
665
666
667
668
669template <>
671 const Order order,
672 const unsigned int i,
673 const unsigned int j,
674 const Point & p,
675 const bool add_p_level)
676{
677 libmesh_assert(elem);
678
679 const ElemType type = elem->type();
680
681 const Order totalorder = order + add_p_level*elem->p_level();
682
683 switch (totalorder)
684 {
685
686 // 1st & 2nd-order Szabo-Babuska.
687 case FIRST:
688 case SECOND:
689 {
690 switch (type)
691 {
692
693 // Szabo-Babuska shape functions on the triangle.
694 case TRI3:
695 case TRI6:
696 case TRI7:
697 {
698 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
699 }
700
701
702 // Szabo-Babuska shape functions on the quadrilateral.
703 case QUAD4:
704 case QUAD8:
705 case QUAD9:
706 {
707 // Compute quad shape functions as a tensor-product
708 const Real xi = p(0);
709 const Real eta = p(1);
710
711 libmesh_assert_less (i, 9);
712
713 // 0 1 2 3 4 5 6 7 8
714 static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
715 static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
716
717 switch (j)
718 {
719 // d()/dxi
720 case 0:
721 return (FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
722 FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
723
724 // d()/deta
725 case 1:
726 return (FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
727 FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
728
729 default:
730 libmesh_error_msg("Invalid j = " << j);
731 }
732 }
733
734 default:
735 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
736 }
737 }
738
739
740
741 // 3rd-order Szabo-Babuska.
742 case THIRD:
743 {
744 switch (type)
745 {
746 // Szabo-Babuska shape functions on the triangle.
747 case TRI6:
748 case TRI7:
749 {
750 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
751 }
752
753
754 // Szabo-Babuska shape functions on the quadrilateral.
755 case QUAD8:
756 case QUAD9:
757 {
758 // Compute quad shape functions as a tensor-product
759 const Real xi = p(0);
760 const Real eta = p(1);
761
762 libmesh_assert_less (i, 16);
763
764 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
765 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
766 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
767
768 const Real f = quad_flip(elem, totalorder, i);
769
770 switch (j)
771 {
772 // d()/dxi
773 case 0:
774 return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
775 FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
776
777 // d()/deta
778 case 1:
779 return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
780 FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
781
782 default:
783 libmesh_error_msg("Invalid j = " << j);
784 }
785 }
786
787 default:
788 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
789 }
790 }
791
792
793
794
795 // 4th-order Szabo-Babuska.
796 case FOURTH:
797 {
798 switch (type)
799 {
800
801 // Szabo-Babuska shape functions on the triangle.
802 case TRI6:
803 case TRI7:
804 {
805 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
806 }
807
808
809 // Szabo-Babuska shape functions on the quadrilateral.
810 case QUAD8:
811 case QUAD9:
812 {
813 // Compute quad shape functions as a tensor-product
814 const Real xi = p(0);
815 const Real eta = p(1);
816
817 libmesh_assert_less (i, 25);
818
819 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
820 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
821 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
822
823 const Real f = quad_flip(elem, totalorder, i);
824
825 switch (j)
826 {
827 // d()/dxi
828 case 0:
829 return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
830 FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
831
832 // d()/deta
833 case 1:
834 return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
835 FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
836
837 default:
838 libmesh_error_msg("Invalid j = " << j);
839 }
840 }
841
842 default:
843 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
844 }
845 }
846
847
848
849
850 // 5th-order Szabo-Babuska.
851 case FIFTH:
852 {
853 // Szabo-Babuska shape functions on the quadrilateral.
854 switch (type)
855 {
856
857 // Szabo-Babuska shape functions on the triangle.
858 case TRI6:
859 case TRI7:
860 {
861 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
862 }
863
864
865 case QUAD8:
866 case QUAD9:
867 {
868 // Compute quad shape functions as a tensor-product
869 const Real xi = p(0);
870 const Real eta = p(1);
871
872 libmesh_assert_less (i, 36);
873
874 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
875 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
876 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
877
878 const Real f = quad_flip(elem, totalorder, i);
879
880 switch (j)
881 {
882 // d()/dxi
883 case 0:
884 return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
885 FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
886
887 // d()/deta
888 case 1:
889 return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
890 FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
891
892 default:
893 libmesh_error_msg("Invalid j = " << j);
894 }
895 }
896
897 default:
898 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
899 }
900 }
901
902
903 // 6th-order Szabo-Babuska.
904 case SIXTH:
905 {
906 // Szabo-Babuska shape functions on the quadrilateral.
907 switch (type)
908 {
909
910 // Szabo-Babuska shape functions on the triangle.
911 case TRI6:
912 case TRI7:
913 {
914 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
915 }
916
917
918 case QUAD8:
919 case QUAD9:
920 {
921 // Compute quad shape functions as a tensor-product
922 const Real xi = p(0);
923 const Real eta = p(1);
924
925 libmesh_assert_less (i, 49);
926
927 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
928 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
929 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
930
931 const Real f = quad_flip(elem, totalorder, i);
932
933 switch (j)
934 {
935 // d()/dxi
936 case 0:
937 return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
938 FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
939
940 // d()/deta
941 case 1:
942 return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
943 FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
944
945 default:
946 libmesh_error_msg("Invalid j = " << j);
947 }
948 }
949
950 default:
951 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
952 }
953 }
954
955
956 // 7th-order Szabo-Babuska.
957 case SEVENTH:
958 {
959 // Szabo-Babuska shape functions on the quadrilateral.
960 switch (type)
961 {
962
963 // Szabo-Babuska shape functions on the triangle.
964 case TRI6:
965 case TRI7:
966 {
967 return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
968 }
969
970
971 case QUAD8:
972 case QUAD9:
973 {
974 // Compute quad shape functions as a tensor-product
975 const Real xi = p(0);
976 const Real eta = p(1);
977
978 libmesh_assert_less (i, 64);
979
980 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
981 static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
982 static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
983
984 const Real f = quad_flip(elem, totalorder, i);
985
986 switch (j)
987 {
988 // d()/dxi
989 case 0:
990 return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
991 FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
992
993 // d()/deta
994 case 1:
995 return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
996 FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
997
998 default:
999 libmesh_error_msg("Invalid j = " << j);
1000 }
1001 }
1002
1003 default:
1004 libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1005 }
1006 }
1007
1008
1009
1010 // by default throw an error;call the orientation-independent shape functions
1011 default:
1012 libmesh_error_msg("ERROR: Unsupported polynomial order!");
1013 }
1014}
1015
1016
1017
1018
1019
1020template <>
1022 const Order,
1023 const unsigned int,
1024 const unsigned int,
1025 const Point &)
1026{
1027 libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
1028 return 0.;
1029}
1030
1031
1032template <>
1034 const Elem * elem,
1035 const unsigned int i,
1036 const unsigned int j,
1037 const Point & p,
1038 const bool add_p_level)
1039{
1040 return FE<2,SZABAB>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1041}
1042
1043
1044
1045
1046#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1047
1048template <>
1050 const Order,
1051 const unsigned int,
1052 const unsigned int,
1053 const Point &)
1054{
1055 static bool warning_given = false;
1056
1057 if (!warning_given)
1058 libMesh::err << "Second derivatives for Szabab elements "
1059 << " are not yet implemented!"
1060 << std::endl;
1061
1062 warning_given = true;
1063 return 0.;
1064}
1065
1066
1067
1068template <>
1070 const Order,
1071 const unsigned int,
1072 const unsigned int,
1073 const Point &,
1074 const bool)
1075{
1076 static bool warning_given = false;
1077
1078 if (!warning_given)
1079 libMesh::err << "Second derivatives for Szabab elements "
1080 << " are not yet implemented!"
1081 << std::endl;
1082
1083 warning_given = true;
1084 return 0.;
1085}
1086
1087
1088template <>
1090 const Elem *,
1091 const unsigned int,
1092 const unsigned int,
1093 const Point &,
1094 const bool)
1095{
1096 static bool warning_given = false;
1097
1098 if (!warning_given)
1099 libMesh::err << "Second derivatives for Szabab elements "
1100 << " are not yet implemented!"
1101 << std::endl;
1102
1103 warning_given = true;
1104 return 0.;
1105}
1106
1107#endif
1108
1109} // namespace libMesh
1110
1111#endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
This is the base class from which all geometric element types are derived.
Definition elem.h:96
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.
OStreamProxy err
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.
libmesh_assert(ctx)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real