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