# Contents of /trunk/finley/src/Quadrature.c

Revision 1342 - (show annotations)
Thu Nov 8 23:56:58 2007 UTC (14 years, 7 months ago) by gross
File MIME type: text/plain
File size: 48957 byte(s)
```improvements on integration schemes for tets and tris: double precision (not all), higher orders
```
 1 2 /* \$Id\$ */ 3 4 /******************************************************* 5 * 6 * Copyright 2003-2007 by ACceSS MNRF 7 * Copyright 2007 by University of Queensland 8 * 9 * http://esscc.uq.edu.au 10 * Primary Business: Queensland, Australia 11 * Licensed under the Open Software License version 3.0 12 * http://www.opensource.org/licenses/osl-3.0.php 13 * 14 *******************************************************/ 15 16 /**************************************************************/ 17 18 /* Finley: */ 19 20 /**************************************************************/ 21 22 #include "Quadrature.h" 23 24 #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)] 25 #define QUADWEIGHTS(_I_) quadWeights[_I_] 26 27 /**************************************************************/ 28 29 /* get a quadrature scheme with numQuadNodes quadrature nodes for the tri */ 30 /* as a queezed scheme on a quad [0,1]^2 */ 31 32 void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) { 33 int i; 34 double Q1,Q2,a,b,c,d,e,f,g,u,v,w; 35 #define DIM 2 36 37 /* the easy cases: */ 38 39 if (numQuadNodes==1) { 40 QUADNODES(0,0)=1./3.; 41 QUADNODES(1,0)=1./3.; 42 QUADWEIGHTS(0)=1./2.; 43 } else if (numQuadNodes==3){ 44 QUADNODES(0,0)=1./2.; 45 QUADNODES(1,0)=0.; 46 QUADWEIGHTS(0)=1./6.; 47 QUADNODES(0,1)=0.; 48 QUADNODES(1,1)=1./2.; 49 QUADWEIGHTS(1)=1./6.; 50 QUADNODES(0,2)=1./2.; 51 QUADNODES(1,2)=1./2.; 52 QUADWEIGHTS(2)=1./6.; 53 } else if (numQuadNodes==4){ 54 QUADNODES(0,0)=1./3.; 55 QUADNODES(1,0)=1./3.; 56 QUADWEIGHTS(0)=-27./96.; 57 QUADNODES(0,1)=0.2; 58 QUADNODES(1,1)=0.2; 59 QUADWEIGHTS(1)=25./96.; 60 QUADNODES(0,2)=0.6; 61 QUADNODES(1,2)=0.2; 62 QUADWEIGHTS(2)=25./96.; 63 QUADNODES(0,3)=0.2; 64 QUADNODES(1,3)=0.6; 65 QUADWEIGHTS(3)=25./96.; 66 } else if (numQuadNodes==6){ 67 QUADWEIGHTS(0) = 0.109951743655322/2.; 68 QUADWEIGHTS(1) = 0.109951743655322/2.; 69 QUADWEIGHTS(2) = 0.109951743655322/2.; 70 QUADWEIGHTS(3) = 0.223381589678011/2.; 71 QUADWEIGHTS(4) = 0.223381589678011/2.; 72 QUADWEIGHTS(5) = 0.223381589678011/2.; 73 74 QUADNODES(0,0) = 0.816847572980459; 75 QUADNODES(0,1) = 0.091576213509771; 76 QUADNODES(0,2) = 0.091576213509771; 77 QUADNODES(0,3) = 0.108103018168070; 78 QUADNODES(0,4) = 0.445948490915965; 79 QUADNODES(0,5) = 0.445948490915965; 80 81 QUADNODES(1,0) = 0.091576213509771; 82 QUADNODES(1,1) = 0.816847572980459; 83 QUADNODES(1,2) = 0.091576213509771; 84 QUADNODES(1,3) = 0.445948490915965; 85 QUADNODES(1,4) = 0.108103018168070; 86 QUADNODES(1,5) = 0.445948490915965; 87 88 } else if (numQuadNodes==7){ 89 QUADNODES(0,0) = 0.33333333333333333; 90 QUADNODES(0,1) = 0.7974269853530872; 91 QUADNODES(0,2) = 0.10128650732345633; 92 QUADNODES(0,3) = 0.10128650732345633; 93 QUADNODES(0,4) = 0.059715871789769809; 94 QUADNODES(0,5) = 0.47014206410511505; 95 QUADNODES(0,6) = 0.47014206410511505; 96 97 QUADNODES(1,0) = 0.33333333333333333; 98 QUADNODES(1,1) = 0.10128650732345633; 99 QUADNODES(1,2) = 0.7974269853530872; 100 QUADNODES(1,3) = 0.10128650732345633; 101 QUADNODES(1,4) = 0.47014206410511505; 102 QUADNODES(1,5) = 0.059715871789769809; 103 QUADNODES(1,6) = 0.47014206410511505; 104 105 QUADWEIGHTS(0) = 0.225/2.; 106 QUADWEIGHTS(1) = 0.12593918054482717/2.; 107 QUADWEIGHTS(2) = 0.12593918054482717/2.; 108 QUADWEIGHTS(3) = 0.12593918054482717/2.; 109 QUADWEIGHTS(4) = 0.13239415278850616/2.; 110 QUADWEIGHTS(5) = 0.13239415278850616/2.; 111 QUADWEIGHTS(6) = 0.13239415278850616/2.; 112 113 } else if (numQuadNodes==12){ 114 a = 0.873821971016996; 115 b = 0.063089014491502; 116 c = 0.501426509658179; 117 d = 0.249286745170910; 118 e = 0.636502499121399; 119 f = 0.310352451033785; 120 g = 0.053145049844816; 121 122 u = 0.050844906370207/2.; 123 v = 0.116786275726379/2.; 124 w = 0.082851075618374/2.; 125 126 QUADNODES(0,0) = a; 127 QUADNODES(0,1) = b; 128 QUADNODES(0,2) = b; 129 QUADNODES(0,3) = c; 130 QUADNODES(0,4) = d; 131 QUADNODES(0,5) = d; 132 QUADNODES(0,6) = e; 133 QUADNODES(0,7) = e; 134 QUADNODES(0,8) = f; 135 QUADNODES(0,9) = f; 136 QUADNODES(0,10) = g; 137 QUADNODES(0,11) = g; 138 139 QUADNODES(1,0) = b; 140 QUADNODES(1,1) = a; 141 QUADNODES(1,2) = b; 142 QUADNODES(1,3) = d; 143 QUADNODES(1,4) = c; 144 QUADNODES(1,5) = d; 145 QUADNODES(1,6) = f; 146 QUADNODES(1,7) = g; 147 QUADNODES(1,8) = e; 148 QUADNODES(1,9) = g; 149 QUADNODES(1,10) = e; 150 QUADNODES(1,11) = f; 151 152 QUADWEIGHTS(0)= u; 153 QUADWEIGHTS(1)= u; 154 QUADWEIGHTS(2)= u; 155 QUADWEIGHTS(3)= v; 156 QUADWEIGHTS(4)= v; 157 QUADWEIGHTS(5)= v; 158 QUADWEIGHTS(6)= w; 159 QUADWEIGHTS(7)= w; 160 QUADWEIGHTS(8)= w; 161 QUADWEIGHTS(9)= w; 162 QUADWEIGHTS(10)= w; 163 QUADWEIGHTS(11)= w; 164 165 } else if (numQuadNodes==13){ 166 QUADWEIGHTS(0) =-0.149570044467670/2.; 167 QUADWEIGHTS(1) = 0.175615257433204/2.; 168 QUADWEIGHTS(2) = 0.175615257433204/2.; 169 QUADWEIGHTS(3) = 0.175615257433204/2.; 170 QUADWEIGHTS(4) = 0.053347235608839/2.; 171 QUADWEIGHTS(5) = 0.053347235608839/2.; 172 QUADWEIGHTS(6) = 0.053347235608839/2.; 173 QUADWEIGHTS(7) = 0.077113760890257/2.; 174 QUADWEIGHTS(8) = 0.077113760890257/2.; 175 QUADWEIGHTS(9) = 0.077113760890257/2.; 176 QUADWEIGHTS(10) = 0.077113760890257/2.; 177 QUADWEIGHTS(11) = 0.077113760890257/2.; 178 QUADWEIGHTS(12) = 0.077113760890257/2.; 179 180 QUADNODES(0,0) = 0.3333333333333333; 181 QUADNODES(0,1) = 0.479308067841923; 182 QUADNODES(0,2) = 0.260345966079038; 183 QUADNODES(0,3) = 0.260345966079038; 184 QUADNODES(0,4) = 0.869739794195568; 185 QUADNODES(0,5) = 0.065130102902216; 186 QUADNODES(0,6) = 0.065130102902216; 187 QUADNODES(0,7) = 0.638444188569809; 188 QUADNODES(0,8) = 0.638444188569809; 189 QUADNODES(0,9) = 0.048690315425316; 190 QUADNODES(0,10) = 0.048690315425316; 191 QUADNODES(0,11) = 0.312865496004875; 192 QUADNODES(0,12) = 0.312865496004875; 193 194 QUADNODES(1,0) = 0.3333333333333333; 195 QUADNODES(1,1) = 0.260345966079038; 196 QUADNODES(1,2) = 0.479308067841923; 197 QUADNODES(1,3) = 0.260345966079038; 198 QUADNODES(1,4) = 0.065130102902216; 199 QUADNODES(1,5) = 0.869739794195568; 200 QUADNODES(1,6) = 0.065130102902216; 201 QUADNODES(1,7) = 0.048690315425316; 202 QUADNODES(1,8) = 0.312865496004875; 203 QUADNODES(1,9) = 0.638444188569809; 204 QUADNODES(1,10) = 0.312865496004875; 205 QUADNODES(1,11) = 0.638444188569809; 206 QUADNODES(1,12) = 0.048690315425316; 207 208 } else if (numQuadNodes==16){ 209 QUADWEIGHTS(0) = 0.07215780; 210 QUADWEIGHTS(1) = 0.04754582; 211 QUADWEIGHTS(2) = 0.04754582; 212 QUADWEIGHTS(3) = 0.04754582; 213 QUADWEIGHTS(4) = 0.01622925; 214 QUADWEIGHTS(5) = 0.01622925; 215 QUADWEIGHTS(6) = 0.01622925; 216 QUADWEIGHTS(7) = 0.05160869; 217 QUADWEIGHTS(8) = 0.05160869; 218 QUADWEIGHTS(9) = 0.05160869; 219 QUADWEIGHTS(10) = 0.01361516; 220 QUADWEIGHTS(11) = 0.01361516; 221 QUADWEIGHTS(12) = 0.01361516; 222 QUADWEIGHTS(13) = 0.01361516; 223 QUADWEIGHTS(14) = 0.01361516; 224 QUADWEIGHTS(15) = 0.01361516; 225 226 QUADNODES(0,0) = 0.3333333; 227 QUADNODES(0,1) = 0.08141482; 228 QUADNODES(0,2) = 0.4592926; 229 QUADNODES(0,3) = 0.4592926; 230 QUADNODES(0,4) = 0.8989055; 231 QUADNODES(0,5) = 0.05054723; 232 QUADNODES(0,6) = 0.05054723; 233 QUADNODES(0,7) = 0.6588614; 234 QUADNODES(0,8) = 0.1705693; 235 QUADNODES(0,9) = 0.1705693; 236 QUADNODES(0,10) = 0.008394777; 237 QUADNODES(0,11) = 0.008394777; 238 QUADNODES(0,12) = 0.7284924; 239 QUADNODES(0,13) = 0.7284924; 240 QUADNODES(0,14) = 0.2631128; 241 QUADNODES(0,15) = 0.2631128; 242 243 QUADNODES(1,0) = 0.3333333; 244 QUADNODES(1,1) = 0.4592926; 245 QUADNODES(1,2) = 0.08141482; 246 QUADNODES(1,3) = 0.4592926; 247 QUADNODES(1,4) = 0.05054723; 248 QUADNODES(1,5) = 0.8989055; 249 QUADNODES(1,6) = 0.05054723; 250 QUADNODES(1,7) = 0.1705693; 251 QUADNODES(1,8) = 0.6588614; 252 QUADNODES(1,9) = 0.1705693; 253 QUADNODES(1,10) = 0.7284924; 254 QUADNODES(1,11) = 0.2631128; 255 QUADNODES(1,12) = 0.008394777; 256 QUADNODES(1,13) = 0.2631128; 257 QUADNODES(1,14) = 0.008394777; 258 QUADNODES(1,15) = 0.7284924; 259 260 } else if (numQuadNodes==19){ 261 QUADWEIGHTS(0) = 0.04856790; 262 QUADWEIGHTS(1) = 0.01566735; 263 QUADWEIGHTS(2) = 0.01566735; 264 QUADWEIGHTS(3) = 0.01566735; 265 QUADWEIGHTS(4) = 0.03891377; 266 QUADWEIGHTS(5) = 0.03891377; 267 QUADWEIGHTS(6) = 0.03891377; 268 QUADWEIGHTS(7) = 0.03982387; 269 QUADWEIGHTS(8) = 0.03982387; 270 QUADWEIGHTS(9) = 0.03982387; 271 QUADWEIGHTS(10) = 0.01278884; 272 QUADWEIGHTS(11) = 0.01278884; 273 QUADWEIGHTS(12) = 0.01278884; 274 QUADWEIGHTS(13) = 0.02164177; 275 QUADWEIGHTS(14) = 0.02164177; 276 QUADWEIGHTS(15) = 0.02164177; 277 QUADWEIGHTS(16) = 0.02164177; 278 QUADWEIGHTS(17) = 0.02164177; 279 QUADWEIGHTS(18) = 0.02164177; 280 281 QUADNODES(0,0) = 0.3333333; 282 QUADNODES(0,1) = 0.02063496; 283 QUADNODES(0,2) = 0.4896825; 284 QUADNODES(0,3) = 0.4896825; 285 QUADNODES(0,4) = 0.1258208; 286 QUADNODES(0,5) = 0.4370896; 287 QUADNODES(0,6) = 0.4370896; 288 QUADNODES(0,7) = 0.6235929; 289 QUADNODES(0,8) = 0.1882035; 290 QUADNODES(0,9) = 0.1882035; 291 QUADNODES(0,10) = 0.9105410; 292 QUADNODES(0,11) = 0.04472951; 293 QUADNODES(0,12) = 0.04472951; 294 QUADNODES(0,13) = 0.03683841; 295 QUADNODES(0,14) = 0.03683841; 296 QUADNODES(0,15) = 0.7411986; 297 QUADNODES(0,16) = 0.7411986; 298 QUADNODES(0,17) = 0.2219630; 299 QUADNODES(0,18) = 0.2219630; 300 301 QUADNODES(1,0) = 0.3333333; 302 QUADNODES(1,1) = 0.4896825; 303 QUADNODES(1,2) = 0.02063496; 304 QUADNODES(1,3) = 0.4896825; 305 QUADNODES(1,4) = 0.4370896; 306 QUADNODES(1,5) = 0.1258208; 307 QUADNODES(1,6) = 0.4370896; 308 QUADNODES(1,7) = 0.1882035; 309 QUADNODES(1,8) = 0.6235929; 310 QUADNODES(1,9) = 0.1882035; 311 QUADNODES(1,10) = 0.04472951; 312 QUADNODES(1,11) = 0.9105410; 313 QUADNODES(1,12) = 0.04472951; 314 QUADNODES(1,13) = 0.7411986; 315 QUADNODES(1,14) = 0.2219630; 316 QUADNODES(1,15) = 0.03683841; 317 QUADNODES(1,16) = 0.2219630; 318 QUADNODES(1,17) = 0.03683841; 319 QUADNODES(1,18) = 0.7411986; 320 } else { 321 322 /* get scheme on [0.1]^2 */ 323 Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights); 324 if (! Finley_noError()) return; 325 326 /* squeeze it: */ 327 328 for (i=0;i 2*MAX_numQuadNodesLine-1) { 1302 sprintf(error_msg,"Finley_Quad_getNumNodesPoint: requested integration order %d on line is too large (>%d).", 1303 order,2*MAX_numQuadNodesLine-1); 1304 Finley_setError(VALUE_ERROR,error_msg); 1305 return -1; 1306 } else { 1307 Finley_resetError(); 1308 return order/2+1; 1309 } 1310 } 1311 } 1312 1313 int Finley_Quad_getNumNodesTri(int order) { 1314 int numQuadNodesLine; 1315 if (order<=1) { 1316 return 1; 1317 } else if (order<=2){ 1318 return 3; 1319 } else if (order<=3){ 1320 return 4; 1321 } else if (order<=4){ 1322 return 6; 1323 } else if (order<=5){ 1324 return 7; 1325 } else if (order<=6){ 1326 return 12; 1327 } else if (order<=7){ 1328 return 13; 1329 } else if (order<=8){ 1330 return 16; 1331 } else if (order<=9){ 1332 return 19; 1333 } else { 1334 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1); 1335 if (Finley_noError()) { 1336 return numQuadNodesLine*numQuadNodesLine; 1337 } else { 1338 return -1; 1339 } 1340 } 1341 } 1342 1343 int Finley_Quad_getNumNodesRec(int order) { 1344 int numQuadNodesLine; 1345 numQuadNodesLine=Finley_Quad_getNumNodesLine(order); 1346 if (Finley_noError()) { 1347 return numQuadNodesLine*numQuadNodesLine; 1348 } else { 1349 return -1; 1350 } 1351 } 1352 1353 int Finley_Quad_getNumNodesTet(int order) { 1354 int numQuadNodesLine; 1355 if (order<=1) { 1356 return 1; 1357 } else if (order<=2){ 1358 return 4; 1359 } else if (order<=3){ 1360 return 5; 1361 } else if (order<=4){ 1362 return 11; 1363 } else if (order<=5){ 1364 return 15; 1365 } else if (order<=6){ 1366 return 24; 1367 } else if (order<=7){ 1368 return 31; 1369 } else if (order<=8){ 1370 return 45; 1371 } else { 1372 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2); 1373 if (Finley_noError()) { 1374 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine; 1375 } else { 1376 return -1; 1377 } 1378 } 1379 } 1380 1381 int Finley_Quad_getNumNodesHex(int order) { 1382 int numQuadNodesLine; 1383 numQuadNodesLine=Finley_Quad_getNumNodesLine(order); 1384 if (Finley_noError()) { 1385 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine; 1386 } else { 1387 return -1; 1388 } 1389 }

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

 ViewVC Help Powered by ViewVC 1.1.26