Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : // Local includes
19 : #include "libmesh/libmesh_config.h"
20 :
21 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22 :
23 : // libmesh includes
24 : #include "libmesh/fe.h"
25 : #include "libmesh/elem.h"
26 : #include "libmesh/utility.h"
27 : #include "libmesh/enum_to_string.h"
28 :
29 : // C++ includes
30 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
31 : #include <cmath> // for std::sqrt
32 :
33 : // Anonymous namespace to hold static std::sqrt values
34 : namespace
35 : {
36 : using namespace libMesh;
37 :
38 : static const Real sqrt2 = std::sqrt(2.);
39 : static const Real sqrt6 = std::sqrt(6.);
40 : static const Real sqrt10 = std::sqrt(10.);
41 : static const Real sqrt14 = std::sqrt(14.);
42 : static const Real sqrt22 = std::sqrt(22.);
43 : static const Real sqrt26 = std::sqrt(26.);
44 :
45 :
46 : // Take care of edge orientation, this is needed at edge shapes with
47 : // (y=0)-asymmetric 1D shapes
48 645612 : Real quad_flip(const Elem * elem,
49 : const unsigned int totalorder,
50 : const unsigned int i)
51 : {
52 53801 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
53 :
54 : // vertex and interior shape functions don't flip
55 645612 : if (i < 4 || i >= 4*totalorder)
56 27597 : return 1;
57 :
58 314448 : const int edge = (i-4) / (totalorder-1);
59 26204 : libmesh_assert_less(edge, 4);
60 26204 : libmesh_assert_greater_equal(edge, 0);
61 :
62 314448 : const int edge_i = i - 4 - edge*(totalorder-1);
63 :
64 : // The "natural" orientation is p1>p0 on edge 0,
65 : // p2>p1 on e1, p2>p3 on e2, p3>p0 on e3
66 324764 : if (edge_i%2 &&
67 123792 : (elem->positive_edge_orientation(edge) ==
68 176200 : (edge < 2)))
69 10316 : return -1;
70 :
71 21046 : return 1;
72 : }
73 :
74 : } // anonymous namespace
75 :
76 :
77 : namespace libMesh
78 : {
79 :
80 :
81 1192112 : LIBMESH_DEFAULT_VECTORIZED_FE(2,SZABAB)
82 :
83 :
84 : template <>
85 2002338 : Real FE<2,SZABAB>::shape(const Elem * elem,
86 : const Order order,
87 : const unsigned int i,
88 : const Point & p,
89 : const bool add_p_level)
90 : {
91 175515 : libmesh_assert(elem);
92 :
93 2002338 : const ElemType type = elem->type();
94 :
95 2353368 : const Order totalorder = order + add_p_level*elem->p_level();
96 :
97 : // Declare that we are using our own special power function
98 : // from the Utility namespace. This saves typing later.
99 : using Utility::pow;
100 :
101 2002338 : switch (totalorder)
102 : {
103 : // 1st & 2nd-order Szabo-Babuska.
104 365742 : case FIRST:
105 : case SECOND:
106 : {
107 333740 : switch (type)
108 : {
109 :
110 : // Szabo-Babuska shape functions on the triangle.
111 302298 : case TRI3:
112 : case TRI6:
113 : case TRI7:
114 : {
115 302298 : const Real l1 = 1-p(0)-p(1);
116 26715 : const Real l2 = p(0);
117 26715 : const Real l3 = p(1);
118 :
119 26715 : libmesh_assert_less (i, 6);
120 :
121 302298 : switch (i)
122 : {
123 5139 : case 0: return l1;
124 58002 : case 1: return l2;
125 58002 : case 2: return l3;
126 :
127 42764 : case 3: return l1*l2*(-4.*sqrt6);
128 42764 : case 4: return l2*l3*(-4.*sqrt6);
129 42764 : case 5: return l3*l1*(-4.*sqrt6);
130 :
131 0 : default:
132 0 : libmesh_error_msg("Invalid i = " << i);
133 : }
134 : }
135 :
136 :
137 : // Szabo-Babuska shape functions on the quadrilateral.
138 63444 : case QUAD4:
139 : case QUAD8:
140 : case QUAD9:
141 : {
142 : // Compute quad shape functions as a tensor-product
143 63444 : const Real xi = p(0);
144 63444 : const Real eta = p(1);
145 :
146 5287 : libmesh_assert_less (i, 9);
147 :
148 : // 0 1 2 3 4 5 6 7 8
149 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
150 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
151 :
152 63444 : return (FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
153 63444 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
154 :
155 : }
156 :
157 0 : default:
158 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
159 : }
160 : }
161 :
162 :
163 : // 3rd-order Szabo-Babuska.
164 578136 : case THIRD:
165 : {
166 578136 : switch (type)
167 : {
168 :
169 : // Szabo-Babuska shape functions on the triangle.
170 463320 : case TRI6:
171 : case TRI7:
172 : {
173 463320 : Real l1 = 1-p(0)-p(1);
174 41220 : Real l2 = p(0);
175 41220 : Real l3 = p(1);
176 :
177 41220 : Real f=1;
178 :
179 41220 : libmesh_assert_less (i, 10);
180 :
181 :
182 463320 : if (i==4 && elem->positive_edge_orientation(0))f=-1;
183 463320 : if (i==6 && elem->positive_edge_orientation(1))f=-1;
184 463320 : if (i==8 && elem->positive_edge_orientation(2))f=-1;
185 :
186 :
187 463320 : switch (i)
188 : {
189 : //nodal modes
190 4122 : case 0: return l1;
191 46332 : case 1: return l2;
192 46332 : case 2: return l3;
193 :
194 : //side modes
195 46332 : case 3: return l1*l2*(-4.*sqrt6);
196 46332 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
197 :
198 46332 : case 5: return l2*l3*(-4.*sqrt6);
199 46332 : case 6: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
200 :
201 46332 : case 7: return l3*l1*(-4.*sqrt6);
202 46332 : case 8: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
203 :
204 : //internal modes
205 46332 : case 9: return l1*l2*l3;
206 :
207 0 : default:
208 0 : libmesh_error_msg("Invalid i = " << i);
209 : }
210 : }
211 :
212 :
213 : // Szabo-Babuska shape functions on the quadrilateral.
214 114816 : case QUAD8:
215 : case QUAD9:
216 : {
217 : // Compute quad shape functions as a tensor-product
218 114816 : const Real xi = p(0);
219 114816 : const Real eta = p(1);
220 :
221 9568 : libmesh_assert_less (i, 16);
222 :
223 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
224 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
225 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
226 :
227 114816 : const Real f = quad_flip(elem, totalorder, i);
228 :
229 114816 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
230 114816 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
231 : }
232 :
233 0 : default:
234 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
235 : }
236 : }
237 :
238 :
239 :
240 :
241 : // 4th-order Szabo-Babuska.
242 1058460 : case FOURTH:
243 : {
244 1058460 : switch (type)
245 : {
246 : // Szabo-Babuska shape functions on the triangle.
247 833160 : case TRI6:
248 : case TRI7:
249 : {
250 833160 : Real l1 = 1-p(0)-p(1);
251 73950 : Real l2 = p(0);
252 73950 : Real l3 = p(1);
253 :
254 73950 : Real f=1;
255 :
256 73950 : libmesh_assert_less (i, 15);
257 :
258 :
259 833160 : if (i== 4 && elem->positive_edge_orientation(0))f=-1;
260 833160 : if (i== 7 && elem->positive_edge_orientation(1))f=-1;
261 833160 : if (i==10 && elem->positive_edge_orientation(2))f=-1;
262 :
263 :
264 833160 : switch (i)
265 : {
266 : //nodal modes
267 4930 : case 0: return l1;
268 55544 : case 1: return l2;
269 55544 : case 2: return l3;
270 :
271 : //side modes
272 55544 : case 3: return l1*l2*(-4.*sqrt6);
273 55544 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
274 60474 : case 5: return l1*l2*(-sqrt14)*(5.*pow<2>(l2-l1)-1);
275 :
276 55544 : case 6: return l2*l3*(-4.*sqrt6);
277 55544 : case 7: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
278 60474 : case 8: return l2*l3*(-sqrt14)*(5.*pow<2>(l3-l2)-1);
279 :
280 55544 : case 9: return l3*l1*(-4.*sqrt6);
281 55544 : case 10: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
282 60474 : case 11: return l3*l1*(-sqrt14)*(5.*pow<2>(l1-l3)-1);
283 :
284 : //internal modes
285 55544 : case 12: return l1*l2*l3;
286 :
287 55544 : case 13: return l1*l2*l3*(l2-l1);
288 55544 : case 14: return l1*l2*l3*(2*l3-1);
289 :
290 0 : default:
291 0 : libmesh_error_msg("Invalid i = " << i);
292 : }
293 : }
294 :
295 :
296 : // Szabo-Babuska shape functions on the quadrilateral.
297 225300 : case QUAD8:
298 : case QUAD9:
299 : {
300 : // Compute quad shape functions as a tensor-product
301 225300 : const Real xi = p(0);
302 225300 : const Real eta = p(1);
303 :
304 18775 : 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 225300 : const Real f = quad_flip(elem, totalorder, i);
311 :
312 225300 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
313 225300 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
314 : }
315 :
316 0 : default:
317 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
318 : }
319 : }
320 :
321 :
322 :
323 :
324 : // 5th-order Szabo-Babuska.
325 0 : case FIFTH:
326 : {
327 0 : switch (type)
328 : {
329 : // Szabo-Babuska shape functions on the triangle.
330 0 : case TRI6:
331 : case TRI7:
332 : {
333 0 : Real l1 = 1-p(0)-p(1);
334 0 : Real l2 = p(0);
335 0 : Real l3 = p(1);
336 :
337 0 : const Real x=l2-l1;
338 0 : const Real y=2.*l3-1;
339 :
340 0 : Real f=1;
341 :
342 0 : libmesh_assert_less (i, 21);
343 :
344 :
345 0 : if ((i== 4||i== 6) && elem->positive_edge_orientation(0))f=-1;
346 0 : if ((i== 8||i==10) && elem->positive_edge_orientation(1))f=-1;
347 0 : if ((i==12||i==14) && elem->positive_edge_orientation(2))f=-1;
348 :
349 :
350 0 : switch (i)
351 : {
352 : //nodal modes
353 0 : case 0: return l1;
354 0 : case 1: return l2;
355 0 : case 2: return l3;
356 :
357 : //side modes
358 0 : case 3: return l1*l2*(-4.*sqrt6);
359 0 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
360 0 : case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
361 0 : case 6: return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
362 :
363 0 : case 7: return l2*l3*(-4.*sqrt6);
364 0 : case 8: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
365 0 : case 9: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
366 0 : case 10: return -f*sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
367 :
368 0 : case 11: return l3*l1*(-4.*sqrt6);
369 0 : case 12: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
370 0 : case 13: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
371 0 : case 14: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
372 :
373 : //internal modes
374 0 : case 15: return l1*l2*l3;
375 :
376 0 : case 16: return l1*l2*l3*x;
377 0 : case 17: return l1*l2*l3*y;
378 :
379 0 : case 18: return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
380 0 : case 19: return l1*l2*l3*x*y;
381 0 : case 20: return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
382 :
383 0 : default:
384 0 : libmesh_error_msg("Invalid i = " << i);
385 : }
386 : } // case TRI6/TRI7
387 :
388 : // Szabo-Babuska shape functions on the quadrilateral.
389 0 : case QUAD8:
390 : case QUAD9:
391 : {
392 : // Compute quad shape functions as a tensor-product
393 0 : const Real xi = p(0);
394 0 : const Real eta = p(1);
395 :
396 0 : libmesh_assert_less (i, 36);
397 :
398 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
399 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
400 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
401 :
402 0 : const Real f = quad_flip(elem, totalorder, i);
403 :
404 0 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
405 0 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
406 :
407 : } // case QUAD8/QUAD9
408 :
409 0 : default:
410 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
411 :
412 : } // switch type
413 :
414 : } // case FIFTH
415 :
416 : // 6th-order Szabo-Babuska.
417 0 : case SIXTH:
418 : {
419 0 : switch (type)
420 : {
421 : // Szabo-Babuska shape functions on the triangle.
422 0 : case TRI6:
423 : case TRI7:
424 : {
425 0 : Real l1 = 1-p(0)-p(1);
426 0 : Real l2 = p(0);
427 0 : Real l3 = p(1);
428 :
429 0 : const Real x=l2-l1;
430 0 : const Real y=2.*l3-1;
431 :
432 0 : Real f=1;
433 :
434 0 : libmesh_assert_less (i, 28);
435 :
436 :
437 0 : if ((i== 4||i== 6) && elem->positive_edge_orientation(0))f=-1;
438 0 : if ((i== 9||i==11) && elem->positive_edge_orientation(1))f=-1;
439 0 : if ((i==14||i==16) && elem->positive_edge_orientation(2))f=-1;
440 :
441 :
442 0 : switch (i)
443 : {
444 : //nodal modes
445 0 : case 0: return l1;
446 0 : case 1: return l2;
447 0 : case 2: return l3;
448 :
449 : //side modes
450 0 : case 3: return l1*l2*(-4.*sqrt6);
451 0 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
452 0 : case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
453 0 : case 6: return f*(-sqrt2)*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
454 0 : case 7: return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
455 :
456 0 : case 8: return l2*l3*(-4.*sqrt6);
457 0 : case 9: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
458 0 : case 10: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
459 0 : case 11: return f*(-sqrt2)*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
460 0 : case 12: return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
461 :
462 0 : case 13: return l3*l1*(-4.*sqrt6);
463 0 : case 14: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
464 0 : case 15: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
465 0 : case 16: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
466 0 : case 17: return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
467 :
468 :
469 :
470 : //internal modes
471 0 : case 18: return l1*l2*l3;
472 :
473 0 : case 19: return l1*l2*l3*x;
474 0 : case 20: return l1*l2*l3*y;
475 :
476 0 : case 21: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2);
477 0 : case 22: return l1*l2*l3*(l2-l1)*(2.0*l3-1.0);
478 0 : case 23: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3);
479 0 : case 24: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
480 0 : case 25: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
481 0 : case 26: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
482 0 : case 27: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
483 :
484 :
485 0 : default:
486 0 : libmesh_error_msg("Invalid i = " << i);
487 : }
488 : } // case TRI6/TRI7
489 :
490 : // Szabo-Babuska shape functions on the quadrilateral.
491 0 : case QUAD8:
492 : case QUAD9:
493 : {
494 : // Compute quad shape functions as a tensor-product
495 0 : const Real xi = p(0);
496 0 : const Real eta = p(1);
497 :
498 0 : libmesh_assert_less (i, 49);
499 :
500 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
501 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
502 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
503 :
504 0 : const Real f = quad_flip(elem, totalorder, i);
505 :
506 0 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
507 0 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
508 :
509 : } // case QUAD8/QUAD9
510 :
511 0 : default:
512 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
513 :
514 : } // switch type
515 :
516 : } // case SIXTH
517 :
518 :
519 : // 7th-order Szabo-Babuska.
520 0 : case SEVENTH:
521 : {
522 0 : switch (type)
523 : {
524 : // Szabo-Babuska shape functions on the triangle.
525 0 : case TRI6:
526 : case TRI7:
527 : {
528 :
529 0 : Real l1 = 1-p(0)-p(1);
530 0 : Real l2 = p(0);
531 0 : Real l3 = p(1);
532 :
533 0 : const Real x=l2-l1;
534 0 : const Real y=2.*l3-1.;
535 :
536 0 : Real f=1;
537 :
538 0 : libmesh_assert_less (i, 36);
539 :
540 :
541 0 : if ((i>= 4&&i<= 8) && elem->positive_edge_orientation(0))f=-1;
542 0 : if ((i>=10&&i<=14) && elem->positive_edge_orientation(1))f=-1;
543 0 : if ((i>=16&&i<=20) && elem->positive_edge_orientation(2))f=-1;
544 :
545 :
546 0 : switch (i)
547 : {
548 : //nodal modes
549 0 : case 0: return l1;
550 0 : case 1: return l2;
551 0 : case 2: return l3;
552 :
553 : //side modes
554 0 : case 3: return l1*l2*(-4.*sqrt6);
555 0 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
556 :
557 0 : case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
558 0 : case 6: return f*-sqrt2*l1*l2*((9.0-21.0*l1*l1)*l1+(-9.0+63.0*l1*l1+(-63.0*l1+21.0*l2)*l2)*l2);
559 0 : case 7: return -sqrt22*l1*l2*(0.5+(-7.0+0.105E2*l1*l1)*l1*l1+((14.0-0.42E2*l1*l1)*l1+(-7.0+0.63E2*l1*l1+(-0.42E2*l1+0.105E2*l2)*l2)*l2)*l2);
560 0 : case 8: return f*-sqrt26*l1*l2*((-0.25E1+(15.0-0.165E2*l1*l1)*l1*l1)*l1+(0.25E1+(-45.0+0.825E2*l1*l1)*l1*l1+((45.0-0.165E3*l1*l1)*l1+(-15.0+0.165E3*l1*l1+(-0.825E2*l1+0.165E2*l2)*l2)*l2)*l2)*l2);
561 :
562 0 : case 9: return l2*l3*(-4.*sqrt6);
563 0 : case 10: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
564 :
565 0 : case 11: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
566 0 : case 12: return f*-sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
567 0 : case 13: return -sqrt22*l2*l3*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l2)*l2)*l2)*l2);
568 0 : case 14: return f*-sqrt26*l2*l3*((0.25E1+(-15.0+0.165E2*l3*l3)*l3*l3)*l3+(-0.25E1+(45.0-0.825E2*l3*l3)*l3*l3+((-45.0+0.165E3*l3*l3)*l3+(15.0-0.165E3*l3*l3+(0.825E2*l3-0.165E2*l2)*l2)*l2)*l2)*l2);
569 :
570 0 : case 15: return l3*l1*(-4.*sqrt6);
571 0 : case 16: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
572 :
573 0 : case 17: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
574 0 : case 18: return -f*sqrt2*l3*l1*((9.-21.*l3*l3)*l3+(-9.+63.*l3*l3+(-63.*l3+21.*l1)*l1)*l1);
575 0 : case 19: return -sqrt22*l3*l1*(0.5+(-7.0+0.105E2*l3*l3)*l3*l3+((14.0-0.42E2*l3*l3)*l3+(-7.0+0.63E2*l3*l3+(-0.42E2*l3+0.105E2*l1)*l1)*l1)*l1);
576 0 : case 20: return f*-sqrt26*l3*l1*((-0.25E1+(15.0-0.165E2*l3*l3)*l3*l3)*l3+(0.25E1+(-45.0+0.825E2*l3*l3)*l3*l3+((45.0-0.165E3*l3*l3)*l3+(-15.0+0.165E3*l3*l3+(-0.825E2*l3+0.165E2*l1)*l1)*l1)*l1)*l1);
577 :
578 :
579 : //internal modes
580 0 : case 21: return l1*l2*l3;
581 :
582 0 : case 22: return l1*l2*l3*x;
583 0 : case 23: return l1*l2*l3*y;
584 :
585 0 : case 24: return l1*l2*l3*0.5*(3.*pow<2>(x)-1.);
586 0 : case 25: return l1*l2*l3*x*y;
587 0 : case 26: return l1*l2*l3*0.5*(3.*pow<2>(y)-1.);
588 :
589 0 : case 27: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2);
590 0 : case 28: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0*l3-1.0);
591 0 : case 29: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
592 0 : case 30: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
593 0 : case 31: return 0.125*l1*l2*l3*((-15.0+(70.0-63.0*l1*l1)*l1*l1)*l1+(15.0+(-210.0+315.0*l1*l1)*l1*l1+((210.0-630.0*l1*l1)*l1+(-70.0+630.0*l1*l1+(-315.0*l1+63.0*l2)*l2)*l2)*l2)*l2);
594 0 : case 32: return 0.5*l1*l2*l3*((3.0-5.0*l1*l1)*l1+(-3.0+15.0*l1*l1+(-15.0*l1+5.0*l2)*l2)*l2)*(2.0*l3-1.0);
595 0 : case 33: return 0.25*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2)*(2.0+(-12.0+12.0*l3)*l3);
596 0 : case 34: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3)*(l2-l1);
597 0 : case 35: return 0.125*l1*l2*l3*(-8.0+(240.0+(-1680.0+(4480.0+(-5040.0+2016.0*l3)*l3)*l3)*l3)*l3);
598 :
599 0 : default:
600 0 : libmesh_error_msg("Invalid i = " << i);
601 : }
602 : } // case TRI6/TRI7
603 :
604 : // Szabo-Babuska shape functions on the quadrilateral.
605 0 : case QUAD8:
606 : case QUAD9:
607 : {
608 : // Compute quad shape functions as a tensor-product
609 0 : const Real xi = p(0);
610 0 : const Real eta = p(1);
611 :
612 0 : libmesh_assert_less (i, 64);
613 :
614 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
615 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
616 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
617 :
618 0 : const Real f = quad_flip(elem, totalorder, i);
619 :
620 0 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
621 0 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
622 :
623 : } // case QUAD8/QUAD9
624 :
625 0 : default:
626 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
627 :
628 : } // switch type
629 :
630 : } // case SEVENTH
631 :
632 :
633 : // by default throw an error
634 0 : default:
635 0 : libmesh_error_msg("ERROR: Unsupported polynomial order!");
636 : } // switch order
637 : }
638 :
639 :
640 :
641 :
642 :
643 : template <>
644 0 : Real FE<2,SZABAB>::shape(const ElemType,
645 : const Order,
646 : const unsigned int,
647 : const Point &)
648 : {
649 0 : libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
650 : return 0.;
651 : }
652 :
653 :
654 :
655 : template <>
656 0 : Real FE<2,SZABAB>::shape(const FEType fet,
657 : const Elem * elem,
658 : const unsigned int i,
659 : const Point & p,
660 : const bool add_p_level)
661 : {
662 0 : return FE<2,SZABAB>::shape(elem, fet.order, i, p, add_p_level);
663 : }
664 :
665 :
666 :
667 :
668 :
669 : template <>
670 868404 : Real FE<2,SZABAB>::shape_deriv(const Elem * elem,
671 : const Order order,
672 : const unsigned int i,
673 : const unsigned int j,
674 : const Point & p,
675 : const bool add_p_level)
676 : {
677 75678 : libmesh_assert(elem);
678 :
679 868404 : const ElemType type = elem->type();
680 :
681 944082 : const Order totalorder = order + add_p_level*elem->p_level();
682 :
683 868404 : switch (totalorder)
684 : {
685 :
686 : // 1st & 2nd-order Szabo-Babuska.
687 179388 : case FIRST:
688 : case SECOND:
689 : {
690 163788 : switch (type)
691 : {
692 :
693 : // Szabo-Babuska shape functions on the triangle.
694 95652 : case TRI3:
695 : case TRI6:
696 : case TRI7:
697 : {
698 95652 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
699 : }
700 :
701 :
702 : // Szabo-Babuska shape functions on the quadrilateral.
703 83736 : case QUAD4:
704 : case QUAD8:
705 : case QUAD9:
706 : {
707 : // Compute quad shape functions as a tensor-product
708 83736 : const Real xi = p(0);
709 83736 : const Real eta = p(1);
710 :
711 6978 : libmesh_assert_less (i, 9);
712 :
713 : // 0 1 2 3 4 5 6 7 8
714 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
715 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
716 :
717 83736 : switch (j)
718 : {
719 : // d()/dxi
720 3489 : case 0:
721 41868 : return (FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
722 41868 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
723 :
724 : // d()/deta
725 3489 : case 1:
726 41868 : return (FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
727 41868 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
728 :
729 0 : default:
730 0 : libmesh_error_msg("Invalid j = " << j);
731 : }
732 : }
733 :
734 0 : default:
735 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
736 : }
737 : }
738 :
739 :
740 :
741 : // 3rd-order Szabo-Babuska.
742 255216 : case THIRD:
743 : {
744 255216 : switch (type)
745 : {
746 : // Szabo-Babuska shape functions on the triangle.
747 142320 : case TRI6:
748 : case TRI7:
749 : {
750 142320 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
751 : }
752 :
753 :
754 : // Szabo-Babuska shape functions on the quadrilateral.
755 112896 : case QUAD8:
756 : case QUAD9:
757 : {
758 : // Compute quad shape functions as a tensor-product
759 112896 : const Real xi = p(0);
760 112896 : const Real eta = p(1);
761 :
762 9408 : libmesh_assert_less (i, 16);
763 :
764 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
765 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
766 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
767 :
768 112896 : const Real f = quad_flip(elem, totalorder, i);
769 :
770 112896 : switch (j)
771 : {
772 : // d()/dxi
773 4704 : case 0:
774 56448 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
775 56448 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
776 :
777 : // d()/deta
778 4704 : case 1:
779 56448 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
780 56448 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
781 :
782 0 : default:
783 0 : libmesh_error_msg("Invalid j = " << j);
784 : }
785 : }
786 :
787 0 : default:
788 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
789 : }
790 : }
791 :
792 :
793 :
794 :
795 : // 4th-order Szabo-Babuska.
796 433800 : case FOURTH:
797 : {
798 433800 : switch (type)
799 : {
800 :
801 : // Szabo-Babuska shape functions on the triangle.
802 241200 : case TRI6:
803 : case TRI7:
804 : {
805 241200 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
806 : }
807 :
808 :
809 : // Szabo-Babuska shape functions on the quadrilateral.
810 192600 : case QUAD8:
811 : case QUAD9:
812 : {
813 : // Compute quad shape functions as a tensor-product
814 192600 : const Real xi = p(0);
815 192600 : const Real eta = p(1);
816 :
817 16050 : libmesh_assert_less (i, 25);
818 :
819 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
820 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
821 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
822 :
823 192600 : const Real f = quad_flip(elem, totalorder, i);
824 :
825 192600 : switch (j)
826 : {
827 : // d()/dxi
828 8025 : case 0:
829 96300 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
830 96300 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
831 :
832 : // d()/deta
833 8025 : case 1:
834 96300 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
835 96300 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
836 :
837 0 : default:
838 0 : libmesh_error_msg("Invalid j = " << j);
839 : }
840 : }
841 :
842 0 : default:
843 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
844 : }
845 : }
846 :
847 :
848 :
849 :
850 : // 5th-order Szabo-Babuska.
851 0 : case FIFTH:
852 : {
853 : // Szabo-Babuska shape functions on the quadrilateral.
854 0 : switch (type)
855 : {
856 :
857 : // Szabo-Babuska shape functions on the triangle.
858 0 : case TRI6:
859 : case TRI7:
860 : {
861 0 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
862 : }
863 :
864 :
865 0 : case QUAD8:
866 : case QUAD9:
867 : {
868 : // Compute quad shape functions as a tensor-product
869 0 : const Real xi = p(0);
870 0 : const Real eta = p(1);
871 :
872 0 : libmesh_assert_less (i, 36);
873 :
874 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
875 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
876 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
877 :
878 0 : const Real f = quad_flip(elem, totalorder, i);
879 :
880 0 : switch (j)
881 : {
882 : // d()/dxi
883 0 : case 0:
884 0 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
885 0 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
886 :
887 : // d()/deta
888 0 : case 1:
889 0 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
890 0 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
891 :
892 0 : default:
893 0 : libmesh_error_msg("Invalid j = " << j);
894 : }
895 : }
896 :
897 0 : default:
898 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
899 : }
900 : }
901 :
902 :
903 : // 6th-order Szabo-Babuska.
904 0 : case SIXTH:
905 : {
906 : // Szabo-Babuska shape functions on the quadrilateral.
907 0 : switch (type)
908 : {
909 :
910 : // Szabo-Babuska shape functions on the triangle.
911 0 : case TRI6:
912 : case TRI7:
913 : {
914 0 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
915 : }
916 :
917 :
918 0 : case QUAD8:
919 : case QUAD9:
920 : {
921 : // Compute quad shape functions as a tensor-product
922 0 : const Real xi = p(0);
923 0 : const Real eta = p(1);
924 :
925 0 : libmesh_assert_less (i, 49);
926 :
927 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
928 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
929 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
930 :
931 0 : const Real f = quad_flip(elem, totalorder, i);
932 :
933 0 : switch (j)
934 : {
935 : // d()/dxi
936 0 : case 0:
937 0 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
938 0 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
939 :
940 : // d()/deta
941 0 : case 1:
942 0 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
943 0 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
944 :
945 0 : default:
946 0 : libmesh_error_msg("Invalid j = " << j);
947 : }
948 : }
949 :
950 0 : default:
951 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
952 : }
953 : }
954 :
955 :
956 : // 7th-order Szabo-Babuska.
957 0 : case SEVENTH:
958 : {
959 : // Szabo-Babuska shape functions on the quadrilateral.
960 0 : switch (type)
961 : {
962 :
963 : // Szabo-Babuska shape functions on the triangle.
964 0 : case TRI6:
965 : case TRI7:
966 : {
967 0 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
968 : }
969 :
970 :
971 0 : case QUAD8:
972 : case QUAD9:
973 : {
974 : // Compute quad shape functions as a tensor-product
975 0 : const Real xi = p(0);
976 0 : const Real eta = p(1);
977 :
978 0 : libmesh_assert_less (i, 64);
979 :
980 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
981 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
982 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
983 :
984 0 : const Real f = quad_flip(elem, totalorder, i);
985 :
986 0 : switch (j)
987 : {
988 : // d()/dxi
989 0 : case 0:
990 0 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
991 0 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
992 :
993 : // d()/deta
994 0 : case 1:
995 0 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
996 0 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
997 :
998 0 : default:
999 0 : libmesh_error_msg("Invalid j = " << j);
1000 : }
1001 : }
1002 :
1003 0 : default:
1004 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1005 : }
1006 : }
1007 :
1008 :
1009 :
1010 : // by default throw an error;call the orientation-independent shape functions
1011 0 : default:
1012 0 : libmesh_error_msg("ERROR: Unsupported polynomial order!");
1013 : }
1014 : }
1015 :
1016 :
1017 :
1018 :
1019 :
1020 : template <>
1021 0 : Real FE<2,SZABAB>::shape_deriv(const ElemType,
1022 : const Order,
1023 : const unsigned int,
1024 : const unsigned int,
1025 : const Point &)
1026 : {
1027 0 : libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
1028 : return 0.;
1029 : }
1030 :
1031 :
1032 : template <>
1033 0 : Real FE<2,SZABAB>::shape_deriv(const FEType fet,
1034 : const Elem * elem,
1035 : const unsigned int i,
1036 : const unsigned int j,
1037 : const Point & p,
1038 : const bool add_p_level)
1039 : {
1040 0 : return FE<2,SZABAB>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1041 : }
1042 :
1043 :
1044 :
1045 :
1046 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1047 :
1048 : template <>
1049 0 : Real FE<2,SZABAB>::shape_second_deriv(const ElemType,
1050 : const Order,
1051 : const unsigned int,
1052 : const unsigned int,
1053 : const Point &)
1054 : {
1055 : static bool warning_given = false;
1056 :
1057 0 : if (!warning_given)
1058 0 : libMesh::err << "Second derivatives for Szabab elements "
1059 0 : << " are not yet implemented!"
1060 0 : << std::endl;
1061 :
1062 0 : warning_given = true;
1063 0 : return 0.;
1064 : }
1065 :
1066 :
1067 :
1068 : template <>
1069 0 : Real FE<2,SZABAB>::shape_second_deriv(const Elem *,
1070 : const Order,
1071 : const unsigned int,
1072 : const unsigned int,
1073 : const Point &,
1074 : const bool)
1075 : {
1076 : static bool warning_given = false;
1077 :
1078 0 : if (!warning_given)
1079 0 : libMesh::err << "Second derivatives for Szabab elements "
1080 0 : << " are not yet implemented!"
1081 0 : << std::endl;
1082 :
1083 0 : warning_given = true;
1084 0 : return 0.;
1085 : }
1086 :
1087 :
1088 : template <>
1089 0 : Real FE<2,SZABAB>::shape_second_deriv(const FEType,
1090 : const Elem *,
1091 : const unsigned int,
1092 : const unsigned int,
1093 : const Point &,
1094 : const bool)
1095 : {
1096 : static bool warning_given = false;
1097 :
1098 0 : if (!warning_given)
1099 0 : libMesh::err << "Second derivatives for Szabab elements "
1100 0 : << " are not yet implemented!"
1101 0 : << std::endl;
1102 :
1103 0 : warning_given = true;
1104 0 : return 0.;
1105 : }
1106 :
1107 : #endif
1108 :
1109 : } // namespace libMesh
1110 :
1111 : #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
