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