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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo