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