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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2034 - (hide 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 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 jgs 82
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 ksteube 1312 /* Finley: */
18 jgs 82
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 gross 1342 double Q1,Q2,a,b,c,d,e,f,g,u,v,w;
34 jgs 82 #define DIM 2
35    
36 btully 1170 /* the easy cases: */
37 jgs 82
38     if (numQuadNodes==1) {
39     QUADNODES(0,0)=1./3.;
40     QUADNODES(1,0)=1./3.;
41 btully 1170 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 gross 1342 } 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 jgs 82 } else {
320    
321     /* get scheme on [0.1]^2 */
322     Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
323 jgs 150 if (! Finley_noError()) return;
324 jgs 82
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 gross 1342
337    
338 jgs 82 }
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 gross 1342 double a,b,c,d,e,f,g,h;
349     double alpha=0.58541019662496852;
350     double beta =0.1381966011250105;
351 jgs 82 #define DIM 3
352    
353 btully 1170 /* the easy cases: */
354 jgs 82 if (numQuadNodes==1) {
355 btully 1170 QUADNODES(0,0)=0.25;
356     QUADNODES(1,0)=0.25;
357     QUADNODES(2,0)=0.25;
358 jgs 82 QUADWEIGHTS(0)=1./6.;
359 btully 1170 } 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 gross 1342
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 jgs 82 } else {
938    
939     /* get scheme on [0.1]^3 */
940    
941     Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
942 jgs 150 if (! Finley_noError()) return;
943 jgs 82
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 jgs 150 char error_msg[LenErrorMsg_MAX];
976 jgs 82 int numQuadNodes1d,i,j,l;
977 gross 1028 double *quadNodes1d=NULL,*quadWeights1d=NULL;
978     bool_t set=FALSE;
979 jgs 82 #define DIM 2
980    
981 gross 1028 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 jgs 82
989 gross 1028 /* get 1D scheme: */
990    
991     Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
992 jgs 82
993 gross 1028 /* make 2D scheme: */
994 jgs 82
995 gross 1028 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 jgs 82 }
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 jgs 150 char error_msg[LenErrorMsg_MAX];
1027 jgs 82 int numQuadNodes1d,i,j,k,l;
1028 gross 1028 double *quadNodes1d=NULL,*quadWeights1d=NULL;
1029 artak 2034 bool_t set=FALSE;
1030 jgs 82 #define DIM 3
1031    
1032     /* find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
1033    
1034 gross 1028 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 jgs 82
1040 gross 1028 /* get 1D scheme: */
1041 jgs 82
1042 gross 1028 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
1043 jgs 82
1044 gross 1028 /* make 3D scheme: */
1045 jgs 82
1046 gross 1028 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 jgs 82 }
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 gross 964 Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: There is no quadrature scheme on points.");
1083 jgs 82 }
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 jgs 150 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
1235 jgs 82 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 gross 1028 double *quadNodesOnFace=NULL;
1262     #define DIM dim
1263     quadNodesOnFace=TMPMEMALLOC(numQuadNodes*(dim-1),double);
1264 jgs 82
1265 gross 1028 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 jgs 82 }
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 gross 964 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
1288 jgs 82 return -1;
1289     } else {
1290     return 0;
1291     }
1292     }
1293    
1294     int Finley_Quad_getNumNodesLine(int order) {
1295 jgs 150 char error_msg[LenErrorMsg_MAX];
1296 jgs 82 if (order <0 ) {
1297 gross 964 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
1298 jgs 82 return -1;
1299     } else {
1300     if (order > 2*MAX_numQuadNodesLine-1) {
1301 gross 964 sprintf(error_msg,"Finley_Quad_getNumNodesPoint: requested integration order %d on line is too large (>%d).",
1302 jgs 82 order,2*MAX_numQuadNodesLine-1);
1303 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
1304 jgs 82 return -1;
1305     } else {
1306 jgs 150 Finley_resetError();
1307 jgs 82 return order/2+1;
1308     }
1309     }
1310     }
1311    
1312     int Finley_Quad_getNumNodesTri(int order) {
1313     int numQuadNodesLine;
1314 gross 1072 if (order<=1) {
1315 jgs 82 return 1;
1316 gross 1342 } else if (order<=2){
1317 btully 1170 return 3;
1318 gross 1342 } else if (order<=3){
1319 btully 1170 return 4;
1320 gross 1342 } 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 jgs 82 } else {
1333     numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
1334 jgs 150 if (Finley_noError()) {
1335 jgs 82 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 jgs 150 if (Finley_noError()) {
1346 jgs 82 return numQuadNodesLine*numQuadNodesLine;
1347     } else {
1348     return -1;
1349     }
1350     }
1351    
1352     int Finley_Quad_getNumNodesTet(int order) {
1353     int numQuadNodesLine;
1354 gross 1072 if (order<=1) {
1355     return 1;
1356 gross 1342 } else if (order<=2){
1357 btully 1170 return 4;
1358 gross 1342 } else if (order<=3){
1359 btully 1170 return 5;
1360 gross 1342 } 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 jgs 82 } else {
1371 gross 1072 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
1372     if (Finley_noError()) {
1373     return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1374     } else {
1375     return -1;
1376     }
1377 jgs 82 }
1378     }
1379    
1380     int Finley_Quad_getNumNodesHex(int order) {
1381     int numQuadNodesLine;
1382     numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
1383 jgs 150 if (Finley_noError()) {
1384 jgs 82 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