Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 818652 : Real quad_flip(const Elem * elem,
49 : const unsigned int totalorder,
50 : const unsigned int i)
51 : {
52 68221 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
53 :
54 : // vertex and interior shape functions don't flip
55 818652 : if (i < 4 || i >= 4*totalorder)
56 34993 : return 1;
57 :
58 398736 : const int edge = (i-4) / (totalorder-1);
59 33228 : libmesh_assert_less(edge, 4);
60 33228 : libmesh_assert_greater_equal(edge, 0);
61 :
62 398736 : 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 411820 : if (edge_i%2 &&
67 157008 : (elem->positive_edge_orientation(edge) ==
68 223464 : (edge < 2)))
69 13084 : return -1;
70 :
71 26686 : return 1;
72 : }
73 :
74 : } // anonymous namespace
75 :
76 :
77 : namespace libMesh
78 : {
79 :
80 :
81 1525144 : LIBMESH_DEFAULT_VECTORIZED_FE(2,SZABAB)
82 :
83 :
84 : template <>
85 2541098 : 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 222539 : libmesh_assert(elem);
92 :
93 2541098 : const ElemType type = elem->type();
94 :
95 2986176 : 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 2541098 : switch (totalorder)
102 : {
103 : // 1st & 2nd-order Szabo-Babuska.
104 454278 : case FIRST:
105 : case SECOND:
106 : {
107 414424 : switch (type)
108 : {
109 :
110 : // Szabo-Babuska shape functions on the triangle.
111 371778 : case TRI3:
112 : case TRISHELL3:
113 : case TRI6:
114 : case TRI7:
115 : {
116 371778 : const Real l1 = 1-p(0)-p(1);
117 32979 : const Real l2 = p(0);
118 32979 : const Real l3 = p(1);
119 :
120 32979 : libmesh_assert_less (i, 6);
121 :
122 371778 : switch (i)
123 : {
124 6363 : case 0: return l1;
125 71562 : case 1: return l2;
126 71562 : case 2: return l3;
127 :
128 52364 : case 3: return l1*l2*(-4.*sqrt6);
129 52364 : case 4: return l2*l3*(-4.*sqrt6);
130 52364 : case 5: return l3*l1*(-4.*sqrt6);
131 :
132 0 : default:
133 0 : libmesh_error_msg("Invalid i = " << i);
134 : }
135 : }
136 :
137 :
138 : // Szabo-Babuska shape functions on the quadrilateral.
139 82500 : case QUAD4:
140 : case QUADSHELL4:
141 : case QUAD8:
142 : case QUADSHELL8:
143 : case QUAD9:
144 : case QUADSHELL9:
145 : {
146 : // Compute quad shape functions as a tensor-product
147 82500 : const Real xi = p(0);
148 82500 : const Real eta = p(1);
149 :
150 6875 : libmesh_assert_less (i, 9);
151 :
152 : // 0 1 2 3 4 5 6 7 8
153 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
154 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
155 :
156 82500 : return (FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
157 82500 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
158 :
159 : }
160 :
161 0 : default:
162 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
163 : }
164 : }
165 :
166 :
167 : // 3rd-order Szabo-Babuska.
168 739640 : case THIRD:
169 : {
170 674740 : switch (type)
171 : {
172 :
173 : // Szabo-Babuska shape functions on the triangle.
174 587960 : case TRI6:
175 : case TRI7:
176 : {
177 587960 : Real l1 = 1-p(0)-p(1);
178 52260 : Real l2 = p(0);
179 52260 : Real l3 = p(1);
180 :
181 52260 : Real f=1;
182 :
183 52260 : libmesh_assert_less (i, 10);
184 :
185 :
186 587960 : if (i==4 && elem->positive_edge_orientation(0))f=-1;
187 587960 : if (i==6 && elem->positive_edge_orientation(1))f=-1;
188 587960 : if (i==8 && elem->positive_edge_orientation(2))f=-1;
189 :
190 :
191 587960 : switch (i)
192 : {
193 : //nodal modes
194 5226 : case 0: return l1;
195 58796 : case 1: return l2;
196 58796 : case 2: return l3;
197 :
198 : //side modes
199 58796 : case 3: return l1*l2*(-4.*sqrt6);
200 58796 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
201 :
202 58796 : case 5: return l2*l3*(-4.*sqrt6);
203 58796 : case 6: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
204 :
205 58796 : case 7: return l3*l1*(-4.*sqrt6);
206 58796 : case 8: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
207 :
208 : //internal modes
209 58796 : case 9: return l1*l2*l3;
210 :
211 0 : default:
212 0 : libmesh_error_msg("Invalid i = " << i);
213 : }
214 : }
215 :
216 :
217 : // Szabo-Babuska shape functions on the quadrilateral.
218 151680 : case QUAD8:
219 : case QUADSHELL8:
220 : case QUAD9:
221 : case QUADSHELL9:
222 : {
223 : // Compute quad shape functions as a tensor-product
224 151680 : const Real xi = p(0);
225 151680 : const Real eta = p(1);
226 :
227 12640 : libmesh_assert_less (i, 16);
228 :
229 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
230 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
231 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
232 :
233 151680 : const Real f = quad_flip(elem, totalorder, i);
234 :
235 151680 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
236 151680 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
237 : }
238 :
239 0 : default:
240 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
241 : }
242 : }
243 :
244 :
245 :
246 :
247 : // 4th-order Szabo-Babuska.
248 1347180 : case FOURTH:
249 : {
250 1229395 : switch (type)
251 : {
252 : // Szabo-Babuska shape functions on the triangle.
253 1048680 : case TRI6:
254 : case TRI7:
255 : {
256 1048680 : Real l1 = 1-p(0)-p(1);
257 92910 : Real l2 = p(0);
258 92910 : Real l3 = p(1);
259 :
260 92910 : Real f=1;
261 :
262 92910 : libmesh_assert_less (i, 15);
263 :
264 :
265 1048680 : if (i== 4 && elem->positive_edge_orientation(0))f=-1;
266 1048680 : if (i== 7 && elem->positive_edge_orientation(1))f=-1;
267 1048680 : if (i==10 && elem->positive_edge_orientation(2))f=-1;
268 :
269 :
270 1048680 : switch (i)
271 : {
272 : //nodal modes
273 6194 : case 0: return l1;
274 69912 : case 1: return l2;
275 69912 : case 2: return l3;
276 :
277 : //side modes
278 69912 : case 3: return l1*l2*(-4.*sqrt6);
279 69912 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
280 76106 : case 5: return l1*l2*(-sqrt14)*(5.*pow<2>(l2-l1)-1);
281 :
282 69912 : case 6: return l2*l3*(-4.*sqrt6);
283 69912 : case 7: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
284 76106 : case 8: return l2*l3*(-sqrt14)*(5.*pow<2>(l3-l2)-1);
285 :
286 69912 : case 9: return l3*l1*(-4.*sqrt6);
287 69912 : case 10: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
288 76106 : case 11: return l3*l1*(-sqrt14)*(5.*pow<2>(l1-l3)-1);
289 :
290 : //internal modes
291 69912 : case 12: return l1*l2*l3;
292 :
293 69912 : case 13: return l1*l2*l3*(l2-l1);
294 69912 : case 14: return l1*l2*l3*(2*l3-1);
295 :
296 0 : default:
297 0 : libmesh_error_msg("Invalid i = " << i);
298 : }
299 : }
300 :
301 :
302 : // Szabo-Babuska shape functions on the quadrilateral.
303 298500 : case QUAD8:
304 : case QUADSHELL8:
305 : case QUAD9:
306 : case QUADSHELL9:
307 : {
308 : // Compute quad shape functions as a tensor-product
309 298500 : const Real xi = p(0);
310 298500 : const Real eta = p(1);
311 :
312 24875 : libmesh_assert_less (i, 25);
313 :
314 : // 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
315 : 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};
316 : 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};
317 :
318 298500 : const Real f = quad_flip(elem, totalorder, i);
319 :
320 298500 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
321 298500 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
322 : }
323 :
324 0 : default:
325 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
326 : }
327 : }
328 :
329 :
330 :
331 :
332 : // 5th-order Szabo-Babuska.
333 0 : case FIFTH:
334 : {
335 0 : switch (type)
336 : {
337 : // Szabo-Babuska shape functions on the triangle.
338 0 : case TRI6:
339 : case TRI7:
340 : {
341 0 : Real l1 = 1-p(0)-p(1);
342 0 : Real l2 = p(0);
343 0 : Real l3 = p(1);
344 :
345 0 : const Real x=l2-l1;
346 0 : const Real y=2.*l3-1;
347 :
348 0 : Real f=1;
349 :
350 0 : libmesh_assert_less (i, 21);
351 :
352 :
353 0 : if ((i== 4||i== 6) && elem->positive_edge_orientation(0))f=-1;
354 0 : if ((i== 8||i==10) && elem->positive_edge_orientation(1))f=-1;
355 0 : if ((i==12||i==14) && elem->positive_edge_orientation(2))f=-1;
356 :
357 :
358 0 : switch (i)
359 : {
360 : //nodal modes
361 0 : case 0: return l1;
362 0 : case 1: return l2;
363 0 : case 2: return l3;
364 :
365 : //side modes
366 0 : case 3: return l1*l2*(-4.*sqrt6);
367 0 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
368 0 : case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
369 0 : case 6: return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
370 :
371 0 : case 7: return l2*l3*(-4.*sqrt6);
372 0 : case 8: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
373 0 : case 9: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
374 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);
375 :
376 0 : case 11: return l3*l1*(-4.*sqrt6);
377 0 : case 12: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
378 0 : case 13: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
379 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);
380 :
381 : //internal modes
382 0 : case 15: return l1*l2*l3;
383 :
384 0 : case 16: return l1*l2*l3*x;
385 0 : case 17: return l1*l2*l3*y;
386 :
387 0 : case 18: return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
388 0 : case 19: return l1*l2*l3*x*y;
389 0 : case 20: return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
390 :
391 0 : default:
392 0 : libmesh_error_msg("Invalid i = " << i);
393 : }
394 : } // case TRI6/TRI7
395 :
396 : // Szabo-Babuska shape functions on the quadrilateral.
397 0 : case QUAD8:
398 : case QUADSHELL8:
399 : case QUAD9:
400 : case QUADSHELL9:
401 : {
402 : // Compute quad shape functions as a tensor-product
403 0 : const Real xi = p(0);
404 0 : const Real eta = p(1);
405 :
406 0 : libmesh_assert_less (i, 36);
407 :
408 : // 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
409 : 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};
410 : 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};
411 :
412 0 : const Real f = quad_flip(elem, totalorder, i);
413 :
414 0 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
415 0 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
416 :
417 : } // case QUAD8/QUAD9
418 :
419 0 : default:
420 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
421 :
422 : } // switch type
423 :
424 : } // case FIFTH
425 :
426 : // 6th-order Szabo-Babuska.
427 0 : case SIXTH:
428 : {
429 0 : switch (type)
430 : {
431 : // Szabo-Babuska shape functions on the triangle.
432 0 : case TRI6:
433 : case TRI7:
434 : {
435 0 : Real l1 = 1-p(0)-p(1);
436 0 : Real l2 = p(0);
437 0 : Real l3 = p(1);
438 :
439 0 : const Real x=l2-l1;
440 0 : const Real y=2.*l3-1;
441 :
442 0 : Real f=1;
443 :
444 0 : libmesh_assert_less (i, 28);
445 :
446 :
447 0 : if ((i== 4||i== 6) && elem->positive_edge_orientation(0))f=-1;
448 0 : if ((i== 9||i==11) && elem->positive_edge_orientation(1))f=-1;
449 0 : if ((i==14||i==16) && elem->positive_edge_orientation(2))f=-1;
450 :
451 :
452 0 : switch (i)
453 : {
454 : //nodal modes
455 0 : case 0: return l1;
456 0 : case 1: return l2;
457 0 : case 2: return l3;
458 :
459 : //side modes
460 0 : case 3: return l1*l2*(-4.*sqrt6);
461 0 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
462 0 : case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
463 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);
464 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);
465 :
466 0 : case 8: return l2*l3*(-4.*sqrt6);
467 0 : case 9: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
468 0 : case 10: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
469 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);
470 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);
471 :
472 0 : case 13: return l3*l1*(-4.*sqrt6);
473 0 : case 14: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
474 0 : case 15: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
475 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);
476 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);
477 :
478 :
479 :
480 : //internal modes
481 0 : case 18: return l1*l2*l3;
482 :
483 0 : case 19: return l1*l2*l3*x;
484 0 : case 20: return l1*l2*l3*y;
485 :
486 0 : case 21: return 0.5*l1*l2*l3*(3.0*l1*l1-1.0+(-6.0*l1+3.0*l2)*l2);
487 0 : case 22: return l1*l2*l3*(l2-l1)*(2.0*l3-1.0);
488 0 : case 23: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3);
489 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);
490 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);
491 0 : case 26: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
492 0 : case 27: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
493 :
494 :
495 0 : default:
496 0 : libmesh_error_msg("Invalid i = " << i);
497 : }
498 : } // case TRI6/TRI7
499 :
500 : // Szabo-Babuska shape functions on the quadrilateral.
501 0 : case QUAD8:
502 : case QUADSHELL8:
503 : case QUAD9:
504 : case QUADSHELL9:
505 : {
506 : // Compute quad shape functions as a tensor-product
507 0 : const Real xi = p(0);
508 0 : const Real eta = p(1);
509 :
510 0 : libmesh_assert_less (i, 49);
511 :
512 : // 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
513 : 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};
514 : 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};
515 :
516 0 : const Real f = quad_flip(elem, totalorder, i);
517 :
518 0 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
519 0 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
520 :
521 : } // case QUAD8/QUAD9
522 :
523 0 : default:
524 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
525 :
526 : } // switch type
527 :
528 : } // case SIXTH
529 :
530 :
531 : // 7th-order Szabo-Babuska.
532 0 : case SEVENTH:
533 : {
534 0 : switch (type)
535 : {
536 : // Szabo-Babuska shape functions on the triangle.
537 0 : case TRI6:
538 : case TRI7:
539 : {
540 :
541 0 : Real l1 = 1-p(0)-p(1);
542 0 : Real l2 = p(0);
543 0 : Real l3 = p(1);
544 :
545 0 : const Real x=l2-l1;
546 0 : const Real y=2.*l3-1.;
547 :
548 0 : Real f=1;
549 :
550 0 : libmesh_assert_less (i, 36);
551 :
552 :
553 0 : if ((i>= 4&&i<= 8) && elem->positive_edge_orientation(0))f=-1;
554 0 : if ((i>=10&&i<=14) && elem->positive_edge_orientation(1))f=-1;
555 0 : if ((i>=16&&i<=20) && elem->positive_edge_orientation(2))f=-1;
556 :
557 :
558 0 : switch (i)
559 : {
560 : //nodal modes
561 0 : case 0: return l1;
562 0 : case 1: return l2;
563 0 : case 2: return l3;
564 :
565 : //side modes
566 0 : case 3: return l1*l2*(-4.*sqrt6);
567 0 : case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
568 :
569 0 : case 5: return -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
570 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);
571 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);
572 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);
573 :
574 0 : case 9: return l2*l3*(-4.*sqrt6);
575 0 : case 10: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
576 :
577 0 : case 11: return -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
578 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);
579 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);
580 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);
581 :
582 0 : case 15: return l3*l1*(-4.*sqrt6);
583 0 : case 16: return f*l3*l1*(-4.*sqrt10)*(l1-l3);
584 :
585 0 : case 17: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
586 0 : case 18: return -f*sqrt2*l3*l1*((9.-21.*l3*l3)*l3+(-9.+63.*l3*l3+(-63.*l3+21.*l1)*l1)*l1);
587 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);
588 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);
589 :
590 :
591 : //internal modes
592 0 : case 21: return l1*l2*l3;
593 :
594 0 : case 22: return l1*l2*l3*x;
595 0 : case 23: return l1*l2*l3*y;
596 :
597 0 : case 24: return l1*l2*l3*0.5*(3.*pow<2>(x)-1.);
598 0 : case 25: return l1*l2*l3*x*y;
599 0 : case 26: return l1*l2*l3*0.5*(3.*pow<2>(y)-1.);
600 :
601 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);
602 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);
603 0 : case 29: return 0.5*l1*l2*l3*(2.0+(-12.0+12.0*l3)*l3)*(l2-l1);
604 0 : case 30: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3);
605 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);
606 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);
607 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);
608 0 : case 34: return 0.5*l1*l2*l3*(-2.0+(24.0+(-60.0+40.0*l3)*l3)*l3)*(l2-l1);
609 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);
610 :
611 0 : default:
612 0 : libmesh_error_msg("Invalid i = " << i);
613 : }
614 : } // case TRI6/TRI7
615 :
616 : // Szabo-Babuska shape functions on the quadrilateral.
617 0 : case QUAD8:
618 : case QUADSHELL8:
619 : case QUAD9:
620 : case QUADSHELL9:
621 : {
622 : // Compute quad shape functions as a tensor-product
623 0 : const Real xi = p(0);
624 0 : const Real eta = p(1);
625 :
626 0 : libmesh_assert_less (i, 64);
627 :
628 : // 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
629 : 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};
630 : 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};
631 :
632 0 : const Real f = quad_flip(elem, totalorder, i);
633 :
634 0 : return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
635 0 : FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
636 :
637 : } // case QUAD8/QUAD9
638 :
639 0 : default:
640 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
641 :
642 : } // switch type
643 :
644 : } // case SEVENTH
645 :
646 :
647 : // by default throw an error
648 0 : default:
649 0 : libmesh_error_msg("ERROR: Unsupported polynomial order!");
650 : } // switch order
651 : }
652 :
653 :
654 :
655 :
656 :
657 : template <>
658 0 : Real FE<2,SZABAB>::shape(const ElemType,
659 : const Order,
660 : const unsigned int,
661 : const Point &)
662 : {
663 0 : libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
664 : return 0.;
665 : }
666 :
667 :
668 :
669 : template <>
670 0 : Real FE<2,SZABAB>::shape(const FEType fet,
671 : const Elem * elem,
672 : const unsigned int i,
673 : const Point & p,
674 : const bool add_p_level)
675 : {
676 0 : return FE<2,SZABAB>::shape(elem, fet.order, i, p, add_p_level);
677 : }
678 :
679 :
680 :
681 :
682 :
683 : template <>
684 1054308 : Real FE<2,SZABAB>::shape_deriv(const Elem * elem,
685 : const Order order,
686 : const unsigned int i,
687 : const unsigned int j,
688 : const Point & p,
689 : const bool add_p_level)
690 : {
691 91950 : libmesh_assert(elem);
692 :
693 1054308 : const ElemType type = elem->type();
694 :
695 1146258 : const Order totalorder = order + add_p_level*elem->p_level();
696 :
697 1054308 : switch (totalorder)
698 : {
699 :
700 : // 1st & 2nd-order Szabo-Babuska.
701 223116 : case FIRST:
702 : case SECOND:
703 : {
704 203692 : switch (type)
705 : {
706 :
707 : // Szabo-Babuska shape functions on the triangle.
708 119412 : case TRI3:
709 : case TRISHELL3:
710 : case TRI6:
711 : case TRI7:
712 : {
713 119412 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
714 : }
715 :
716 :
717 : // Szabo-Babuska shape functions on the quadrilateral.
718 103704 : case QUAD4:
719 : case QUADSHELL4:
720 : case QUAD8:
721 : case QUADSHELL8:
722 : case QUAD9:
723 : case QUADSHELL9:
724 : {
725 : // Compute quad shape functions as a tensor-product
726 103704 : const Real xi = p(0);
727 103704 : const Real eta = p(1);
728 :
729 8642 : libmesh_assert_less (i, 9);
730 :
731 : // 0 1 2 3 4 5 6 7 8
732 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
733 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
734 :
735 103704 : switch (j)
736 : {
737 : // d()/dxi
738 4321 : case 0:
739 51852 : return (FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
740 51852 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
741 :
742 : // d()/deta
743 4321 : case 1:
744 51852 : return (FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
745 51852 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
746 :
747 0 : default:
748 0 : libmesh_error_msg("Invalid j = " << j);
749 : }
750 : }
751 :
752 0 : default:
753 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
754 : }
755 : }
756 :
757 :
758 :
759 : // 3rd-order Szabo-Babuska.
760 311472 : case THIRD:
761 : {
762 284296 : switch (type)
763 : {
764 : // Szabo-Babuska shape functions on the triangle.
765 174000 : case TRI6:
766 : case TRI7:
767 : {
768 174000 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
769 : }
770 :
771 :
772 : // Szabo-Babuska shape functions on the quadrilateral.
773 137472 : case QUAD8:
774 : case QUADSHELL8:
775 : case QUAD9:
776 : case QUADSHELL9:
777 : {
778 : // Compute quad shape functions as a tensor-product
779 137472 : const Real xi = p(0);
780 137472 : const Real eta = p(1);
781 :
782 11456 : libmesh_assert_less (i, 16);
783 :
784 : // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
785 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
786 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
787 :
788 137472 : const Real f = quad_flip(elem, totalorder, i);
789 :
790 137472 : switch (j)
791 : {
792 : // d()/dxi
793 5728 : case 0:
794 68736 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
795 68736 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
796 :
797 : // d()/deta
798 5728 : case 1:
799 68736 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
800 68736 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
801 :
802 0 : default:
803 0 : libmesh_error_msg("Invalid j = " << j);
804 : }
805 : }
806 :
807 0 : default:
808 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
809 : }
810 : }
811 :
812 :
813 :
814 :
815 : // 4th-order Szabo-Babuska.
816 519720 : case FOURTH:
817 : {
818 474370 : switch (type)
819 : {
820 :
821 : // Szabo-Babuska shape functions on the triangle.
822 288720 : case TRI6:
823 : case TRI7:
824 : {
825 288720 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
826 : }
827 :
828 :
829 : // Szabo-Babuska shape functions on the quadrilateral.
830 231000 : case QUAD8:
831 : case QUADSHELL8:
832 : case QUAD9:
833 : case QUADSHELL9:
834 : {
835 : // Compute quad shape functions as a tensor-product
836 231000 : const Real xi = p(0);
837 231000 : const Real eta = p(1);
838 :
839 19250 : libmesh_assert_less (i, 25);
840 :
841 : // 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
842 : 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};
843 : 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};
844 :
845 231000 : const Real f = quad_flip(elem, totalorder, i);
846 :
847 231000 : switch (j)
848 : {
849 : // d()/dxi
850 9625 : case 0:
851 115500 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
852 115500 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
853 :
854 : // d()/deta
855 9625 : case 1:
856 115500 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
857 115500 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
858 :
859 0 : default:
860 0 : libmesh_error_msg("Invalid j = " << j);
861 : }
862 : }
863 :
864 0 : default:
865 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
866 : }
867 : }
868 :
869 :
870 :
871 :
872 : // 5th-order Szabo-Babuska.
873 0 : case FIFTH:
874 : {
875 : // Szabo-Babuska shape functions on the quadrilateral.
876 0 : switch (type)
877 : {
878 :
879 : // Szabo-Babuska shape functions on the triangle.
880 0 : case TRI6:
881 : case TRI7:
882 : {
883 0 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
884 : }
885 :
886 :
887 0 : case QUAD8:
888 : case QUADSHELL8:
889 : case QUAD9:
890 : case QUADSHELL9:
891 : {
892 : // Compute quad shape functions as a tensor-product
893 0 : const Real xi = p(0);
894 0 : const Real eta = p(1);
895 :
896 0 : libmesh_assert_less (i, 36);
897 :
898 : // 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
899 : 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};
900 : 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};
901 :
902 0 : const Real f = quad_flip(elem, totalorder, i);
903 :
904 0 : switch (j)
905 : {
906 : // d()/dxi
907 0 : case 0:
908 0 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
909 0 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
910 :
911 : // d()/deta
912 0 : case 1:
913 0 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
914 0 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
915 :
916 0 : default:
917 0 : libmesh_error_msg("Invalid j = " << j);
918 : }
919 : }
920 :
921 0 : default:
922 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
923 : }
924 : }
925 :
926 :
927 : // 6th-order Szabo-Babuska.
928 0 : case SIXTH:
929 : {
930 : // Szabo-Babuska shape functions on the quadrilateral.
931 0 : switch (type)
932 : {
933 :
934 : // Szabo-Babuska shape functions on the triangle.
935 0 : case TRI6:
936 : case TRI7:
937 : {
938 0 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
939 : }
940 :
941 :
942 0 : case QUAD8:
943 : case QUADSHELL8:
944 : case QUAD9:
945 : case QUADSHELL9:
946 : {
947 : // Compute quad shape functions as a tensor-product
948 0 : const Real xi = p(0);
949 0 : const Real eta = p(1);
950 :
951 0 : libmesh_assert_less (i, 49);
952 :
953 : // 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
954 : 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};
955 : 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};
956 :
957 0 : const Real f = quad_flip(elem, totalorder, i);
958 :
959 0 : switch (j)
960 : {
961 : // d()/dxi
962 0 : case 0:
963 0 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
964 0 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
965 :
966 : // d()/deta
967 0 : case 1:
968 0 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
969 0 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
970 :
971 0 : default:
972 0 : libmesh_error_msg("Invalid j = " << j);
973 : }
974 : }
975 :
976 0 : default:
977 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
978 : }
979 : }
980 :
981 :
982 : // 7th-order Szabo-Babuska.
983 0 : case SEVENTH:
984 : {
985 : // Szabo-Babuska shape functions on the quadrilateral.
986 0 : switch (type)
987 : {
988 :
989 : // Szabo-Babuska shape functions on the triangle.
990 0 : case TRI6:
991 : case TRI7:
992 : {
993 0 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,SZABAB>::shape);
994 : }
995 :
996 :
997 0 : case QUAD8:
998 : case QUADSHELL8:
999 : case QUAD9:
1000 : case QUADSHELL9:
1001 : {
1002 : // Compute quad shape functions as a tensor-product
1003 0 : const Real xi = p(0);
1004 0 : const Real eta = p(1);
1005 :
1006 0 : libmesh_assert_less (i, 64);
1007 :
1008 : // 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
1009 : 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};
1010 : 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};
1011 :
1012 0 : const Real f = quad_flip(elem, totalorder, i);
1013 :
1014 0 : switch (j)
1015 : {
1016 : // d()/dxi
1017 0 : case 0:
1018 0 : return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1019 0 : FE<1,SZABAB>::shape (EDGE3, totalorder, i1[i], eta));
1020 :
1021 : // d()/deta
1022 0 : case 1:
1023 0 : return f*(FE<1,SZABAB>::shape (EDGE3, totalorder, i0[i], xi)*
1024 0 : FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1025 :
1026 0 : default:
1027 0 : libmesh_error_msg("Invalid j = " << j);
1028 : }
1029 : }
1030 :
1031 0 : default:
1032 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1033 : }
1034 : }
1035 :
1036 :
1037 :
1038 : // by default throw an error;call the orientation-independent shape functions
1039 0 : default:
1040 0 : libmesh_error_msg("ERROR: Unsupported polynomial order!");
1041 : }
1042 : }
1043 :
1044 :
1045 :
1046 :
1047 :
1048 : template <>
1049 0 : Real FE<2,SZABAB>::shape_deriv(const ElemType,
1050 : const Order,
1051 : const unsigned int,
1052 : const unsigned int,
1053 : const Point &)
1054 : {
1055 0 : libmesh_error_msg("Szabo-Babuska polynomials require the element type \nbecause edge orientation is needed.");
1056 : return 0.;
1057 : }
1058 :
1059 :
1060 : template <>
1061 0 : Real FE<2,SZABAB>::shape_deriv(const FEType fet,
1062 : const Elem * elem,
1063 : const unsigned int i,
1064 : const unsigned int j,
1065 : const Point & p,
1066 : const bool add_p_level)
1067 : {
1068 0 : return FE<2,SZABAB>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1069 : }
1070 :
1071 :
1072 :
1073 :
1074 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1075 :
1076 : template <>
1077 0 : Real FE<2,SZABAB>::shape_second_deriv(const ElemType,
1078 : const Order,
1079 : const unsigned int,
1080 : const unsigned int,
1081 : const Point &)
1082 : {
1083 : static bool warning_given = false;
1084 :
1085 0 : if (!warning_given)
1086 0 : libMesh::err << "Second derivatives for Szabab elements "
1087 0 : << " are not yet implemented!"
1088 0 : << std::endl;
1089 :
1090 0 : warning_given = true;
1091 0 : return 0.;
1092 : }
1093 :
1094 :
1095 :
1096 : template <>
1097 0 : Real FE<2,SZABAB>::shape_second_deriv(const Elem *,
1098 : const Order,
1099 : const unsigned int,
1100 : const unsigned int,
1101 : const Point &,
1102 : const bool)
1103 : {
1104 : static bool warning_given = false;
1105 :
1106 0 : if (!warning_given)
1107 0 : libMesh::err << "Second derivatives for Szabab elements "
1108 0 : << " are not yet implemented!"
1109 0 : << std::endl;
1110 :
1111 0 : warning_given = true;
1112 0 : return 0.;
1113 : }
1114 :
1115 :
1116 : template <>
1117 0 : Real FE<2,SZABAB>::shape_second_deriv(const FEType,
1118 : const Elem *,
1119 : const unsigned int,
1120 : const unsigned int,
1121 : const Point &,
1122 : const bool)
1123 : {
1124 : static bool warning_given = false;
1125 :
1126 0 : if (!warning_given)
1127 0 : libMesh::err << "Second derivatives for Szabab elements "
1128 0 : << " are not yet implemented!"
1129 0 : << std::endl;
1130 :
1131 0 : warning_given = true;
1132 0 : return 0.;
1133 : }
1134 :
1135 : #endif
1136 :
1137 : } // namespace libMesh
1138 :
1139 : #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
