libMesh
fe_szabab_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::sqrt
23 
24 
25 // Local includes
26 #include "libmesh/libmesh_config.h"
27 
28 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
29 
30 #include "libmesh/fe.h"
31 #include "libmesh/elem.h"
32 #include "libmesh/utility.h"
33 
34 
35 // Anonymous namespace to hold static std::sqrt values
36 namespace
37 {
38 using libMesh::Real;
39 
40 static const Real sqrt2 = std::sqrt(2.);
41 static const Real sqrt6 = std::sqrt(6.);
42 static const Real sqrt10 = std::sqrt(10.);
43 static const Real sqrt14 = std::sqrt(14.);
44 static const Real sqrt22 = std::sqrt(22.);
45 static const Real sqrt26 = std::sqrt(26.);
46 }
47 
48 
49 namespace libMesh
50 {
51 
52 template <>
54  const Order,
55  const unsigned int,
56  const Point &)
57 {
58  libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
59  return 0.;
60 }
61 
62 
63 
64 template <>
66  const Order order,
67  const unsigned int i,
68  const Point & p)
69 {
70  libmesh_assert(elem);
71 
72  const ElemType type = elem->type();
73 
74  const Order totalorder = static_cast<Order>(order + elem->p_level());
75 
76  // Declare that we are using our own special power function
77  // from the Utility namespace. This saves typing later.
78  using Utility::pow;
79 
80  switch (totalorder)
81  {
82  // 1st & 2nd-order Szabo-Babuska.
83  case FIRST:
84  case SECOND:
85  {
86  switch (type)
87  {
88 
89  // Szabo-Babuska shape functions on the triangle.
90  case TRI3:
91  case TRI6:
92  {
93  const Real l1 = 1-p(0)-p(1);
94  const Real l2 = p(0);
95  const Real l3 = p(1);
96 
97  libmesh_assert_less (i, 6);
98 
99  switch (i)
100  {
101  case 0: return l1;
102  case 1: return l2;
103  case 2: return l3;
104 
105  case 3: return l1*l2*(-4.*sqrt6);
106  case 4: return l2*l3*(-4.*sqrt6);
107  case 5: return l3*l1*(-4.*sqrt6);
108 
109  default:
110  libmesh_error_msg("Invalid i = " << i);
111  }
112  }
113 
114 
115  // Szabo-Babuska shape functions on the quadrilateral.
116  case QUAD4:
117  case QUAD8:
118  case QUAD9:
119  {
120  // Compute quad shape functions as a tensor-product
121  const Real xi = p(0);
122  const Real eta = p(1);
123 
124  libmesh_assert_less (i, 9);
125 
126  // 0 1 2 3 4 5 6 7 8
127  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
128  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
129 
130  return (FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
131  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
132 
133  }
134 
135  default:
136  libmesh_error_msg("Invalid element type = " << type);
137  }
138  }
139 
140 
141  // 3rd-order Szabo-Babuska.
142  case THIRD:
143  {
144  switch (type)
145  {
146 
147  // Szabo-Babuska shape functions on the triangle.
148  case TRI6:
149  {
150  Real l1 = 1-p(0)-p(1);
151  Real l2 = p(0);
152  Real l3 = p(1);
153 
154  Real f=1;
155 
156  libmesh_assert_less (i, 10);
157 
158 
159  if (i==4 && (elem->point(0) > elem->point(1)))f=-1;
160  if (i==6 && (elem->point(1) > elem->point(2)))f=-1;
161  if (i==8 && (elem->point(2) > elem->point(0)))f=-1;
162 
163 
164  switch (i)
165  {
166  //nodal modes
167  case 0: return l1;
168  case 1: return l2;
169  case 2: return l3;
170 
171  //side modes
172  case 3: return l1*l2*(-4.*sqrt6);
173  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
174 
175  case 5: return l2*l3*(-4.*sqrt6);
176  case 6: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
177 
178  case 7: return l3*l1*(-4.*sqrt6);
179  case 8: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
180 
181  //internal modes
182  case 9: return l1*l2*l3;
183 
184  default:
185  libmesh_error_msg("Invalid i = " << i);
186  }
187  }
188 
189 
190  // Szabo-Babuska shape functions on the quadrilateral.
191  case QUAD8:
192  case QUAD9:
193  {
194  // Compute quad shape functions as a tensor-product
195  const Real xi = p(0);
196  const Real eta = p(1);
197 
198  libmesh_assert_less (i, 16);
199 
200  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
201  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
202  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
203 
204  Real f=1.;
205 
206  // take care of edge orientation, this is needed at
207  // edge shapes with (y=0)-asymmetric 1D shapes, these have
208  // one 1D shape index being 0 or 1, the other one being odd and >=3
209 
210  switch(i)
211  {
212  case 5: // edge 0 points
213  if (elem->point(0) > elem->point(1))f = -1.;
214  break;
215  case 7: // edge 1 points
216  if (elem->point(1) > elem->point(2))f = -1.;
217  break;
218  case 9: // edge 2 points
219  if (elem->point(3) > elem->point(2))f = -1.;
220  break;
221  case 11: // edge 3 points
222  if (elem->point(0) > elem->point(3))f = -1.;
223  break;
224 
225  default:
226  // Everything else keeps f=1
227  break;
228  }
229 
230  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
231  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
232  }
233 
234  default:
235  libmesh_error_msg("Invalid element type = " << type);
236  }
237  }
238 
239 
240 
241 
242  // 4th-order Szabo-Babuska.
243  case FOURTH:
244  {
245  switch (type)
246  {
247  // Szabo-Babuska shape functions on the triangle.
248  case TRI6:
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->point(0) > elem->point(1)))f=-1;
260  if (i== 7 && (elem->point(1) > elem->point(2)))f=-1;
261  if (i==10 && (elem->point(2) > elem->point(0)))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  Real f=1.;
311 
312  switch(i)
313  {
314  case 5: // edge 0 points
315  if (elem->point(0) > elem->point(1))f = -1.;
316  break;
317  case 8: // edge 1 points
318  if (elem->point(1) > elem->point(2))f = -1.;
319  break;
320  case 11: // edge 2 points
321  if (elem->point(3) > elem->point(2))f = -1.;
322  break;
323  case 14: // edge 3 points
324  if (elem->point(0) > elem->point(3))f = -1.;
325  break;
326 
327  default:
328  // Everything else keeps f=1
329  break;
330  }
331 
332  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
333  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
334  }
335 
336  default:
337  libmesh_error_msg("Invalid element type = " << type);
338  }
339  }
340 
341 
342 
343 
344  // 5th-order Szabo-Babuska.
345  case FIFTH:
346  {
347  switch (type)
348  {
349  // Szabo-Babuska shape functions on the triangle.
350  case TRI6:
351  {
352  Real l1 = 1-p(0)-p(1);
353  Real l2 = p(0);
354  Real l3 = p(1);
355 
356  const Real x=l2-l1;
357  const Real y=2.*l3-1;
358 
359  Real f=1;
360 
361  libmesh_assert_less (i, 21);
362 
363 
364  if ((i== 4||i== 6) && (elem->point(0) > elem->point(1)))f=-1;
365  if ((i== 8||i==10) && (elem->point(1) > elem->point(2)))f=-1;
366  if ((i==12||i==14) && (elem->point(2) > elem->point(0)))f=-1;
367 
368 
369  switch (i)
370  {
371  //nodal modes
372  case 0: return l1;
373  case 1: return l2;
374  case 2: return l3;
375 
376  //side modes
377  case 3: return l1*l2*(-4.*sqrt6);
378  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
379  case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
380  case 6: return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
381 
382  case 7: return l2*l3*(-4.*sqrt6);
383  case 8: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
384  case 9: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
385  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);
386 
387  case 11: return l3*l1*(-4.*sqrt6);
388  case 12: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
389  case 13: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
390  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);
391 
392  //internal modes
393  case 15: return l1*l2*l3;
394 
395  case 16: return l1*l2*l3*x;
396  case 17: return l1*l2*l3*y;
397 
398  case 18: return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
399  case 19: return l1*l2*l3*x*y;
400  case 20: return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
401 
402  default:
403  libmesh_error_msg("Invalid i = " << i);
404  }
405  } // case TRI6
406 
407  // Szabo-Babuska shape functions on the quadrilateral.
408  case QUAD8:
409  case QUAD9:
410  {
411  // Compute quad shape functions as a tensor-product
412  const Real xi = p(0);
413  const Real eta = p(1);
414 
415  libmesh_assert_less (i, 36);
416 
417  // 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
418  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};
419  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};
420 
421  Real f=1.;
422 
423  switch(i)
424  {
425  case 5: // edge 0 points
426  case 7:
427  if (elem->point(0) > elem->point(1))f = -1.;
428  break;
429  case 9: // edge 1 points
430  case 11:
431  if (elem->point(1) > elem->point(2))f = -1.;
432  break;
433  case 13: // edge 2 points
434  case 15:
435  if (elem->point(3) > elem->point(2))f = -1.;
436  break;
437  case 14: // edge 3 points
438  case 19:
439  if (elem->point(0) > elem->point(3))f = -1.;
440  break;
441 
442  default:
443  // Everything else keeps f=1
444  break;
445  }
446 
447  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
448  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
449 
450  } // case QUAD8/QUAD9
451 
452  default:
453  libmesh_error_msg("Invalid element type = " << type);
454 
455  } // switch type
456 
457  } // case FIFTH
458 
459  // 6th-order Szabo-Babuska.
460  case SIXTH:
461  {
462  switch (type)
463  {
464  // Szabo-Babuska shape functions on the triangle.
465  case TRI6:
466  {
467  Real l1 = 1-p(0)-p(1);
468  Real l2 = p(0);
469  Real l3 = p(1);
470 
471  const Real x=l2-l1;
472  const Real y=2.*l3-1;
473 
474  Real f=1;
475 
476  libmesh_assert_less (i, 28);
477 
478 
479  if ((i== 4||i== 6) && (elem->point(0) > elem->point(1)))f=-1;
480  if ((i== 9||i==11) && (elem->point(1) > elem->point(2)))f=-1;
481  if ((i==14||i==16) && (elem->point(2) > elem->point(0)))f=-1;
482 
483 
484  switch (i)
485  {
486  //nodal modes
487  case 0: return l1;
488  case 1: return l2;
489  case 2: return l3;
490 
491  //side modes
492  case 3: return l1*l2*(-4.*sqrt6);
493  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
494  case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
495  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);
496  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);
497 
498  case 8: return l2*l3*(-4.*sqrt6);
499  case 9: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
500  case 10: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
501  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);
502  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);
503 
504  case 13: return l3*l1*(-4.*sqrt6);
505  case 14: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
506  case 15: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
507  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);
508  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);
509 
510 
511 
512  //internal modes
513  case 18: return l1*l2*l3;
514 
515  case 19: return l1*l2*l3*x;
516  case 20: return l1*l2*l3*y;
517 
518  case 21: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2);
519  case 22: return l1*l2*l3*(l2-l1)*(2.0*l3-1.0);
520  case 23: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3);
521  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);
522  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);
523  case 26: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
524  case 27: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
525 
526 
527  default:
528  libmesh_error_msg("Invalid i = " << i);
529  }
530  } // case TRI6
531 
532  // Szabo-Babuska shape functions on the quadrilateral.
533  case QUAD8:
534  case QUAD9:
535  {
536  // Compute quad shape functions as a tensor-product
537  const Real xi = p(0);
538  const Real eta = p(1);
539 
540  libmesh_assert_less (i, 49);
541 
542  // 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
543  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};
544  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};
545 
546  Real f=1.;
547 
548  switch(i)
549  {
550  case 5: // edge 0 points
551  case 7:
552  if (elem->point(0) > elem->point(1))f = -1.;
553  break;
554  case 10: // edge 1 points
555  case 12:
556  if (elem->point(1) > elem->point(2))f = -1.;
557  break;
558  case 15: // edge 2 points
559  case 17:
560  if (elem->point(3) > elem->point(2))f = -1.;
561  break;
562  case 20: // edge 3 points
563  case 22:
564  if (elem->point(0) > elem->point(3))f = -1.;
565  break;
566 
567  default:
568  // Everything else keeps f=1
569  break;
570  }
571 
572  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
573  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
574 
575  } // case QUAD8/QUAD9
576 
577  default:
578  libmesh_error_msg("Invalid element type = " << type);
579 
580  } // switch type
581 
582  } // case SIXTH
583 
584 
585  // 7th-order Szabo-Babuska.
586  case SEVENTH:
587  {
588  switch (type)
589  {
590  // Szabo-Babuska shape functions on the triangle.
591  case TRI6:
592  {
593 
594  Real l1 = 1-p(0)-p(1);
595  Real l2 = p(0);
596  Real l3 = p(1);
597 
598  const Real x=l2-l1;
599  const Real y=2.*l3-1.;
600 
601  Real f=1;
602 
603  libmesh_assert_less (i, 36);
604 
605 
606  if ((i>= 4&&i<= 8) && (elem->point(0) > elem->point(1)))f=-1;
607  if ((i>=10&&i<=14) && (elem->point(1) > elem->point(2)))f=-1;
608  if ((i>=16&&i<=20) && (elem->point(2) > elem->point(0)))f=-1;
609 
610 
611  switch (i)
612  {
613  //nodal modes
614  case 0: return l1;
615  case 1: return l2;
616  case 2: return l3;
617 
618  //side modes
619  case 3: return l1*l2*(-4.*sqrt6);
620  case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
621 
622  case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
623  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);
624  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);
625  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);
626 
627  case 9: return l2*l3*(-4.*sqrt6);
628  case 10: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
629 
630  case 11: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
631  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);
632  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);
633  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);
634 
635  case 15: return l3*l1*(-4.*sqrt6);
636  case 16: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
637 
638  case 17: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
639  case 18: return -f*sqrt2*l3*l1*((9.-21.*l3*l3)*l3+(-9.+63.*l3*l3+(-63.*l3+21.*l1)*l1)*l1);
640  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);
641  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);
642 
643 
644  //internal modes
645  case 21: return l1*l2*l3;
646 
647  case 22: return l1*l2*l3*x;
648  case 23: return l1*l2*l3*y;
649 
650  case 24: return l1*l2*l3*0.5*(3.*pow<2>(x)-1.);
651  case 25: return l1*l2*l3*x*y;
652  case 26: return l1*l2*l3*0.5*(3.*pow<2>(y)-1.);
653 
654  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);
655  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);
656  case 29: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
657  case 30: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
658  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);
659  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);
660  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);
661  case 34: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3)*(l2-l1);
662  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);
663 
664  default:
665  libmesh_error_msg("Invalid i = " << i);
666  }
667  } // case TRI6
668 
669  // Szabo-Babuska shape functions on the quadrilateral.
670  case QUAD8:
671  case QUAD9:
672  {
673  // Compute quad shape functions as a tensor-product
674  const Real xi = p(0);
675  const Real eta = p(1);
676 
677  libmesh_assert_less (i, 64);
678 
679  // 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
680  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};
681  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};
682 
683  Real f=1.;
684 
685  switch(i)
686  {
687  case 5: // edge 0 points
688  case 7:
689  case 9:
690  if (elem->point(0) > elem->point(1))f = -1.;
691  break;
692  case 11: // edge 1 points
693  case 13:
694  case 15:
695  if (elem->point(1) > elem->point(2))f = -1.;
696  break;
697  case 17: // edge 2 points
698  case 19:
699  case 21:
700  if (elem->point(3) > elem->point(2))f = -1.;
701  break;
702  case 23: // edge 3 points
703  case 25:
704  case 27:
705  if (elem->point(0) > elem->point(3))f = -1.;
706  break;
707 
708  default:
709  // Everything else keeps f=1
710  break;
711  }
712 
713  return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
714  FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
715 
716  } // case QUAD8/QUAD9
717 
718  default:
719  libmesh_error_msg("Invalid element type = " << type);
720 
721  } // switch type
722 
723  } // case SEVENTH
724 
725 
726  // by default throw an error
727  default:
728  libmesh_error_msg("ERROR: Unsupported polynomial order!");
729  } // switch order
730 
731  libmesh_error_msg("We'll never get here!");
732  return 0.;
733 }
734 
735 
736 
737 
738 
739 template <>
741  const Order,
742  const unsigned int,
743  const unsigned int,
744  const Point &)
745 {
746  libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
747  return 0.;
748 }
749 
750 
751 
752 template <>
754  const Order order,
755  const unsigned int i,
756  const unsigned int j,
757  const Point & p)
758 {
759  libmesh_assert(elem);
760 
761  const ElemType type = elem->type();
762 
763  const Order totalorder = static_cast<Order>(order + elem->p_level());
764 
765  switch (totalorder)
766  {
767 
768  // 1st & 2nd-order Szabo-Babuska.
769  case FIRST:
770  case SECOND:
771  {
772  switch (type)
773  {
774 
775  // Szabo-Babuska shape functions on the triangle.
776  case TRI3:
777  case TRI6:
778  {
779  // Here we use finite differences to compute the derivatives!
780  const Real eps = 1.e-6;
781 
782  libmesh_assert_less (i, 6);
783  libmesh_assert_less (j, 2);
784 
785  switch (j)
786  {
787  // d()/dxi
788  case 0:
789  {
790  const Point pp(p(0)+eps, p(1));
791  const Point pm(p(0)-eps, p(1));
792 
793  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
794  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
795  }
796 
797  // d()/deta
798  case 1:
799  {
800  const Point pp(p(0), p(1)+eps);
801  const Point pm(p(0), p(1)-eps);
802 
803  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
804  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
805  }
806 
807  default:
808  libmesh_error_msg("Invalid j = " << j);
809  }
810  }
811 
812 
813 
814  // Szabo-Babuska shape functions on the quadrilateral.
815  case QUAD4:
816  case QUAD8:
817  case QUAD9:
818  {
819  // Compute quad shape functions as a tensor-product
820  const Real xi = p(0);
821  const Real eta = p(1);
822 
823  libmesh_assert_less (i, 9);
824 
825  // 0 1 2 3 4 5 6 7 8
826  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
827  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
828 
829  switch (j)
830  {
831  // d()/dxi
832  case 0:
833  return (FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
834  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
835 
836  // d()/deta
837  case 1:
838  return (FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
839  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
840 
841  default:
842  libmesh_error_msg("Invalid j = " << j);
843  }
844  }
845 
846  default:
847  libmesh_error_msg("Invalid element type = " << type);
848  }
849  }
850 
851 
852 
853  // 3rd-order Szabo-Babuska.
854  case THIRD:
855  {
856  switch (type)
857  {
858  // Szabo-Babuska shape functions on the triangle.
859  case TRI6:
860  {
861  // Here we use finite differences to compute the derivatives!
862  const Real eps = 1.e-6;
863 
864  libmesh_assert_less (i, 10);
865  libmesh_assert_less (j, 2);
866 
867  switch (j)
868  {
869  // d()/dxi
870  case 0:
871  {
872  const Point pp(p(0)+eps, p(1));
873  const Point pm(p(0)-eps, p(1));
874 
875  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
876  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
877  }
878 
879  // d()/deta
880  case 1:
881  {
882  const Point pp(p(0), p(1)+eps);
883  const Point pm(p(0), p(1)-eps);
884 
885  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
886  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
887  }
888 
889 
890  default:
891  libmesh_error_msg("Invalid j = " << j);
892  }
893  }
894 
895 
896  // Szabo-Babuska shape functions on the quadrilateral.
897  case QUAD8:
898  case QUAD9:
899  {
900  // Compute quad shape functions as a tensor-product
901  const Real xi = p(0);
902  const Real eta = p(1);
903 
904  libmesh_assert_less (i, 16);
905 
906  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
907  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
908  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
909 
910  Real f=1.;
911 
912  switch(i)
913  {
914  case 5: // edge 0 points
915  if (elem->point(0) > elem->point(1))f = -1.;
916  break;
917  case 7: // edge 1 points
918  if (elem->point(1) > elem->point(2))f = -1.;
919  break;
920  case 9: // edge 2 points
921  if (elem->point(3) > elem->point(2))f = -1.;
922  break;
923  case 11: // edge 3 points
924  if (elem->point(0) > elem->point(3))f = -1.;
925  break;
926 
927  default:
928  // Everything else keeps f=1
929  break;
930  }
931 
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 = " << type);
952  }
953  }
954 
955 
956 
957 
958  // 4th-order Szabo-Babuska.
959  case FOURTH:
960  {
961  switch (type)
962  {
963 
964  // Szabo-Babuska shape functions on the triangle.
965  case TRI6:
966  {
967  // Here we use finite differences to compute the derivatives!
968  const Real eps = 1.e-6;
969 
970  libmesh_assert_less (i, 15);
971  libmesh_assert_less (j, 2);
972 
973  switch (j)
974  {
975  // d()/dxi
976  case 0:
977  {
978  const Point pp(p(0)+eps, p(1));
979  const Point pm(p(0)-eps, p(1));
980 
981  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
982  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
983  }
984 
985  // d()/deta
986  case 1:
987  {
988  const Point pp(p(0), p(1)+eps);
989  const Point pm(p(0), p(1)-eps);
990 
991  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
992  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
993  }
994 
995 
996  default:
997  libmesh_error_msg("Invalid j = " << j);
998  }
999  }
1000 
1001 
1002 
1003  // Szabo-Babuska shape functions on the quadrilateral.
1004  case QUAD8:
1005  case QUAD9:
1006  {
1007  // Compute quad shape functions as a tensor-product
1008  const Real xi = p(0);
1009  const Real eta = p(1);
1010 
1011  libmesh_assert_less (i, 25);
1012 
1013  // 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
1014  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};
1015  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};
1016 
1017  Real f=1.;
1018 
1019  switch(i)
1020  {
1021  case 5: // edge 0 points
1022  if (elem->point(0) > elem->point(1))f = -1.;
1023  break;
1024  case 8: // edge 1 points
1025  if (elem->point(1) > elem->point(2))f = -1.;
1026  break;
1027  case 11: // edge 2 points
1028  if (elem->point(3) > elem->point(2))f = -1.;
1029  break;
1030  case 14: // edge 3 points
1031  if (elem->point(0) > elem->point(3))f = -1.;
1032  break;
1033 
1034  default:
1035  // Everything else keeps f=1
1036  break;
1037  }
1038 
1039 
1040  switch (j)
1041  {
1042  // d()/dxi
1043  case 0:
1044  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1045  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1046 
1047  // d()/deta
1048  case 1:
1049  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1050  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1051 
1052  default:
1053  libmesh_error_msg("Invalid j = " << j);
1054  }
1055  }
1056 
1057  default:
1058  libmesh_error_msg("Invalid element type = " << type);
1059  }
1060  }
1061 
1062 
1063 
1064 
1065  // 5th-order Szabo-Babuska.
1066  case FIFTH:
1067  {
1068  // Szabo-Babuska shape functions on the quadrilateral.
1069  switch (type)
1070  {
1071 
1072  // Szabo-Babuska shape functions on the triangle.
1073  case TRI6:
1074  {
1075  // Here we use finite differences to compute the derivatives!
1076  const Real eps = 1.e-6;
1077 
1078  libmesh_assert_less (i, 21);
1079  libmesh_assert_less (j, 2);
1080 
1081  switch (j)
1082  {
1083  // d()/dxi
1084  case 0:
1085  {
1086  const Point pp(p(0)+eps, p(1));
1087  const Point pm(p(0)-eps, p(1));
1088 
1089  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1090  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1091  }
1092 
1093  // d()/deta
1094  case 1:
1095  {
1096  const Point pp(p(0), p(1)+eps);
1097  const Point pm(p(0), p(1)-eps);
1098 
1099  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1100  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1101  }
1102 
1103  default:
1104  libmesh_error_msg("Invalid j = " << j);
1105  }
1106  }
1107 
1108 
1109 
1110  case QUAD8:
1111  case QUAD9:
1112  {
1113  // Compute quad shape functions as a tensor-product
1114  const Real xi = p(0);
1115  const Real eta = p(1);
1116 
1117  libmesh_assert_less (i, 36);
1118 
1119  // 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
1120  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};
1121  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};
1122 
1123  Real f=1.;
1124 
1125  switch(i)
1126  {
1127  case 5: // edge 0 points
1128  case 7:
1129  if (elem->point(0) > elem->point(1))f = -1.;
1130  break;
1131  case 9: // edge 1 points
1132  case 11:
1133  if (elem->point(1) > elem->point(2))f = -1.;
1134  break;
1135  case 13: // edge 2 points
1136  case 15:
1137  if (elem->point(3) > elem->point(2))f = -1.;
1138  break;
1139  case 14: // edge 3 points
1140  case 19:
1141  if (elem->point(0) > elem->point(3))f = -1.;
1142  break;
1143 
1144  default:
1145  // Everything else keeps f=1
1146  break;
1147  }
1148 
1149 
1150  switch (j)
1151  {
1152  // d()/dxi
1153  case 0:
1154  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1155  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1156 
1157  // d()/deta
1158  case 1:
1159  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1160  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1161 
1162  default:
1163  libmesh_error_msg("Invalid j = " << j);
1164  }
1165  }
1166 
1167  default:
1168  libmesh_error_msg("Invalid element type = " << type);
1169  }
1170  }
1171 
1172 
1173  // 6th-order Szabo-Babuska.
1174  case SIXTH:
1175  {
1176  // Szabo-Babuska shape functions on the quadrilateral.
1177  switch (type)
1178  {
1179 
1180  // Szabo-Babuska shape functions on the triangle.
1181  case TRI6:
1182  {
1183  // Here we use finite differences to compute the derivatives!
1184  const Real eps = 1.e-6;
1185 
1186  libmesh_assert_less (i, 28);
1187  libmesh_assert_less (j, 2);
1188 
1189  switch (j)
1190  {
1191  // d()/dxi
1192  case 0:
1193  {
1194  const Point pp(p(0)+eps, p(1));
1195  const Point pm(p(0)-eps, p(1));
1196 
1197  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1198  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1199  }
1200 
1201  // d()/deta
1202  case 1:
1203  {
1204  const Point pp(p(0), p(1)+eps);
1205  const Point pm(p(0), p(1)-eps);
1206 
1207  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1208  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1209  }
1210 
1211  default:
1212  libmesh_error_msg("Invalid j = " << j);
1213  }
1214  }
1215 
1216 
1217 
1218  case QUAD8:
1219  case QUAD9:
1220  {
1221  // Compute quad shape functions as a tensor-product
1222  const Real xi = p(0);
1223  const Real eta = p(1);
1224 
1225  libmesh_assert_less (i, 49);
1226 
1227  // 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
1228  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};
1229  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};
1230 
1231  Real f=1.;
1232 
1233  switch(i)
1234  {
1235  case 5: // edge 0 points
1236  case 7:
1237  if (elem->point(0) > elem->point(1))f = -1.;
1238  break;
1239  case 10: // edge 1 points
1240  case 12:
1241  if (elem->point(1) > elem->point(2))f = -1.;
1242  break;
1243  case 15: // edge 2 points
1244  case 17:
1245  if (elem->point(3) > elem->point(2))f = -1.;
1246  break;
1247  case 20: // edge 3 points
1248  case 22:
1249  if (elem->point(0) > elem->point(3))f = -1.;
1250  break;
1251 
1252  default:
1253  // Everything else keeps f=1
1254  break;
1255  }
1256 
1257 
1258  switch (j)
1259  {
1260  // d()/dxi
1261  case 0:
1262  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1263  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1264 
1265  // d()/deta
1266  case 1:
1267  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1268  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1269 
1270  default:
1271  libmesh_error_msg("Invalid j = " << j);
1272  }
1273  }
1274 
1275  default:
1276  libmesh_error_msg("Invalid element type = " << type);
1277  }
1278  }
1279 
1280 
1281  // 7th-order Szabo-Babuska.
1282  case SEVENTH:
1283  {
1284  // Szabo-Babuska shape functions on the quadrilateral.
1285  switch (type)
1286  {
1287 
1288  // Szabo-Babuska shape functions on the triangle.
1289  case TRI6:
1290  {
1291  // Here we use finite differences to compute the derivatives!
1292  const Real eps = 1.e-6;
1293 
1294  libmesh_assert_less (i, 36);
1295  libmesh_assert_less (j, 2);
1296 
1297  switch (j)
1298  {
1299  // d()/dxi
1300  case 0:
1301  {
1302  const Point pp(p(0)+eps, p(1));
1303  const Point pm(p(0)-eps, p(1));
1304 
1305  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1306  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1307  }
1308 
1309  // d()/deta
1310  case 1:
1311  {
1312  const Point pp(p(0), p(1)+eps);
1313  const Point pm(p(0), p(1)-eps);
1314 
1315  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1316  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1317  }
1318 
1319  default:
1320  libmesh_error_msg("Invalid j = " << j);
1321  }
1322  }
1323 
1324 
1325 
1326  case QUAD8:
1327  case QUAD9:
1328  {
1329  // Compute quad shape functions as a tensor-product
1330  const Real xi = p(0);
1331  const Real eta = p(1);
1332 
1333  libmesh_assert_less (i, 64);
1334 
1335  // 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
1336  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};
1337  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};
1338 
1339  Real f=1.;
1340 
1341  switch(i)
1342  {
1343  case 5: // edge 0 points
1344  case 7:
1345  case 9:
1346  if (elem->point(0) > elem->point(1))f = -1.;
1347  break;
1348  case 11: // edge 1 points
1349  case 13:
1350  case 15:
1351  if (elem->point(1) > elem->point(2))f = -1.;
1352  break;
1353  case 17: // edge 2 points
1354  case 19:
1355  case 21:
1356  if (elem->point(3) > elem->point(2))f = -1.;
1357  break;
1358  case 23: // edge 3 points
1359  case 25:
1360  case 27:
1361  if (elem->point(0) > elem->point(3))f = -1.;
1362  break;
1363 
1364  default:
1365  // Everything else keeps f=1
1366  break;
1367  }
1368 
1369 
1370  switch (j)
1371  {
1372  // d()/dxi
1373  case 0:
1374  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1375  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1376 
1377  // d()/deta
1378  case 1:
1379  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1380  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1381 
1382  default:
1383  libmesh_error_msg("Invalid j = " << j);
1384  }
1385  }
1386 
1387  default:
1388  libmesh_error_msg("Invalid element type = " << type);
1389  }
1390  }
1391 
1392 
1393 
1394  // by default throw an error;call the orientation-independent shape functions
1395  default:
1396  libmesh_error_msg("ERROR: Unsupported polynomial order!");
1397  }
1398 
1399  libmesh_error_msg("We'll never get here!");
1400  return 0.;
1401 }
1402 
1403 
1404 
1405 template <>
1407  const Order,
1408  const unsigned int,
1409  const unsigned int,
1410  const Point &)
1411 {
1412  static bool warning_given = false;
1413 
1414  if (!warning_given)
1415  libMesh::err << "Second derivatives for Szabab elements "
1416  << " are not yet implemented!"
1417  << std::endl;
1418 
1419  warning_given = true;
1420  return 0.;
1421 }
1422 
1423 
1424 
1425 template <>
1427  const Order,
1428  const unsigned int,
1429  const unsigned int,
1430  const Point &)
1431 {
1432  static bool warning_given = false;
1433 
1434  if (!warning_given)
1435  libMesh::err << "Second derivatives for Szabab elements "
1436  << " are not yet implemented!"
1437  << std::endl;
1438 
1439  warning_given = true;
1440  return 0.;
1441 }
1442 
1443 } // namespace libMesh
1444 
1445 
1446 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
OStreamProxy err
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
A specific instantiation of the FEBase class.
Definition: fe.h:89
T pow(const T &x)
Definition: utility.h:195
PetscErrorCode Vec x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)