libMesh
fe_xyz_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 
24 #include "libmesh/elem.h"
25 #include "libmesh/int_range.h"
26 
27 
28 namespace libMesh
29 {
30 
31 
32 LIBMESH_DEFAULT_VECTORIZED_FE(3,XYZ)
33 
34 
35 template <>
36 Real 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 
215 template <>
216 Real 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 
227 template <>
228 Real 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 
241 template <>
242 Real 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 
703 template <>
704 Real 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 
716 template <>
717 Real 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 
731 template <>
732 Real 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 
1411 template <>
1412 Real 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 
1423 template <>
1424 Real 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
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3109
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
distance
Real distance(const Point &p)
Definition: subdomains_ex3.C:51
libMesh::LIBMESH_DEFAULT_VECTORIZED_FE
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
Definition: fe_bernstein_shape_0D.C:30
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:328
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::libmesh_assert
libmesh_assert(ctx)
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:2666
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:140
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2454
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::Elem::vertex_average
Point vertex_average() const
Definition: elem.C:686
Generated on Thu Oct 30 2025 06:31:51 for libMesh by   doxygen 1.8.14