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