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