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