1 #include <stdlib.h>
2
3 #include <FCMangle.h>
4 #define sympyr FortranCInterface_GLOBAL_(sympyr, SYMPYR)
5
6 typedef double DARR4[4];
7
8 int pyrIntPnt(int,DARR4**,double **);
9
sympyr(int * n1,double pt[][4],double wt[],int * err)10 void sympyr(int *n1, double pt[][4], double wt[], int *err)
11 {
12 double *lwt;
13 DARR4 *lpt;
14 int i,j;
15 *err = pyrIntPnt(*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 mpt455 -0.455341801261479548
29 #define ppt455 0.455341801261479548
30 #define mpt577 -0.577350269189625764
31 #define ppt577 0.577350269189625764
32 #define mpt122 -0.122008467928146216
33 #define ppt122 0.122008467928146216
34 #define ppt622 0.622008467928146212
35 #define ppt044 0.0446581987385204510
36
37 /* these are the rule 3 int points and weights */
38
39 #define zero 0.0000000000000000
40 #define mpt687 -0.6872983346207417
41 #define ppt687 0.6872983346207417
42 #define mpt774 -0.7745966692414834
43 #define ppt774 0.7745966692414834
44 #define mpt387 -0.3872983346207416
45 #define ppt387 0.3872983346207416
46 #define mpt087 -0.08729833462074170
47 #define ppt087 0.08729833462074170
48 #define ppt134 0.1349962850858610
49 #define ppt215 0.2159940561373775
50 #define ppt345 0.3455904898198040
51 #define ppt068 0.06858710562414265
52 #define ppt109 0.1097393689986282
53 #define ppt175 0.1755829903978052
54 #define ppt002 0.002177926162424265
55 #define ppt003 0.003484681859878818
56 #define ppt005 0.005575490975806120
57
58 /* these are the rule 4 integration points */
59 /*NOT FIXED YET*/
60 #define Qp41 -0.801346029378026
61 #define Qp42 -0.316375532749811
62 #define Qp43 0.316375532749811
63 #define Qp44 0.801346029378026
64 #define Qp45 -0.8611363116
65
66 #define Qp46 -0.576953166749811
67 #define Qp47 -0.227784076803672
68 #define Qp48 0.227784076803672
69 #define Qp49 0.576953166749811
70 #define Qp410 -0.3399810436
71
72 #define Qp411 -0.284183144850188
73 #define Qp412 -0.112196966796327
74 #define Qp413 0.112196966796327
75 #define Qp414 0.284183144850188
76 #define Qp415 0.3399810436
77
78 #define Qp416 -0.0597902822219738
79 #define Qp417 -0.0236055108501886
80 #define Qp418 0.0236055108501886
81 #define Qp419 0.0597902822219738
82 #define Qp420 0.8611363116
83
84 /* Rule 1*/
85
86 static double rstw1[][4] = {
87 { 0.0, 0.0, 0.0, 0.0 }
88 };
89
90 static double twt1[] = { 8.0000000 };
91
92
93 /* Rule 2*/
94
95
96 static double rstw8[][4] = {
97 {mpt455,mpt455,mpt577,0.0},
98 {ppt455,mpt455,mpt577,0.0},
99 {mpt455,ppt455,mpt577,0.0},
100 {ppt455,ppt455,mpt577,0.0},
101 {mpt122,mpt122,ppt577,0.0},
102 {ppt122,mpt122,ppt577,0.0},
103 {mpt122,ppt122,ppt577,0.0},
104 {ppt122,ppt122,ppt577,0.0}
105 };
106
107 static double twt8[] = {ppt622,ppt622,ppt622,ppt622,
108 ppt044,ppt044,ppt044,ppt044};
109
110 /* Rule 3*/
111
112 static double rstw27[][4] = {
113 { mpt687, mpt687, mpt774, 0.0 },
114 { zero, mpt687, mpt774, 0.0 },
115 { ppt687, mpt687, mpt774, 0.0 },
116 { mpt687, zero, mpt774, 0.0 },
117 { zero, zero, mpt774, 0.0 },
118 { ppt687, zero, mpt774, 0.0 },
119 { mpt687, ppt687, mpt774, 0.0 },
120 { zero, ppt687, mpt774, 0.0 },
121 { ppt687, ppt687, mpt774, 0.0 },
122 { mpt387, mpt387, zero, 0.0 },
123 { zero, mpt387, zero, 0.0 },
124 { ppt387, mpt387, zero, 0.0 },
125 { mpt387, zero, zero, 0.0 },
126 { zero, zero, zero, 0.0 },
127 { ppt387, zero, zero, 0.0 },
128 { mpt387, ppt387, zero, 0.0 },
129 { zero, ppt387, zero, 0.0 },
130 { ppt387, ppt387, zero, 0.0 },
131 { mpt087, mpt087, ppt774, 0.0 },
132 { zero, mpt087, ppt774, 0.0 },
133 { ppt087, mpt087, ppt774, 0.0 },
134 { mpt087, zero, ppt774, 0.0 },
135 { zero, zero, ppt774, 0.0 },
136 { ppt087, zero, ppt774, 0.0 },
137 { mpt087, ppt087, ppt774, 0.0 },
138 { zero, ppt087, ppt774, 0.0 },
139 { ppt087, ppt087, ppt774, 0.0 }
140 };
141
142
143 static double twt27[] =
144 {0.1349962850858610, 0.2159940561373775, 0.1349962850858610,
145 0.2159940561373775, 0.3455904898198040, 0.2159940561373775,
146 0.1349962850858610, 0.2159940561373775, 0.1349962850858610,
147 0.06858710562414265, 0.1097393689986282, 0.06858710562414265,
148 0.1097393689986282, 0.1755829903978052, 0.1097393689986282,
149 0.06858710562414265, 0.1097393689986282, 0.06858710562414265,
150 0.002177926162424265, 0.003484681859878818, 0.002177926162424265,
151 0.003484681859878822, 0.005575490975806120, 0.003484681859878825,
152 0.002177926162424265, 0.003484681859878822, 0.002177926162424265};
153
154 /* Rule 4 */
155
156 static double rstw64[][4] = {
157 { Qp41, Qp41, Qp45, 0.0 },
158 { Qp42, Qp41, Qp45, 0.0 },
159 { Qp43, Qp41, Qp45, 0.0 },
160 { Qp44, Qp41, Qp45, 0.0 },
161 { Qp41, Qp42, Qp45, 0.0 },
162 { Qp42, Qp42, Qp45, 0.0 },
163 { Qp43, Qp42, Qp45, 0.0 },
164 { Qp44, Qp42, Qp45, 0.0 },
165 { Qp41, Qp43, Qp45, 0.0 },
166 { Qp42, Qp43, Qp45, 0.0 },
167 { Qp43, Qp43, Qp45, 0.0 },
168 { Qp44, Qp43, Qp45, 0.0 },
169 { Qp41, Qp44, Qp45, 0.0 },
170 { Qp42, Qp44, Qp45, 0.0 },
171 { Qp43, Qp44, Qp45, 0.0 },
172 { Qp44, Qp44, Qp45, 0.0 },
173 { Qp46, Qp46, Qp410, 0.0 },
174 { Qp47, Qp46, Qp410, 0.0 },
175 { Qp48, Qp46, Qp410, 0.0 },
176 { Qp49, Qp46, Qp410, 0.0 },
177 { Qp46, Qp47, Qp410, 0.0 },
178 { Qp47, Qp47, Qp410, 0.0 },
179 { Qp48, Qp47, Qp410, 0.0 },
180 { Qp49, Qp47, Qp410, 0.0 },
181 { Qp46, Qp48, Qp410, 0.0 },
182 { Qp47, Qp48, Qp410, 0.0 },
183 { Qp48, Qp48, Qp410, 0.0 },
184 { Qp49, Qp48, Qp410, 0.0 },
185 { Qp46, Qp49, Qp410, 0.0 },
186 { Qp47, Qp49, Qp410, 0.0 },
187 { Qp48, Qp49, Qp410, 0.0 },
188 { Qp49, Qp49, Qp410, 0.0 },
189 { Qp411, Qp411, Qp415, 0.0 },
190 { Qp412, Qp411, Qp415, 0.0 },
191 { Qp413, Qp411, Qp415, 0.0 },
192 { Qp414, Qp411, Qp415, 0.0 },
193 { Qp411, Qp412, Qp415, 0.0 },
194 { Qp412, Qp412, Qp415, 0.0 },
195 { Qp413, Qp412, Qp415, 0.0 },
196 { Qp414, Qp412, Qp415, 0.0 },
197 { Qp411, Qp413, Qp415, 0.0 },
198 { Qp412, Qp413, Qp415, 0.0 },
199 { Qp413, Qp413, Qp415, 0.0 },
200 { Qp414, Qp413, Qp415, 0.0 },
201 { Qp411, Qp414, Qp415, 0.0 },
202 { Qp412, Qp414, Qp415, 0.0 },
203 { Qp413, Qp414, Qp415, 0.0 },
204 { Qp414, Qp414, Qp415, 0.0 },
205 { Qp416, Qp416, Qp420, 0.0 },
206 { Qp417, Qp416, Qp420, 0.0 },
207 { Qp418, Qp416, Qp420, 0.0 },
208 { Qp419, Qp416, Qp420, 0.0 },
209 { Qp416, Qp417, Qp420, 0.0 },
210 { Qp417, Qp417, Qp420, 0.0 },
211 { Qp418, Qp417, Qp420, 0.0 },
212 { Qp419, Qp417, Qp420, 0.0 },
213 { Qp416, Qp418, Qp420, 0.0 },
214 { Qp417, Qp418, Qp420, 0.0 },
215 { Qp418, Qp418, Qp420, 0.0 },
216 { Qp419, Qp418, Qp420, 0.0 },
217 { Qp416, Qp419, Qp420, 0.0 },
218 { Qp417, Qp419, Qp420, 0.0 },
219 { Qp418, Qp419, Qp420, 0.0 },
220 { Qp419, Qp419, Qp420, 0.0 }
221 };
222
223 #ifdef __cplusplus
224 extern "C" {
225 #endif
226
pyrIntPnt(int nint,DARR4 ** bcord,double ** wt)227 int pyrIntPnt(int nint, DARR4 **bcord, double **wt)
228 {
229 int retval = 1;
230 int i,j,k,l;
231 DARR4 *rstw;
232 double *twt;
233
234 /* Rule 4 & 5*/
235
236 double bp5[]={-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459};
237 double bw4[]={0.3478548446,0.6521451548,0.6521451548,0.3478548446};
238 double bw5[]={0.2369268850,0.4786286708,0.5019607843,0.4786286708,0.2369268850};
239
240 if( nint == 1){ *bcord = rstw1; *wt = twt1;}
241
242 else if( nint == 8 ){*bcord = rstw8 ; *wt = twt8; }
243
244 else if( nint == 27 ) {*bcord = rstw27 ; *wt = twt27; }
245
246 else if( nint == 64 ) {
247
248 /* rstw = (DARR4 *)malloc(64*sizeof(DARR4)); */
249 twt = (double *)malloc(64*sizeof(double));
250
251 i=j=k=0;
252
253 for(l=0;l<64;l++){
254 twt[l]=bw4[i]*bw4[j]*bw4[k];
255 i++;
256 if(i == 4){
257 i=0 ; j++ ;
258 if(j == 4) { j=0; k++;}
259 }
260 }
261 *bcord = rstw64; *wt = twt;
262
263 } else if( nint == 125) {
264
265 rstw = (DARR4 *)malloc(125*sizeof(DARR4));
266 twt = (double *)malloc(125*sizeof(double));
267
268 i=j=k=0;
269
270 for(l=0;l<125;l++){
271 rstw[l][0]=bp5[i];
272 rstw[l][1]=bp5[j];
273 rstw[l][2]=bp5[k];
274 rstw[l][3]=0.0;
275 twt[l]=bw5[i]*bw5[j]*bw5[k];
276 i++;
277 if(i == 5){
278 i=0 ; j++ ;
279 if(j == 5) { j=0; k++;}
280 }
281 }
282 *bcord = rstw; *wt = twt;
283
284 } else {
285
286 fprintf(stderr,"\n%d integration points unsupported; give {8,27,64,125,.}\n",nint);
287 retval = 0;
288 }
289 return retval ;
290 }
291
292 #ifdef __cplusplus
293 }
294 #endif
295