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 39043963 : LIBMESH_DEFAULT_VECTORIZED_FE(3,MONOMIAL)
31 :
32 :
33 : template <>
34 40879683 : 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 40879683 : const Real xi = p(0);
42 40879683 : const Real eta = p(1);
43 40879683 : const Real zeta = p(2);
44 :
45 10428747 : 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 40879683 : switch (i)
51 : {
52 : // constant
53 1005665 : case 0:
54 1005665 : return 1.;
55 :
56 : // linears
57 2557736 : case 1:
58 2557736 : return xi;
59 :
60 2557736 : case 2:
61 2557736 : return eta;
62 :
63 2557736 : case 3:
64 2557736 : return zeta;
65 :
66 : // quadratics
67 928728 : case 4:
68 928728 : return xi*xi;
69 :
70 928728 : case 5:
71 928728 : return xi*eta;
72 :
73 928728 : case 6:
74 928728 : return eta*eta;
75 :
76 928728 : case 7:
77 928728 : return xi*zeta;
78 :
79 928728 : case 8:
80 928728 : return zeta*eta;
81 :
82 928728 : case 9:
83 928728 : return zeta*zeta;
84 :
85 : // cubics
86 843665 : case 10:
87 843665 : return xi*xi*xi;
88 :
89 843665 : case 11:
90 843665 : return xi*xi*eta;
91 :
92 843665 : case 12:
93 843665 : return xi*eta*eta;
94 :
95 843665 : case 13:
96 843665 : return eta*eta*eta;
97 :
98 843665 : case 14:
99 843665 : return xi*xi*zeta;
100 :
101 843665 : case 15:
102 843665 : return xi*eta*zeta;
103 :
104 843665 : case 16:
105 843665 : return eta*eta*zeta;
106 :
107 843665 : case 17:
108 843665 : return xi*zeta*zeta;
109 :
110 843665 : case 18:
111 843665 : return eta*zeta*zeta;
112 :
113 843665 : case 19:
114 843665 : return zeta*zeta*zeta;
115 :
116 : // quartics
117 548002 : case 20:
118 548002 : return xi*xi*xi*xi;
119 :
120 548002 : case 21:
121 548002 : return xi*xi*xi*eta;
122 :
123 548002 : case 22:
124 548002 : return xi*xi*eta*eta;
125 :
126 548002 : case 23:
127 548002 : return xi*eta*eta*eta;
128 :
129 548002 : case 24:
130 548002 : return eta*eta*eta*eta;
131 :
132 548002 : case 25:
133 548002 : return xi*xi*xi*zeta;
134 :
135 548002 : case 26:
136 548002 : return xi*xi*eta*zeta;
137 :
138 548002 : case 27:
139 548002 : return xi*eta*eta*zeta;
140 :
141 548002 : case 28:
142 548002 : return eta*eta*eta*zeta;
143 :
144 548002 : case 29:
145 548002 : return xi*xi*zeta*zeta;
146 :
147 548002 : case 30:
148 548002 : return xi*eta*zeta*zeta;
149 :
150 548002 : case 31:
151 548002 : return eta*eta*zeta*zeta;
152 :
153 548002 : case 32:
154 548002 : return xi*zeta*zeta*zeta;
155 :
156 548002 : case 33:
157 548002 : return eta*zeta*zeta*zeta;
158 :
159 548002 : case 34:
160 548002 : return zeta*zeta*zeta*zeta;
161 :
162 1785735 : default:
163 1785735 : unsigned int o = 0;
164 41926122 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
165 6987687 : const int i2 = i - (o*(o+1)*(o+2)/6);
166 6987687 : int block=o, nz = 0;
167 18633832 : for (; block < i2; block += (o-nz+1)) { nz++; }
168 6987687 : const int nx = block - i2;
169 6987687 : const int ny = o - nx - nz;
170 1785735 : Real val = 1.;
171 18633832 : for (int index=0; index != nx; index++)
172 15217615 : val *= xi;
173 18633832 : for (int index=0; index != ny; index++)
174 15217615 : val *= eta;
175 18633832 : for (int index=0; index != nz; index++)
176 15217615 : val *= zeta;
177 1785735 : 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 40879683 : 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 10428747 : libmesh_assert(elem);
197 :
198 : // call the orientation-independent shape functions
199 61737177 : 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 35371728 : 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 9395241 : libmesh_assert_less (j, 3);
228 :
229 9395241 : 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 35371728 : const Real xi = p(0);
235 35371728 : const Real eta = p(1);
236 35371728 : const Real zeta = p(2);
237 :
238 : // monomials. since they are hierarchic we only need one case block.
239 35371728 : switch (j)
240 : {
241 : // d()/dxi
242 11790576 : case 0:
243 : {
244 11790576 : switch (i)
245 : {
246 : // constant
247 188114 : case 0:
248 188114 : return 0.;
249 :
250 : // linear
251 376052 : case 1:
252 376052 : return 1.;
253 :
254 188026 : case 2:
255 188026 : return 0.;
256 :
257 188026 : case 3:
258 188026 : return 0.;
259 :
260 : // quadratic
261 311914 : case 4:
262 311914 : return 2.*xi;
263 :
264 166794 : case 5:
265 166794 : return eta;
266 :
267 83397 : case 6:
268 83397 : return 0.;
269 :
270 166794 : case 7:
271 166794 : return zeta;
272 :
273 83397 : case 8:
274 83397 : return 0.;
275 :
276 83397 : case 9:
277 83397 : return 0.;
278 :
279 : // cubic
280 252107 : case 10:
281 252107 : return 3.*xi*xi;
282 :
283 252107 : case 11:
284 252107 : return 2.*xi*eta;
285 :
286 252107 : case 12:
287 252107 : return eta*eta;
288 :
289 67590 : case 13:
290 67590 : return 0.;
291 :
292 252107 : case 14:
293 252107 : return 2.*xi*zeta;
294 :
295 252107 : case 15:
296 252107 : return eta*zeta;
297 :
298 67590 : case 16:
299 67590 : return 0.;
300 :
301 252107 : case 17:
302 252107 : return zeta*zeta;
303 :
304 67590 : case 18:
305 67590 : return 0.;
306 :
307 67590 : case 19:
308 67590 : return 0.;
309 :
310 : // quartics
311 172255 : case 20:
312 172255 : return 4.*xi*xi*xi;
313 :
314 172255 : case 21:
315 172255 : return 3.*xi*xi*eta;
316 :
317 172255 : case 22:
318 172255 : return 2.*xi*eta*eta;
319 :
320 172255 : case 23:
321 172255 : return eta*eta*eta;
322 :
323 46194 : case 24:
324 46194 : return 0.;
325 :
326 172255 : case 25:
327 172255 : return 3.*xi*xi*zeta;
328 :
329 172255 : case 26:
330 172255 : return 2.*xi*eta*zeta;
331 :
332 172255 : case 27:
333 172255 : return eta*eta*zeta;
334 :
335 46194 : case 28:
336 46194 : return 0.;
337 :
338 172255 : case 29:
339 172255 : return 2.*xi*zeta*zeta;
340 :
341 172255 : case 30:
342 172255 : return eta*zeta*zeta;
343 :
344 46194 : case 31:
345 46194 : return 0.;
346 :
347 172255 : case 32:
348 172255 : return zeta*zeta*zeta;
349 :
350 46194 : case 33:
351 46194 : return 0.;
352 :
353 46194 : case 34:
354 46194 : return 0.;
355 :
356 510363 : default:
357 510363 : unsigned int o = 0;
358 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
359 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
360 1900185 : int block=o, nz = 0;
361 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
362 1900185 : const int nx = block - i2;
363 1900185 : const int ny = o - nx - nz;
364 1900185 : Real val = nx;
365 3709885 : for (int index=1; index < nx; index++)
366 2830426 : val *= xi;
367 5067160 : for (int index=0; index != ny; index++)
368 4187701 : val *= eta;
369 5067160 : for (int index=0; index != nz; index++)
370 4187701 : val *= zeta;
371 510363 : return val;
372 : }
373 : }
374 :
375 :
376 : // d()/deta
377 11790576 : case 1:
378 : {
379 11790576 : switch (i)
380 : {
381 : // constant
382 188114 : case 0:
383 188114 : return 0.;
384 :
385 : // linear
386 188026 : case 1:
387 188026 : return 0.;
388 :
389 376052 : case 2:
390 376052 : return 1.;
391 :
392 188026 : case 3:
393 188026 : return 0.;
394 :
395 : // quadratic
396 83397 : case 4:
397 83397 : return 0.;
398 :
399 166794 : case 5:
400 166794 : return xi;
401 :
402 311914 : case 6:
403 311914 : return 2.*eta;
404 :
405 83397 : case 7:
406 83397 : return 0.;
407 :
408 166794 : case 8:
409 166794 : return zeta;
410 :
411 83397 : case 9:
412 83397 : return 0.;
413 :
414 : // cubic
415 67590 : case 10:
416 67590 : return 0.;
417 :
418 252107 : case 11:
419 252107 : return xi*xi;
420 :
421 252107 : case 12:
422 252107 : return 2.*xi*eta;
423 :
424 252107 : case 13:
425 252107 : return 3.*eta*eta;
426 :
427 67590 : case 14:
428 67590 : return 0.;
429 :
430 252107 : case 15:
431 252107 : return xi*zeta;
432 :
433 252107 : case 16:
434 252107 : return 2.*eta*zeta;
435 :
436 67590 : case 17:
437 67590 : return 0.;
438 :
439 252107 : case 18:
440 252107 : return zeta*zeta;
441 :
442 67590 : case 19:
443 67590 : return 0.;
444 :
445 : // quartics
446 46194 : case 20:
447 46194 : return 0.;
448 :
449 172255 : case 21:
450 172255 : return xi*xi*xi;
451 :
452 172255 : case 22:
453 172255 : return 2.*xi*xi*eta;
454 :
455 172255 : case 23:
456 172255 : return 3.*xi*eta*eta;
457 :
458 172255 : case 24:
459 172255 : return 4.*eta*eta*eta;
460 :
461 46194 : case 25:
462 46194 : return 0.;
463 :
464 172255 : case 26:
465 172255 : return xi*xi*zeta;
466 :
467 172255 : case 27:
468 172255 : return 2.*xi*eta*zeta;
469 :
470 172255 : case 28:
471 172255 : return 3.*eta*eta*zeta;
472 :
473 46194 : case 29:
474 46194 : return 0.;
475 :
476 172255 : case 30:
477 172255 : return xi*zeta*zeta;
478 :
479 172255 : case 31:
480 172255 : return 2.*eta*zeta*zeta;
481 :
482 46194 : case 32:
483 46194 : return 0.;
484 :
485 172255 : case 33:
486 172255 : return zeta*zeta*zeta;
487 :
488 46194 : case 34:
489 46194 : return 0.;
490 :
491 510363 : default:
492 510363 : unsigned int o = 0;
493 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
494 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
495 1900185 : int block=o, nz = 0;
496 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
497 1900185 : const int nx = block - i2;
498 1900185 : const int ny = o - nx - nz;
499 1900185 : Real val = ny;
500 5067160 : for (int index=0; index != nx; index++)
501 4187701 : val *= xi;
502 3709885 : for (int index=1; index < ny; index++)
503 2830426 : val *= eta;
504 5067160 : for (int index=0; index != nz; index++)
505 4187701 : val *= zeta;
506 510363 : return val;
507 : }
508 : }
509 :
510 :
511 : // d()/dzeta
512 11790576 : case 2:
513 : {
514 11790576 : switch (i)
515 : {
516 : // constant
517 188114 : case 0:
518 188114 : return 0.;
519 :
520 : // linear
521 188026 : case 1:
522 188026 : return 0.;
523 :
524 188026 : case 2:
525 188026 : return 0.;
526 :
527 376052 : case 3:
528 376052 : return 1.;
529 :
530 : // quadratic
531 83397 : case 4:
532 83397 : return 0.;
533 :
534 83397 : case 5:
535 83397 : return 0.;
536 :
537 83397 : case 6:
538 83397 : return 0.;
539 :
540 166794 : case 7:
541 166794 : return xi;
542 :
543 166794 : case 8:
544 166794 : return eta;
545 :
546 311914 : case 9:
547 311914 : return 2.*zeta;
548 :
549 : // cubic
550 67590 : case 10:
551 67590 : return 0.;
552 :
553 67590 : case 11:
554 67590 : return 0.;
555 :
556 67590 : case 12:
557 67590 : return 0.;
558 :
559 67590 : case 13:
560 67590 : return 0.;
561 :
562 252107 : case 14:
563 252107 : return xi*xi;
564 :
565 252107 : case 15:
566 252107 : return xi*eta;
567 :
568 252107 : case 16:
569 252107 : return eta*eta;
570 :
571 252107 : case 17:
572 252107 : return 2.*xi*zeta;
573 :
574 252107 : case 18:
575 252107 : return 2.*eta*zeta;
576 :
577 252107 : case 19:
578 252107 : return 3.*zeta*zeta;
579 :
580 : // quartics
581 46194 : case 20:
582 46194 : return 0.;
583 :
584 46194 : case 21:
585 46194 : return 0.;
586 :
587 46194 : case 22:
588 46194 : return 0.;
589 :
590 46194 : case 23:
591 46194 : return 0.;
592 :
593 46194 : case 24:
594 46194 : return 0.;
595 :
596 172255 : case 25:
597 172255 : return xi*xi*xi;
598 :
599 172255 : case 26:
600 172255 : return xi*xi*eta;
601 :
602 172255 : case 27:
603 172255 : return xi*eta*eta;
604 :
605 172255 : case 28:
606 172255 : return eta*eta*eta;
607 :
608 172255 : case 29:
609 172255 : return 2.*xi*xi*zeta;
610 :
611 172255 : case 30:
612 172255 : return 2.*xi*eta*zeta;
613 :
614 172255 : case 31:
615 172255 : return 2.*eta*eta*zeta;
616 :
617 172255 : case 32:
618 172255 : return 3.*xi*zeta*zeta;
619 :
620 172255 : case 33:
621 172255 : return 3.*eta*zeta*zeta;
622 :
623 172255 : case 34:
624 172255 : return 4.*zeta*zeta*zeta;
625 :
626 510363 : default:
627 510363 : unsigned int o = 0;
628 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
629 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
630 1900185 : int block=o, nz = 0;
631 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
632 1900185 : const int nx = block - i2;
633 1900185 : const int ny = o - nx - nz;
634 1900185 : Real val = nz;
635 5067160 : for (int index=0; index != nx; index++)
636 4187701 : val *= xi;
637 5067160 : for (int index=0; index != ny; index++)
638 4187701 : val *= eta;
639 3709885 : for (int index=1; index < nz; index++)
640 2830426 : val *= zeta;
641 510363 : 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 35371728 : 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 9395241 : libmesh_assert(elem);
667 :
668 : // call the orientation-independent shape function derivatives
669 54162210 : 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 62150112 : 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 16642002 : libmesh_assert_less (j, 6);
699 :
700 16642002 : 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 62150112 : const Real xi = p(0);
705 62150112 : const Real eta = p(1);
706 62150112 : const Real zeta = p(2);
707 :
708 : // monomials. since they are hierarchic we only need one case block.
709 62150112 : switch (j)
710 : {
711 : // d^2()/dxi^2
712 10358352 : case 0:
713 : {
714 10358352 : switch (i)
715 : {
716 : // constant
717 394112 : case 0:
718 :
719 : // linear
720 : case 1:
721 : case 2:
722 : case 3:
723 394112 : return 0.;
724 :
725 : // quadratic
726 166794 : case 4:
727 166794 : return 2.;
728 :
729 416985 : case 5:
730 : case 6:
731 : case 7:
732 : case 8:
733 : case 9:
734 416985 : return 0.;
735 :
736 : // cubic
737 252107 : case 10:
738 252107 : return 6.*xi;
739 :
740 252107 : case 11:
741 252107 : return 2.*eta;
742 :
743 135180 : case 12:
744 : case 13:
745 135180 : return 0.;
746 :
747 252107 : case 14:
748 252107 : return 2.*zeta;
749 :
750 337950 : case 15:
751 : case 16:
752 : case 17:
753 : case 18:
754 : case 19:
755 337950 : return 0.;
756 :
757 : // quartics
758 172255 : case 20:
759 172255 : return 12.*xi*xi;
760 :
761 172255 : case 21:
762 172255 : return 6.*xi*eta;
763 :
764 172255 : case 22:
765 172255 : return 2.*eta*eta;
766 :
767 92388 : case 23:
768 : case 24:
769 92388 : return 0.;
770 :
771 172255 : case 25:
772 172255 : return 6.*xi*zeta;
773 :
774 172255 : case 26:
775 172255 : return 2.*eta*zeta;
776 :
777 92388 : case 27:
778 : case 28:
779 92388 : return 0.;
780 :
781 172255 : case 29:
782 172255 : return 2.*zeta*zeta;
783 :
784 230970 : case 30:
785 : case 31:
786 : case 32:
787 : case 33:
788 : case 34:
789 230970 : return 0.;
790 :
791 510363 : default:
792 510363 : unsigned int o = 0;
793 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
794 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
795 1900185 : int block=o, nz = 0;
796 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
797 1900185 : const int nx = block - i2;
798 1900185 : const int ny = o - nx - nz;
799 1900185 : Real val = nx * (nx - 1);
800 2805035 : for (int index=2; index < nx; index++)
801 1925576 : val *= xi;
802 5067160 : for (int index=0; index != ny; index++)
803 4187701 : val *= eta;
804 5067160 : for (int index=0; index != nz; index++)
805 4187701 : val *= zeta;
806 510363 : return val;
807 : }
808 : }
809 :
810 :
811 : // d^2()/dxideta
812 10358352 : case 1:
813 : {
814 10358352 : switch (i)
815 : {
816 : // constant
817 394112 : case 0:
818 :
819 : // linear
820 : case 1:
821 : case 2:
822 : case 3:
823 394112 : return 0.;
824 :
825 : // quadratic
826 83397 : case 4:
827 83397 : return 0.;
828 :
829 166794 : case 5:
830 166794 : return 1.;
831 :
832 333588 : case 6:
833 : case 7:
834 : case 8:
835 : case 9:
836 333588 : return 0.;
837 :
838 : // cubic
839 67590 : case 10:
840 67590 : return 0.;
841 :
842 252107 : case 11:
843 252107 : return 2.*xi;
844 :
845 252107 : case 12:
846 252107 : return 2.*eta;
847 :
848 135180 : case 13:
849 : case 14:
850 135180 : return 0.;
851 :
852 135180 : case 15:
853 135180 : return zeta;
854 :
855 270360 : case 16:
856 : case 17:
857 : case 18:
858 : case 19:
859 270360 : return 0.;
860 :
861 : // quartics
862 46194 : case 20:
863 46194 : return 0.;
864 :
865 172255 : case 21:
866 172255 : return 3.*xi*xi;
867 :
868 172255 : case 22:
869 172255 : return 4.*xi*eta;
870 :
871 172255 : case 23:
872 172255 : return 3.*eta*eta;
873 :
874 92388 : case 24:
875 : case 25:
876 92388 : return 0.;
877 :
878 172255 : case 26:
879 172255 : return 2.*xi*zeta;
880 :
881 172255 : case 27:
882 172255 : return 2.*eta*zeta;
883 :
884 92388 : case 28:
885 : case 29:
886 92388 : return 0.;
887 :
888 172255 : case 30:
889 172255 : return zeta*zeta;
890 :
891 184776 : case 31:
892 : case 32:
893 : case 33:
894 : case 34:
895 184776 : return 0.;
896 :
897 510363 : default:
898 510363 : unsigned int o = 0;
899 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
900 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
901 1900185 : int block=o, nz = 0;
902 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
903 1900185 : const int nx = block - i2;
904 1900185 : const int ny = o - nx - nz;
905 1900185 : Real val = nx * ny;
906 3709885 : for (int index=1; index < nx; index++)
907 2830426 : val *= xi;
908 3709885 : for (int index=1; index < ny; index++)
909 2830426 : val *= eta;
910 5067160 : for (int index=0; index != nz; index++)
911 4187701 : val *= zeta;
912 510363 : return val;
913 : }
914 : }
915 :
916 :
917 : // d^2()/deta^2
918 10358352 : case 2:
919 : {
920 10358352 : switch (i)
921 : {
922 : // constant
923 394112 : case 0:
924 :
925 : // linear
926 : case 1:
927 : case 2:
928 : case 3:
929 394112 : return 0.;
930 :
931 : // quadratic
932 166794 : case 4:
933 : case 5:
934 166794 : return 0.;
935 :
936 166794 : case 6:
937 166794 : return 2.;
938 :
939 250191 : case 7:
940 : case 8:
941 : case 9:
942 250191 : return 0.;
943 :
944 : // cubic
945 135180 : case 10:
946 : case 11:
947 135180 : return 0.;
948 :
949 252107 : case 12:
950 252107 : return 2.*xi;
951 252107 : case 13:
952 252107 : return 6.*eta;
953 :
954 135180 : case 14:
955 : case 15:
956 135180 : return 0.;
957 :
958 252107 : case 16:
959 252107 : return 2.*zeta;
960 :
961 202770 : case 17:
962 : case 18:
963 : case 19:
964 202770 : return 0.;
965 :
966 : // quartics
967 92388 : case 20:
968 : case 21:
969 92388 : return 0.;
970 :
971 172255 : case 22:
972 172255 : return 2.*xi*xi;
973 :
974 172255 : case 23:
975 172255 : return 6.*xi*eta;
976 :
977 172255 : case 24:
978 172255 : return 12.*eta*eta;
979 :
980 92388 : case 25:
981 : case 26:
982 92388 : return 0.;
983 :
984 172255 : case 27:
985 172255 : return 2.*xi*zeta;
986 :
987 172255 : case 28:
988 172255 : return 6.*eta*zeta;
989 :
990 92388 : case 29:
991 : case 30:
992 92388 : return 0.;
993 :
994 172255 : case 31:
995 172255 : return 2.*zeta*zeta;
996 :
997 138582 : case 32:
998 : case 33:
999 : case 34:
1000 138582 : return 0.;
1001 :
1002 510363 : default:
1003 510363 : unsigned int o = 0;
1004 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1005 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
1006 1900185 : int block=o, nz = 0;
1007 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
1008 1900185 : const int nx = block - i2;
1009 1900185 : const int ny = o - nx - nz;
1010 1900185 : Real val = ny * (ny - 1);
1011 5067160 : for (int index=0; index != nx; index++)
1012 4187701 : val *= xi;
1013 2805035 : for (int index=2; index < ny; index++)
1014 1925576 : val *= eta;
1015 5067160 : for (int index=0; index != nz; index++)
1016 4187701 : val *= zeta;
1017 510363 : return val;
1018 : }
1019 : }
1020 :
1021 :
1022 : // d^2()/dxidzeta
1023 10358352 : case 3:
1024 : {
1025 10358352 : switch (i)
1026 : {
1027 : // constant
1028 394112 : case 0:
1029 :
1030 : // linear
1031 : case 1:
1032 : case 2:
1033 : case 3:
1034 394112 : return 0.;
1035 :
1036 : // quadratic
1037 250191 : case 4:
1038 : case 5:
1039 : case 6:
1040 250191 : return 0.;
1041 :
1042 166794 : case 7:
1043 166794 : return 1.;
1044 :
1045 166794 : case 8:
1046 : case 9:
1047 166794 : return 0.;
1048 :
1049 : // cubic
1050 270360 : case 10:
1051 : case 11:
1052 : case 12:
1053 : case 13:
1054 270360 : return 0.;
1055 :
1056 252107 : case 14:
1057 252107 : return 2.*xi;
1058 :
1059 135180 : case 15:
1060 135180 : return eta;
1061 :
1062 67590 : case 16:
1063 67590 : return 0.;
1064 :
1065 252107 : case 17:
1066 252107 : return 2.*zeta;
1067 :
1068 135180 : case 18:
1069 : case 19:
1070 135180 : return 0.;
1071 :
1072 : // quartics
1073 230970 : case 20:
1074 : case 21:
1075 : case 22:
1076 : case 23:
1077 : case 24:
1078 230970 : return 0.;
1079 :
1080 172255 : case 25:
1081 172255 : return 3.*xi*xi;
1082 :
1083 172255 : case 26:
1084 172255 : return 2.*xi*eta;
1085 :
1086 172255 : case 27:
1087 172255 : return eta*eta;
1088 :
1089 46194 : case 28:
1090 46194 : return 0.;
1091 :
1092 172255 : case 29:
1093 172255 : return 4.*xi*zeta;
1094 :
1095 172255 : case 30:
1096 172255 : return 2.*eta*zeta;
1097 :
1098 46194 : case 31:
1099 46194 : return 0.;
1100 :
1101 172255 : case 32:
1102 172255 : return 3.*zeta*zeta;
1103 :
1104 92388 : case 33:
1105 : case 34:
1106 92388 : return 0.;
1107 :
1108 510363 : default:
1109 510363 : unsigned int o = 0;
1110 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1111 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
1112 1900185 : int block=o, nz = 0;
1113 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
1114 1900185 : const int nx = block - i2;
1115 1900185 : const int ny = o - nx - nz;
1116 1900185 : Real val = nx * nz;
1117 3709885 : for (int index=1; index < nx; index++)
1118 2830426 : val *= xi;
1119 5067160 : for (int index=0; index != ny; index++)
1120 4187701 : val *= eta;
1121 3709885 : for (int index=1; index < nz; index++)
1122 2830426 : val *= zeta;
1123 510363 : return val;
1124 : }
1125 : }
1126 :
1127 : // d^2()/detadzeta
1128 10358352 : case 4:
1129 : {
1130 10358352 : switch (i)
1131 : {
1132 : // constant
1133 394112 : case 0:
1134 :
1135 : // linear
1136 : case 1:
1137 : case 2:
1138 : case 3:
1139 394112 : return 0.;
1140 :
1141 : // quadratic
1142 333588 : case 4:
1143 : case 5:
1144 : case 6:
1145 : case 7:
1146 333588 : return 0.;
1147 :
1148 166794 : case 8:
1149 166794 : return 1.;
1150 :
1151 83397 : case 9:
1152 83397 : return 0.;
1153 :
1154 : // cubic
1155 337950 : case 10:
1156 : case 11:
1157 : case 12:
1158 : case 13:
1159 : case 14:
1160 337950 : return 0.;
1161 :
1162 135180 : case 15:
1163 135180 : return xi;
1164 :
1165 252107 : case 16:
1166 252107 : return 2.*eta;
1167 :
1168 67590 : case 17:
1169 67590 : return 0.;
1170 :
1171 252107 : case 18:
1172 252107 : return 2.*zeta;
1173 :
1174 67590 : case 19:
1175 67590 : return 0.;
1176 :
1177 : // quartics
1178 277164 : case 20:
1179 : case 21:
1180 : case 22:
1181 : case 23:
1182 : case 24:
1183 : case 25:
1184 277164 : return 0.;
1185 :
1186 172255 : case 26:
1187 172255 : return xi*xi;
1188 :
1189 172255 : case 27:
1190 172255 : return 2.*xi*eta;
1191 :
1192 172255 : case 28:
1193 172255 : return 3.*eta*eta;
1194 :
1195 46194 : case 29:
1196 46194 : return 0.;
1197 :
1198 172255 : case 30:
1199 172255 : return 2.*xi*zeta;
1200 :
1201 172255 : case 31:
1202 172255 : return 4.*eta*zeta;
1203 :
1204 46194 : case 32:
1205 46194 : return 0.;
1206 :
1207 172255 : case 33:
1208 172255 : return 3.*zeta*zeta;
1209 :
1210 46194 : case 34:
1211 46194 : return 0.;
1212 :
1213 510363 : default:
1214 510363 : unsigned int o = 0;
1215 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1216 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
1217 1900185 : int block=o, nz = 0;
1218 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
1219 1900185 : const int nx = block - i2;
1220 1900185 : const int ny = o - nx - nz;
1221 1900185 : Real val = ny * nz;
1222 5067160 : for (int index=0; index != nx; index++)
1223 4187701 : val *= xi;
1224 3709885 : for (int index=1; index < ny; index++)
1225 2830426 : val *= eta;
1226 3709885 : for (int index=1; index < nz; index++)
1227 2830426 : val *= zeta;
1228 510363 : return val;
1229 : }
1230 : }
1231 :
1232 :
1233 : // d^2()/dzeta^2
1234 10358352 : case 5:
1235 : {
1236 10358352 : switch (i)
1237 : {
1238 : // constant
1239 394112 : case 0:
1240 :
1241 : // linear
1242 : case 1:
1243 : case 2:
1244 : case 3:
1245 394112 : return 0.;
1246 :
1247 : // quadratic
1248 416985 : case 4:
1249 : case 5:
1250 : case 6:
1251 : case 7:
1252 : case 8:
1253 416985 : return 0.;
1254 :
1255 166794 : case 9:
1256 166794 : return 2.;
1257 :
1258 : // cubic
1259 473130 : case 10:
1260 : case 11:
1261 : case 12:
1262 : case 13:
1263 : case 14:
1264 : case 15:
1265 : case 16:
1266 473130 : return 0.;
1267 :
1268 252107 : case 17:
1269 252107 : return 2.*xi;
1270 :
1271 252107 : case 18:
1272 252107 : return 2.*eta;
1273 :
1274 252107 : case 19:
1275 252107 : return 6.*zeta;
1276 :
1277 : // quartics
1278 415746 : 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 415746 : return 0.;
1288 :
1289 172255 : case 29:
1290 172255 : return 2.*xi*xi;
1291 :
1292 172255 : case 30:
1293 172255 : return 2.*xi*eta;
1294 :
1295 172255 : case 31:
1296 172255 : return 2.*eta*eta;
1297 :
1298 172255 : case 32:
1299 172255 : return 6.*xi*zeta;
1300 :
1301 172255 : case 33:
1302 172255 : return 6.*eta*zeta;
1303 :
1304 172255 : case 34:
1305 172255 : return 12.*zeta*zeta;
1306 :
1307 510363 : default:
1308 510363 : unsigned int o = 0;
1309 11401110 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1310 1900185 : const int i2 = i - (o*(o+1)*(o+2)/6);
1311 1900185 : int block=o, nz = 0;
1312 5067160 : for (; block < i2; block += (o-nz+1)) { nz++; }
1313 1900185 : const int nx = block - i2;
1314 1900185 : const int ny = o - nx - nz;
1315 1900185 : Real val = nz * (nz - 1);
1316 5067160 : for (int index=0; index != nx; index++)
1317 4187701 : val *= xi;
1318 5067160 : for (int index=0; index != ny; index++)
1319 4187701 : val *= eta;
1320 2805035 : for (int index=2; index < nz; index++)
1321 1925576 : val *= zeta;
1322 510363 : 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 62150112 : 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 16642002 : libmesh_assert(elem);
1348 :
1349 : // call the orientation-independent shape function derivatives
1350 95434116 : 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
|
