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 :
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 34837611 : LIBMESH_DEFAULT_VECTORIZED_FE(3,MONOMIAL)
31 :
32 :
33 : template <>
34 33485291 : 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 33485291 : const Real xi = p(0);
42 33485291 : const Real eta = p(1);
43 33485291 : const Real zeta = p(2);
44 :
45 8557371 : 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 33485291 : switch (i)
51 : {
52 : // constant
53 954533 : case 0:
54 954533 : return 1.;
55 :
56 : // linears
57 2356696 : case 1:
58 2356696 : return xi;
59 :
60 2356696 : case 2:
61 2356696 : return eta;
62 :
63 2356696 : case 3:
64 2356696 : return zeta;
65 :
66 : // quadratics
67 735072 : case 4:
68 735072 : return xi*xi;
69 :
70 735072 : case 5:
71 735072 : return xi*eta;
72 :
73 735072 : case 6:
74 735072 : return eta*eta;
75 :
76 735072 : case 7:
77 735072 : return xi*zeta;
78 :
79 735072 : case 8:
80 735072 : return zeta*eta;
81 :
82 735072 : case 9:
83 735072 : return zeta*zeta;
84 :
85 : // cubics
86 664361 : case 10:
87 664361 : return xi*xi*xi;
88 :
89 664361 : case 11:
90 664361 : return xi*xi*eta;
91 :
92 664361 : case 12:
93 664361 : return xi*eta*eta;
94 :
95 664361 : case 13:
96 664361 : return eta*eta*eta;
97 :
98 664361 : case 14:
99 664361 : return xi*xi*zeta;
100 :
101 664361 : case 15:
102 664361 : return xi*eta*zeta;
103 :
104 664361 : case 16:
105 664361 : return eta*eta*zeta;
106 :
107 664361 : case 17:
108 664361 : return xi*zeta*zeta;
109 :
110 664361 : case 18:
111 664361 : return eta*zeta*zeta;
112 :
113 664361 : case 19:
114 664361 : return zeta*zeta*zeta;
115 :
116 : // quartics
117 418866 : case 20:
118 418866 : return xi*xi*xi*xi;
119 :
120 418866 : case 21:
121 418866 : return xi*xi*xi*eta;
122 :
123 418866 : case 22:
124 418866 : return xi*xi*eta*eta;
125 :
126 418866 : case 23:
127 418866 : return xi*eta*eta*eta;
128 :
129 418866 : case 24:
130 418866 : return eta*eta*eta*eta;
131 :
132 418866 : case 25:
133 418866 : return xi*xi*xi*zeta;
134 :
135 418866 : case 26:
136 418866 : return xi*xi*eta*zeta;
137 :
138 418866 : case 27:
139 418866 : return xi*eta*eta*zeta;
140 :
141 418866 : case 28:
142 418866 : return eta*eta*eta*zeta;
143 :
144 418866 : case 29:
145 418866 : return xi*xi*zeta*zeta;
146 :
147 418866 : case 30:
148 418866 : return xi*eta*zeta*zeta;
149 :
150 418866 : case 31:
151 418866 : return eta*eta*zeta*zeta;
152 :
153 418866 : case 32:
154 418866 : return xi*zeta*zeta*zeta;
155 :
156 418866 : case 33:
157 418866 : return eta*zeta*zeta*zeta;
158 :
159 418866 : case 34:
160 418866 : return zeta*zeta*zeta*zeta;
161 :
162 1357335 : default:
163 1357335 : unsigned int o = 0;
164 31738266 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
165 5289711 : const int i2 = i - (o*(o+1)*(o+2)/6);
166 5289711 : int block=o, nz = 0;
167 14105896 : for (; block < i2; block += (o-nz+1)) { nz++; }
168 5289711 : const int nx = block - i2;
169 5289711 : const int ny = o - nx - nz;
170 1357335 : Real val = 1.;
171 14105896 : for (int index=0; index != nx; index++)
172 11530855 : val *= xi;
173 14105896 : for (int index=0; index != ny; index++)
174 11530855 : val *= eta;
175 14105896 : for (int index=0; index != nz; index++)
176 11530855 : 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 33485291 : 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 8557371 : libmesh_assert(elem);
197 :
198 : // call the orientation-independent shape functions
199 50600033 : 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 32326488 : 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 8565657 : libmesh_assert_less (j, 3);
228 :
229 8565657 : 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 32326488 : const Real xi = p(0);
235 32326488 : const Real eta = p(1);
236 32326488 : const Real zeta = p(2);
237 :
238 : // monomials. since they are hierarchic we only need one case block.
239 32326488 : switch (j)
240 : {
241 : // d()/dxi
242 10775496 : case 0:
243 : {
244 10775496 : switch (i)
245 : {
246 : // constant
247 178322 : case 0:
248 178322 : return 0.;
249 :
250 : // linear
251 356468 : case 1:
252 356468 : return 1.;
253 :
254 178234 : case 2:
255 178234 : return 0.;
256 :
257 178234 : case 3:
258 178234 : return 0.;
259 :
260 : // quadratic
261 279634 : case 4:
262 279634 : 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 226595 : case 10:
281 226595 : return 3.*xi*xi;
282 :
283 226595 : case 11:
284 226595 : return 2.*xi*eta;
285 :
286 226595 : case 12:
287 226595 : return eta*eta;
288 :
289 60630 : case 13:
290 60630 : return 0.;
291 :
292 226595 : case 14:
293 226595 : return 2.*xi*zeta;
294 :
295 226595 : case 15:
296 226595 : return eta*zeta;
297 :
298 60630 : case 16:
299 60630 : return 0.;
300 :
301 226595 : case 17:
302 226595 : 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 155695 : case 20:
312 155695 : return 4.*xi*xi*xi;
313 :
314 155695 : case 21:
315 155695 : return 3.*xi*xi*eta;
316 :
317 155695 : case 22:
318 155695 : return 2.*xi*eta*eta;
319 :
320 155695 : case 23:
321 155695 : return eta*eta*eta;
322 :
323 41682 : case 24:
324 41682 : return 0.;
325 :
326 155695 : case 25:
327 155695 : return 3.*xi*xi*zeta;
328 :
329 155695 : case 26:
330 155695 : return 2.*xi*eta*zeta;
331 :
332 155695 : case 27:
333 155695 : return eta*eta*zeta;
334 :
335 41682 : case 28:
336 41682 : return 0.;
337 :
338 155695 : case 29:
339 155695 : return 2.*xi*zeta*zeta;
340 :
341 155695 : case 30:
342 155695 : return eta*zeta*zeta;
343 :
344 41682 : case 31:
345 41682 : return 0.;
346 :
347 155695 : case 32:
348 155695 : 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 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
359 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
360 1726305 : int block=o, nz = 0;
361 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
362 1726305 : const int nx = block - i2;
363 1726305 : const int ny = o - nx - nz;
364 1726305 : Real val = nx;
365 3370405 : for (int index=1; index < nx; index++)
366 2570074 : val *= xi;
367 4603480 : for (int index=0; index != ny; index++)
368 3803149 : val *= eta;
369 4603480 : for (int index=0; index != nz; index++)
370 3803149 : val *= zeta;
371 462987 : return val;
372 : }
373 : }
374 :
375 :
376 : // d()/deta
377 10775496 : case 1:
378 : {
379 10775496 : switch (i)
380 : {
381 : // constant
382 178322 : case 0:
383 178322 : return 0.;
384 :
385 : // linear
386 178234 : case 1:
387 178234 : return 0.;
388 :
389 356468 : case 2:
390 356468 : return 1.;
391 :
392 178234 : case 3:
393 178234 : return 0.;
394 :
395 : // quadratic
396 74613 : case 4:
397 74613 : return 0.;
398 :
399 149226 : case 5:
400 149226 : return xi;
401 :
402 279634 : case 6:
403 279634 : 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 226595 : case 11:
419 226595 : return xi*xi;
420 :
421 226595 : case 12:
422 226595 : return 2.*xi*eta;
423 :
424 226595 : case 13:
425 226595 : return 3.*eta*eta;
426 :
427 60630 : case 14:
428 60630 : return 0.;
429 :
430 226595 : case 15:
431 226595 : return xi*zeta;
432 :
433 226595 : case 16:
434 226595 : return 2.*eta*zeta;
435 :
436 60630 : case 17:
437 60630 : return 0.;
438 :
439 226595 : case 18:
440 226595 : 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 155695 : case 21:
450 155695 : return xi*xi*xi;
451 :
452 155695 : case 22:
453 155695 : return 2.*xi*xi*eta;
454 :
455 155695 : case 23:
456 155695 : return 3.*xi*eta*eta;
457 :
458 155695 : case 24:
459 155695 : return 4.*eta*eta*eta;
460 :
461 41682 : case 25:
462 41682 : return 0.;
463 :
464 155695 : case 26:
465 155695 : return xi*xi*zeta;
466 :
467 155695 : case 27:
468 155695 : return 2.*xi*eta*zeta;
469 :
470 155695 : case 28:
471 155695 : return 3.*eta*eta*zeta;
472 :
473 41682 : case 29:
474 41682 : return 0.;
475 :
476 155695 : case 30:
477 155695 : return xi*zeta*zeta;
478 :
479 155695 : case 31:
480 155695 : return 2.*eta*zeta*zeta;
481 :
482 41682 : case 32:
483 41682 : return 0.;
484 :
485 155695 : case 33:
486 155695 : return zeta*zeta*zeta;
487 :
488 41682 : case 34:
489 41682 : return 0.;
490 :
491 462987 : default:
492 462987 : unsigned int o = 0;
493 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
494 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
495 1726305 : int block=o, nz = 0;
496 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
497 1726305 : const int nx = block - i2;
498 1726305 : const int ny = o - nx - nz;
499 1726305 : Real val = ny;
500 4603480 : for (int index=0; index != nx; index++)
501 3803149 : val *= xi;
502 3370405 : for (int index=1; index < ny; index++)
503 2570074 : val *= eta;
504 4603480 : for (int index=0; index != nz; index++)
505 3803149 : val *= zeta;
506 462987 : return val;
507 : }
508 : }
509 :
510 :
511 : // d()/dzeta
512 10775496 : case 2:
513 : {
514 10775496 : switch (i)
515 : {
516 : // constant
517 178322 : case 0:
518 178322 : return 0.;
519 :
520 : // linear
521 178234 : case 1:
522 178234 : return 0.;
523 :
524 178234 : case 2:
525 178234 : return 0.;
526 :
527 356468 : case 3:
528 356468 : 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 279634 : case 9:
547 279634 : 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 226595 : case 14:
563 226595 : return xi*xi;
564 :
565 226595 : case 15:
566 226595 : return xi*eta;
567 :
568 226595 : case 16:
569 226595 : return eta*eta;
570 :
571 226595 : case 17:
572 226595 : return 2.*xi*zeta;
573 :
574 226595 : case 18:
575 226595 : return 2.*eta*zeta;
576 :
577 226595 : case 19:
578 226595 : 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 155695 : case 25:
597 155695 : return xi*xi*xi;
598 :
599 155695 : case 26:
600 155695 : return xi*xi*eta;
601 :
602 155695 : case 27:
603 155695 : return xi*eta*eta;
604 :
605 155695 : case 28:
606 155695 : return eta*eta*eta;
607 :
608 155695 : case 29:
609 155695 : return 2.*xi*xi*zeta;
610 :
611 155695 : case 30:
612 155695 : return 2.*xi*eta*zeta;
613 :
614 155695 : case 31:
615 155695 : return 2.*eta*eta*zeta;
616 :
617 155695 : case 32:
618 155695 : return 3.*xi*zeta*zeta;
619 :
620 155695 : case 33:
621 155695 : return 3.*eta*zeta*zeta;
622 :
623 155695 : case 34:
624 155695 : return 4.*zeta*zeta*zeta;
625 :
626 462987 : default:
627 462987 : unsigned int o = 0;
628 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
629 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
630 1726305 : int block=o, nz = 0;
631 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
632 1726305 : const int nx = block - i2;
633 1726305 : const int ny = o - nx - nz;
634 1726305 : Real val = nz;
635 4603480 : for (int index=0; index != nx; index++)
636 3803149 : val *= xi;
637 4603480 : for (int index=0; index != ny; index++)
638 3803149 : val *= eta;
639 3370405 : for (int index=1; index < nz; index++)
640 2570074 : 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 32326488 : 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 8565657 : libmesh_assert(elem);
667 :
668 : // call the orientation-independent shape function derivatives
669 49457802 : 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 56059632 : 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 14982834 : libmesh_assert_less (j, 6);
699 :
700 14982834 : 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 56059632 : const Real xi = p(0);
705 56059632 : const Real eta = p(1);
706 56059632 : const Real zeta = p(2);
707 :
708 : // monomials. since they are hierarchic we only need one case block.
709 56059632 : switch (j)
710 : {
711 : // d^2()/dxi^2
712 9343272 : case 0:
713 : {
714 9343272 : switch (i)
715 : {
716 : // constant
717 354944 : case 0:
718 :
719 : // linear
720 : case 1:
721 : case 2:
722 : case 3:
723 354944 : 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 226595 : case 10:
738 226595 : return 6.*xi;
739 :
740 226595 : case 11:
741 226595 : return 2.*eta;
742 :
743 121260 : case 12:
744 : case 13:
745 121260 : return 0.;
746 :
747 226595 : case 14:
748 226595 : 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 155695 : case 20:
759 155695 : return 12.*xi*xi;
760 :
761 155695 : case 21:
762 155695 : return 6.*xi*eta;
763 :
764 155695 : case 22:
765 155695 : return 2.*eta*eta;
766 :
767 83364 : case 23:
768 : case 24:
769 83364 : return 0.;
770 :
771 155695 : case 25:
772 155695 : return 6.*xi*zeta;
773 :
774 155695 : case 26:
775 155695 : return 2.*eta*zeta;
776 :
777 83364 : case 27:
778 : case 28:
779 83364 : return 0.;
780 :
781 155695 : case 29:
782 155695 : 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 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
794 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
795 1726305 : int block=o, nz = 0;
796 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
797 1726305 : const int nx = block - i2;
798 1726305 : const int ny = o - nx - nz;
799 1726305 : Real val = nx * (nx - 1);
800 2548355 : for (int index=2; index < nx; index++)
801 1748024 : val *= xi;
802 4603480 : for (int index=0; index != ny; index++)
803 3803149 : val *= eta;
804 4603480 : for (int index=0; index != nz; index++)
805 3803149 : val *= zeta;
806 462987 : return val;
807 : }
808 : }
809 :
810 :
811 : // d^2()/dxideta
812 9343272 : case 1:
813 : {
814 9343272 : switch (i)
815 : {
816 : // constant
817 354944 : case 0:
818 :
819 : // linear
820 : case 1:
821 : case 2:
822 : case 3:
823 354944 : 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 226595 : case 11:
843 226595 : return 2.*xi;
844 :
845 226595 : case 12:
846 226595 : 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 155695 : case 21:
866 155695 : return 3.*xi*xi;
867 :
868 155695 : case 22:
869 155695 : return 4.*xi*eta;
870 :
871 155695 : case 23:
872 155695 : return 3.*eta*eta;
873 :
874 83364 : case 24:
875 : case 25:
876 83364 : return 0.;
877 :
878 155695 : case 26:
879 155695 : return 2.*xi*zeta;
880 :
881 155695 : case 27:
882 155695 : return 2.*eta*zeta;
883 :
884 83364 : case 28:
885 : case 29:
886 83364 : return 0.;
887 :
888 155695 : case 30:
889 155695 : 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 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
900 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
901 1726305 : int block=o, nz = 0;
902 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
903 1726305 : const int nx = block - i2;
904 1726305 : const int ny = o - nx - nz;
905 1726305 : Real val = nx * ny;
906 3370405 : for (int index=1; index < nx; index++)
907 2570074 : val *= xi;
908 3370405 : for (int index=1; index < ny; index++)
909 2570074 : val *= eta;
910 4603480 : for (int index=0; index != nz; index++)
911 3803149 : val *= zeta;
912 462987 : return val;
913 : }
914 : }
915 :
916 :
917 : // d^2()/deta^2
918 9343272 : case 2:
919 : {
920 9343272 : switch (i)
921 : {
922 : // constant
923 354944 : case 0:
924 :
925 : // linear
926 : case 1:
927 : case 2:
928 : case 3:
929 354944 : 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 226595 : case 12:
950 226595 : return 2.*xi;
951 226595 : case 13:
952 226595 : return 6.*eta;
953 :
954 121260 : case 14:
955 : case 15:
956 121260 : return 0.;
957 :
958 226595 : case 16:
959 226595 : 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 155695 : case 22:
972 155695 : return 2.*xi*xi;
973 :
974 155695 : case 23:
975 155695 : return 6.*xi*eta;
976 :
977 155695 : case 24:
978 155695 : return 12.*eta*eta;
979 :
980 83364 : case 25:
981 : case 26:
982 83364 : return 0.;
983 :
984 155695 : case 27:
985 155695 : return 2.*xi*zeta;
986 :
987 155695 : case 28:
988 155695 : return 6.*eta*zeta;
989 :
990 83364 : case 29:
991 : case 30:
992 83364 : return 0.;
993 :
994 155695 : case 31:
995 155695 : 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 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1005 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
1006 1726305 : int block=o, nz = 0;
1007 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
1008 1726305 : const int nx = block - i2;
1009 1726305 : const int ny = o - nx - nz;
1010 1726305 : Real val = ny * (ny - 1);
1011 4603480 : for (int index=0; index != nx; index++)
1012 3803149 : val *= xi;
1013 2548355 : for (int index=2; index < ny; index++)
1014 1748024 : val *= eta;
1015 4603480 : for (int index=0; index != nz; index++)
1016 3803149 : val *= zeta;
1017 462987 : return val;
1018 : }
1019 : }
1020 :
1021 :
1022 : // d^2()/dxidzeta
1023 9343272 : case 3:
1024 : {
1025 9343272 : switch (i)
1026 : {
1027 : // constant
1028 354944 : case 0:
1029 :
1030 : // linear
1031 : case 1:
1032 : case 2:
1033 : case 3:
1034 354944 : 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 226595 : case 14:
1057 226595 : return 2.*xi;
1058 :
1059 121260 : case 15:
1060 121260 : return eta;
1061 :
1062 60630 : case 16:
1063 60630 : return 0.;
1064 :
1065 226595 : case 17:
1066 226595 : 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 155695 : case 25:
1081 155695 : return 3.*xi*xi;
1082 :
1083 155695 : case 26:
1084 155695 : return 2.*xi*eta;
1085 :
1086 155695 : case 27:
1087 155695 : return eta*eta;
1088 :
1089 41682 : case 28:
1090 41682 : return 0.;
1091 :
1092 155695 : case 29:
1093 155695 : return 4.*xi*zeta;
1094 :
1095 155695 : case 30:
1096 155695 : return 2.*eta*zeta;
1097 :
1098 41682 : case 31:
1099 41682 : return 0.;
1100 :
1101 155695 : case 32:
1102 155695 : 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 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1111 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
1112 1726305 : int block=o, nz = 0;
1113 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
1114 1726305 : const int nx = block - i2;
1115 1726305 : const int ny = o - nx - nz;
1116 1726305 : Real val = nx * nz;
1117 3370405 : for (int index=1; index < nx; index++)
1118 2570074 : val *= xi;
1119 4603480 : for (int index=0; index != ny; index++)
1120 3803149 : val *= eta;
1121 3370405 : for (int index=1; index < nz; index++)
1122 2570074 : val *= zeta;
1123 462987 : return val;
1124 : }
1125 : }
1126 :
1127 : // d^2()/detadzeta
1128 9343272 : case 4:
1129 : {
1130 9343272 : switch (i)
1131 : {
1132 : // constant
1133 354944 : case 0:
1134 :
1135 : // linear
1136 : case 1:
1137 : case 2:
1138 : case 3:
1139 354944 : 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 226595 : case 16:
1166 226595 : return 2.*eta;
1167 :
1168 60630 : case 17:
1169 60630 : return 0.;
1170 :
1171 226595 : case 18:
1172 226595 : 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 155695 : case 26:
1187 155695 : return xi*xi;
1188 :
1189 155695 : case 27:
1190 155695 : return 2.*xi*eta;
1191 :
1192 155695 : case 28:
1193 155695 : return 3.*eta*eta;
1194 :
1195 41682 : case 29:
1196 41682 : return 0.;
1197 :
1198 155695 : case 30:
1199 155695 : return 2.*xi*zeta;
1200 :
1201 155695 : case 31:
1202 155695 : return 4.*eta*zeta;
1203 :
1204 41682 : case 32:
1205 41682 : return 0.;
1206 :
1207 155695 : case 33:
1208 155695 : return 3.*zeta*zeta;
1209 :
1210 41682 : case 34:
1211 41682 : return 0.;
1212 :
1213 462987 : default:
1214 462987 : unsigned int o = 0;
1215 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1216 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
1217 1726305 : int block=o, nz = 0;
1218 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
1219 1726305 : const int nx = block - i2;
1220 1726305 : const int ny = o - nx - nz;
1221 1726305 : Real val = ny * nz;
1222 4603480 : for (int index=0; index != nx; index++)
1223 3803149 : val *= xi;
1224 3370405 : for (int index=1; index < ny; index++)
1225 2570074 : val *= eta;
1226 3370405 : for (int index=1; index < nz; index++)
1227 2570074 : val *= zeta;
1228 462987 : return val;
1229 : }
1230 : }
1231 :
1232 :
1233 : // d^2()/dzeta^2
1234 9343272 : case 5:
1235 : {
1236 9343272 : switch (i)
1237 : {
1238 : // constant
1239 354944 : case 0:
1240 :
1241 : // linear
1242 : case 1:
1243 : case 2:
1244 : case 3:
1245 354944 : 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 226595 : case 17:
1269 226595 : return 2.*xi;
1270 :
1271 226595 : case 18:
1272 226595 : return 2.*eta;
1273 :
1274 226595 : case 19:
1275 226595 : 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 155695 : case 29:
1290 155695 : return 2.*xi*xi;
1291 :
1292 155695 : case 30:
1293 155695 : return 2.*xi*eta;
1294 :
1295 155695 : case 31:
1296 155695 : return 2.*eta*eta;
1297 :
1298 155695 : case 32:
1299 155695 : return 6.*xi*zeta;
1300 :
1301 155695 : case 33:
1302 155695 : return 6.*eta*zeta;
1303 :
1304 155695 : case 34:
1305 155695 : return 12.*zeta*zeta;
1306 :
1307 462987 : default:
1308 462987 : unsigned int o = 0;
1309 10357830 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1310 1726305 : const int i2 = i - (o*(o+1)*(o+2)/6);
1311 1726305 : int block=o, nz = 0;
1312 4603480 : for (; block < i2; block += (o-nz+1)) { nz++; }
1313 1726305 : const int nx = block - i2;
1314 1726305 : const int ny = o - nx - nz;
1315 1726305 : Real val = nz * (nz - 1);
1316 4603480 : for (int index=0; index != nx; index++)
1317 3803149 : val *= xi;
1318 4603480 : for (int index=0; index != ny; index++)
1319 3803149 : val *= eta;
1320 2548355 : for (int index=2; index < nz; index++)
1321 1748024 : 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 56059632 : 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 14982834 : libmesh_assert(elem);
1348 :
1349 : // call the orientation-independent shape function derivatives
1350 86025300 : 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
|
