1 /*$Id$*/ 2 #include <math.h> 3 4 #ifdef __cplusplus 5 extern "C" { 6 #endif 7 8 void brickToTet(double xi,double eta, double zeta, 9 double *r, double *s, double *t, double *J); 10 11 void quadToTri(double xi,double eta, 12 double *r, double *s,double *J); 13 double quadToTriJac(double xi,double eta); 14 int GaussLegendre1D(int, double **,double **); 15 16 int GaussLegendreTet(int n1,int n2,int n3,double GLrstw[][4],double *GLwt) { 17 /* get degenerate n1Xn2Xn3 Gauss-Legendre scheme to integrate over a tet */ 18 int i,j,k,index=0; 19 double *pt1,*pt2,*pt3,*wt1,*wt2,*wt3,dJ; 20 const double six=6.000000000000000; 21 22 GaussLegendre1D(n1,&pt1,&wt1); 23 GaussLegendre1D(n2,&pt2,&wt2); 24 GaussLegendre1D(n3,&pt3,&wt3); 25 for(i=0; i < n1; i++) { 26 for(j=0; j < n2; j++) { 27 for(k=0; k < n3; k++) { 28 brickToTet(pt1[i],pt2[j],pt3[k],&GLrstw[index][0], 29 &GLrstw[index][1],&GLrstw[index][2],&dJ); 30 GLrstw[index][3] = 1.0e0-GLrstw[index][0]-GLrstw[index][1] 31 -GLrstw[index][2]; 32 GLwt[index++] = dJ*wt1[i]*wt2[j]*wt3[k]*six; 33 } 34 } 35 } 36 return index; 37 } 38 int GaussLegendreTri(int n1,int n2,double GLr[][3],double *GLwt) { 39 /* get degenerate n1Xn2 Gauss-Legendre scheme to integrate over a tet */ 40 int i,j,index=0; 41 double *pt1,*pt2,*wt1,*wt2,dJ; 42 const double two = 2.0000000000000000; 43 44 GaussLegendre1D(n1,&pt1,&wt1); 45 GaussLegendre1D(n2,&pt2,&wt2); 46 for(i=0; i < n1; i++) { 47 for(j=0; j < n2; j++) { 48 quadToTri(pt1[i],pt2[j],&GLr[index][0],&GLr[index][1],&dJ); 49 GLr[index][2] = 1.0e0-GLr[index][0]-GLr[index][1]; 50 GLwt[index++] = dJ*wt1[i]*wt2[j]*two; 51 } 52 } 53 return index; 54 } 55 56 void brickToTet(double xi,double eta, double zeta, 57 double *r, double *s, double *t, double *J) { 58 double r1,rs1; 59 *r = 0.5e0*(1.0e0+xi); 60 r1 = 1.0e0-(*r); 61 *s = 0.5e0*(1.0e0+eta)*r1; 62 rs1 = 1.0e0-(*r)-(*s); 63 *t = 0.5e0*(1.0e0+zeta)*rs1; 64 *J = 0.125e0*r1*rs1; 65 } 66 67 void quadToTri(double xi,double eta,double *r, double *s, double *J) { 68 double r1; 69 *r = 0.5e0*(1.0e0+xi); 70 r1 = 1.0e0-(*r); 71 *s = 0.5e0*(1.0e0+eta)*r1; 72 *J = 0.25e0*r1; 73 } 74 75 double quadToTriJac(double xi,double eta) { 76 return 0.125e0*(1.0e0-eta); 77 } 78 79 #ifdef __cplusplus 80 } 81 #endif 82 83