1 #include <FCMangle.h> 2 #define symwdg FortranCInterface_GLOBAL_(symwdg,SYMWDG) 3 4 typedef double DARR4[4]; 5 6 int wdgIntPnt(int,DARR4**,double **); 7 8 void symwdg(int *n1, double pt[][4], double wt[], int *err) 9 { 10 double *lwt; 11 DARR4 *lpt; 12 int i,j; 13 *err = wdgIntPnt(*n1, &lpt, &lwt); 14 for(i=0; i < *n1; i++) { 15 wt[i] = lwt[i]; 16 for(j=0; j <4; j++) 17 pt[i][j]=lpt[i][j]; 18 } 19 } 20 21 /*$Id$*/ 22 #include <stdio.h> 23 24 /* constants for 3D - 2pt rule */ 25 26 #define Qp21 -0.577350269189626 27 #define Qp22 0.577350269189626 28 #define Pp23 0.166666666666667 29 #define Pp24 0.666666666666667 30 #define Ww32 0.666666666666667 /* Pw22 * Qw2 */ 31 32 /* constants for 3D - rule 3 (18 point) */ 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 41 #define Pp34 0.816847572980459 42 #define Pp35 0.091576213509771 43 #define Pp36 0.108103018168070 44 #define Pp37 0.445948490915965 45 #define Pw33 0.219903487310644 46 #define Pw34 0.446763179356022 47 48 #define Ww3331 0.122168604061469 /* Pw33 * Qw31 */ 49 #define Ww3332 0.195469766498350 /* Pw33 * Qw32 */ 50 #define Ww3333 0.122168604061469 /* Pw33 * Qw33 */ 51 #define Ww3431 0.248201766308901 /* Pw34 * Qw31 */ 52 #define Ww3432 0.397122826094242 /* Pw34 * Qw32 */ 53 #define Ww3433 0.248201766308901 /* Pw34 * Qw33 */ 54 55 /* constants for 3D rule 4 (48 point) */ 56 57 #define Qp41 -0.861136311594053 58 #define Qp42 -0.339981043584856 59 #define Qp43 0.339981043584856 60 #define Qp44 0.861136311594053 61 #define Qw41 0.347854845137454 62 #define Qw42 0.652145154862544 63 #define Qw43 0.652145154862544 64 #define Qw44 0.347854845137544 65 66 #define Pp41 0.873821971016996 67 #define Pp42 0.063089014491502 68 #define Pp43 0.501426509658179 69 #define Pp44 0.249286745170910 70 #define Pp45 0.636502499121399 71 #define Pp46 0.310352451033785 72 #define Pp47 0.053145049844816 73 #define Pw41 0.101689812740414 /* S in symtri.c times 2 */ 74 #define Pw42 0.233572551452758 /* T in symtri.c times 2 */ 75 #define Pw43 0.165702151236748 /* U in symtri.c times 2 */ 76 77 /* the weights should go */ 78 /*(Pw41,Pw41,Pw41,Pw42,Pw42,Pw42,Pw43,Pw43,Pw43,Pw43,Pw43,Pw43)*Qw41 */ 79 /* >> *Qw42,Qw43,Qw44 */ 80 81 #define Ww4141 0.0353732940628734 /* Pw41 * Qw41 */ 82 #define Ww4142 0.0663165186775404 /* Pw41 * Qw42 */ 83 #define Ww4143 0.0663165186775404 /* Pw41 * Qw43 */ 84 #define Ww4144 0.0353732940628734 /* Pw41 * Qw44 */ 85 #define Ww4241 0.0812493437139591 /* Pw42 * Qw41 */ 86 #define Ww4242 0.152323207738798 /* Pw42 * Qw42 */ 87 #define Ww4243 0.152323207738798 /* Pw42 * Qw43 */ 88 #define Ww4244 0.0812493437139591 /* Pw42 * Qw44 */ 89 #define Ww4341 0.057640296157402 /* Pw43 * Qw41 */ 90 #define Ww4342 0.108061855079346 /* Pw43 * Qw42 */ 91 #define Ww4343 0.108061855079346 /* Pw43 * Qw43 */ 92 #define Ww4344 0.057640296157402 /* Pw43 * Qw44 */ 93 94 /* typedef double DARR4[4] ; */ 95 96 /* Rule 1*/ 97 98 static double rstw1[][4] = { 99 { 0.333333333333, 0.333333333333, 0.0, 0.0 } 100 }; 101 102 static double twt1[] = { 4.0000000 }; 103 104 105 106 /* Rule 2 */ 107 108 static double rstw6[][4] = { 109 {Pp23,Pp23,Qp21,0.0}, 110 {Pp24,Pp23,Qp21,0.0}, 111 {Pp23,Pp24,Qp21,0.0}, 112 {Pp23,Pp23,Qp22,0.0}, 113 {Pp24,Pp23,Qp22,0.0}, 114 {Pp23,Pp24,Qp22,0.0} 115 }; 116 117 static double twt6[] = {Ww32,Ww32,Ww32,Ww32,Ww32,Ww32}; 118 119 /* Rule 3 */ 120 121 static double rstw18[][4] = { 122 {Pp34,Pp35,Qp31,0.0}, 123 {Pp35,Pp34,Qp31,0.0}, 124 {Pp35,Pp35,Qp31,0.0}, 125 {Pp36,Pp37,Qp31,0.0}, 126 {Pp37,Pp36,Qp31,0.0}, 127 {Pp37,Pp37,Qp31,0.0}, 128 {Pp34,Pp35,Qp32,0.0}, 129 {Pp35,Pp34,Qp32,0.0}, 130 {Pp35,Pp35,Qp32,0.0}, 131 {Pp36,Pp37,Qp32,0.0}, 132 {Pp37,Pp36,Qp32,0.0}, 133 {Pp37,Pp37,Qp32,0.0}, 134 {Pp34,Pp35,Qp33,0.0}, 135 {Pp35,Pp34,Qp33,0.0}, 136 {Pp35,Pp35,Qp33,0.0}, 137 {Pp36,Pp37,Qp33,0.0}, 138 {Pp37,Pp36,Qp33,0.0}, 139 {Pp37,Pp37,Qp33,0.0} 140 }; 141 142 static double twt18[] = 143 { Ww3331, Ww3331, Ww3331, Ww3431, Ww3431, Ww3431, Ww3332, Ww3332, Ww3332, 144 Ww3432, Ww3432, Ww3432, Ww3333, Ww3333, Ww3333, Ww3433, Ww3433, Ww3433}; 145 146 /* Rule 4 */ 147 148 static double rstw48[][4] = { 149 {Pp41,Pp42,Qp41,0.0}, 150 {Pp42,Pp41,Qp41,0.0}, 151 {Pp42,Pp42,Qp41,0.0}, 152 {Pp43,Pp44,Qp41,0.0}, 153 {Pp44,Pp43,Qp41,0.0}, 154 {Pp44,Pp44,Qp41,0.0}, 155 {Pp45,Pp46,Qp41,0.0}, 156 {Pp46,Pp45,Qp41,0.0}, 157 {Pp45,Pp47,Qp41,0.0}, 158 {Pp46,Pp47,Qp41,0.0}, 159 {Pp47,Pp46,Qp41,0.0}, 160 {Pp47,Pp45,Qp41,0.0}, 161 {Pp41,Pp42,Qp42,0.0}, 162 {Pp42,Pp41,Qp42,0.0}, 163 {Pp42,Pp42,Qp42,0.0}, 164 {Pp43,Pp44,Qp42,0.0}, 165 {Pp44,Pp43,Qp42,0.0}, 166 {Pp44,Pp44,Qp42,0.0}, 167 {Pp45,Pp46,Qp42,0.0}, 168 {Pp46,Pp45,Qp42,0.0}, 169 {Pp45,Pp47,Qp42,0.0}, 170 {Pp46,Pp47,Qp42,0.0}, 171 {Pp47,Pp46,Qp42,0.0}, 172 {Pp47,Pp45,Qp42,0.0}, 173 {Pp41,Pp42,Qp43,0.0}, 174 {Pp42,Pp41,Qp43,0.0}, 175 {Pp42,Pp42,Qp43,0.0}, 176 {Pp43,Pp44,Qp43,0.0}, 177 {Pp44,Pp43,Qp43,0.0}, 178 {Pp44,Pp44,Qp43,0.0}, 179 {Pp45,Pp46,Qp43,0.0}, 180 {Pp46,Pp45,Qp43,0.0}, 181 {Pp45,Pp47,Qp43,0.0}, 182 {Pp46,Pp47,Qp43,0.0}, 183 {Pp47,Pp46,Qp43,0.0}, 184 {Pp47,Pp45,Qp43,0.0}, 185 {Pp41,Pp42,Qp44,0.0}, 186 {Pp42,Pp41,Qp44,0.0}, 187 {Pp42,Pp42,Qp44,0.0}, 188 {Pp43,Pp44,Qp44,0.0}, 189 {Pp44,Pp43,Qp44,0.0}, 190 {Pp44,Pp44,Qp44,0.0}, 191 {Pp45,Pp46,Qp44,0.0}, 192 {Pp46,Pp45,Qp44,0.0}, 193 {Pp45,Pp47,Qp44,0.0}, 194 {Pp46,Pp47,Qp44,0.0}, 195 {Pp47,Pp46,Qp44,0.0}, 196 {Pp47,Pp45,Qp44,0.0} 197 }; 198 199 static double twt48[] = 200 { Ww4141, Ww4141, Ww4141, Ww4241, Ww4241, Ww4241, Ww4341, Ww4341, Ww4341, 201 Ww4341, Ww4341, Ww4341, Ww4142, Ww4142, Ww4142, Ww4242, Ww4242, Ww4242, 202 Ww4342, Ww4342, Ww4342, Ww4342, Ww4342, Ww4342, Ww4143, Ww4143, Ww4143, 203 Ww4243, Ww4243, Ww4243, Ww4343, Ww4343, Ww4343, Ww4343, Ww4343, Ww4343, 204 Ww4144, Ww4144, Ww4144, Ww4244, Ww4244, Ww4244, Ww4344, Ww4344, Ww4344, 205 Ww4344, Ww4344, Ww4344}; 206 207 #ifdef __cplusplus 208 extern "C" { 209 #endif 210 211 int wdgIntPnt(int nint, DARR4 **bcord, double **wt) 212 { 213 int retval = 1; 214 215 if( nint == 6 ){*bcord = rstw6 ; *wt = twt6; } 216 else if( nint == 1 ){*bcord = rstw1 ; *wt = twt1; } 217 else if( nint == 18 ){*bcord = rstw18 ; *wt = twt18; } 218 else if( nint == 48 ){*bcord = rstw48 ; *wt = twt48; } 219 else 220 { 221 fprintf(stderr,"\n%d integration points unsupported symwdg.c; give {6,18,48}\n",nint); 222 retval = 0; 223 } 224 return retval ; 225 } 226 227 #ifdef __cplusplus 228 } 229 #endif 230