1 #include <stdlib.h>
2
3 #include <FCMangle.h>
4 #define symhex FortranCInterface_GLOBAL_(symhex, SYMHEX)
5
6 typedef double DARR4[4];
7
8 int hexIntPnt(int,DARR4**,double **);
9
symhex(int * n1,double pt[][4],double wt[],int * err)10 void symhex(int *n1, double pt[][4], double wt[], int *err)
11 {
12 double *lwt;
13 DARR4 *lpt;
14 int i,j;
15 *err = hexIntPnt(*n1, &lpt, &lwt);
16 for(i=0; i < *n1; i++) {
17 wt[i] = lwt[i];
18 for(j=0; j <4; j++)
19 pt[i][j]=lpt[i][j];
20 }
21 }
22
23 /*$Id$*/
24 #include <stdio.h>
25
26 /* these are the rule 2 int points and weights */
27
28 #define Qp21 -0.577350269189626
29 #define Qp22 0.577350269189626
30 #define Qw2 1.000000000000000
31
32 /* these are the rule 3 int points and weights */
33
34 #define Qp31 -0.774596669241483
35 #define Qp32 0.000000000000000
36 #define Qp33 0.774596669241483
37 #define Qw31 0.555555555555556
38 #define Qw32 0.888888888888889
39 #define Qw33 0.555555555555556
40 #define Qw311 0.308641975308642
41 #define Qw321 0.493827160493831
42 #define Qw322 0.790123456790124
43
44
45 #define Qw3111 0.171467764060361
46 #define Qw3112 0.274348422496571
47 #define Qw3113 0.171467764060361
48 #define Qw3121 0.274348422496571
49 #define Qw3122 0.438957475994513
50 #define Qw3123 0.274348422496571
51 #define Qw3131 0.171467764060361
52 #define Qw3132 0.274348422496571
53 #define Qw3133 0.171467764060361
54
55 #define Qw3211 0.274348422496571
56 #define Qw3212 0.438957475994513
57 #define Qw3213 0.274348422496571
58 #define Qw3221 0.438957475994513
59 #define Qw3222 0.702331961591221
60 #define Qw3223 0.438957475994513
61 #define Qw3231 0.274348422496571
62 #define Qw3232 0.438957475994513
63 #define Qw3233 0.274348422496571
64
65 #define Qw3311 0.171467764060361
66 #define Qw3312 0.274348422496571
67 #define Qw3313 0.171467764060361
68 #define Qw3321 0.274348422496571
69 #define Qw3322 0.438957475994513
70 #define Qw3323 0.274348422496571
71 #define Qw3331 0.171467764060361
72 #define Qw3332 0.274348422496571
73 #define Qw3333 0.171467764060357
74
75 /* Rule 1*/
76
77 static double rstw1[][4] = {
78 { 0.0, 0.0, 0.0, 0.0 }
79 };
80
81 static double twt1[] = { 8.0000000 };
82
83
84 /* Rule 2*/
85
86
87 static double rstw8[][4] = {
88 {Qp21,Qp21,Qp21,0.0},
89 {Qp22,Qp21,Qp21,0.0},
90 {Qp21,Qp22,Qp21,0.0},
91 {Qp22,Qp22,Qp21,0.0},
92 {Qp21,Qp21,Qp22,0.0},
93 {Qp22,Qp21,Qp22,0.0},
94 {Qp21,Qp22,Qp22,0.0},
95 {Qp22,Qp22,Qp22,0.0}
96 };
97
98 static double twt8[] = {Qw2,Qw2,Qw2,Qw2,Qw2,Qw2,Qw2,Qw2};
99
100 /* Rule 3*/
101
102 static double rstw27[][4] = {
103 { Qp31, Qp31, Qp31, 0.0 },
104 { Qp32, Qp31, Qp31, 0.0 },
105 { Qp33, Qp31, Qp31, 0.0 },
106 { Qp31, Qp32, Qp31, 0.0 },
107 { Qp32, Qp32, Qp31, 0.0 },
108 { Qp33, Qp32, Qp31, 0.0 },
109 { Qp31, Qp33, Qp31, 0.0 },
110 { Qp32, Qp33, Qp31, 0.0 },
111 { Qp33, Qp33, Qp31, 0.0 },
112 { Qp31, Qp31, Qp32, 0.0 },
113 { Qp32, Qp31, Qp32, 0.0 },
114 { Qp33, Qp31, Qp32, 0.0 },
115 { Qp31, Qp32, Qp32, 0.0 },
116 { Qp32, Qp32, Qp32, 0.0 },
117 { Qp33, Qp32, Qp32, 0.0 },
118 { Qp31, Qp33, Qp32, 0.0 },
119 { Qp32, Qp33, Qp32, 0.0 },
120 { Qp33, Qp33, Qp32, 0.0 },
121 { Qp31, Qp31, Qp33, 0.0 },
122 { Qp32, Qp31, Qp33, 0.0 },
123 { Qp33, Qp31, Qp33, 0.0 },
124 { Qp31, Qp32, Qp33, 0.0 },
125 { Qp32, Qp32, Qp33, 0.0 },
126 { Qp33, Qp32, Qp33, 0.0 },
127 { Qp31, Qp33, Qp33, 0.0 },
128 { Qp32, Qp33, Qp33, 0.0 },
129 { Qp33, Qp33, Qp33, 0.0 }
130 };
131
132
133 static double twt27[] =
134 { Qw3111, Qw3211, Qw3311, Qw3121, Qw3221, Qw3321, Qw3131, Qw3231, Qw3331,
135 Qw3112, Qw3212, Qw3312, Qw3122, Qw3222, Qw3322, Qw3132, Qw3232, Qw3332,
136 Qw3113, Qw3213, Qw3313, Qw3123, Qw3223, Qw3323, Qw3133, Qw3233, Qw3333 };
137
138
139
140 #ifdef __cplusplus
141 extern "C" {
142 #endif
143
hexIntPnt(int nint,DARR4 ** bcord,double ** wt)144 int hexIntPnt(int nint, DARR4 **bcord, double **wt)
145 {
146 int retval = 1;
147 int i,j,k,l;
148 DARR4 *rstw;
149 double *twt;
150
151 /* Rule 4 & 5*/
152 /* rule 5 has not been tested */
153 double bp4[]={-0.8611363115940526,-0.3399810435848563,0.339981043584856,0.8611363115940526};
154 double bp5[]={-0.9061798459386640,-0.5384693101056831,0,0.5384693101056831,0.9061798459386640};
155 double bw4[]={0.3478548451374539,0.6521451548625461,0.6521451548625461,0.3478548451374539};
156 double bw5[]={0.2369268850561891,0.4786286704993665,0.5688888888888889,0.4786286704993665,0.2369268850561891};
157
158
159 if( nint == 1){ *bcord = rstw1; *wt = twt1;}
160
161 else if( nint == 8 ){*bcord = rstw8 ; *wt = twt8; }
162
163 else if( nint == 27 ) {*bcord = rstw27 ; *wt = twt27; }
164
165 else if( nint == 64 ) {
166
167 rstw = (DARR4 *)malloc(64*sizeof(DARR4));
168 twt = (double *)malloc(64*sizeof(double));
169
170 i=j=k=0;
171
172 for(l=0;l<64;l++){
173 rstw[l][0]=bp4[i];
174 rstw[l][1]=bp4[j];
175 rstw[l][2]=bp4[k];
176 rstw[l][3]=0.0;
177 twt[l]=bw4[i]*bw4[j]*bw4[k];
178 i++;
179 if(i == 4){
180 i=0 ; j++ ;
181 if(j == 4) { j=0; k++;}
182 }
183 }
184 *bcord = rstw; *wt = twt;
185
186 } else if( nint == 125) {
187
188 rstw = (DARR4 *)malloc(125*sizeof(DARR4));
189 twt = (double *)malloc(125*sizeof(double));
190
191 i=j=k=0;
192
193 for(l=0;l<125;l++){
194 rstw[l][0]=bp5[i];
195 rstw[l][1]=bp5[j];
196 rstw[l][2]=bp5[k];
197 rstw[l][3]=0.0;
198 twt[l]=bw5[i]*bw5[j]*bw5[k];
199 i++;
200 if(i == 5){
201 i=0 ; j++ ;
202 if(j == 5) { j=0; k++;}
203 }
204 }
205 *bcord = rstw; *wt = twt;
206
207 } else {
208
209 fprintf(stderr,"\n%d integration points unsupported; give {8,27,64,125,.}\n",nint);
210 retval = 0;
211 }
212 return retval ;
213 }
214
215 #ifdef __cplusplus
216 }
217 #endif
218