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 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 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