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 125692439 : LIBMESH_DEFAULT_VECTORIZED_FE(3,MONOMIAL)
31 :
32 :
33 : template <>
34 122207771 : 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 122207771 : const Real xi = p(0);
42 122207771 : const Real eta = p(1);
43 122207771 : const Real zeta = p(2);
44 :
45 10030055 : 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 122207771 : switch (i)
51 : {
52 : // constant
53 995350 : case 0:
54 995350 : return 1.;
55 :
56 : // linears
57 7608693 : case 1:
58 7608693 : return xi;
59 :
60 7608693 : case 2:
61 7608693 : return eta;
62 :
63 7608693 : case 3:
64 7608693 : return zeta;
65 :
66 : // quadratics
67 2811176 : case 4:
68 2811176 : return xi*xi;
69 :
70 2811176 : case 5:
71 2811176 : return xi*eta;
72 :
73 2811176 : case 6:
74 2811176 : return eta*eta;
75 :
76 2811176 : case 7:
77 2811176 : return xi*zeta;
78 :
79 2811176 : case 8:
80 2811176 : return zeta*eta;
81 :
82 2811176 : case 9:
83 2811176 : return zeta*zeta;
84 :
85 : // cubics
86 2545166 : case 10:
87 2545166 : return xi*xi*xi;
88 :
89 2545166 : case 11:
90 2545166 : return xi*xi*eta;
91 :
92 2545166 : case 12:
93 2545166 : return xi*eta*eta;
94 :
95 2545166 : case 13:
96 2545166 : return eta*eta*eta;
97 :
98 2545166 : case 14:
99 2545166 : return xi*xi*zeta;
100 :
101 2545166 : case 15:
102 2545166 : return xi*eta*zeta;
103 :
104 2545166 : case 16:
105 2545166 : return eta*eta*zeta;
106 :
107 2545166 : case 17:
108 2545166 : return xi*zeta*zeta;
109 :
110 2545166 : case 18:
111 2545166 : return eta*zeta*zeta;
112 :
113 2545166 : case 19:
114 2545166 : return zeta*zeta*zeta;
115 :
116 : // quartics
117 1630092 : case 20:
118 1630092 : return xi*xi*xi*xi;
119 :
120 1630092 : case 21:
121 1630092 : return xi*xi*xi*eta;
122 :
123 1630092 : case 22:
124 1630092 : return xi*xi*eta*eta;
125 :
126 1630092 : case 23:
127 1630092 : return xi*eta*eta*eta;
128 :
129 1630092 : case 24:
130 1630092 : return eta*eta*eta*eta;
131 :
132 1630092 : case 25:
133 1630092 : return xi*xi*xi*zeta;
134 :
135 1630092 : case 26:
136 1630092 : return xi*xi*eta*zeta;
137 :
138 1630092 : case 27:
139 1630092 : return xi*eta*eta*zeta;
140 :
141 1630092 : case 28:
142 1630092 : return eta*eta*eta*zeta;
143 :
144 1630092 : case 29:
145 1630092 : return xi*xi*zeta*zeta;
146 :
147 1630092 : case 30:
148 1630092 : return xi*eta*zeta*zeta;
149 :
150 1630092 : case 31:
151 1630092 : return eta*eta*zeta*zeta;
152 :
153 1630092 : case 32:
154 1630092 : return xi*zeta*zeta*zeta;
155 :
156 1630092 : case 33:
157 1630092 : return eta*zeta*zeta*zeta;
158 :
159 1630092 : case 34:
160 1630092 : return zeta*zeta*zeta*zeta;
161 :
162 1690479 : default:
163 1690479 : unsigned int o = 0;
164 124090092 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
165 20681682 : const int i2 = i - (o*(o+1)*(o+2)/6);
166 20681682 : int block=o, nz = 0;
167 55151152 : for (; block < i2; block += (o-nz+1)) { nz++; }
168 20681682 : const int nx = block - i2;
169 20681682 : const int ny = o - nx - nz;
170 1690479 : Real val = 1.;
171 55151152 : for (int index=0; index != nx; index++)
172 37850428 : val *= xi;
173 55151152 : for (int index=0; index != ny; index++)
174 37850428 : val *= eta;
175 55151152 : for (int index=0; index != nz; index++)
176 37850428 : val *= zeta;
177 1690479 : 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 122207771 : 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 10030055 : libmesh_assert(elem);
197 :
198 : // call the orientation-independent shape functions
199 142267881 : 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 117706860 : 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 117706860 : const Real xi = p(0);
235 117706860 : const Real eta = p(1);
236 117706860 : const Real zeta = p(2);
237 :
238 : // monomials. since they are hierarchic we only need one case block.
239 117706860 : switch (j)
240 : {
241 : // d()/dxi
242 39235620 : case 0:
243 : {
244 39235620 : 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 1051096 : case 4:
262 1051096 : 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 855866 : case 10:
281 855866 : return 3.*xi*xi;
282 :
283 855866 : case 11:
284 855866 : return 2.*xi*eta;
285 :
286 855866 : case 12:
287 855866 : return eta*eta;
288 :
289 67590 : case 13:
290 67590 : return 0.;
291 :
292 855866 : case 14:
293 855866 : return 2.*xi*zeta;
294 :
295 855866 : case 15:
296 855866 : return eta*zeta;
297 :
298 67590 : case 16:
299 67590 : return 0.;
300 :
301 855866 : case 17:
302 855866 : 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 577246 : case 20:
312 577246 : return 4.*xi*xi*xi;
313 :
314 577246 : case 21:
315 577246 : return 3.*xi*xi*eta;
316 :
317 577246 : case 22:
318 577246 : return 2.*xi*eta*eta;
319 :
320 577246 : case 23:
321 577246 : return eta*eta*eta;
322 :
323 46194 : case 24:
324 46194 : return 0.;
325 :
326 577246 : case 25:
327 577246 : return 3.*xi*xi*zeta;
328 :
329 577246 : case 26:
330 577246 : return 2.*xi*eta*zeta;
331 :
332 577246 : case 27:
333 577246 : return eta*eta*zeta;
334 :
335 46194 : case 28:
336 46194 : return 0.;
337 :
338 577246 : case 29:
339 577246 : return 2.*xi*zeta*zeta;
340 :
341 577246 : case 30:
342 577246 : return eta*zeta*zeta;
343 :
344 46194 : case 31:
345 46194 : return 0.;
346 :
347 577246 : case 32:
348 577246 : 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 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
359 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
360 6426546 : int block=o, nz = 0;
361 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
362 6426546 : const int nx = block - i2;
363 6426546 : const int ny = o - nx - nz;
364 6426546 : Real val = nx;
365 12547066 : for (int index=1; index < nx; index++)
366 7141246 : val *= xi;
367 17137456 : for (int index=0; index != ny; index++)
368 11731636 : val *= eta;
369 17137456 : for (int index=0; index != nz; index++)
370 11731636 : val *= zeta;
371 510363 : return val;
372 : }
373 : }
374 :
375 :
376 : // d()/deta
377 39235620 : case 1:
378 : {
379 39235620 : 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 1051096 : case 6:
403 1051096 : 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 855866 : case 11:
419 855866 : return xi*xi;
420 :
421 855866 : case 12:
422 855866 : return 2.*xi*eta;
423 :
424 855866 : case 13:
425 855866 : return 3.*eta*eta;
426 :
427 67590 : case 14:
428 67590 : return 0.;
429 :
430 855866 : case 15:
431 855866 : return xi*zeta;
432 :
433 855866 : case 16:
434 855866 : return 2.*eta*zeta;
435 :
436 67590 : case 17:
437 67590 : return 0.;
438 :
439 855866 : case 18:
440 855866 : 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 577246 : case 21:
450 577246 : return xi*xi*xi;
451 :
452 577246 : case 22:
453 577246 : return 2.*xi*xi*eta;
454 :
455 577246 : case 23:
456 577246 : return 3.*xi*eta*eta;
457 :
458 577246 : case 24:
459 577246 : return 4.*eta*eta*eta;
460 :
461 46194 : case 25:
462 46194 : return 0.;
463 :
464 577246 : case 26:
465 577246 : return xi*xi*zeta;
466 :
467 577246 : case 27:
468 577246 : return 2.*xi*eta*zeta;
469 :
470 577246 : case 28:
471 577246 : return 3.*eta*eta*zeta;
472 :
473 46194 : case 29:
474 46194 : return 0.;
475 :
476 577246 : case 30:
477 577246 : return xi*zeta*zeta;
478 :
479 577246 : case 31:
480 577246 : return 2.*eta*zeta*zeta;
481 :
482 46194 : case 32:
483 46194 : return 0.;
484 :
485 577246 : case 33:
486 577246 : return zeta*zeta*zeta;
487 :
488 46194 : case 34:
489 46194 : return 0.;
490 :
491 510363 : default:
492 510363 : unsigned int o = 0;
493 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
494 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
495 6426546 : int block=o, nz = 0;
496 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
497 6426546 : const int nx = block - i2;
498 6426546 : const int ny = o - nx - nz;
499 6426546 : Real val = ny;
500 17137456 : for (int index=0; index != nx; index++)
501 11731636 : val *= xi;
502 12547066 : for (int index=1; index < ny; index++)
503 7141246 : val *= eta;
504 17137456 : for (int index=0; index != nz; index++)
505 11731636 : val *= zeta;
506 510363 : return val;
507 : }
508 : }
509 :
510 :
511 : // d()/dzeta
512 39235620 : case 2:
513 : {
514 39235620 : 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 1051096 : case 9:
547 1051096 : 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 855866 : case 14:
563 855866 : return xi*xi;
564 :
565 855866 : case 15:
566 855866 : return xi*eta;
567 :
568 855866 : case 16:
569 855866 : return eta*eta;
570 :
571 855866 : case 17:
572 855866 : return 2.*xi*zeta;
573 :
574 855866 : case 18:
575 855866 : return 2.*eta*zeta;
576 :
577 855866 : case 19:
578 855866 : 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 577246 : case 25:
597 577246 : return xi*xi*xi;
598 :
599 577246 : case 26:
600 577246 : return xi*xi*eta;
601 :
602 577246 : case 27:
603 577246 : return xi*eta*eta;
604 :
605 577246 : case 28:
606 577246 : return eta*eta*eta;
607 :
608 577246 : case 29:
609 577246 : return 2.*xi*xi*zeta;
610 :
611 577246 : case 30:
612 577246 : return 2.*xi*eta*zeta;
613 :
614 577246 : case 31:
615 577246 : return 2.*eta*eta*zeta;
616 :
617 577246 : case 32:
618 577246 : return 3.*xi*zeta*zeta;
619 :
620 577246 : case 33:
621 577246 : return 3.*eta*zeta*zeta;
622 :
623 577246 : case 34:
624 577246 : return 4.*zeta*zeta*zeta;
625 :
626 510363 : default:
627 510363 : unsigned int o = 0;
628 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
629 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
630 6426546 : int block=o, nz = 0;
631 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
632 6426546 : const int nx = block - i2;
633 6426546 : const int ny = o - nx - nz;
634 6426546 : Real val = nz;
635 17137456 : for (int index=0; index != nx; index++)
636 11731636 : val *= xi;
637 17137456 : for (int index=0; index != ny; index++)
638 11731636 : val *= eta;
639 12547066 : for (int index=1; index < nz; index++)
640 7141246 : 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 117706860 : 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 136497342 : 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 209604888 : 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 209604888 : const Real xi = p(0);
705 209604888 : const Real eta = p(1);
706 209604888 : const Real zeta = p(2);
707 :
708 : // monomials. since they are hierarchic we only need one case block.
709 209604888 : switch (j)
710 : {
711 : // d^2()/dxi^2
712 34934148 : case 0:
713 : {
714 34934148 : 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 855866 : case 10:
738 855866 : return 6.*xi;
739 :
740 855866 : case 11:
741 855866 : return 2.*eta;
742 :
743 135180 : case 12:
744 : case 13:
745 135180 : return 0.;
746 :
747 855866 : case 14:
748 855866 : 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 577246 : case 20:
759 577246 : return 12.*xi*xi;
760 :
761 577246 : case 21:
762 577246 : return 6.*xi*eta;
763 :
764 577246 : case 22:
765 577246 : return 2.*eta*eta;
766 :
767 92388 : case 23:
768 : case 24:
769 92388 : return 0.;
770 :
771 577246 : case 25:
772 577246 : return 6.*xi*zeta;
773 :
774 577246 : case 26:
775 577246 : return 2.*eta*zeta;
776 :
777 92388 : case 27:
778 : case 28:
779 92388 : return 0.;
780 :
781 577246 : case 29:
782 577246 : 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 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
794 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
795 6426546 : int block=o, nz = 0;
796 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
797 6426546 : const int nx = block - i2;
798 6426546 : const int ny = o - nx - nz;
799 6426546 : Real val = nx * (nx - 1);
800 9486806 : for (int index=2; index < nx; index++)
801 4080986 : val *= xi;
802 17137456 : for (int index=0; index != ny; index++)
803 11731636 : val *= eta;
804 17137456 : for (int index=0; index != nz; index++)
805 11731636 : val *= zeta;
806 510363 : return val;
807 : }
808 : }
809 :
810 :
811 : // d^2()/dxideta
812 34934148 : case 1:
813 : {
814 34934148 : 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 855866 : case 11:
843 855866 : return 2.*xi;
844 :
845 855866 : case 12:
846 855866 : 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 577246 : case 21:
866 577246 : return 3.*xi*xi;
867 :
868 577246 : case 22:
869 577246 : return 4.*xi*eta;
870 :
871 577246 : case 23:
872 577246 : return 3.*eta*eta;
873 :
874 92388 : case 24:
875 : case 25:
876 92388 : return 0.;
877 :
878 577246 : case 26:
879 577246 : return 2.*xi*zeta;
880 :
881 577246 : case 27:
882 577246 : return 2.*eta*zeta;
883 :
884 92388 : case 28:
885 : case 29:
886 92388 : return 0.;
887 :
888 577246 : case 30:
889 577246 : 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 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
900 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
901 6426546 : int block=o, nz = 0;
902 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
903 6426546 : const int nx = block - i2;
904 6426546 : const int ny = o - nx - nz;
905 6426546 : Real val = nx * ny;
906 12547066 : for (int index=1; index < nx; index++)
907 7141246 : val *= xi;
908 12547066 : for (int index=1; index < ny; index++)
909 7141246 : val *= eta;
910 17137456 : for (int index=0; index != nz; index++)
911 11731636 : val *= zeta;
912 510363 : return val;
913 : }
914 : }
915 :
916 :
917 : // d^2()/deta^2
918 34934148 : case 2:
919 : {
920 34934148 : 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 855866 : case 12:
950 855866 : return 2.*xi;
951 855866 : case 13:
952 855866 : return 6.*eta;
953 :
954 135180 : case 14:
955 : case 15:
956 135180 : return 0.;
957 :
958 855866 : case 16:
959 855866 : 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 577246 : case 22:
972 577246 : return 2.*xi*xi;
973 :
974 577246 : case 23:
975 577246 : return 6.*xi*eta;
976 :
977 577246 : case 24:
978 577246 : return 12.*eta*eta;
979 :
980 92388 : case 25:
981 : case 26:
982 92388 : return 0.;
983 :
984 577246 : case 27:
985 577246 : return 2.*xi*zeta;
986 :
987 577246 : case 28:
988 577246 : return 6.*eta*zeta;
989 :
990 92388 : case 29:
991 : case 30:
992 92388 : return 0.;
993 :
994 577246 : case 31:
995 577246 : 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 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1005 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
1006 6426546 : int block=o, nz = 0;
1007 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
1008 6426546 : const int nx = block - i2;
1009 6426546 : const int ny = o - nx - nz;
1010 6426546 : Real val = ny * (ny - 1);
1011 17137456 : for (int index=0; index != nx; index++)
1012 11731636 : val *= xi;
1013 9486806 : for (int index=2; index < ny; index++)
1014 4080986 : val *= eta;
1015 17137456 : for (int index=0; index != nz; index++)
1016 11731636 : val *= zeta;
1017 510363 : return val;
1018 : }
1019 : }
1020 :
1021 :
1022 : // d^2()/dxidzeta
1023 34934148 : case 3:
1024 : {
1025 34934148 : 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 855866 : case 14:
1057 855866 : return 2.*xi;
1058 :
1059 135180 : case 15:
1060 135180 : return eta;
1061 :
1062 67590 : case 16:
1063 67590 : return 0.;
1064 :
1065 855866 : case 17:
1066 855866 : 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 577246 : case 25:
1081 577246 : return 3.*xi*xi;
1082 :
1083 577246 : case 26:
1084 577246 : return 2.*xi*eta;
1085 :
1086 577246 : case 27:
1087 577246 : return eta*eta;
1088 :
1089 46194 : case 28:
1090 46194 : return 0.;
1091 :
1092 577246 : case 29:
1093 577246 : return 4.*xi*zeta;
1094 :
1095 577246 : case 30:
1096 577246 : return 2.*eta*zeta;
1097 :
1098 46194 : case 31:
1099 46194 : return 0.;
1100 :
1101 577246 : case 32:
1102 577246 : 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 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1111 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
1112 6426546 : int block=o, nz = 0;
1113 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
1114 6426546 : const int nx = block - i2;
1115 6426546 : const int ny = o - nx - nz;
1116 6426546 : Real val = nx * nz;
1117 12547066 : for (int index=1; index < nx; index++)
1118 7141246 : val *= xi;
1119 17137456 : for (int index=0; index != ny; index++)
1120 11731636 : val *= eta;
1121 12547066 : for (int index=1; index < nz; index++)
1122 7141246 : val *= zeta;
1123 510363 : return val;
1124 : }
1125 : }
1126 :
1127 : // d^2()/detadzeta
1128 34934148 : case 4:
1129 : {
1130 34934148 : 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 855866 : case 16:
1166 855866 : return 2.*eta;
1167 :
1168 67590 : case 17:
1169 67590 : return 0.;
1170 :
1171 855866 : case 18:
1172 855866 : 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 577246 : case 26:
1187 577246 : return xi*xi;
1188 :
1189 577246 : case 27:
1190 577246 : return 2.*xi*eta;
1191 :
1192 577246 : case 28:
1193 577246 : return 3.*eta*eta;
1194 :
1195 46194 : case 29:
1196 46194 : return 0.;
1197 :
1198 577246 : case 30:
1199 577246 : return 2.*xi*zeta;
1200 :
1201 577246 : case 31:
1202 577246 : return 4.*eta*zeta;
1203 :
1204 46194 : case 32:
1205 46194 : return 0.;
1206 :
1207 577246 : case 33:
1208 577246 : return 3.*zeta*zeta;
1209 :
1210 46194 : case 34:
1211 46194 : return 0.;
1212 :
1213 510363 : default:
1214 510363 : unsigned int o = 0;
1215 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1216 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
1217 6426546 : int block=o, nz = 0;
1218 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
1219 6426546 : const int nx = block - i2;
1220 6426546 : const int ny = o - nx - nz;
1221 6426546 : Real val = ny * nz;
1222 17137456 : for (int index=0; index != nx; index++)
1223 11731636 : val *= xi;
1224 12547066 : for (int index=1; index < ny; index++)
1225 7141246 : val *= eta;
1226 12547066 : for (int index=1; index < nz; index++)
1227 7141246 : val *= zeta;
1228 510363 : return val;
1229 : }
1230 : }
1231 :
1232 :
1233 : // d^2()/dzeta^2
1234 34934148 : case 5:
1235 : {
1236 34934148 : 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 855866 : case 17:
1269 855866 : return 2.*xi;
1270 :
1271 855866 : case 18:
1272 855866 : return 2.*eta;
1273 :
1274 855866 : case 19:
1275 855866 : 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 577246 : case 29:
1290 577246 : return 2.*xi*xi;
1291 :
1292 577246 : case 30:
1293 577246 : return 2.*xi*eta;
1294 :
1295 577246 : case 31:
1296 577246 : return 2.*eta*eta;
1297 :
1298 577246 : case 32:
1299 577246 : return 6.*xi*zeta;
1300 :
1301 577246 : case 33:
1302 577246 : return 6.*eta*zeta;
1303 :
1304 577246 : case 34:
1305 577246 : return 12.*zeta*zeta;
1306 :
1307 510363 : default:
1308 510363 : unsigned int o = 0;
1309 38559276 : for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { }
1310 6426546 : const int i2 = i - (o*(o+1)*(o+2)/6);
1311 6426546 : int block=o, nz = 0;
1312 17137456 : for (; block < i2; block += (o-nz+1)) { nz++; }
1313 6426546 : const int nx = block - i2;
1314 6426546 : const int ny = o - nx - nz;
1315 6426546 : Real val = nz * (nz - 1);
1316 17137456 : for (int index=0; index != nx; index++)
1317 11731636 : val *= xi;
1318 17137456 : for (int index=0; index != ny; index++)
1319 11731636 : val *= eta;
1320 9486806 : for (int index=2; index < nz; index++)
1321 4080986 : 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 209604888 : 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 242888892 : 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
|
