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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo