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