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
GaussLegendreTet(int n1,int n2,int n3,double GLrstw[][4],double * GLwt)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 }
GaussLegendreTri(int n1,int n2,double GLr[][3],double * GLwt)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
brickToTet(double xi,double eta,double zeta,double * r,double * s,double * t,double * J)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
quadToTri(double xi,double eta,double * r,double * s,double * J)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
quadToTriJac(double xi,double eta)75 double quadToTriJac(double xi,double eta) {
76 return 0.125e0*(1.0e0-eta);
77 }
78
79 #ifdef __cplusplus
80 }
81 #endif
82
83