libMesh
fe_szabab_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // 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
34 namespace
35 {
36 using namespace libMesh;
37 
38 static const Real sqrt2 = std::sqrt(2.);
39 static const Real sqrt6 = std::sqrt(6.);
40 static const Real sqrt10 = std::sqrt(10.);
41 static const Real sqrt14 = std::sqrt(14.);
42 static const Real sqrt22 = std::sqrt(22.);
43 static 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
48 Real 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 
77 namespace libMesh
78 {
79 
80 
81 LIBMESH_DEFAULT_VECTORIZED_FE(2,SZABAB)
82 
83 
84 template <>
85 Real 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 
643 template <>
644 Real FE<2,SZABAB>::shape(const ElemType,
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 
655 template <>
656 Real FE<2,SZABAB>::shape(const FEType fet,
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 
669 template <>
670 Real FE<2,SZABAB>::shape_deriv(const Elem * elem,
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 
1020 template <>
1021 Real FE<2,SZABAB>::shape_deriv(const ElemType,
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 
1032 template <>
1033 Real FE<2,SZABAB>::shape_deriv(const FEType fet,
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 
1048 template <>
1049 Real FE<2,SZABAB>::shape_second_deriv(const ElemType,
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 
1068 template <>
1069 Real FE<2,SZABAB>::shape_second_deriv(const Elem *,
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 
1088 template <>
1089 Real FE<2,SZABAB>::shape_second_deriv(const FEType,
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
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::err
OStreamProxy err
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3109
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
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&#39;t be...
Definition: fe.C:911
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition: fe_bernstein_shape_0D.C:30
libMesh::FE
A specific instantiation of the FEBase class.
Definition: fe.h:127
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:328
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition: elem.C:3604
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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)
Generated on Thu Oct 30 2025 20:01:53 for libMesh by   doxygen 1.8.14