fe_szabab_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 <>
53 Real FE<2,SZABAB>::shape(const ElemType,
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 <>
65 Real FE<2,SZABAB>::shape(const Elem * elem,
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 
732 
733 
734 
735 
736 template <>
737 Real FE<2,SZABAB>::shape_deriv(const ElemType,
738  const Order,
739  const unsigned int,
740  const unsigned int,
741  const Point &)
742 {
743  libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
744  return 0.;
745 }
746 
747 
748 
749 template <>
750 Real FE<2,SZABAB>::shape_deriv(const Elem * elem,
751  const Order order,
752  const unsigned int i,
753  const unsigned int j,
754  const Point & p)
755 {
756  libmesh_assert(elem);
757 
758  const ElemType type = elem->type();
759 
760  const Order totalorder = static_cast<Order>(order + elem->p_level());
761 
762  switch (totalorder)
763  {
764 
765  // 1st & 2nd-order Szabo-Babuska.
766  case FIRST:
767  case SECOND:
768  {
769  switch (type)
770  {
771 
772  // Szabo-Babuska shape functions on the triangle.
773  case TRI3:
774  case TRI6:
775  {
776  // Here we use finite differences to compute the derivatives!
777  const Real eps = 1.e-6;
778 
779  libmesh_assert_less (i, 6);
780  libmesh_assert_less (j, 2);
781 
782  switch (j)
783  {
784  // d()/dxi
785  case 0:
786  {
787  const Point pp(p(0)+eps, p(1));
788  const Point pm(p(0)-eps, p(1));
789 
790  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
791  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
792  }
793 
794  // d()/deta
795  case 1:
796  {
797  const Point pp(p(0), p(1)+eps);
798  const Point pm(p(0), p(1)-eps);
799 
800  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
801  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
802  }
803 
804  default:
805  libmesh_error_msg("Invalid j = " << j);
806  }
807  }
808 
809 
810 
811  // Szabo-Babuska shape functions on the quadrilateral.
812  case QUAD4:
813  case QUAD8:
814  case QUAD9:
815  {
816  // Compute quad shape functions as a tensor-product
817  const Real xi = p(0);
818  const Real eta = p(1);
819 
820  libmesh_assert_less (i, 9);
821 
822  // 0 1 2 3 4 5 6 7 8
823  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
824  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
825 
826  switch (j)
827  {
828  // d()/dxi
829  case 0:
830  return (FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
831  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
832 
833  // d()/deta
834  case 1:
835  return (FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
836  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
837 
838  default:
839  libmesh_error_msg("Invalid j = " << j);
840  }
841  }
842 
843  default:
844  libmesh_error_msg("Invalid element type = " << type);
845  }
846  }
847 
848 
849 
850  // 3rd-order Szabo-Babuska.
851  case THIRD:
852  {
853  switch (type)
854  {
855  // Szabo-Babuska shape functions on the triangle.
856  case TRI6:
857  {
858  // Here we use finite differences to compute the derivatives!
859  const Real eps = 1.e-6;
860 
861  libmesh_assert_less (i, 10);
862  libmesh_assert_less (j, 2);
863 
864  switch (j)
865  {
866  // d()/dxi
867  case 0:
868  {
869  const Point pp(p(0)+eps, p(1));
870  const Point pm(p(0)-eps, p(1));
871 
872  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
873  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
874  }
875 
876  // d()/deta
877  case 1:
878  {
879  const Point pp(p(0), p(1)+eps);
880  const Point pm(p(0), p(1)-eps);
881 
882  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
883  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
884  }
885 
886 
887  default:
888  libmesh_error_msg("Invalid j = " << j);
889  }
890  }
891 
892 
893  // Szabo-Babuska shape functions on the quadrilateral.
894  case QUAD8:
895  case QUAD9:
896  {
897  // Compute quad shape functions as a tensor-product
898  const Real xi = p(0);
899  const Real eta = p(1);
900 
901  libmesh_assert_less (i, 16);
902 
903  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
904  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
905  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
906 
907  Real f=1.;
908 
909  switch(i)
910  {
911  case 5: // edge 0 points
912  if (elem->point(0) > elem->point(1))f = -1.;
913  break;
914  case 7: // edge 1 points
915  if (elem->point(1) > elem->point(2))f = -1.;
916  break;
917  case 9: // edge 2 points
918  if (elem->point(3) > elem->point(2))f = -1.;
919  break;
920  case 11: // edge 3 points
921  if (elem->point(0) > elem->point(3))f = -1.;
922  break;
923 
924  default:
925  // Everything else keeps f=1
926  break;
927  }
928 
929 
930  switch (j)
931  {
932  // d()/dxi
933  case 0:
934  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
935  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
936 
937  // d()/deta
938  case 1:
939  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
940  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
941 
942  default:
943  libmesh_error_msg("Invalid j = " << j);
944  }
945  }
946 
947  default:
948  libmesh_error_msg("Invalid element type = " << type);
949  }
950  }
951 
952 
953 
954 
955  // 4th-order Szabo-Babuska.
956  case FOURTH:
957  {
958  switch (type)
959  {
960 
961  // Szabo-Babuska shape functions on the triangle.
962  case TRI6:
963  {
964  // Here we use finite differences to compute the derivatives!
965  const Real eps = 1.e-6;
966 
967  libmesh_assert_less (i, 15);
968  libmesh_assert_less (j, 2);
969 
970  switch (j)
971  {
972  // d()/dxi
973  case 0:
974  {
975  const Point pp(p(0)+eps, p(1));
976  const Point pm(p(0)-eps, p(1));
977 
978  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
979  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
980  }
981 
982  // d()/deta
983  case 1:
984  {
985  const Point pp(p(0), p(1)+eps);
986  const Point pm(p(0), p(1)-eps);
987 
988  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
989  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
990  }
991 
992 
993  default:
994  libmesh_error_msg("Invalid j = " << j);
995  }
996  }
997 
998 
999 
1000  // Szabo-Babuska shape functions on the quadrilateral.
1001  case QUAD8:
1002  case QUAD9:
1003  {
1004  // Compute quad shape functions as a tensor-product
1005  const Real xi = p(0);
1006  const Real eta = p(1);
1007 
1008  libmesh_assert_less (i, 25);
1009 
1010  // 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
1011  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};
1012  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};
1013 
1014  Real f=1.;
1015 
1016  switch(i)
1017  {
1018  case 5: // edge 0 points
1019  if (elem->point(0) > elem->point(1))f = -1.;
1020  break;
1021  case 8: // edge 1 points
1022  if (elem->point(1) > elem->point(2))f = -1.;
1023  break;
1024  case 11: // edge 2 points
1025  if (elem->point(3) > elem->point(2))f = -1.;
1026  break;
1027  case 14: // edge 3 points
1028  if (elem->point(0) > elem->point(3))f = -1.;
1029  break;
1030 
1031  default:
1032  // Everything else keeps f=1
1033  break;
1034  }
1035 
1036 
1037  switch (j)
1038  {
1039  // d()/dxi
1040  case 0:
1041  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1042  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1043 
1044  // d()/deta
1045  case 1:
1046  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1047  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1048 
1049  default:
1050  libmesh_error_msg("Invalid j = " << j);
1051  }
1052  }
1053 
1054  default:
1055  libmesh_error_msg("Invalid element type = " << type);
1056  }
1057  }
1058 
1059 
1060 
1061 
1062  // 5th-order Szabo-Babuska.
1063  case FIFTH:
1064  {
1065  // Szabo-Babuska shape functions on the quadrilateral.
1066  switch (type)
1067  {
1068 
1069  // Szabo-Babuska shape functions on the triangle.
1070  case TRI6:
1071  {
1072  // Here we use finite differences to compute the derivatives!
1073  const Real eps = 1.e-6;
1074 
1075  libmesh_assert_less (i, 21);
1076  libmesh_assert_less (j, 2);
1077 
1078  switch (j)
1079  {
1080  // d()/dxi
1081  case 0:
1082  {
1083  const Point pp(p(0)+eps, p(1));
1084  const Point pm(p(0)-eps, p(1));
1085 
1086  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1087  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1088  }
1089 
1090  // d()/deta
1091  case 1:
1092  {
1093  const Point pp(p(0), p(1)+eps);
1094  const Point pm(p(0), p(1)-eps);
1095 
1096  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1097  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1098  }
1099 
1100  default:
1101  libmesh_error_msg("Invalid j = " << j);
1102  }
1103  }
1104 
1105 
1106 
1107  case QUAD8:
1108  case QUAD9:
1109  {
1110  // Compute quad shape functions as a tensor-product
1111  const Real xi = p(0);
1112  const Real eta = p(1);
1113 
1114  libmesh_assert_less (i, 36);
1115 
1116  // 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
1117  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};
1118  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};
1119 
1120  Real f=1.;
1121 
1122  switch(i)
1123  {
1124  case 5: // edge 0 points
1125  case 7:
1126  if (elem->point(0) > elem->point(1))f = -1.;
1127  break;
1128  case 9: // edge 1 points
1129  case 11:
1130  if (elem->point(1) > elem->point(2))f = -1.;
1131  break;
1132  case 13: // edge 2 points
1133  case 15:
1134  if (elem->point(3) > elem->point(2))f = -1.;
1135  break;
1136  case 14: // edge 3 points
1137  case 19:
1138  if (elem->point(0) > elem->point(3))f = -1.;
1139  break;
1140 
1141  default:
1142  // Everything else keeps f=1
1143  break;
1144  }
1145 
1146 
1147  switch (j)
1148  {
1149  // d()/dxi
1150  case 0:
1151  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1152  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1153 
1154  // d()/deta
1155  case 1:
1156  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1157  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1158 
1159  default:
1160  libmesh_error_msg("Invalid j = " << j);
1161  }
1162  }
1163 
1164  default:
1165  libmesh_error_msg("Invalid element type = " << type);
1166  }
1167  }
1168 
1169 
1170  // 6th-order Szabo-Babuska.
1171  case SIXTH:
1172  {
1173  // Szabo-Babuska shape functions on the quadrilateral.
1174  switch (type)
1175  {
1176 
1177  // Szabo-Babuska shape functions on the triangle.
1178  case TRI6:
1179  {
1180  // Here we use finite differences to compute the derivatives!
1181  const Real eps = 1.e-6;
1182 
1183  libmesh_assert_less (i, 28);
1184  libmesh_assert_less (j, 2);
1185 
1186  switch (j)
1187  {
1188  // d()/dxi
1189  case 0:
1190  {
1191  const Point pp(p(0)+eps, p(1));
1192  const Point pm(p(0)-eps, p(1));
1193 
1194  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1195  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1196  }
1197 
1198  // d()/deta
1199  case 1:
1200  {
1201  const Point pp(p(0), p(1)+eps);
1202  const Point pm(p(0), p(1)-eps);
1203 
1204  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1205  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1206  }
1207 
1208  default:
1209  libmesh_error_msg("Invalid j = " << j);
1210  }
1211  }
1212 
1213 
1214 
1215  case QUAD8:
1216  case QUAD9:
1217  {
1218  // Compute quad shape functions as a tensor-product
1219  const Real xi = p(0);
1220  const Real eta = p(1);
1221 
1222  libmesh_assert_less (i, 49);
1223 
1224  // 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
1225  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};
1226  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};
1227 
1228  Real f=1.;
1229 
1230  switch(i)
1231  {
1232  case 5: // edge 0 points
1233  case 7:
1234  if (elem->point(0) > elem->point(1))f = -1.;
1235  break;
1236  case 10: // edge 1 points
1237  case 12:
1238  if (elem->point(1) > elem->point(2))f = -1.;
1239  break;
1240  case 15: // edge 2 points
1241  case 17:
1242  if (elem->point(3) > elem->point(2))f = -1.;
1243  break;
1244  case 20: // edge 3 points
1245  case 22:
1246  if (elem->point(0) > elem->point(3))f = -1.;
1247  break;
1248 
1249  default:
1250  // Everything else keeps f=1
1251  break;
1252  }
1253 
1254 
1255  switch (j)
1256  {
1257  // d()/dxi
1258  case 0:
1259  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1260  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1261 
1262  // d()/deta
1263  case 1:
1264  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1265  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1266 
1267  default:
1268  libmesh_error_msg("Invalid j = " << j);
1269  }
1270  }
1271 
1272  default:
1273  libmesh_error_msg("Invalid element type = " << type);
1274  }
1275  }
1276 
1277 
1278  // 7th-order Szabo-Babuska.
1279  case SEVENTH:
1280  {
1281  // Szabo-Babuska shape functions on the quadrilateral.
1282  switch (type)
1283  {
1284 
1285  // Szabo-Babuska shape functions on the triangle.
1286  case TRI6:
1287  {
1288  // Here we use finite differences to compute the derivatives!
1289  const Real eps = 1.e-6;
1290 
1291  libmesh_assert_less (i, 36);
1292  libmesh_assert_less (j, 2);
1293 
1294  switch (j)
1295  {
1296  // d()/dxi
1297  case 0:
1298  {
1299  const Point pp(p(0)+eps, p(1));
1300  const Point pm(p(0)-eps, p(1));
1301 
1302  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1303  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1304  }
1305 
1306  // d()/deta
1307  case 1:
1308  {
1309  const Point pp(p(0), p(1)+eps);
1310  const Point pm(p(0), p(1)-eps);
1311 
1312  return (FE<2,SZABAB>::shape(elem, order, i, pp) -
1313  FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
1314  }
1315 
1316  default:
1317  libmesh_error_msg("Invalid j = " << j);
1318  }
1319  }
1320 
1321 
1322 
1323  case QUAD8:
1324  case QUAD9:
1325  {
1326  // Compute quad shape functions as a tensor-product
1327  const Real xi = p(0);
1328  const Real eta = p(1);
1329 
1330  libmesh_assert_less (i, 64);
1331 
1332  // 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
1333  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};
1334  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};
1335 
1336  Real f=1.;
1337 
1338  switch(i)
1339  {
1340  case 5: // edge 0 points
1341  case 7:
1342  case 9:
1343  if (elem->point(0) > elem->point(1))f = -1.;
1344  break;
1345  case 11: // edge 1 points
1346  case 13:
1347  case 15:
1348  if (elem->point(1) > elem->point(2))f = -1.;
1349  break;
1350  case 17: // edge 2 points
1351  case 19:
1352  case 21:
1353  if (elem->point(3) > elem->point(2))f = -1.;
1354  break;
1355  case 23: // edge 3 points
1356  case 25:
1357  case 27:
1358  if (elem->point(0) > elem->point(3))f = -1.;
1359  break;
1360 
1361  default:
1362  // Everything else keeps f=1
1363  break;
1364  }
1365 
1366 
1367  switch (j)
1368  {
1369  // d()/dxi
1370  case 0:
1371  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1372  FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1373 
1374  // d()/deta
1375  case 1:
1376  return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1377  FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1378 
1379  default:
1380  libmesh_error_msg("Invalid j = " << j);
1381  }
1382  }
1383 
1384  default:
1385  libmesh_error_msg("Invalid element type = " << type);
1386  }
1387  }
1388 
1389 
1390 
1391  // by default throw an error;call the orientation-independent shape functions
1392  default:
1393  libmesh_error_msg("ERROR: Unsupported polynomial order!");
1394  }
1395 }
1396 
1397 
1398 
1399 template <>
1400 Real FE<2,SZABAB>::shape_second_deriv(const ElemType,
1401  const Order,
1402  const unsigned int,
1403  const unsigned int,
1404  const Point &)
1405 {
1406  static bool warning_given = false;
1407 
1408  if (!warning_given)
1409  libMesh::err << "Second derivatives for Szabab elements "
1410  << " are not yet implemented!"
1411  << std::endl;
1412 
1413  warning_given = true;
1414  return 0.;
1415 }
1416 
1417 
1418 
1419 template <>
1420 Real FE<2,SZABAB>::shape_second_deriv(const Elem *,
1421  const Order,
1422  const unsigned int,
1423  const unsigned int,
1424  const Point &)
1425 {
1426  static bool warning_given = false;
1427 
1428  if (!warning_given)
1429  libMesh::err << "Second derivatives for Szabab elements "
1430  << " are not yet implemented!"
1431  << std::endl;
1432 
1433  warning_given = true;
1434  return 0.;
1435 }
1436 
1437 } // namespace libMesh
1438 
1439 
1440 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
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:2555
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:198
libMesh::err
OStreamProxy err(std::cerr)
Definition: libmesh_common.h:229
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libmesh_config.h
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
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 by doxygen