/[escript]/trunk/finley/src/Quadrature.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2034 - (show annotations)
Thu Nov 13 03:03:42 2008 UTC (10 years, 10 months ago) by artak
File MIME type: text/plain
File size: 48928 byte(s)
variable <set> is initilized to get rid of compiler warnning errors. finley/src/Quadrature.c
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Finley: */
18
19 /**************************************************************/
20
21 #include "Quadrature.h"
22
23 #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
24 #define QUADWEIGHTS(_I_) quadWeights[_I_]
25
26 /**************************************************************/
27
28 /* get a quadrature scheme with numQuadNodes quadrature nodes for the tri */
29 /* as a queezed scheme on a quad [0,1]^2 */
30
31 void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
32 int i;
33 double Q1,Q2,a,b,c,d,e,f,g,u,v,w;
34 #define DIM 2
35
36 /* the easy cases: */
37
38 if (numQuadNodes==1) {
39 QUADNODES(0,0)=1./3.;
40 QUADNODES(1,0)=1./3.;
41 QUADWEIGHTS(0)=1./2.;
42 } else if (numQuadNodes==3){
43 QUADNODES(0,0)=1./2.;
44 QUADNODES(1,0)=0.;
45 QUADWEIGHTS(0)=1./6.;
46 QUADNODES(0,1)=0.;
47 QUADNODES(1,1)=1./2.;
48 QUADWEIGHTS(1)=1./6.;
49 QUADNODES(0,2)=1./2.;
50 QUADNODES(1,2)=1./2.;
51 QUADWEIGHTS(2)=1./6.;
52 } else if (numQuadNodes==4){
53 QUADNODES(0,0)=1./3.;
54 QUADNODES(1,0)=1./3.;
55 QUADWEIGHTS(0)=-27./96.;
56 QUADNODES(0,1)=0.2;
57 QUADNODES(1,1)=0.2;
58 QUADWEIGHTS(1)=25./96.;
59 QUADNODES(0,2)=0.6;
60 QUADNODES(1,2)=0.2;
61 QUADWEIGHTS(2)=25./96.;
62 QUADNODES(0,3)=0.2;
63 QUADNODES(1,3)=0.6;
64 QUADWEIGHTS(3)=25./96.;
65 } else if (numQuadNodes==6){
66 QUADWEIGHTS(0) = 0.109951743655322/2.;
67 QUADWEIGHTS(1) = 0.109951743655322/2.;
68 QUADWEIGHTS(2) = 0.109951743655322/2.;
69 QUADWEIGHTS(3) = 0.223381589678011/2.;
70 QUADWEIGHTS(4) = 0.223381589678011/2.;
71 QUADWEIGHTS(5) = 0.223381589678011/2.;
72
73 QUADNODES(0,0) = 0.816847572980459;
74 QUADNODES(0,1) = 0.091576213509771;
75 QUADNODES(0,2) = 0.091576213509771;
76 QUADNODES(0,3) = 0.108103018168070;
77 QUADNODES(0,4) = 0.445948490915965;
78 QUADNODES(0,5) = 0.445948490915965;
79
80 QUADNODES(1,0) = 0.091576213509771;
81 QUADNODES(1,1) = 0.816847572980459;
82 QUADNODES(1,2) = 0.091576213509771;
83 QUADNODES(1,3) = 0.445948490915965;
84 QUADNODES(1,4) = 0.108103018168070;
85 QUADNODES(1,5) = 0.445948490915965;
86
87 } else if (numQuadNodes==7){
88 QUADNODES(0,0) = 0.33333333333333333;
89 QUADNODES(0,1) = 0.7974269853530872;
90 QUADNODES(0,2) = 0.10128650732345633;
91 QUADNODES(0,3) = 0.10128650732345633;
92 QUADNODES(0,4) = 0.059715871789769809;
93 QUADNODES(0,5) = 0.47014206410511505;
94 QUADNODES(0,6) = 0.47014206410511505;
95
96 QUADNODES(1,0) = 0.33333333333333333;
97 QUADNODES(1,1) = 0.10128650732345633;
98 QUADNODES(1,2) = 0.7974269853530872;
99 QUADNODES(1,3) = 0.10128650732345633;
100 QUADNODES(1,4) = 0.47014206410511505;
101 QUADNODES(1,5) = 0.059715871789769809;
102 QUADNODES(1,6) = 0.47014206410511505;
103
104 QUADWEIGHTS(0) = 0.225/2.;
105 QUADWEIGHTS(1) = 0.12593918054482717/2.;
106 QUADWEIGHTS(2) = 0.12593918054482717/2.;
107 QUADWEIGHTS(3) = 0.12593918054482717/2.;
108 QUADWEIGHTS(4) = 0.13239415278850616/2.;
109 QUADWEIGHTS(5) = 0.13239415278850616/2.;
110 QUADWEIGHTS(6) = 0.13239415278850616/2.;
111
112 } else if (numQuadNodes==12){
113 a = 0.873821971016996;
114 b = 0.063089014491502;
115 c = 0.501426509658179;
116 d = 0.249286745170910;
117 e = 0.636502499121399;
118 f = 0.310352451033785;
119 g = 0.053145049844816;
120
121 u = 0.050844906370207/2.;
122 v = 0.116786275726379/2.;
123 w = 0.082851075618374/2.;
124
125 QUADNODES(0,0) = a;
126 QUADNODES(0,1) = b;
127 QUADNODES(0,2) = b;
128 QUADNODES(0,3) = c;
129 QUADNODES(0,4) = d;
130 QUADNODES(0,5) = d;
131 QUADNODES(0,6) = e;
132 QUADNODES(0,7) = e;
133 QUADNODES(0,8) = f;
134 QUADNODES(0,9) = f;
135 QUADNODES(0,10) = g;
136 QUADNODES(0,11) = g;
137
138 QUADNODES(1,0) = b;
139 QUADNODES(1,1) = a;
140 QUADNODES(1,2) = b;
141 QUADNODES(1,3) = d;
142 QUADNODES(1,4) = c;
143 QUADNODES(1,5) = d;
144 QUADNODES(1,6) = f;
145 QUADNODES(1,7) = g;
146 QUADNODES(1,8) = e;
147 QUADNODES(1,9) = g;
148 QUADNODES(1,10) = e;
149 QUADNODES(1,11) = f;
150
151 QUADWEIGHTS(0)= u;
152 QUADWEIGHTS(1)= u;
153 QUADWEIGHTS(2)= u;
154 QUADWEIGHTS(3)= v;
155 QUADWEIGHTS(4)= v;
156 QUADWEIGHTS(5)= v;
157 QUADWEIGHTS(6)= w;
158 QUADWEIGHTS(7)= w;
159 QUADWEIGHTS(8)= w;
160 QUADWEIGHTS(9)= w;
161 QUADWEIGHTS(10)= w;
162 QUADWEIGHTS(11)= w;
163
164 } else if (numQuadNodes==13){
165 QUADWEIGHTS(0) =-0.149570044467670/2.;
166 QUADWEIGHTS(1) = 0.175615257433204/2.;
167 QUADWEIGHTS(2) = 0.175615257433204/2.;
168 QUADWEIGHTS(3) = 0.175615257433204/2.;
169 QUADWEIGHTS(4) = 0.053347235608839/2.;
170 QUADWEIGHTS(5) = 0.053347235608839/2.;
171 QUADWEIGHTS(6) = 0.053347235608839/2.;
172 QUADWEIGHTS(7) = 0.077113760890257/2.;
173 QUADWEIGHTS(8) = 0.077113760890257/2.;
174 QUADWEIGHTS(9) = 0.077113760890257/2.;
175 QUADWEIGHTS(10) = 0.077113760890257/2.;
176 QUADWEIGHTS(11) = 0.077113760890257/2.;
177 QUADWEIGHTS(12) = 0.077113760890257/2.;
178
179 QUADNODES(0,0) = 0.3333333333333333;
180 QUADNODES(0,1) = 0.479308067841923;
181 QUADNODES(0,2) = 0.260345966079038;
182 QUADNODES(0,3) = 0.260345966079038;
183 QUADNODES(0,4) = 0.869739794195568;
184 QUADNODES(0,5) = 0.065130102902216;
185 QUADNODES(0,6) = 0.065130102902216;
186 QUADNODES(0,7) = 0.638444188569809;
187 QUADNODES(0,8) = 0.638444188569809;
188 QUADNODES(0,9) = 0.048690315425316;
189 QUADNODES(0,10) = 0.048690315425316;
190 QUADNODES(0,11) = 0.312865496004875;
191 QUADNODES(0,12) = 0.312865496004875;
192
193 QUADNODES(1,0) = 0.3333333333333333;
194 QUADNODES(1,1) = 0.260345966079038;
195 QUADNODES(1,2) = 0.479308067841923;
196 QUADNODES(1,3) = 0.260345966079038;
197 QUADNODES(1,4) = 0.065130102902216;
198 QUADNODES(1,5) = 0.869739794195568;
199 QUADNODES(1,6) = 0.065130102902216;
200 QUADNODES(1,7) = 0.048690315425316;
201 QUADNODES(1,8) = 0.312865496004875;
202 QUADNODES(1,9) = 0.638444188569809;
203 QUADNODES(1,10) = 0.312865496004875;
204 QUADNODES(1,11) = 0.638444188569809;
205 QUADNODES(1,12) = 0.048690315425316;
206
207 } else if (numQuadNodes==16){
208 QUADWEIGHTS(0) = 0.07215780;
209 QUADWEIGHTS(1) = 0.04754582;
210 QUADWEIGHTS(2) = 0.04754582;
211 QUADWEIGHTS(3) = 0.04754582;
212 QUADWEIGHTS(4) = 0.01622925;
213 QUADWEIGHTS(5) = 0.01622925;
214 QUADWEIGHTS(6) = 0.01622925;
215 QUADWEIGHTS(7) = 0.05160869;
216 QUADWEIGHTS(8) = 0.05160869;
217 QUADWEIGHTS(9) = 0.05160869;
218 QUADWEIGHTS(10) = 0.01361516;
219 QUADWEIGHTS(11) = 0.01361516;
220 QUADWEIGHTS(12) = 0.01361516;
221 QUADWEIGHTS(13) = 0.01361516;
222 QUADWEIGHTS(14) = 0.01361516;
223 QUADWEIGHTS(15) = 0.01361516;
224
225 QUADNODES(0,0) = 0.3333333;
226 QUADNODES(0,1) = 0.08141482;
227 QUADNODES(0,2) = 0.4592926;
228 QUADNODES(0,3) = 0.4592926;
229 QUADNODES(0,4) = 0.8989055;
230 QUADNODES(0,5) = 0.05054723;
231 QUADNODES(0,6) = 0.05054723;
232 QUADNODES(0,7) = 0.6588614;
233 QUADNODES(0,8) = 0.1705693;
234 QUADNODES(0,9) = 0.1705693;
235 QUADNODES(0,10) = 0.008394777;
236 QUADNODES(0,11) = 0.008394777;
237 QUADNODES(0,12) = 0.7284924;
238 QUADNODES(0,13) = 0.7284924;
239 QUADNODES(0,14) = 0.2631128;
240 QUADNODES(0,15) = 0.2631128;
241
242 QUADNODES(1,0) = 0.3333333;
243 QUADNODES(1,1) = 0.4592926;
244 QUADNODES(1,2) = 0.08141482;
245 QUADNODES(1,3) = 0.4592926;
246 QUADNODES(1,4) = 0.05054723;
247 QUADNODES(1,5) = 0.8989055;
248 QUADNODES(1,6) = 0.05054723;
249 QUADNODES(1,7) = 0.1705693;
250 QUADNODES(1,8) = 0.6588614;
251 QUADNODES(1,9) = 0.1705693;
252 QUADNODES(1,10) = 0.7284924;
253 QUADNODES(1,11) = 0.2631128;
254 QUADNODES(1,12) = 0.008394777;
255 QUADNODES(1,13) = 0.2631128;
256 QUADNODES(1,14) = 0.008394777;
257 QUADNODES(1,15) = 0.7284924;
258
259 } else if (numQuadNodes==19){
260 QUADWEIGHTS(0) = 0.04856790;
261 QUADWEIGHTS(1) = 0.01566735;
262 QUADWEIGHTS(2) = 0.01566735;
263 QUADWEIGHTS(3) = 0.01566735;
264 QUADWEIGHTS(4) = 0.03891377;
265 QUADWEIGHTS(5) = 0.03891377;
266 QUADWEIGHTS(6) = 0.03891377;
267 QUADWEIGHTS(7) = 0.03982387;
268 QUADWEIGHTS(8) = 0.03982387;
269 QUADWEIGHTS(9) = 0.03982387;
270 QUADWEIGHTS(10) = 0.01278884;
271 QUADWEIGHTS(11) = 0.01278884;
272 QUADWEIGHTS(12) = 0.01278884;
273 QUADWEIGHTS(13) = 0.02164177;
274 QUADWEIGHTS(14) = 0.02164177;
275 QUADWEIGHTS(15) = 0.02164177;
276 QUADWEIGHTS(16) = 0.02164177;
277 QUADWEIGHTS(17) = 0.02164177;
278 QUADWEIGHTS(18) = 0.02164177;
279
280 QUADNODES(0,0) = 0.3333333;
281 QUADNODES(0,1) = 0.02063496;
282 QUADNODES(0,2) = 0.4896825;
283 QUADNODES(0,3) = 0.4896825;
284 QUADNODES(0,4) = 0.1258208;
285 QUADNODES(0,5) = 0.4370896;
286 QUADNODES(0,6) = 0.4370896;
287 QUADNODES(0,7) = 0.6235929;
288 QUADNODES(0,8) = 0.1882035;
289 QUADNODES(0,9) = 0.1882035;
290 QUADNODES(0,10) = 0.9105410;
291 QUADNODES(0,11) = 0.04472951;
292 QUADNODES(0,12) = 0.04472951;
293 QUADNODES(0,13) = 0.03683841;
294 QUADNODES(0,14) = 0.03683841;
295 QUADNODES(0,15) = 0.7411986;
296 QUADNODES(0,16) = 0.7411986;
297 QUADNODES(0,17) = 0.2219630;
298 QUADNODES(0,18) = 0.2219630;
299
300 QUADNODES(1,0) = 0.3333333;
301 QUADNODES(1,1) = 0.4896825;
302 QUADNODES(1,2) = 0.02063496;
303 QUADNODES(1,3) = 0.4896825;
304 QUADNODES(1,4) = 0.4370896;
305 QUADNODES(1,5) = 0.1258208;
306 QUADNODES(1,6) = 0.4370896;
307 QUADNODES(1,7) = 0.1882035;
308 QUADNODES(1,8) = 0.6235929;
309 QUADNODES(1,9) = 0.1882035;
310 QUADNODES(1,10) = 0.04472951;
311 QUADNODES(1,11) = 0.9105410;
312 QUADNODES(1,12) = 0.04472951;
313 QUADNODES(1,13) = 0.7411986;
314 QUADNODES(1,14) = 0.2219630;
315 QUADNODES(1,15) = 0.03683841;
316 QUADNODES(1,16) = 0.2219630;
317 QUADNODES(1,17) = 0.03683841;
318 QUADNODES(1,18) = 0.7411986;
319 } else {
320
321 /* get scheme on [0.1]^2 */
322 Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
323 if (! Finley_noError()) return;
324
325 /* squeeze it: */
326
327 for (i=0;i<numQuadNodes;i++) {
328 Q1=QUADNODES(0,i);
329 Q2=QUADNODES(1,i);
330 QUADWEIGHTS(i)=QUADWEIGHTS(i)*(1.-(1./2.)*(Q1+Q2));
331 QUADNODES(0,i)=Q1*(1.-(1./2.)*Q2);
332 QUADNODES(1,i)=Q2*(1.-(1./2.)*Q1);
333 }
334 }
335 #undef DIM
336
337
338 }
339
340 /**************************************************************/
341
342 /* get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
343 /* as a queezed scheme on a hex [0,1]^3 */
344
345 void Finley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
346 int i;
347 double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;
348 double a,b,c,d,e,f,g,h;
349 double alpha=0.58541019662496852;
350 double beta =0.1381966011250105;
351 #define DIM 3
352
353 /* the easy cases: */
354 if (numQuadNodes==1) {
355 QUADNODES(0,0)=0.25;
356 QUADNODES(1,0)=0.25;
357 QUADNODES(2,0)=0.25;
358 QUADWEIGHTS(0)=1./6.;
359 } else if (numQuadNodes==4){
360 QUADNODES(0,0)=beta;
361 QUADNODES(1,0)=beta;
362 QUADNODES(2,0)=beta;
363 QUADWEIGHTS(0)=1./24.;
364 QUADNODES(0,1)=alpha;
365 QUADNODES(1,1)=beta;
366 QUADNODES(2,1)=beta;
367 QUADWEIGHTS(1)=1./24.;
368 QUADNODES(0,2)=beta;
369 QUADNODES(1,2)=alpha;
370 QUADNODES(2,2)=beta;
371 QUADWEIGHTS(2)=1./24.;
372 QUADNODES(0,3)=beta;
373 QUADNODES(1,3)=beta;
374 QUADNODES(2,3)=alpha;
375 QUADWEIGHTS(3)=1./24.;
376 } else if (numQuadNodes==5){
377 QUADNODES(0,0)=1./4.;
378 QUADNODES(1,0)=1./4.;
379 QUADNODES(2,0)=1./4.;
380 QUADWEIGHTS(0)=-2./15.;
381 QUADNODES(0,1)=1./6.;
382 QUADNODES(1,1)=1./6.;
383 QUADNODES(2,1)=1./6.;
384 QUADWEIGHTS(1)=3./40.;
385 QUADNODES(0,2)=1./2.;
386 QUADNODES(1,2)=1./6.;
387 QUADNODES(2,2)=1./6.;
388 QUADWEIGHTS(2)=3./40.;
389 QUADNODES(0,3)=1./6.;
390 QUADNODES(1,3)=1./2.;
391 QUADNODES(2,3)=1./6.;
392 QUADWEIGHTS(3)=3./40.;
393 QUADNODES(0,4)=1./6.;
394 QUADNODES(1,4)=1./6.;
395 QUADNODES(2,4)=1./2.;
396 QUADWEIGHTS(4)=3./40.;
397
398 } else if (numQuadNodes==11){
399
400 a = 0.25;
401 b = 11.0/14.0;
402 c = 1.0/14.0;
403 d = 0.25 * (1.0 + sqrt ( 5.0 / 14.0 ) );
404 e = 0.25 * (1.0 - sqrt ( 5.0 / 14.0 ) );
405 f = -74.0 / 5625.0;
406 g = 343.0 / 45000.0;
407 h = 56.0 / 2250.0;
408
409 QUADWEIGHTS(401-401) = f;
410 QUADWEIGHTS(402-401) = g;
411 QUADWEIGHTS(403-401) = g;
412 QUADWEIGHTS(404-401) = g;
413 QUADWEIGHTS(405-401) = g;
414 QUADWEIGHTS(406-401) = h;
415 QUADWEIGHTS(407-401) = h;
416 QUADWEIGHTS(408-401) = h;
417 QUADWEIGHTS(409-401) = h;
418 QUADWEIGHTS(410-401) = h;
419 QUADWEIGHTS(411-401) = h;
420
421 QUADNODES(0,401-401) = a;
422 QUADNODES(0,402-401) = b;
423 QUADNODES(0,403-401) = c;
424 QUADNODES(0,404-401) = c;
425 QUADNODES(0,405-401) = c;
426 QUADNODES(0,406-401) = d;
427 QUADNODES(0,407-401) = d;
428 QUADNODES(0,408-401) = d;
429 QUADNODES(0,409-401) = e;
430 QUADNODES(0,410-401) = e;
431 QUADNODES(0,411-401) = e;
432
433 QUADNODES(1,401-401) = a;
434 QUADNODES(1,402-401) = c;
435 QUADNODES(1,403-401) = b;
436 QUADNODES(1,404-401) = c;
437 QUADNODES(1,405-401) = c;
438 QUADNODES(1,406-401) = d;
439 QUADNODES(1,407-401) = e;
440 QUADNODES(1,408-401) = e;
441 QUADNODES(1,409-401) = d;
442 QUADNODES(1,410-401) = d;
443 QUADNODES(1,411-401) = e;
444
445 QUADNODES(2,401-401) = a;
446 QUADNODES(2,402-401) = c;
447 QUADNODES(2,403-401) = c;
448 QUADNODES(2,404-401) = b;
449 QUADNODES(2,405-401) = c;
450 QUADNODES(2,406-401) = e;
451 QUADNODES(2,407-401) = d;
452 QUADNODES(2,408-401) = e;
453 QUADNODES(2,409-401) = d;
454 QUADNODES(2,410-401) = e;
455 QUADNODES(2,411-401) = d;
456
457 } else if (numQuadNodes==15){
458 QUADWEIGHTS(412-412) = 0.019753086419753086;
459 QUADWEIGHTS(413-412) = 0.01198951396316977;
460 QUADWEIGHTS(414-412) = 0.01198951396316977;
461 QUADWEIGHTS(415-412) = 0.01198951396316977;
462 QUADWEIGHTS(416-412) = 0.01198951396316977;
463 QUADWEIGHTS(417-412) = 0.011511367871045397;
464 QUADWEIGHTS(418-412) = 0.011511367871045397;
465 QUADWEIGHTS(419-412) = 0.011511367871045397;
466 QUADWEIGHTS(420-412) = 0.011511367871045397;
467 QUADWEIGHTS(421-412) = 0.0088183421516754845;
468 QUADWEIGHTS(422-412) = 0.0088183421516754845;
469 QUADWEIGHTS(423-412) = 0.0088183421516754845;
470 QUADWEIGHTS(424-412) = 0.0088183421516754845;
471 QUADWEIGHTS(425-412) = 0.0088183421516754845;
472 QUADWEIGHTS(426-412) = 0.0088183421516754845;
473
474 QUADNODES(0,412-412) = 0.2500000;
475 QUADNODES(0,413-412) = 0.091971078052723032;
476 QUADNODES(0,414-412) = 0.72408676584183096;
477 QUADNODES(0,415-412) = 0.091971078052723032;
478 QUADNODES(0,416-412) = 0.091971078052723032;
479 QUADNODES(0,417-412) = 0.31979362782962989;
480 QUADNODES(0,418-412) = 0.040619116511110234;
481 QUADNODES(0,419-412) = 0.31979362782962989;
482 QUADNODES(0,420-412) = 0.31979362782962989;
483 QUADNODES(0,421-412) = 0.056350832689629149;
484 QUADNODES(0,422-412) = 0.056350832689629149;
485 QUADNODES(0,423-412) = 0.056350832689629149;
486 QUADNODES(0,424-412) = 0.4436491673103708;
487 QUADNODES(0,425-412) = 0.4436491673103708;
488 QUADNODES(0,426-412) = 0.4436491673103708;
489
490 QUADNODES(1,412-412) = 0.2500000;
491 QUADNODES(1,413-412) = 0.091971078052723032;
492 QUADNODES(1,414-412) = 0.091971078052723032;
493 QUADNODES(1,415-412) = 0.72408676584183096;
494 QUADNODES(1,416-412) = 0.091971078052723032;
495 QUADNODES(1,417-412) = 0.31979362782962989;
496 QUADNODES(1,418-412) = 0.31979362782962989;
497 QUADNODES(1,419-412) = 0.040619116511110234;
498 QUADNODES(1,420-412) = 0.31979362782962989;
499 QUADNODES(1,421-412) = 0.056350832689629149;
500 QUADNODES(1,422-412) = 0.4436491673103708;
501 QUADNODES(1,423-412) = 0.4436491673103708;
502 QUADNODES(1,424-412) = 0.056350832689629149;
503 QUADNODES(1,425-412) = 0.056350832689629149;
504 QUADNODES(1,426-412) = 0.4436491673103708;
505
506 QUADNODES(2,412-412) = 0.2500000;
507 QUADNODES(2,413-412) = 0.091971078052723032;
508 QUADNODES(2,414-412) = 0.091971078052723032;
509 QUADNODES(2,415-412) = 0.091971078052723032;
510 QUADNODES(2,416-412) = 0.72408676584183096;
511 QUADNODES(2,417-412) = 0.31979362782962989;
512 QUADNODES(2,418-412) = 0.31979362782962989;
513 QUADNODES(2,419-412) = 0.31979362782962989;
514 QUADNODES(2,420-412) = 0.040619116511110234;
515 QUADNODES(2,421-412) = 0.4436491673103708;
516 QUADNODES(2,422-412) = 0.056350832689629149;
517 QUADNODES(2,423-412) = 0.4436491673103708;
518 QUADNODES(2,424-412) = 0.056350832689629149;
519 QUADNODES(2,425-412) = 0.4436491673103708;
520 QUADNODES(2,426-412) = 0.056350832689629149;
521
522 } else if (numQuadNodes==24){
523 QUADWEIGHTS(427-427) = 0.006653792;
524 QUADWEIGHTS(428-427) = 0.006653792;
525 QUADWEIGHTS(429-427) = 0.006653792;
526 QUADWEIGHTS(430-427) = 0.006653792;
527 QUADWEIGHTS(431-427) = 0.001679535;
528 QUADWEIGHTS(432-427) = 0.001679535;
529 QUADWEIGHTS(433-427) = 0.001679535;
530 QUADWEIGHTS(434-427) = 0.001679535;
531 QUADWEIGHTS(435-427) = 0.009226197;
532 QUADWEIGHTS(436-427) = 0.009226197;
533 QUADWEIGHTS(437-427) = 0.009226197;
534 QUADWEIGHTS(438-427) = 0.009226197;
535 QUADWEIGHTS(439-427) = 0.008035714;
536 QUADWEIGHTS(440-427) = 0.008035714;
537 QUADWEIGHTS(441-427) = 0.008035714;
538 QUADWEIGHTS(442-427) = 0.008035714;
539 QUADWEIGHTS(443-427) = 0.008035714;
540 QUADWEIGHTS(444-427) = 0.008035714;
541 QUADWEIGHTS(445-427) = 0.008035714;
542 QUADWEIGHTS(446-427) = 0.008035714;
543 QUADWEIGHTS(447-427) = 0.008035714;
544 QUADWEIGHTS(448-427) = 0.008035714;
545 QUADWEIGHTS(449-427) = 0.008035714;
546 QUADWEIGHTS(450-427) = 0.008035714;
547
548 QUADNODES(0,427-427) = 0.3561914;
549 QUADNODES(0,428-427) = 0.2146029;
550 QUADNODES(0,429-427) = 0.2146029;
551 QUADNODES(0,430-427) = 0.2146029;
552 QUADNODES(0,431-427) = 0.8779781;
553 QUADNODES(0,432-427) = 0.04067396;
554 QUADNODES(0,433-427) = 0.04067396;
555 QUADNODES(0,434-427) = 0.04067396;
556 QUADNODES(0,435-427) = 0.03298633;
557 QUADNODES(0,436-427) = 0.3223379;
558 QUADNODES(0,437-427) = 0.3223379;
559 QUADNODES(0,438-427) = 0.3223379;
560 QUADNODES(0,439-427) = 0.6030057;
561 QUADNODES(0,440-427) = 0.6030057;
562 QUADNODES(0,441-427) = 0.6030057;
563 QUADNODES(0,442-427) = 0.2696723;
564 QUADNODES(0,443-427) = 0.2696723;
565 QUADNODES(0,444-427) = 0.2696723;
566 QUADNODES(0,445-427) = 0.06366100;
567 QUADNODES(0,446-427) = 0.06366100;
568 QUADNODES(0,447-427) = 0.06366100;
569 QUADNODES(0,448-427) = 0.06366100;
570 QUADNODES(0,449-427) = 0.06366100;
571 QUADNODES(0,450-427) = 0.06366100;
572
573 QUADNODES(1,427-427) = 0.2146029;
574 QUADNODES(1,428-427) = 0.3561914;
575 QUADNODES(1,429-427) = 0.2146029;
576 QUADNODES(1,430-427) = 0.2146029;
577 QUADNODES(1,431-427) = 0.04067396;
578 QUADNODES(1,432-427) = 0.8779781;
579 QUADNODES(1,433-427) = 0.04067396;
580 QUADNODES(1,434-427) = 0.04067396;
581 QUADNODES(1,435-427) = 0.3223379;
582 QUADNODES(1,436-427) = 0.03298633;
583 QUADNODES(1,437-427) = 0.3223379;
584 QUADNODES(1,438-427) = 0.3223379;
585 QUADNODES(1,439-427) = 0.2696723;
586 QUADNODES(1,440-427) = 0.06366100;
587 QUADNODES(1,441-427) = 0.06366100;
588 QUADNODES(1,442-427) = 0.6030057;
589 QUADNODES(1,443-427) = 0.06366100;
590 QUADNODES(1,444-427) = 0.06366100;
591 QUADNODES(1,445-427) = 0.6030057;
592 QUADNODES(1,446-427) = 0.6030057;
593 QUADNODES(1,447-427) = 0.2696723;
594 QUADNODES(1,448-427) = 0.2696723;
595 QUADNODES(1,449-427) = 0.06366100;
596 QUADNODES(1,450-427) = 0.06366100;
597
598 QUADNODES(2,427-427) = 0.2146029;
599 QUADNODES(2,428-427) = 0.2146029;
600 QUADNODES(2,429-427) = 0.3561914;
601 QUADNODES(2,430-427) = 0.2146029;
602 QUADNODES(2,431-427) = 0.04067396;
603 QUADNODES(2,432-427) = 0.04067396;
604 QUADNODES(2,433-427) = 0.8779781;
605 QUADNODES(2,434-427) = 0.04067396;
606 QUADNODES(2,435-427) = 0.3223379;
607 QUADNODES(2,436-427) = 0.3223379;
608 QUADNODES(2,437-427) = 0.03298633;
609 QUADNODES(2,438-427) = 0.3223379;
610 QUADNODES(2,439-427) = 0.06366100;
611 QUADNODES(2,440-427) = 0.2696723;
612 QUADNODES(2,441-427) = 0.06366100;
613 QUADNODES(2,442-427) = 0.06366100;
614 QUADNODES(2,443-427) = 0.6030057;
615 QUADNODES(2,444-427) = 0.06366100;
616 QUADNODES(2,445-427) = 0.2696723;
617 QUADNODES(2,446-427) = 0.06366100;
618 QUADNODES(2,447-427) = 0.6030057;
619 QUADNODES(2,448-427) = 0.06366100;
620 QUADNODES(2,449-427) = 0.6030057;
621 QUADNODES(2,450-427) = 0.2696723;
622
623 } else if (numQuadNodes==31){
624 QUADWEIGHTS(451-451) = 0.01826422;
625 QUADWEIGHTS(452-451) = 0.01059994;
626 QUADWEIGHTS(453-451) = 0.01059994;
627 QUADWEIGHTS(454-451) = 0.01059994;
628 QUADWEIGHTS(455-451) = 0.01059994;
629 QUADWEIGHTS(456-451) =-0.06251774;
630 QUADWEIGHTS(457-451) =-0.06251774;
631 QUADWEIGHTS(458-451) =-0.06251774;
632 QUADWEIGHTS(459-451) =-0.06251774;
633 QUADWEIGHTS(460-451) = 0.004891425;
634 QUADWEIGHTS(461-451) = 0.004891425;
635 QUADWEIGHTS(462-451) = 0.004891425;
636 QUADWEIGHTS(463-451) = 0.004891425;
637 QUADWEIGHTS(464-451) = 0.0009700176;
638 QUADWEIGHTS(465-451) = 0.0009700176;
639 QUADWEIGHTS(466-451) = 0.0009700176;
640 QUADWEIGHTS(467-451) = 0.0009700176;
641 QUADWEIGHTS(468-451) = 0.0009700176;
642 QUADWEIGHTS(469-451) = 0.0009700176;
643 QUADWEIGHTS(470-451) = 0.02755732;
644 QUADWEIGHTS(471-451) = 0.02755732;
645 QUADWEIGHTS(472-451) = 0.02755732;
646 QUADWEIGHTS(473-451) = 0.02755732;
647 QUADWEIGHTS(474-451) = 0.02755732;
648 QUADWEIGHTS(475-451) = 0.02755732;
649 QUADWEIGHTS(476-451) = 0.02755732;
650 QUADWEIGHTS(477-451) = 0.02755732;
651 QUADWEIGHTS(478-451) = 0.02755732;
652 QUADWEIGHTS(479-451) = 0.02755732;
653 QUADWEIGHTS(480-451) = 0.02755732;
654 QUADWEIGHTS(481-451) = 0.02755732;
655
656 QUADNODES(0,451-451) = 0.2500000;
657 QUADNODES(0,452-451) = 0.7653604;
658 QUADNODES(0,453-451) = 0.07821319;
659 QUADNODES(0,454-451) = 0.07821319;
660 QUADNODES(0,455-451) = 0.07821319;
661 QUADNODES(0,456-451) = 0.6344704;
662 QUADNODES(0,457-451) = 0.1218432;
663 QUADNODES(0,458-451) = 0.1218432;
664 QUADNODES(0,459-451) = 0.1218432;
665 QUADNODES(0,460-451) = 0.002382507;
666 QUADNODES(0,461-451) = 0.3325392;
667 QUADNODES(0,462-451) = 0.3325392;
668 QUADNODES(0,463-451) = 0.3325392;
669 QUADNODES(0,464-451) = 0.0000000;
670 QUADNODES(0,465-451) = 0.0000000;
671 QUADNODES(0,466-451) = 0.0000000;
672 QUADNODES(0,467-451) = 0.5000000;
673 QUADNODES(0,468-451) = 0.5000000;
674 QUADNODES(0,469-451) = 0.5000000;
675 QUADNODES(0,470-451) = 0.6000000;
676 QUADNODES(0,471-451) = 0.6000000;
677 QUADNODES(0,472-451) = 0.6000000;
678 QUADNODES(0,473-451) = 0.2000000;
679 QUADNODES(0,474-451) = 0.2000000;
680 QUADNODES(0,475-451) = 0.2000000;
681 QUADNODES(0,476-451) = 0.1000000;
682 QUADNODES(0,477-451) = 0.1000000;
683 QUADNODES(0,478-451) = 0.1000000;
684 QUADNODES(0,479-451) = 0.1000000;
685 QUADNODES(0,480-451) = 0.1000000;
686 QUADNODES(0,481-451) = 0.1000000;
687
688 QUADNODES(1,451-451) = 0.2500000;
689 QUADNODES(1,452-451) = 0.07821319;
690 QUADNODES(1,453-451) = 0.7653604;
691 QUADNODES(1,454-451) = 0.07821319;
692 QUADNODES(1,455-451) = 0.07821319;
693 QUADNODES(1,456-451) = 0.1218432;
694 QUADNODES(1,457-451) = 0.6344704;
695 QUADNODES(1,458-451) = 0.1218432;
696 QUADNODES(1,459-451) = 0.1218432;
697 QUADNODES(1,460-451) = 0.3325392;
698 QUADNODES(1,461-451) = 0.002382507;
699 QUADNODES(1,462-451) = 0.3325392;
700 QUADNODES(1,463-451) = 0.3325392;
701 QUADNODES(1,464-451) = 0.0000000;
702 QUADNODES(1,465-451) = 0.5000000;
703 QUADNODES(1,466-451) = 0.5000000;
704 QUADNODES(1,467-451) = 0.0000000;
705 QUADNODES(1,468-451) = 0.0000000;
706 QUADNODES(1,469-451) = 0.5000000;
707 QUADNODES(1,470-451) = 0.2000000;
708 QUADNODES(1,471-451) = 0.1000000;
709 QUADNODES(1,472-451) = 0.1000000;
710 QUADNODES(1,473-451) = 0.6000000;
711 QUADNODES(1,474-451) = 0.1000000;
712 QUADNODES(1,475-451) = 0.1000000;
713 QUADNODES(1,476-451) = 0.6000000;
714 QUADNODES(1,477-451) = 0.6000000;
715 QUADNODES(1,478-451) = 0.2000000;
716 QUADNODES(1,479-451) = 0.2000000;
717 QUADNODES(1,480-451) = 0.1000000;
718 QUADNODES(1,481-451) = 0.1000000;
719
720 QUADNODES(2,451-451) = 0.2500000;
721 QUADNODES(2,452-451) = 0.07821319;
722 QUADNODES(2,453-451) = 0.07821319;
723 QUADNODES(2,454-451) = 0.7653604;
724 QUADNODES(2,455-451) = 0.07821319;
725 QUADNODES(2,456-451) = 0.1218432;
726 QUADNODES(2,457-451) = 0.1218432;
727 QUADNODES(2,458-451) = 0.6344704;
728 QUADNODES(2,459-451) = 0.1218432;
729 QUADNODES(2,460-451) = 0.3325392;
730 QUADNODES(2,461-451) = 0.3325392;
731 QUADNODES(2,462-451) = 0.002382507;
732 QUADNODES(2,463-451) = 0.3325392;
733 QUADNODES(2,464-451) = 0.5000000;
734 QUADNODES(2,465-451) = 0.0000000;
735 QUADNODES(2,466-451) = 0.5000000;
736 QUADNODES(2,467-451) = 0.0000000;
737 QUADNODES(2,468-451) = 0.5000000;
738 QUADNODES(2,469-451) = 0.0000000;
739 QUADNODES(2,470-451) = 0.1000000;
740 QUADNODES(2,471-451) = 0.2000000;
741 QUADNODES(2,472-451) = 0.1000000;
742 QUADNODES(2,473-451) = 0.1000000;
743 QUADNODES(2,474-451) = 0.6000000;
744 QUADNODES(2,475-451) = 0.1000000;
745 QUADNODES(2,476-451) = 0.2000000;
746 QUADNODES(2,477-451) = 0.1000000;
747 QUADNODES(2,478-451) = 0.6000000;
748 QUADNODES(2,479-451) = 0.1000000;
749 QUADNODES(2,480-451) = 0.6000000;
750 QUADNODES(2,481-451) = 0.2000000;
751
752 } else if (numQuadNodes==45){
753 QUADWEIGHTS(482-482) =-0.03932701;
754 QUADWEIGHTS(483-482) = 0.004081316;
755 QUADWEIGHTS(484-482) = 0.004081316;
756 QUADWEIGHTS(485-482) = 0.004081316;
757 QUADWEIGHTS(486-482) = 0.004081316;
758 QUADWEIGHTS(487-482) = 0.0006580868;
759 QUADWEIGHTS(488-482) = 0.0006580868;
760 QUADWEIGHTS(489-482) = 0.0006580868;
761 QUADWEIGHTS(490-482) = 0.0006580868;
762 QUADWEIGHTS(491-482) = 0.004384259;
763 QUADWEIGHTS(492-482) = 0.004384259;
764 QUADWEIGHTS(493-482) = 0.004384259;
765 QUADWEIGHTS(494-482) = 0.004384259;
766 QUADWEIGHTS(495-482) = 0.004384259;
767 QUADWEIGHTS(496-482) = 0.004384259;
768 QUADWEIGHTS(497-482) = 0.01383006;
769 QUADWEIGHTS(498-482) = 0.01383006;
770 QUADWEIGHTS(499-482) = 0.01383006;
771 QUADWEIGHTS(500-482) = 0.01383006;
772 QUADWEIGHTS(501-482) = 0.01383006;
773 QUADWEIGHTS(502-482) = 0.01383006;
774 QUADWEIGHTS(503-482) = 0.004240437;
775 QUADWEIGHTS(504-482) = 0.004240437;
776 QUADWEIGHTS(505-482) = 0.004240437;
777 QUADWEIGHTS(506-482) = 0.004240437;
778 QUADWEIGHTS(507-482) = 0.004240437;
779 QUADWEIGHTS(508-482) = 0.004240437;
780 QUADWEIGHTS(509-482) = 0.004240437;
781 QUADWEIGHTS(510-482) = 0.004240437;
782 QUADWEIGHTS(511-482) = 0.004240437;
783 QUADWEIGHTS(512-482) = 0.004240437;
784 QUADWEIGHTS(513-482) = 0.004240437;
785 QUADWEIGHTS(514-482) = 0.004240437;
786 QUADWEIGHTS(515-482) = 0.002238740;
787 QUADWEIGHTS(516-482) = 0.002238740;
788 QUADWEIGHTS(517-482) = 0.002238740;
789 QUADWEIGHTS(518-482) = 0.002238740;
790 QUADWEIGHTS(519-482) = 0.002238740;
791 QUADWEIGHTS(520-482) = 0.002238740;
792 QUADWEIGHTS(521-482) = 0.002238740;
793 QUADWEIGHTS(522-482) = 0.002238740;
794 QUADWEIGHTS(523-482) = 0.002238740;
795 QUADWEIGHTS(524-482) = 0.002238740;
796 QUADWEIGHTS(525-482) = 0.002238740;
797 QUADWEIGHTS(526-482) = 0.002238740;
798
799 QUADNODES(0,482-482) = 0.2500000;
800 QUADNODES(0,483-482) = 0.6175872;
801 QUADNODES(0,484-482) = 0.1274709;
802 QUADNODES(0,485-482) = 0.1274709;
803 QUADNODES(0,486-482) = 0.1274709;
804 QUADNODES(0,487-482) = 0.9037635;
805 QUADNODES(0,488-482) = 0.03207883;
806 QUADNODES(0,489-482) = 0.03207883;
807 QUADNODES(0,490-482) = 0.03207883;
808 QUADNODES(0,491-482) = 0.4502229;
809 QUADNODES(0,492-482) = 0.4502229;
810 QUADNODES(0,493-482) = 0.4502229;
811 QUADNODES(0,494-482) = 0.04977710;
812 QUADNODES(0,495-482) = 0.04977710;
813 QUADNODES(0,496-482) = 0.04977710;
814 QUADNODES(0,497-482) = 0.3162696;
815 QUADNODES(0,498-482) = 0.3162696;
816 QUADNODES(0,499-482) = 0.3162696;
817 QUADNODES(0,500-482) = 0.1837304;
818 QUADNODES(0,501-482) = 0.1837304;
819 QUADNODES(0,502-482) = 0.1837304;
820 QUADNODES(0,503-482) = 0.5132800;
821 QUADNODES(0,504-482) = 0.5132800;
822 QUADNODES(0,505-482) = 0.5132800;
823 QUADNODES(0,506-482) = 0.02291779;
824 QUADNODES(0,507-482) = 0.02291779;
825 QUADNODES(0,508-482) = 0.02291779;
826 QUADNODES(0,509-482) = 0.2319011;
827 QUADNODES(0,510-482) = 0.2319011;
828 QUADNODES(0,511-482) = 0.2319011;
829 QUADNODES(0,512-482) = 0.2319011;
830 QUADNODES(0,513-482) = 0.2319011;
831 QUADNODES(0,514-482) = 0.2319011;
832 QUADNODES(0,515-482) = 0.1937465;
833 QUADNODES(0,516-482) = 0.1937465;
834 QUADNODES(0,517-482) = 0.1937465;
835 QUADNODES(0,518-482) = 0.7303134;
836 QUADNODES(0,519-482) = 0.7303134;
837 QUADNODES(0,520-482) = 0.7303134;
838 QUADNODES(0,521-482) = 0.03797005;
839 QUADNODES(0,522-482) = 0.03797005;
840 QUADNODES(0,523-482) = 0.03797005;
841 QUADNODES(0,524-482) = 0.03797005;
842 QUADNODES(0,525-482) = 0.03797005;
843 QUADNODES(0,526-482) = 0.03797005;
844
845 QUADNODES(1,482-482) = 0.2500000;
846 QUADNODES(1,483-482) = 0.1274709;
847 QUADNODES(1,484-482) = 0.6175872;
848 QUADNODES(1,485-482) = 0.1274709;
849 QUADNODES(1,486-482) = 0.1274709;
850 QUADNODES(1,487-482) = 0.03207883;
851 QUADNODES(1,488-482) = 0.9037635;
852 QUADNODES(1,489-482) = 0.03207883;
853 QUADNODES(1,490-482) = 0.03207883;
854 QUADNODES(1,491-482) = 0.4502229;
855 QUADNODES(1,492-482) = 0.04977710;
856 QUADNODES(1,493-482) = 0.04977710;
857 QUADNODES(1,494-482) = 0.4502229;
858 QUADNODES(1,495-482) = 0.4502229;
859 QUADNODES(1,496-482) = 0.04977710;
860 QUADNODES(1,497-482) = 0.3162696;
861 QUADNODES(1,498-482) = 0.1837304;
862 QUADNODES(1,499-482) = 0.1837304;
863 QUADNODES(1,500-482) = 0.3162696;
864 QUADNODES(1,501-482) = 0.3162696;
865 QUADNODES(1,502-482) = 0.1837304;
866 QUADNODES(1,503-482) = 0.02291779;
867 QUADNODES(1,504-482) = 0.2319011;
868 QUADNODES(1,505-482) = 0.2319011;
869 QUADNODES(1,506-482) = 0.5132800;
870 QUADNODES(1,507-482) = 0.2319011;
871 QUADNODES(1,508-482) = 0.2319011;
872 QUADNODES(1,509-482) = 0.5132800;
873 QUADNODES(1,510-482) = 0.5132800;
874 QUADNODES(1,511-482) = 0.02291779;
875 QUADNODES(1,512-482) = 0.02291779;
876 QUADNODES(1,513-482) = 0.2319011;
877 QUADNODES(1,514-482) = 0.2319011;
878 QUADNODES(1,515-482) = 0.7303134;
879 QUADNODES(1,516-482) = 0.03797005;
880 QUADNODES(1,517-482) = 0.03797005;
881 QUADNODES(1,518-482) = 0.1937465;
882 QUADNODES(1,519-482) = 0.03797005;
883 QUADNODES(1,520-482) = 0.03797005;
884 QUADNODES(1,521-482) = 0.1937465;
885 QUADNODES(1,522-482) = 0.1937465;
886 QUADNODES(1,523-482) = 0.7303134;
887 QUADNODES(1,524-482) = 0.7303134;
888 QUADNODES(1,525-482) = 0.03797005;
889 QUADNODES(1,526-482) = 0.03797005;
890
891 QUADNODES(2,482-482) = 0.2500000;
892 QUADNODES(2,483-482) = 0.1274709;
893 QUADNODES(2,484-482) = 0.1274709;
894 QUADNODES(2,485-482) = 0.6175872;
895 QUADNODES(2,486-482) = 0.1274709;
896 QUADNODES(2,487-482) = 0.03207883;
897 QUADNODES(2,488-482) = 0.03207883;
898 QUADNODES(2,489-482) = 0.9037635;
899 QUADNODES(2,490-482) = 0.03207883;
900 QUADNODES(2,491-482) = 0.04977710;
901 QUADNODES(2,492-482) = 0.4502229;
902 QUADNODES(2,493-482) = 0.04977710;
903 QUADNODES(2,494-482) = 0.4502229;
904 QUADNODES(2,495-482) = 0.04977710;
905 QUADNODES(2,496-482) = 0.4502229;
906 QUADNODES(2,497-482) = 0.1837304;
907 QUADNODES(2,498-482) = 0.3162696;
908 QUADNODES(2,499-482) = 0.1837304;
909 QUADNODES(2,500-482) = 0.3162696;
910 QUADNODES(2,501-482) = 0.1837304;
911 QUADNODES(2,502-482) = 0.3162696;
912 QUADNODES(2,503-482) = 0.2319011;
913 QUADNODES(2,504-482) = 0.02291779;
914 QUADNODES(2,505-482) = 0.2319011;
915 QUADNODES(2,506-482) = 0.2319011;
916 QUADNODES(2,507-482) = 0.5132800;
917 QUADNODES(2,508-482) = 0.2319011;
918 QUADNODES(2,509-482) = 0.02291779;
919 QUADNODES(2,510-482) = 0.2319011;
920 QUADNODES(2,511-482) = 0.5132800;
921 QUADNODES(2,512-482) = 0.2319011;
922 QUADNODES(2,513-482) = 0.5132800;
923 QUADNODES(2,514-482) = 0.02291779;
924 QUADNODES(2,515-482) = 0.03797005;
925 QUADNODES(2,516-482) = 0.7303134;
926 QUADNODES(2,517-482) = 0.03797005;
927 QUADNODES(2,518-482) = 0.03797005;
928 QUADNODES(2,519-482) = 0.1937465;
929 QUADNODES(2,520-482) = 0.03797005;
930 QUADNODES(2,521-482) = 0.7303134;
931 QUADNODES(2,522-482) = 0.03797005;
932 QUADNODES(2,523-482) = 0.1937465;
933 QUADNODES(2,524-482) = 0.03797005;
934 QUADNODES(2,525-482) = 0.1937465;
935 QUADNODES(2,526-482) = 0.7303134;
936
937 } else {
938
939 /* get scheme on [0.1]^3 */
940
941 Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
942 if (! Finley_noError()) return;
943
944 /* squeeze it: */
945 for (i=0;i<numQuadNodes;i++) {
946 Q1=QUADNODES(0,i);
947 Q2=QUADNODES(1,i);
948 Q3=QUADNODES(2,i);
949
950 JA11= (1./3.)*Q2*Q3-(1./2.)*(Q2+Q3) +1.;
951 JA12= (1./3.)*Q1*Q3-(1./2.)*Q1;
952 JA13= (1./3.)*Q1*Q2-(1./2.)*Q1;
953 JA21= (1./3.)*Q2*Q3-(1./2.)*Q2;
954 JA22= (1./3.)*Q1*Q3-(1./2.)*(Q1+Q3) +1.;
955 JA23= (1./3.)*Q1*Q2-(1./2.)*Q2;
956 JA31= (1./3.)*Q2*Q3-(1./2.)*Q3;
957 JA32= (1./3.)*Q1*Q3-(1./2.)*Q3;
958 JA33= (1./3.)*Q1*Q2-(1./2.)*(Q1+Q2) +1.;
959 DET=JA11*JA22*JA33+JA12*JA23*JA31+JA13*JA21*JA32-JA13*JA22*JA31-JA11*JA23*JA32-JA12*JA21*JA33;
960 quadWeights[i]=quadWeights[i]*ABS(DET);
961 QUADNODES(0,i)=Q1*((1./3.)*Q2*Q3-(1./2.)*(Q2+Q3)+1.);
962 QUADNODES(1,i)=Q2*((1./3.)*Q1*Q3-(1./2.)*(Q1+Q3)+1.);
963 QUADNODES(2,i)=Q3*((1./3.)*Q1*Q2-(1./2.)*(Q1+Q2)+1.);
964 }
965 }
966 #undef DIM
967 }
968
969 /**************************************************************/
970
971 /* get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */
972 /* as a X-product of a 1D scheme. */
973
974 void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
975 char error_msg[LenErrorMsg_MAX];
976 int numQuadNodes1d,i,j,l;
977 double *quadNodes1d=NULL,*quadWeights1d=NULL;
978 bool_t set=FALSE;
979 #define DIM 2
980
981 quadNodes1d=TMPMEMALLOC(numQuadNodes,double);
982 quadWeights1d=TMPMEMALLOC(numQuadNodes,double);
983 if (! ( Finley_checkPtr(quadNodes1d) || Finley_checkPtr(quadWeights1d) ) ) {
984 /* find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */
985
986 for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
987 if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
988
989 /* get 1D scheme: */
990
991 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
992
993 /* make 2D scheme: */
994
995 if (Finley_noError()) {
996 l=0;
997 for (i=0;i<numQuadNodes1d;i++) {
998 for (j=0;j<numQuadNodes1d;j++) {
999 QUADNODES(0,l)=quadNodes1d[i];
1000 QUADNODES(1,l)=quadNodes1d[j];
1001 QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];
1002 l++;
1003 }
1004 }
1005 set=TRUE;
1006 break;
1007 }
1008 }
1009 }
1010 if (!set) {
1011 sprintf(error_msg,"Finley_Quad_getNodesRec: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
1012 Finley_setError(VALUE_ERROR,error_msg);
1013 }
1014 TMPMEMFREE(quadNodes1d);
1015 TMPMEMFREE(quadWeights1d);
1016 }
1017 #undef DIM
1018 }
1019
1020 /**************************************************************/
1021
1022 /* get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
1023 /* as a X-product of a 1D scheme. */
1024
1025 void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
1026 char error_msg[LenErrorMsg_MAX];
1027 int numQuadNodes1d,i,j,k,l;
1028 double *quadNodes1d=NULL,*quadWeights1d=NULL;
1029 bool_t set=FALSE;
1030 #define DIM 3
1031
1032 /* find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
1033
1034 quadNodes1d=TMPMEMALLOC(numQuadNodes,double);
1035 quadWeights1d=TMPMEMALLOC(numQuadNodes,double);
1036 if (! ( Finley_checkPtr(quadNodes1d) || Finley_checkPtr(quadWeights1d) ) ) {
1037 for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
1038 if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
1039
1040 /* get 1D scheme: */
1041
1042 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
1043
1044 /* make 3D scheme: */
1045
1046 if (Finley_noError()) {
1047 l=0;
1048 for (i=0;i<numQuadNodes1d;i++) {
1049 for (j=0;j<numQuadNodes1d;j++) {
1050 for (k=0;k<numQuadNodes1d;k++) {
1051 QUADNODES(0,l)=quadNodes1d[i];
1052 QUADNODES(1,l)=quadNodes1d[j];
1053 QUADNODES(2,l)=quadNodes1d[k];
1054 QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];
1055 l++;
1056 }
1057 }
1058 }
1059 set=TRUE;
1060 break;
1061 }
1062 }
1063 }
1064 if (!set) {
1065 sprintf(error_msg,"Finley_Quad_getNodesHex: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
1066 Finley_setError(VALUE_ERROR,error_msg);
1067 }
1068 TMPMEMFREE(quadNodes1d);
1069 TMPMEMFREE(quadWeights1d);
1070 }
1071 #undef DIM
1072 }
1073
1074 /**************************************************************/
1075
1076 /* get a quadrature scheme with numQuadNodes quadrature nodes for a point. As there */
1077 /* in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
1078 /* an error. */
1079
1080 void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
1081 if (numQuadNodes!=0) {
1082 Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: There is no quadrature scheme on points.");
1083 }
1084 }
1085
1086 /**************************************************************/
1087
1088 /* get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
1089 /* The nodes and weights are set from a table. */
1090
1091 void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
1092 switch(numQuadNodes) {
1093 case 1:
1094 quadNodes[0]=0.5;
1095 quadWeights[0]=1.;
1096 break;
1097
1098 case 2:
1099 quadNodes[0]=(1.-.577350269189626)/2.;
1100 quadNodes[1]=(1.+.577350269189626)/2.;
1101 quadWeights[0]=.5;
1102 quadWeights[1]=.5;
1103 break;
1104
1105 case 3:
1106 quadNodes[0]=(1.-.774596669241483)/2.;
1107 quadNodes[1]=.5;
1108 quadNodes[2]=(1.+.774596669241483)/2.;
1109 quadWeights[0]=5./18.;
1110 quadWeights[1]=4./ 9.;
1111 quadWeights[2]=5./18.;
1112 break;
1113
1114 case 4:
1115 quadNodes[0]=(1.-.861136311594053)/2.;
1116 quadNodes[1]=(1.-.339981043584856)/2.;
1117 quadNodes[2]=(1.+.339981043584856)/2.;
1118 quadNodes[3]=(1.+.861136311594053)/2.;
1119 quadWeights[0]=.347854845137454/2.;
1120 quadWeights[1]=.652145154862546/2.;
1121 quadWeights[2]=.652145154862546/2.;
1122 quadWeights[3]=.347854845137454/2.;
1123 break;
1124
1125 case 5:
1126 quadNodes[0]=(1.-.906179845938664)/2.;
1127 quadNodes[1]=(1.-.538469310105683)/2.;
1128 quadNodes[2]= .5;
1129 quadNodes[3]=(1.+.538469310105683)/2.;
1130 quadNodes[4]=(1.+.906179845938664)/2.;
1131 quadWeights[0]=.236926885056189/2.;
1132 quadWeights[1]=.478628670499366/2.;
1133 quadWeights[2]=.568888888888889/2.;
1134 quadWeights[3]=.478628670499366/2.;
1135 quadWeights[4]=.236926885056189/2.;
1136 break;
1137
1138 case 6:
1139 quadNodes[0]=(1.-.932469514203152)/2.;
1140 quadNodes[1]=(1.-.661209386466265)/2.;
1141 quadNodes[2]=(1.-.238619186083197)/2.;
1142 quadNodes[3]=(1.+.238619186083197)/2.;
1143 quadNodes[4]=(1.+.661209386466265)/2.;
1144 quadNodes[5]=(1.+.932469514203152)/2.;
1145 quadWeights[0]=.171324492379170/2.;
1146 quadWeights[1]=.360761573048139/2.;
1147 quadWeights[2]=.467913934572691/2.;
1148 quadWeights[3]=.467913934572691/2.;
1149 quadWeights[4]=.360761573048139/2.;
1150 quadWeights[5]=.171324492379170/2.;
1151 break;
1152
1153 case 7:
1154 quadNodes[0]=(1.-.949107912342759)/2.;
1155 quadNodes[1]=(1.-.741531185599394)/2.;
1156 quadNodes[2]=(1.-.405845151377397)/2.;
1157 quadNodes[3]=0.5;
1158 quadNodes[4]=(1.+.405845151377397)/2.;
1159 quadNodes[5]=(1.+.741531185599394)/2.;
1160 quadNodes[6]=(1.+.949107912342759)/2.;
1161 quadWeights[0]= .129484966168870/2.;
1162 quadWeights[1]= .279705391489277/2.;
1163 quadWeights[2]= .381830050505119/2.;
1164 quadWeights[3]= .417959183673469/2.;
1165 quadWeights[4]= .381830050505119/2.;
1166 quadWeights[5]= .279705391489277/2.;
1167 quadWeights[6]= .129484966168870/2.;
1168 break;
1169
1170 case 8:
1171 quadNodes[0]=(1.-.960289856497536)/2.;
1172 quadNodes[1]=(1.-.796666477413627)/2.;
1173 quadNodes[2]=(1.-.525532409916329)/2.;
1174 quadNodes[3]=(1.-.183434642495650)/2.;
1175 quadNodes[4]=(1.+.183434642495650)/2.;
1176 quadNodes[5]=(1.+.525532409916329)/2.;
1177 quadNodes[6]=(1.+.796666477413627)/2.;
1178 quadNodes[7]=(1.+.960289856497536)/2.;
1179 quadWeights[0]= .101228536290376/2.;
1180 quadWeights[1]= .222381034453374/2.;
1181 quadWeights[2]= .313706645877887/2.;
1182 quadWeights[3]= .362683783378362/2.;
1183 quadWeights[4]= .362683783378362/2.;
1184 quadWeights[5]= .313706645877887/2.;
1185 quadWeights[6]= .222381034453374/2.;
1186 quadWeights[7]= .101228536290376/2.;
1187 break;
1188
1189 case 9:
1190 quadNodes[0]=(1.-.968160239507626)/2.;
1191 quadNodes[1]=(1.-.836031107326636)/2.;
1192 quadNodes[2]=(1.-.613371432700590)/2.;
1193 quadNodes[3]=(1.-.324253423403809)/2.;
1194 quadNodes[4]= .5;
1195 quadNodes[5]=(1.+.324253423403809)/2.;
1196 quadNodes[6]=(1.+.613371432700590)/2.;
1197 quadNodes[7]=(1.+.836031107326636)/2.;
1198 quadNodes[8]=(1.+.968160239507626)/2.;
1199 quadWeights[0]= .081274388361574/2.;
1200 quadWeights[1]= .180648160694857/2.;
1201 quadWeights[2]= .260610696402935/2.;
1202 quadWeights[3]= .312347077040003/2.;
1203 quadWeights[4]= .330239355001260/2.;
1204 quadWeights[5]= .312347077040003/2.;
1205 quadWeights[6]= .260610696402935/2.;
1206 quadWeights[7]= .180648160694857/2.;
1207 quadWeights[8]= .081274388361574/2.;
1208 break;
1209
1210 case 10:
1211 quadNodes[0]=(1.-.973906528517172)/2.;
1212 quadNodes[1]=(1.-.865063366688985)/2.;
1213 quadNodes[2]=(1.-.679409568299024)/2.;
1214 quadNodes[3]=(1.-.433395394129247)/2.;
1215 quadNodes[4]=(1.-.148874338981631)/2.;
1216 quadNodes[5]=(1.+.148874338981631)/2.;
1217 quadNodes[6]=(1.+.433395394129247)/2.;
1218 quadNodes[7]=(1.+.679409568299024)/2.;
1219 quadNodes[8]=(1.+.865063366688985)/2.;
1220 quadNodes[9]=(1.+.973906528517172)/2.;
1221 quadWeights[0]= .066671344308688/2.;
1222 quadWeights[1]= .149451349150581/2.;
1223 quadWeights[2]= .219086362515982/2.;
1224 quadWeights[3]= .269266719309996/2.;
1225 quadWeights[4]= .295524224714753/2.;
1226 quadWeights[5]= .295524224714753/2.;
1227 quadWeights[6]= .269266719309996/2.;
1228 quadWeights[7]= .219086362515982/2.;
1229 quadWeights[8]= .149451349150581/2.;
1230 quadWeights[9]= .066671344308688/2.;
1231 break;
1232
1233 default:
1234 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
1235 break;
1236 }
1237 }
1238 /**************************************************************/
1239 /* */
1240 /* the following function are used define the meshes on the surface in the xy-plane */
1241
1242 /* triangle surface on a tetrahedron */
1243 void Finley_Quad_getNodesTriOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
1244 Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesTri);
1245 }
1246 /* rectangular surface on a hexahedron */
1247 void Finley_Quad_getNodesRecOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
1248 Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesRec);
1249 }
1250 /* line surface on a triangle or rectangle */
1251 void Finley_Quad_getNodesLineOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
1252 Finley_Quad_makeNodesOnFace(2,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesLine);
1253 }
1254 /* point surface on a line */
1255 void Finley_Quad_getNodesPointOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
1256 Finley_Quad_makeNodesOnFace(1,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesPoint);
1257 }
1258
1259 void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {
1260 int q,i;
1261 double *quadNodesOnFace=NULL;
1262 #define DIM dim
1263 quadNodesOnFace=TMPMEMALLOC(numQuadNodes*(dim-1),double);
1264
1265 if (! Finley_checkPtr(quadNodesOnFace) ) {
1266 getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);
1267 if (Finley_noError()) {
1268 for (q=0;q<numQuadNodes;q++) {
1269 for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];
1270 QUADNODES(dim-1,q)=0;
1271 }
1272 }
1273 TMPMEMFREE(quadNodesOnFace);
1274 }
1275 #undef DIM
1276 }
1277
1278 /**************************************************************/
1279
1280 /* The following functions Finley_Quad_getNumNodes* return the nmber of quadrature points needed to */
1281 /* achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
1282 /* to consider the accuracy reduction through the construction process. */
1283
1284
1285 int Finley_Quad_getNumNodesPoint(int order) {
1286 if (order <0 ) {
1287 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
1288 return -1;
1289 } else {
1290 return 0;
1291 }
1292 }
1293
1294 int Finley_Quad_getNumNodesLine(int order) {
1295 char error_msg[LenErrorMsg_MAX];
1296 if (order <0 ) {
1297 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
1298 return -1;
1299 } else {
1300 if (order > 2*MAX_numQuadNodesLine-1) {
1301 sprintf(error_msg,"Finley_Quad_getNumNodesPoint: requested integration order %d on line is too large (>%d).",
1302 order,2*MAX_numQuadNodesLine-1);
1303 Finley_setError(VALUE_ERROR,error_msg);
1304 return -1;
1305 } else {
1306 Finley_resetError();
1307 return order/2+1;
1308 }
1309 }
1310 }
1311
1312 int Finley_Quad_getNumNodesTri(int order) {
1313 int numQuadNodesLine;
1314 if (order<=1) {
1315 return 1;
1316 } else if (order<=2){
1317 return 3;
1318 } else if (order<=3){
1319 return 4;
1320 } else if (order<=4){
1321 return 6;
1322 } else if (order<=5){
1323 return 7;
1324 } else if (order<=6){
1325 return 12;
1326 } else if (order<=7){
1327 return 13;
1328 } else if (order<=8){
1329 return 16;
1330 } else if (order<=9){
1331 return 19;
1332 } else {
1333 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
1334 if (Finley_noError()) {
1335 return numQuadNodesLine*numQuadNodesLine;
1336 } else {
1337 return -1;
1338 }
1339 }
1340 }
1341
1342 int Finley_Quad_getNumNodesRec(int order) {
1343 int numQuadNodesLine;
1344 numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
1345 if (Finley_noError()) {
1346 return numQuadNodesLine*numQuadNodesLine;
1347 } else {
1348 return -1;
1349 }
1350 }
1351
1352 int Finley_Quad_getNumNodesTet(int order) {
1353 int numQuadNodesLine;
1354 if (order<=1) {
1355 return 1;
1356 } else if (order<=2){
1357 return 4;
1358 } else if (order<=3){
1359 return 5;
1360 } else if (order<=4){
1361 return 11;
1362 } else if (order<=5){
1363 return 15;
1364 } else if (order<=6){
1365 return 24;
1366 } else if (order<=7){
1367 return 31;
1368 } else if (order<=8){
1369 return 45;
1370 } else {
1371 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
1372 if (Finley_noError()) {
1373 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1374 } else {
1375 return -1;
1376 }
1377 }
1378 }
1379
1380 int Finley_Quad_getNumNodesHex(int order) {
1381 int numQuadNodesLine;
1382 numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
1383 if (Finley_noError()) {
1384 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1385 } else {
1386 return -1;
1387 }
1388 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26