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