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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 68083 byte(s)
Merging dudley and scons updates from branches

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 gross 2748 /* Finley: quadrature schemes */
18 jgs 82
19     /**************************************************************/
20    
21     #include "Quadrature.h"
22 jfenwick 3259 #include "esysUtils/index.h"
23     #include "esysUtils/mem.h"
24 jgs 82
25 gross 2748
26 jgs 82 #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
27     #define QUADWEIGHTS(_I_) quadWeights[_I_]
28    
29     /**************************************************************/
30    
31 gross 2748 Finley_QuadInfo Finley_QuadInfoList[]={
32     {PointQuad, "Point", 0, 1, Finley_Quad_getNodesPoint, Finley_Quad_getNumNodesPoint, Finley_Quad_MacroPoint} ,
33 gross 2850 {LineQuad, "Line", 1, 2, Finley_Quad_getNodesLine, Finley_Quad_getNumNodesLine, Finley_Quad_MacroLine} ,
34     {TriQuad, "Tri", 2, 3, Finley_Quad_getNodesTri, Finley_Quad_getNumNodesTri, Finley_Quad_MacroTri},
35 gross 2748 {RecQuad, "Rec", 2, 4, Finley_Quad_getNodesRec, Finley_Quad_getNumNodesRec, Finley_Quad_MacroRec},
36     {TetQuad, "Tet", 3, 4, Finley_Quad_getNodesTet, Finley_Quad_getNumNodesTet, Finley_Quad_MacroTet},
37     {HexQuad, "Hex", 3, 8, Finley_Quad_getNodesHex, Finley_Quad_getNumNodesHex, Finley_Quad_MacroHex},
38     {NoQuad, "NoType", 0, 1, Finley_Quad_getNodesPoint, Finley_Quad_getNumNodesPoint, Finley_Quad_MacroPoint}
39     };
40    
41     Finley_QuadInfo* Finley_QuadInfo_getInfo(Finley_QuadTypeId id)
42     {
43     int ptr=0;
44     Finley_QuadInfo* out=NULL;
45     while (Finley_QuadInfoList[ptr].TypeId!=NoQuad && out==NULL) {
46     if (Finley_QuadInfoList[ptr].TypeId==id) out=&(Finley_QuadInfoList[ptr]);
47     ptr++;
48     }
49     if (out==NULL) {
50     Finley_setError(VALUE_ERROR,"Finley_QuadInfo_getInfo: canot find requested quadrature scheme.");
51     }
52     return out;
53     }
54    
55     /**************************************************************/
56    
57 jgs 82 /* get a quadrature scheme with numQuadNodes quadrature nodes for the tri */
58     /* as a queezed scheme on a quad [0,1]^2 */
59    
60     void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
61     int i;
62 gross 1342 double Q1,Q2,a,b,c,d,e,f,g,u,v,w;
63 jgs 82 #define DIM 2
64    
65 btully 1170 /* the easy cases: */
66 jgs 82
67     if (numQuadNodes==1) {
68     QUADNODES(0,0)=1./3.;
69     QUADNODES(1,0)=1./3.;
70 btully 1170 QUADWEIGHTS(0)=1./2.;
71     } else if (numQuadNodes==3){
72     QUADNODES(0,0)=1./2.;
73     QUADNODES(1,0)=0.;
74     QUADWEIGHTS(0)=1./6.;
75     QUADNODES(0,1)=0.;
76     QUADNODES(1,1)=1./2.;
77     QUADWEIGHTS(1)=1./6.;
78     QUADNODES(0,2)=1./2.;
79     QUADNODES(1,2)=1./2.;
80     QUADWEIGHTS(2)=1./6.;
81     } else if (numQuadNodes==4){
82     QUADNODES(0,0)=1./3.;
83     QUADNODES(1,0)=1./3.;
84     QUADWEIGHTS(0)=-27./96.;
85     QUADNODES(0,1)=0.2;
86     QUADNODES(1,1)=0.2;
87     QUADWEIGHTS(1)=25./96.;
88     QUADNODES(0,2)=0.6;
89     QUADNODES(1,2)=0.2;
90     QUADWEIGHTS(2)=25./96.;
91     QUADNODES(0,3)=0.2;
92     QUADNODES(1,3)=0.6;
93     QUADWEIGHTS(3)=25./96.;
94 gross 1342 } else if (numQuadNodes==6){
95     QUADWEIGHTS(0) = 0.109951743655322/2.;
96     QUADWEIGHTS(1) = 0.109951743655322/2.;
97     QUADWEIGHTS(2) = 0.109951743655322/2.;
98     QUADWEIGHTS(3) = 0.223381589678011/2.;
99     QUADWEIGHTS(4) = 0.223381589678011/2.;
100     QUADWEIGHTS(5) = 0.223381589678011/2.;
101    
102     QUADNODES(0,0) = 0.816847572980459;
103     QUADNODES(0,1) = 0.091576213509771;
104     QUADNODES(0,2) = 0.091576213509771;
105     QUADNODES(0,3) = 0.108103018168070;
106     QUADNODES(0,4) = 0.445948490915965;
107     QUADNODES(0,5) = 0.445948490915965;
108    
109     QUADNODES(1,0) = 0.091576213509771;
110     QUADNODES(1,1) = 0.816847572980459;
111     QUADNODES(1,2) = 0.091576213509771;
112     QUADNODES(1,3) = 0.445948490915965;
113     QUADNODES(1,4) = 0.108103018168070;
114     QUADNODES(1,5) = 0.445948490915965;
115    
116     } else if (numQuadNodes==7){
117     QUADNODES(0,0) = 0.33333333333333333;
118     QUADNODES(0,1) = 0.7974269853530872;
119     QUADNODES(0,2) = 0.10128650732345633;
120     QUADNODES(0,3) = 0.10128650732345633;
121     QUADNODES(0,4) = 0.059715871789769809;
122     QUADNODES(0,5) = 0.47014206410511505;
123     QUADNODES(0,6) = 0.47014206410511505;
124    
125     QUADNODES(1,0) = 0.33333333333333333;
126     QUADNODES(1,1) = 0.10128650732345633;
127     QUADNODES(1,2) = 0.7974269853530872;
128     QUADNODES(1,3) = 0.10128650732345633;
129     QUADNODES(1,4) = 0.47014206410511505;
130     QUADNODES(1,5) = 0.059715871789769809;
131     QUADNODES(1,6) = 0.47014206410511505;
132    
133     QUADWEIGHTS(0) = 0.225/2.;
134     QUADWEIGHTS(1) = 0.12593918054482717/2.;
135     QUADWEIGHTS(2) = 0.12593918054482717/2.;
136     QUADWEIGHTS(3) = 0.12593918054482717/2.;
137     QUADWEIGHTS(4) = 0.13239415278850616/2.;
138     QUADWEIGHTS(5) = 0.13239415278850616/2.;
139     QUADWEIGHTS(6) = 0.13239415278850616/2.;
140    
141     } else if (numQuadNodes==12){
142     a = 0.873821971016996;
143     b = 0.063089014491502;
144     c = 0.501426509658179;
145     d = 0.249286745170910;
146     e = 0.636502499121399;
147     f = 0.310352451033785;
148     g = 0.053145049844816;
149    
150     u = 0.050844906370207/2.;
151     v = 0.116786275726379/2.;
152     w = 0.082851075618374/2.;
153    
154     QUADNODES(0,0) = a;
155     QUADNODES(0,1) = b;
156     QUADNODES(0,2) = b;
157     QUADNODES(0,3) = c;
158     QUADNODES(0,4) = d;
159     QUADNODES(0,5) = d;
160     QUADNODES(0,6) = e;
161     QUADNODES(0,7) = e;
162     QUADNODES(0,8) = f;
163     QUADNODES(0,9) = f;
164     QUADNODES(0,10) = g;
165     QUADNODES(0,11) = g;
166    
167     QUADNODES(1,0) = b;
168     QUADNODES(1,1) = a;
169     QUADNODES(1,2) = b;
170     QUADNODES(1,3) = d;
171     QUADNODES(1,4) = c;
172     QUADNODES(1,5) = d;
173     QUADNODES(1,6) = f;
174     QUADNODES(1,7) = g;
175     QUADNODES(1,8) = e;
176     QUADNODES(1,9) = g;
177     QUADNODES(1,10) = e;
178     QUADNODES(1,11) = f;
179    
180     QUADWEIGHTS(0)= u;
181     QUADWEIGHTS(1)= u;
182     QUADWEIGHTS(2)= u;
183     QUADWEIGHTS(3)= v;
184     QUADWEIGHTS(4)= v;
185     QUADWEIGHTS(5)= v;
186     QUADWEIGHTS(6)= w;
187     QUADWEIGHTS(7)= w;
188     QUADWEIGHTS(8)= w;
189     QUADWEIGHTS(9)= w;
190     QUADWEIGHTS(10)= w;
191     QUADWEIGHTS(11)= w;
192    
193     } else if (numQuadNodes==13){
194     QUADWEIGHTS(0) =-0.149570044467670/2.;
195     QUADWEIGHTS(1) = 0.175615257433204/2.;
196     QUADWEIGHTS(2) = 0.175615257433204/2.;
197     QUADWEIGHTS(3) = 0.175615257433204/2.;
198     QUADWEIGHTS(4) = 0.053347235608839/2.;
199     QUADWEIGHTS(5) = 0.053347235608839/2.;
200     QUADWEIGHTS(6) = 0.053347235608839/2.;
201     QUADWEIGHTS(7) = 0.077113760890257/2.;
202     QUADWEIGHTS(8) = 0.077113760890257/2.;
203     QUADWEIGHTS(9) = 0.077113760890257/2.;
204     QUADWEIGHTS(10) = 0.077113760890257/2.;
205     QUADWEIGHTS(11) = 0.077113760890257/2.;
206     QUADWEIGHTS(12) = 0.077113760890257/2.;
207    
208     QUADNODES(0,0) = 0.3333333333333333;
209     QUADNODES(0,1) = 0.479308067841923;
210     QUADNODES(0,2) = 0.260345966079038;
211     QUADNODES(0,3) = 0.260345966079038;
212     QUADNODES(0,4) = 0.869739794195568;
213     QUADNODES(0,5) = 0.065130102902216;
214     QUADNODES(0,6) = 0.065130102902216;
215     QUADNODES(0,7) = 0.638444188569809;
216     QUADNODES(0,8) = 0.638444188569809;
217     QUADNODES(0,9) = 0.048690315425316;
218     QUADNODES(0,10) = 0.048690315425316;
219     QUADNODES(0,11) = 0.312865496004875;
220     QUADNODES(0,12) = 0.312865496004875;
221    
222     QUADNODES(1,0) = 0.3333333333333333;
223     QUADNODES(1,1) = 0.260345966079038;
224     QUADNODES(1,2) = 0.479308067841923;
225     QUADNODES(1,3) = 0.260345966079038;
226     QUADNODES(1,4) = 0.065130102902216;
227     QUADNODES(1,5) = 0.869739794195568;
228     QUADNODES(1,6) = 0.065130102902216;
229     QUADNODES(1,7) = 0.048690315425316;
230     QUADNODES(1,8) = 0.312865496004875;
231     QUADNODES(1,9) = 0.638444188569809;
232     QUADNODES(1,10) = 0.312865496004875;
233     QUADNODES(1,11) = 0.638444188569809;
234     QUADNODES(1,12) = 0.048690315425316;
235    
236     } else if (numQuadNodes==16){
237     QUADWEIGHTS(0) = 0.07215780;
238     QUADWEIGHTS(1) = 0.04754582;
239     QUADWEIGHTS(2) = 0.04754582;
240     QUADWEIGHTS(3) = 0.04754582;
241     QUADWEIGHTS(4) = 0.01622925;
242     QUADWEIGHTS(5) = 0.01622925;
243     QUADWEIGHTS(6) = 0.01622925;
244     QUADWEIGHTS(7) = 0.05160869;
245     QUADWEIGHTS(8) = 0.05160869;
246     QUADWEIGHTS(9) = 0.05160869;
247     QUADWEIGHTS(10) = 0.01361516;
248     QUADWEIGHTS(11) = 0.01361516;
249     QUADWEIGHTS(12) = 0.01361516;
250     QUADWEIGHTS(13) = 0.01361516;
251     QUADWEIGHTS(14) = 0.01361516;
252     QUADWEIGHTS(15) = 0.01361516;
253    
254     QUADNODES(0,0) = 0.3333333;
255     QUADNODES(0,1) = 0.08141482;
256     QUADNODES(0,2) = 0.4592926;
257     QUADNODES(0,3) = 0.4592926;
258     QUADNODES(0,4) = 0.8989055;
259     QUADNODES(0,5) = 0.05054723;
260     QUADNODES(0,6) = 0.05054723;
261     QUADNODES(0,7) = 0.6588614;
262     QUADNODES(0,8) = 0.1705693;
263     QUADNODES(0,9) = 0.1705693;
264     QUADNODES(0,10) = 0.008394777;
265     QUADNODES(0,11) = 0.008394777;
266     QUADNODES(0,12) = 0.7284924;
267     QUADNODES(0,13) = 0.7284924;
268     QUADNODES(0,14) = 0.2631128;
269     QUADNODES(0,15) = 0.2631128;
270    
271     QUADNODES(1,0) = 0.3333333;
272     QUADNODES(1,1) = 0.4592926;
273     QUADNODES(1,2) = 0.08141482;
274     QUADNODES(1,3) = 0.4592926;
275     QUADNODES(1,4) = 0.05054723;
276     QUADNODES(1,5) = 0.8989055;
277     QUADNODES(1,6) = 0.05054723;
278     QUADNODES(1,7) = 0.1705693;
279     QUADNODES(1,8) = 0.6588614;
280     QUADNODES(1,9) = 0.1705693;
281     QUADNODES(1,10) = 0.7284924;
282     QUADNODES(1,11) = 0.2631128;
283     QUADNODES(1,12) = 0.008394777;
284     QUADNODES(1,13) = 0.2631128;
285     QUADNODES(1,14) = 0.008394777;
286     QUADNODES(1,15) = 0.7284924;
287    
288     } else if (numQuadNodes==19){
289     QUADWEIGHTS(0) = 0.04856790;
290     QUADWEIGHTS(1) = 0.01566735;
291     QUADWEIGHTS(2) = 0.01566735;
292     QUADWEIGHTS(3) = 0.01566735;
293     QUADWEIGHTS(4) = 0.03891377;
294     QUADWEIGHTS(5) = 0.03891377;
295     QUADWEIGHTS(6) = 0.03891377;
296     QUADWEIGHTS(7) = 0.03982387;
297     QUADWEIGHTS(8) = 0.03982387;
298     QUADWEIGHTS(9) = 0.03982387;
299     QUADWEIGHTS(10) = 0.01278884;
300     QUADWEIGHTS(11) = 0.01278884;
301     QUADWEIGHTS(12) = 0.01278884;
302     QUADWEIGHTS(13) = 0.02164177;
303     QUADWEIGHTS(14) = 0.02164177;
304     QUADWEIGHTS(15) = 0.02164177;
305     QUADWEIGHTS(16) = 0.02164177;
306     QUADWEIGHTS(17) = 0.02164177;
307     QUADWEIGHTS(18) = 0.02164177;
308    
309     QUADNODES(0,0) = 0.3333333;
310     QUADNODES(0,1) = 0.02063496;
311     QUADNODES(0,2) = 0.4896825;
312     QUADNODES(0,3) = 0.4896825;
313     QUADNODES(0,4) = 0.1258208;
314     QUADNODES(0,5) = 0.4370896;
315     QUADNODES(0,6) = 0.4370896;
316     QUADNODES(0,7) = 0.6235929;
317     QUADNODES(0,8) = 0.1882035;
318     QUADNODES(0,9) = 0.1882035;
319     QUADNODES(0,10) = 0.9105410;
320     QUADNODES(0,11) = 0.04472951;
321     QUADNODES(0,12) = 0.04472951;
322     QUADNODES(0,13) = 0.03683841;
323     QUADNODES(0,14) = 0.03683841;
324     QUADNODES(0,15) = 0.7411986;
325     QUADNODES(0,16) = 0.7411986;
326     QUADNODES(0,17) = 0.2219630;
327     QUADNODES(0,18) = 0.2219630;
328    
329     QUADNODES(1,0) = 0.3333333;
330     QUADNODES(1,1) = 0.4896825;
331     QUADNODES(1,2) = 0.02063496;
332     QUADNODES(1,3) = 0.4896825;
333     QUADNODES(1,4) = 0.4370896;
334     QUADNODES(1,5) = 0.1258208;
335     QUADNODES(1,6) = 0.4370896;
336     QUADNODES(1,7) = 0.1882035;
337     QUADNODES(1,8) = 0.6235929;
338     QUADNODES(1,9) = 0.1882035;
339     QUADNODES(1,10) = 0.04472951;
340     QUADNODES(1,11) = 0.9105410;
341     QUADNODES(1,12) = 0.04472951;
342     QUADNODES(1,13) = 0.7411986;
343     QUADNODES(1,14) = 0.2219630;
344     QUADNODES(1,15) = 0.03683841;
345     QUADNODES(1,16) = 0.2219630;
346     QUADNODES(1,17) = 0.03683841;
347     QUADNODES(1,18) = 0.7411986;
348 jgs 82 } else {
349    
350     /* get scheme on [0.1]^2 */
351     Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
352 jgs 150 if (! Finley_noError()) return;
353 jgs 82
354     /* squeeze it: */
355    
356     for (i=0;i<numQuadNodes;i++) {
357     Q1=QUADNODES(0,i);
358     Q2=QUADNODES(1,i);
359     QUADWEIGHTS(i)=QUADWEIGHTS(i)*(1.-(1./2.)*(Q1+Q2));
360     QUADNODES(0,i)=Q1*(1.-(1./2.)*Q2);
361     QUADNODES(1,i)=Q2*(1.-(1./2.)*Q1);
362     }
363     }
364     #undef DIM
365 gross 1342
366    
367 jgs 82 }
368    
369     /**************************************************************/
370    
371     /* get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
372     /* as a queezed scheme on a hex [0,1]^3 */
373    
374     void Finley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
375     int i;
376     double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;
377 gross 1342 double a,b,c,d,e,f,g,h;
378     double alpha=0.58541019662496852;
379     double beta =0.1381966011250105;
380 jgs 82 #define DIM 3
381    
382 btully 1170 /* the easy cases: */
383 jgs 82 if (numQuadNodes==1) {
384 btully 1170 QUADNODES(0,0)=0.25;
385     QUADNODES(1,0)=0.25;
386     QUADNODES(2,0)=0.25;
387 jgs 82 QUADWEIGHTS(0)=1./6.;
388 btully 1170 } else if (numQuadNodes==4){
389     QUADNODES(0,0)=beta;
390     QUADNODES(1,0)=beta;
391     QUADNODES(2,0)=beta;
392     QUADWEIGHTS(0)=1./24.;
393     QUADNODES(0,1)=alpha;
394     QUADNODES(1,1)=beta;
395     QUADNODES(2,1)=beta;
396     QUADWEIGHTS(1)=1./24.;
397     QUADNODES(0,2)=beta;
398     QUADNODES(1,2)=alpha;
399     QUADNODES(2,2)=beta;
400     QUADWEIGHTS(2)=1./24.;
401     QUADNODES(0,3)=beta;
402     QUADNODES(1,3)=beta;
403     QUADNODES(2,3)=alpha;
404     QUADWEIGHTS(3)=1./24.;
405     } else if (numQuadNodes==5){
406     QUADNODES(0,0)=1./4.;
407     QUADNODES(1,0)=1./4.;
408     QUADNODES(2,0)=1./4.;
409     QUADWEIGHTS(0)=-2./15.;
410     QUADNODES(0,1)=1./6.;
411     QUADNODES(1,1)=1./6.;
412     QUADNODES(2,1)=1./6.;
413     QUADWEIGHTS(1)=3./40.;
414     QUADNODES(0,2)=1./2.;
415     QUADNODES(1,2)=1./6.;
416     QUADNODES(2,2)=1./6.;
417     QUADWEIGHTS(2)=3./40.;
418     QUADNODES(0,3)=1./6.;
419     QUADNODES(1,3)=1./2.;
420     QUADNODES(2,3)=1./6.;
421     QUADWEIGHTS(3)=3./40.;
422     QUADNODES(0,4)=1./6.;
423     QUADNODES(1,4)=1./6.;
424     QUADNODES(2,4)=1./2.;
425     QUADWEIGHTS(4)=3./40.;
426 gross 1342
427     } else if (numQuadNodes==11){
428    
429     a = 0.25;
430     b = 11.0/14.0;
431     c = 1.0/14.0;
432     d = 0.25 * (1.0 + sqrt ( 5.0 / 14.0 ) );
433     e = 0.25 * (1.0 - sqrt ( 5.0 / 14.0 ) );
434     f = -74.0 / 5625.0;
435     g = 343.0 / 45000.0;
436     h = 56.0 / 2250.0;
437    
438     QUADWEIGHTS(401-401) = f;
439     QUADWEIGHTS(402-401) = g;
440     QUADWEIGHTS(403-401) = g;
441     QUADWEIGHTS(404-401) = g;
442     QUADWEIGHTS(405-401) = g;
443     QUADWEIGHTS(406-401) = h;
444     QUADWEIGHTS(407-401) = h;
445     QUADWEIGHTS(408-401) = h;
446     QUADWEIGHTS(409-401) = h;
447     QUADWEIGHTS(410-401) = h;
448     QUADWEIGHTS(411-401) = h;
449    
450     QUADNODES(0,401-401) = a;
451     QUADNODES(0,402-401) = b;
452     QUADNODES(0,403-401) = c;
453     QUADNODES(0,404-401) = c;
454     QUADNODES(0,405-401) = c;
455     QUADNODES(0,406-401) = d;
456     QUADNODES(0,407-401) = d;
457     QUADNODES(0,408-401) = d;
458     QUADNODES(0,409-401) = e;
459     QUADNODES(0,410-401) = e;
460     QUADNODES(0,411-401) = e;
461    
462     QUADNODES(1,401-401) = a;
463     QUADNODES(1,402-401) = c;
464     QUADNODES(1,403-401) = b;
465     QUADNODES(1,404-401) = c;
466     QUADNODES(1,405-401) = c;
467     QUADNODES(1,406-401) = d;
468     QUADNODES(1,407-401) = e;
469     QUADNODES(1,408-401) = e;
470     QUADNODES(1,409-401) = d;
471     QUADNODES(1,410-401) = d;
472     QUADNODES(1,411-401) = e;
473    
474     QUADNODES(2,401-401) = a;
475     QUADNODES(2,402-401) = c;
476     QUADNODES(2,403-401) = c;
477     QUADNODES(2,404-401) = b;
478     QUADNODES(2,405-401) = c;
479     QUADNODES(2,406-401) = e;
480     QUADNODES(2,407-401) = d;
481     QUADNODES(2,408-401) = e;
482     QUADNODES(2,409-401) = d;
483     QUADNODES(2,410-401) = e;
484     QUADNODES(2,411-401) = d;
485    
486     } else if (numQuadNodes==15){
487     QUADWEIGHTS(412-412) = 0.019753086419753086;
488     QUADWEIGHTS(413-412) = 0.01198951396316977;
489     QUADWEIGHTS(414-412) = 0.01198951396316977;
490     QUADWEIGHTS(415-412) = 0.01198951396316977;
491     QUADWEIGHTS(416-412) = 0.01198951396316977;
492     QUADWEIGHTS(417-412) = 0.011511367871045397;
493     QUADWEIGHTS(418-412) = 0.011511367871045397;
494     QUADWEIGHTS(419-412) = 0.011511367871045397;
495     QUADWEIGHTS(420-412) = 0.011511367871045397;
496     QUADWEIGHTS(421-412) = 0.0088183421516754845;
497     QUADWEIGHTS(422-412) = 0.0088183421516754845;
498     QUADWEIGHTS(423-412) = 0.0088183421516754845;
499     QUADWEIGHTS(424-412) = 0.0088183421516754845;
500     QUADWEIGHTS(425-412) = 0.0088183421516754845;
501     QUADWEIGHTS(426-412) = 0.0088183421516754845;
502    
503     QUADNODES(0,412-412) = 0.2500000;
504     QUADNODES(0,413-412) = 0.091971078052723032;
505     QUADNODES(0,414-412) = 0.72408676584183096;
506     QUADNODES(0,415-412) = 0.091971078052723032;
507     QUADNODES(0,416-412) = 0.091971078052723032;
508     QUADNODES(0,417-412) = 0.31979362782962989;
509     QUADNODES(0,418-412) = 0.040619116511110234;
510     QUADNODES(0,419-412) = 0.31979362782962989;
511     QUADNODES(0,420-412) = 0.31979362782962989;
512     QUADNODES(0,421-412) = 0.056350832689629149;
513     QUADNODES(0,422-412) = 0.056350832689629149;
514     QUADNODES(0,423-412) = 0.056350832689629149;
515     QUADNODES(0,424-412) = 0.4436491673103708;
516     QUADNODES(0,425-412) = 0.4436491673103708;
517     QUADNODES(0,426-412) = 0.4436491673103708;
518    
519     QUADNODES(1,412-412) = 0.2500000;
520     QUADNODES(1,413-412) = 0.091971078052723032;
521     QUADNODES(1,414-412) = 0.091971078052723032;
522     QUADNODES(1,415-412) = 0.72408676584183096;
523     QUADNODES(1,416-412) = 0.091971078052723032;
524     QUADNODES(1,417-412) = 0.31979362782962989;
525     QUADNODES(1,418-412) = 0.31979362782962989;
526     QUADNODES(1,419-412) = 0.040619116511110234;
527     QUADNODES(1,420-412) = 0.31979362782962989;
528     QUADNODES(1,421-412) = 0.056350832689629149;
529     QUADNODES(1,422-412) = 0.4436491673103708;
530     QUADNODES(1,423-412) = 0.4436491673103708;
531     QUADNODES(1,424-412) = 0.056350832689629149;
532     QUADNODES(1,425-412) = 0.056350832689629149;
533     QUADNODES(1,426-412) = 0.4436491673103708;
534    
535     QUADNODES(2,412-412) = 0.2500000;
536     QUADNODES(2,413-412) = 0.091971078052723032;
537     QUADNODES(2,414-412) = 0.091971078052723032;
538     QUADNODES(2,415-412) = 0.091971078052723032;
539     QUADNODES(2,416-412) = 0.72408676584183096;
540     QUADNODES(2,417-412) = 0.31979362782962989;
541     QUADNODES(2,418-412) = 0.31979362782962989;
542     QUADNODES(2,419-412) = 0.31979362782962989;
543     QUADNODES(2,420-412) = 0.040619116511110234;
544     QUADNODES(2,421-412) = 0.4436491673103708;
545     QUADNODES(2,422-412) = 0.056350832689629149;
546     QUADNODES(2,423-412) = 0.4436491673103708;
547     QUADNODES(2,424-412) = 0.056350832689629149;
548     QUADNODES(2,425-412) = 0.4436491673103708;
549     QUADNODES(2,426-412) = 0.056350832689629149;
550    
551     } else if (numQuadNodes==24){
552     QUADWEIGHTS(427-427) = 0.006653792;
553     QUADWEIGHTS(428-427) = 0.006653792;
554     QUADWEIGHTS(429-427) = 0.006653792;
555     QUADWEIGHTS(430-427) = 0.006653792;
556     QUADWEIGHTS(431-427) = 0.001679535;
557     QUADWEIGHTS(432-427) = 0.001679535;
558     QUADWEIGHTS(433-427) = 0.001679535;
559     QUADWEIGHTS(434-427) = 0.001679535;
560     QUADWEIGHTS(435-427) = 0.009226197;
561     QUADWEIGHTS(436-427) = 0.009226197;
562     QUADWEIGHTS(437-427) = 0.009226197;
563     QUADWEIGHTS(438-427) = 0.009226197;
564     QUADWEIGHTS(439-427) = 0.008035714;
565     QUADWEIGHTS(440-427) = 0.008035714;
566     QUADWEIGHTS(441-427) = 0.008035714;
567     QUADWEIGHTS(442-427) = 0.008035714;
568     QUADWEIGHTS(443-427) = 0.008035714;
569     QUADWEIGHTS(444-427) = 0.008035714;
570     QUADWEIGHTS(445-427) = 0.008035714;
571     QUADWEIGHTS(446-427) = 0.008035714;
572     QUADWEIGHTS(447-427) = 0.008035714;
573     QUADWEIGHTS(448-427) = 0.008035714;
574     QUADWEIGHTS(449-427) = 0.008035714;
575     QUADWEIGHTS(450-427) = 0.008035714;
576    
577     QUADNODES(0,427-427) = 0.3561914;
578     QUADNODES(0,428-427) = 0.2146029;
579     QUADNODES(0,429-427) = 0.2146029;
580     QUADNODES(0,430-427) = 0.2146029;
581     QUADNODES(0,431-427) = 0.8779781;
582     QUADNODES(0,432-427) = 0.04067396;
583     QUADNODES(0,433-427) = 0.04067396;
584     QUADNODES(0,434-427) = 0.04067396;
585     QUADNODES(0,435-427) = 0.03298633;
586     QUADNODES(0,436-427) = 0.3223379;
587     QUADNODES(0,437-427) = 0.3223379;
588     QUADNODES(0,438-427) = 0.3223379;
589     QUADNODES(0,439-427) = 0.6030057;
590     QUADNODES(0,440-427) = 0.6030057;
591     QUADNODES(0,441-427) = 0.6030057;
592     QUADNODES(0,442-427) = 0.2696723;
593     QUADNODES(0,443-427) = 0.2696723;
594     QUADNODES(0,444-427) = 0.2696723;
595     QUADNODES(0,445-427) = 0.06366100;
596     QUADNODES(0,446-427) = 0.06366100;
597     QUADNODES(0,447-427) = 0.06366100;
598     QUADNODES(0,448-427) = 0.06366100;
599     QUADNODES(0,449-427) = 0.06366100;
600     QUADNODES(0,450-427) = 0.06366100;
601    
602     QUADNODES(1,427-427) = 0.2146029;
603     QUADNODES(1,428-427) = 0.3561914;
604     QUADNODES(1,429-427) = 0.2146029;
605     QUADNODES(1,430-427) = 0.2146029;
606     QUADNODES(1,431-427) = 0.04067396;
607     QUADNODES(1,432-427) = 0.8779781;
608     QUADNODES(1,433-427) = 0.04067396;
609     QUADNODES(1,434-427) = 0.04067396;
610     QUADNODES(1,435-427) = 0.3223379;
611     QUADNODES(1,436-427) = 0.03298633;
612     QUADNODES(1,437-427) = 0.3223379;
613     QUADNODES(1,438-427) = 0.3223379;
614     QUADNODES(1,439-427) = 0.2696723;
615     QUADNODES(1,440-427) = 0.06366100;
616     QUADNODES(1,441-427) = 0.06366100;
617     QUADNODES(1,442-427) = 0.6030057;
618     QUADNODES(1,443-427) = 0.06366100;
619     QUADNODES(1,444-427) = 0.06366100;
620     QUADNODES(1,445-427) = 0.6030057;
621     QUADNODES(1,446-427) = 0.6030057;
622     QUADNODES(1,447-427) = 0.2696723;
623     QUADNODES(1,448-427) = 0.2696723;
624     QUADNODES(1,449-427) = 0.06366100;
625     QUADNODES(1,450-427) = 0.06366100;
626    
627     QUADNODES(2,427-427) = 0.2146029;
628     QUADNODES(2,428-427) = 0.2146029;
629     QUADNODES(2,429-427) = 0.3561914;
630     QUADNODES(2,430-427) = 0.2146029;
631     QUADNODES(2,431-427) = 0.04067396;
632     QUADNODES(2,432-427) = 0.04067396;
633     QUADNODES(2,433-427) = 0.8779781;
634     QUADNODES(2,434-427) = 0.04067396;
635     QUADNODES(2,435-427) = 0.3223379;
636     QUADNODES(2,436-427) = 0.3223379;
637     QUADNODES(2,437-427) = 0.03298633;
638     QUADNODES(2,438-427) = 0.3223379;
639     QUADNODES(2,439-427) = 0.06366100;
640     QUADNODES(2,440-427) = 0.2696723;
641     QUADNODES(2,441-427) = 0.06366100;
642     QUADNODES(2,442-427) = 0.06366100;
643     QUADNODES(2,443-427) = 0.6030057;
644     QUADNODES(2,444-427) = 0.06366100;
645     QUADNODES(2,445-427) = 0.2696723;
646     QUADNODES(2,446-427) = 0.06366100;
647     QUADNODES(2,447-427) = 0.6030057;
648     QUADNODES(2,448-427) = 0.06366100;
649     QUADNODES(2,449-427) = 0.6030057;
650     QUADNODES(2,450-427) = 0.2696723;
651    
652     } else if (numQuadNodes==31){
653     QUADWEIGHTS(451-451) = 0.01826422;
654     QUADWEIGHTS(452-451) = 0.01059994;
655     QUADWEIGHTS(453-451) = 0.01059994;
656     QUADWEIGHTS(454-451) = 0.01059994;
657     QUADWEIGHTS(455-451) = 0.01059994;
658     QUADWEIGHTS(456-451) =-0.06251774;
659     QUADWEIGHTS(457-451) =-0.06251774;
660     QUADWEIGHTS(458-451) =-0.06251774;
661     QUADWEIGHTS(459-451) =-0.06251774;
662     QUADWEIGHTS(460-451) = 0.004891425;
663     QUADWEIGHTS(461-451) = 0.004891425;
664     QUADWEIGHTS(462-451) = 0.004891425;
665     QUADWEIGHTS(463-451) = 0.004891425;
666     QUADWEIGHTS(464-451) = 0.0009700176;
667     QUADWEIGHTS(465-451) = 0.0009700176;
668     QUADWEIGHTS(466-451) = 0.0009700176;
669     QUADWEIGHTS(467-451) = 0.0009700176;
670     QUADWEIGHTS(468-451) = 0.0009700176;
671     QUADWEIGHTS(469-451) = 0.0009700176;
672     QUADWEIGHTS(470-451) = 0.02755732;
673     QUADWEIGHTS(471-451) = 0.02755732;
674     QUADWEIGHTS(472-451) = 0.02755732;
675     QUADWEIGHTS(473-451) = 0.02755732;
676     QUADWEIGHTS(474-451) = 0.02755732;
677     QUADWEIGHTS(475-451) = 0.02755732;
678     QUADWEIGHTS(476-451) = 0.02755732;
679     QUADWEIGHTS(477-451) = 0.02755732;
680     QUADWEIGHTS(478-451) = 0.02755732;
681     QUADWEIGHTS(479-451) = 0.02755732;
682     QUADWEIGHTS(480-451) = 0.02755732;
683     QUADWEIGHTS(481-451) = 0.02755732;
684    
685     QUADNODES(0,451-451) = 0.2500000;
686     QUADNODES(0,452-451) = 0.7653604;
687     QUADNODES(0,453-451) = 0.07821319;
688     QUADNODES(0,454-451) = 0.07821319;
689     QUADNODES(0,455-451) = 0.07821319;
690     QUADNODES(0,456-451) = 0.6344704;
691     QUADNODES(0,457-451) = 0.1218432;
692     QUADNODES(0,458-451) = 0.1218432;
693     QUADNODES(0,459-451) = 0.1218432;
694     QUADNODES(0,460-451) = 0.002382507;
695     QUADNODES(0,461-451) = 0.3325392;
696     QUADNODES(0,462-451) = 0.3325392;
697     QUADNODES(0,463-451) = 0.3325392;
698     QUADNODES(0,464-451) = 0.0000000;
699     QUADNODES(0,465-451) = 0.0000000;
700     QUADNODES(0,466-451) = 0.0000000;
701     QUADNODES(0,467-451) = 0.5000000;
702     QUADNODES(0,468-451) = 0.5000000;
703     QUADNODES(0,469-451) = 0.5000000;
704     QUADNODES(0,470-451) = 0.6000000;
705     QUADNODES(0,471-451) = 0.6000000;
706     QUADNODES(0,472-451) = 0.6000000;
707     QUADNODES(0,473-451) = 0.2000000;
708     QUADNODES(0,474-451) = 0.2000000;
709     QUADNODES(0,475-451) = 0.2000000;
710     QUADNODES(0,476-451) = 0.1000000;
711     QUADNODES(0,477-451) = 0.1000000;
712     QUADNODES(0,478-451) = 0.1000000;
713     QUADNODES(0,479-451) = 0.1000000;
714     QUADNODES(0,480-451) = 0.1000000;
715     QUADNODES(0,481-451) = 0.1000000;
716    
717     QUADNODES(1,451-451) = 0.2500000;
718     QUADNODES(1,452-451) = 0.07821319;
719     QUADNODES(1,453-451) = 0.7653604;
720     QUADNODES(1,454-451) = 0.07821319;
721     QUADNODES(1,455-451) = 0.07821319;
722     QUADNODES(1,456-451) = 0.1218432;
723     QUADNODES(1,457-451) = 0.6344704;
724     QUADNODES(1,458-451) = 0.1218432;
725     QUADNODES(1,459-451) = 0.1218432;
726     QUADNODES(1,460-451) = 0.3325392;
727     QUADNODES(1,461-451) = 0.002382507;
728     QUADNODES(1,462-451) = 0.3325392;
729     QUADNODES(1,463-451) = 0.3325392;
730     QUADNODES(1,464-451) = 0.0000000;
731     QUADNODES(1,465-451) = 0.5000000;
732     QUADNODES(1,466-451) = 0.5000000;
733     QUADNODES(1,467-451) = 0.0000000;
734     QUADNODES(1,468-451) = 0.0000000;
735     QUADNODES(1,469-451) = 0.5000000;
736     QUADNODES(1,470-451) = 0.2000000;
737     QUADNODES(1,471-451) = 0.1000000;
738     QUADNODES(1,472-451) = 0.1000000;
739     QUADNODES(1,473-451) = 0.6000000;
740     QUADNODES(1,474-451) = 0.1000000;
741     QUADNODES(1,475-451) = 0.1000000;
742     QUADNODES(1,476-451) = 0.6000000;
743     QUADNODES(1,477-451) = 0.6000000;
744     QUADNODES(1,478-451) = 0.2000000;
745     QUADNODES(1,479-451) = 0.2000000;
746     QUADNODES(1,480-451) = 0.1000000;
747     QUADNODES(1,481-451) = 0.1000000;
748    
749     QUADNODES(2,451-451) = 0.2500000;
750     QUADNODES(2,452-451) = 0.07821319;
751     QUADNODES(2,453-451) = 0.07821319;
752     QUADNODES(2,454-451) = 0.7653604;
753     QUADNODES(2,455-451) = 0.07821319;
754     QUADNODES(2,456-451) = 0.1218432;
755     QUADNODES(2,457-451) = 0.1218432;
756     QUADNODES(2,458-451) = 0.6344704;
757     QUADNODES(2,459-451) = 0.1218432;
758     QUADNODES(2,460-451) = 0.3325392;
759     QUADNODES(2,461-451) = 0.3325392;
760     QUADNODES(2,462-451) = 0.002382507;
761     QUADNODES(2,463-451) = 0.3325392;
762     QUADNODES(2,464-451) = 0.5000000;
763     QUADNODES(2,465-451) = 0.0000000;
764     QUADNODES(2,466-451) = 0.5000000;
765     QUADNODES(2,467-451) = 0.0000000;
766     QUADNODES(2,468-451) = 0.5000000;
767     QUADNODES(2,469-451) = 0.0000000;
768     QUADNODES(2,470-451) = 0.1000000;
769     QUADNODES(2,471-451) = 0.2000000;
770     QUADNODES(2,472-451) = 0.1000000;
771     QUADNODES(2,473-451) = 0.1000000;
772     QUADNODES(2,474-451) = 0.6000000;
773     QUADNODES(2,475-451) = 0.1000000;
774     QUADNODES(2,476-451) = 0.2000000;
775     QUADNODES(2,477-451) = 0.1000000;
776     QUADNODES(2,478-451) = 0.6000000;
777     QUADNODES(2,479-451) = 0.1000000;
778     QUADNODES(2,480-451) = 0.6000000;
779     QUADNODES(2,481-451) = 0.2000000;
780    
781     } else if (numQuadNodes==45){
782     QUADWEIGHTS(482-482) =-0.03932701;
783     QUADWEIGHTS(483-482) = 0.004081316;
784     QUADWEIGHTS(484-482) = 0.004081316;
785     QUADWEIGHTS(485-482) = 0.004081316;
786     QUADWEIGHTS(486-482) = 0.004081316;
787     QUADWEIGHTS(487-482) = 0.0006580868;
788     QUADWEIGHTS(488-482) = 0.0006580868;
789     QUADWEIGHTS(489-482) = 0.0006580868;
790     QUADWEIGHTS(490-482) = 0.0006580868;
791     QUADWEIGHTS(491-482) = 0.004384259;
792     QUADWEIGHTS(492-482) = 0.004384259;
793     QUADWEIGHTS(493-482) = 0.004384259;
794     QUADWEIGHTS(494-482) = 0.004384259;
795     QUADWEIGHTS(495-482) = 0.004384259;
796     QUADWEIGHTS(496-482) = 0.004384259;
797     QUADWEIGHTS(497-482) = 0.01383006;
798     QUADWEIGHTS(498-482) = 0.01383006;
799     QUADWEIGHTS(499-482) = 0.01383006;
800     QUADWEIGHTS(500-482) = 0.01383006;
801     QUADWEIGHTS(501-482) = 0.01383006;
802     QUADWEIGHTS(502-482) = 0.01383006;
803     QUADWEIGHTS(503-482) = 0.004240437;
804     QUADWEIGHTS(504-482) = 0.004240437;
805     QUADWEIGHTS(505-482) = 0.004240437;
806     QUADWEIGHTS(506-482) = 0.004240437;
807     QUADWEIGHTS(507-482) = 0.004240437;
808     QUADWEIGHTS(508-482) = 0.004240437;
809     QUADWEIGHTS(509-482) = 0.004240437;
810     QUADWEIGHTS(510-482) = 0.004240437;
811     QUADWEIGHTS(511-482) = 0.004240437;
812     QUADWEIGHTS(512-482) = 0.004240437;
813     QUADWEIGHTS(513-482) = 0.004240437;
814     QUADWEIGHTS(514-482) = 0.004240437;
815     QUADWEIGHTS(515-482) = 0.002238740;
816     QUADWEIGHTS(516-482) = 0.002238740;
817     QUADWEIGHTS(517-482) = 0.002238740;
818     QUADWEIGHTS(518-482) = 0.002238740;
819     QUADWEIGHTS(519-482) = 0.002238740;
820     QUADWEIGHTS(520-482) = 0.002238740;
821     QUADWEIGHTS(521-482) = 0.002238740;
822     QUADWEIGHTS(522-482) = 0.002238740;
823     QUADWEIGHTS(523-482) = 0.002238740;
824     QUADWEIGHTS(524-482) = 0.002238740;
825     QUADWEIGHTS(525-482) = 0.002238740;
826     QUADWEIGHTS(526-482) = 0.002238740;
827    
828     QUADNODES(0,482-482) = 0.2500000;
829     QUADNODES(0,483-482) = 0.6175872;
830     QUADNODES(0,484-482) = 0.1274709;
831     QUADNODES(0,485-482) = 0.1274709;
832     QUADNODES(0,486-482) = 0.1274709;
833     QUADNODES(0,487-482) = 0.9037635;
834     QUADNODES(0,488-482) = 0.03207883;
835     QUADNODES(0,489-482) = 0.03207883;
836     QUADNODES(0,490-482) = 0.03207883;
837     QUADNODES(0,491-482) = 0.4502229;
838     QUADNODES(0,492-482) = 0.4502229;
839     QUADNODES(0,493-482) = 0.4502229;
840     QUADNODES(0,494-482) = 0.04977710;
841     QUADNODES(0,495-482) = 0.04977710;
842     QUADNODES(0,496-482) = 0.04977710;
843     QUADNODES(0,497-482) = 0.3162696;
844     QUADNODES(0,498-482) = 0.3162696;
845     QUADNODES(0,499-482) = 0.3162696;
846     QUADNODES(0,500-482) = 0.1837304;
847     QUADNODES(0,501-482) = 0.1837304;
848     QUADNODES(0,502-482) = 0.1837304;
849     QUADNODES(0,503-482) = 0.5132800;
850     QUADNODES(0,504-482) = 0.5132800;
851     QUADNODES(0,505-482) = 0.5132800;
852     QUADNODES(0,506-482) = 0.02291779;
853     QUADNODES(0,507-482) = 0.02291779;
854     QUADNODES(0,508-482) = 0.02291779;
855     QUADNODES(0,509-482) = 0.2319011;
856     QUADNODES(0,510-482) = 0.2319011;
857     QUADNODES(0,511-482) = 0.2319011;
858     QUADNODES(0,512-482) = 0.2319011;
859     QUADNODES(0,513-482) = 0.2319011;
860     QUADNODES(0,514-482) = 0.2319011;
861     QUADNODES(0,515-482) = 0.1937465;
862     QUADNODES(0,516-482) = 0.1937465;
863     QUADNODES(0,517-482) = 0.1937465;
864     QUADNODES(0,518-482) = 0.7303134;
865     QUADNODES(0,519-482) = 0.7303134;
866     QUADNODES(0,520-482) = 0.7303134;
867     QUADNODES(0,521-482) = 0.03797005;
868     QUADNODES(0,522-482) = 0.03797005;
869     QUADNODES(0,523-482) = 0.03797005;
870     QUADNODES(0,524-482) = 0.03797005;
871     QUADNODES(0,525-482) = 0.03797005;
872     QUADNODES(0,526-482) = 0.03797005;
873    
874     QUADNODES(1,482-482) = 0.2500000;
875     QUADNODES(1,483-482) = 0.1274709;
876     QUADNODES(1,484-482) = 0.6175872;
877     QUADNODES(1,485-482) = 0.1274709;
878     QUADNODES(1,486-482) = 0.1274709;
879     QUADNODES(1,487-482) = 0.03207883;
880     QUADNODES(1,488-482) = 0.9037635;
881     QUADNODES(1,489-482) = 0.03207883;
882     QUADNODES(1,490-482) = 0.03207883;
883     QUADNODES(1,491-482) = 0.4502229;
884     QUADNODES(1,492-482) = 0.04977710;
885     QUADNODES(1,493-482) = 0.04977710;
886     QUADNODES(1,494-482) = 0.4502229;
887     QUADNODES(1,495-482) = 0.4502229;
888     QUADNODES(1,496-482) = 0.04977710;
889     QUADNODES(1,497-482) = 0.3162696;
890     QUADNODES(1,498-482) = 0.1837304;
891     QUADNODES(1,499-482) = 0.1837304;
892     QUADNODES(1,500-482) = 0.3162696;
893     QUADNODES(1,501-482) = 0.3162696;
894     QUADNODES(1,502-482) = 0.1837304;
895     QUADNODES(1,503-482) = 0.02291779;
896     QUADNODES(1,504-482) = 0.2319011;
897     QUADNODES(1,505-482) = 0.2319011;
898     QUADNODES(1,506-482) = 0.5132800;
899     QUADNODES(1,507-482) = 0.2319011;
900     QUADNODES(1,508-482) = 0.2319011;
901     QUADNODES(1,509-482) = 0.5132800;
902     QUADNODES(1,510-482) = 0.5132800;
903     QUADNODES(1,511-482) = 0.02291779;
904     QUADNODES(1,512-482) = 0.02291779;
905     QUADNODES(1,513-482) = 0.2319011;
906     QUADNODES(1,514-482) = 0.2319011;
907     QUADNODES(1,515-482) = 0.7303134;
908     QUADNODES(1,516-482) = 0.03797005;
909     QUADNODES(1,517-482) = 0.03797005;
910     QUADNODES(1,518-482) = 0.1937465;
911     QUADNODES(1,519-482) = 0.03797005;
912     QUADNODES(1,520-482) = 0.03797005;
913     QUADNODES(1,521-482) = 0.1937465;
914     QUADNODES(1,522-482) = 0.1937465;
915     QUADNODES(1,523-482) = 0.7303134;
916     QUADNODES(1,524-482) = 0.7303134;
917     QUADNODES(1,525-482) = 0.03797005;
918     QUADNODES(1,526-482) = 0.03797005;
919    
920     QUADNODES(2,482-482) = 0.2500000;
921     QUADNODES(2,483-482) = 0.1274709;
922     QUADNODES(2,484-482) = 0.1274709;
923     QUADNODES(2,485-482) = 0.6175872;
924     QUADNODES(2,486-482) = 0.1274709;
925     QUADNODES(2,487-482) = 0.03207883;
926     QUADNODES(2,488-482) = 0.03207883;
927     QUADNODES(2,489-482) = 0.9037635;
928     QUADNODES(2,490-482) = 0.03207883;
929     QUADNODES(2,491-482) = 0.04977710;
930     QUADNODES(2,492-482) = 0.4502229;
931     QUADNODES(2,493-482) = 0.04977710;
932     QUADNODES(2,494-482) = 0.4502229;
933     QUADNODES(2,495-482) = 0.04977710;
934     QUADNODES(2,496-482) = 0.4502229;
935     QUADNODES(2,497-482) = 0.1837304;
936     QUADNODES(2,498-482) = 0.3162696;
937     QUADNODES(2,499-482) = 0.1837304;
938     QUADNODES(2,500-482) = 0.3162696;
939     QUADNODES(2,501-482) = 0.1837304;
940     QUADNODES(2,502-482) = 0.3162696;
941     QUADNODES(2,503-482) = 0.2319011;
942     QUADNODES(2,504-482) = 0.02291779;
943     QUADNODES(2,505-482) = 0.2319011;
944     QUADNODES(2,506-482) = 0.2319011;
945     QUADNODES(2,507-482) = 0.5132800;
946     QUADNODES(2,508-482) = 0.2319011;
947     QUADNODES(2,509-482) = 0.02291779;
948     QUADNODES(2,510-482) = 0.2319011;
949     QUADNODES(2,511-482) = 0.5132800;
950     QUADNODES(2,512-482) = 0.2319011;
951     QUADNODES(2,513-482) = 0.5132800;
952     QUADNODES(2,514-482) = 0.02291779;
953     QUADNODES(2,515-482) = 0.03797005;
954     QUADNODES(2,516-482) = 0.7303134;
955     QUADNODES(2,517-482) = 0.03797005;
956     QUADNODES(2,518-482) = 0.03797005;
957     QUADNODES(2,519-482) = 0.1937465;
958     QUADNODES(2,520-482) = 0.03797005;
959     QUADNODES(2,521-482) = 0.7303134;
960     QUADNODES(2,522-482) = 0.03797005;
961     QUADNODES(2,523-482) = 0.1937465;
962     QUADNODES(2,524-482) = 0.03797005;
963     QUADNODES(2,525-482) = 0.1937465;
964     QUADNODES(2,526-482) = 0.7303134;
965    
966 jgs 82 } else {
967    
968     /* get scheme on [0.1]^3 */
969    
970     Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
971 jgs 150 if (! Finley_noError()) return;
972 jgs 82
973     /* squeeze it: */
974     for (i=0;i<numQuadNodes;i++) {
975     Q1=QUADNODES(0,i);
976     Q2=QUADNODES(1,i);
977     Q3=QUADNODES(2,i);
978    
979     JA11= (1./3.)*Q2*Q3-(1./2.)*(Q2+Q3) +1.;
980     JA12= (1./3.)*Q1*Q3-(1./2.)*Q1;
981     JA13= (1./3.)*Q1*Q2-(1./2.)*Q1;
982     JA21= (1./3.)*Q2*Q3-(1./2.)*Q2;
983     JA22= (1./3.)*Q1*Q3-(1./2.)*(Q1+Q3) +1.;
984     JA23= (1./3.)*Q1*Q2-(1./2.)*Q2;
985     JA31= (1./3.)*Q2*Q3-(1./2.)*Q3;
986     JA32= (1./3.)*Q1*Q3-(1./2.)*Q3;
987     JA33= (1./3.)*Q1*Q2-(1./2.)*(Q1+Q2) +1.;
988     DET=JA11*JA22*JA33+JA12*JA23*JA31+JA13*JA21*JA32-JA13*JA22*JA31-JA11*JA23*JA32-JA12*JA21*JA33;
989     quadWeights[i]=quadWeights[i]*ABS(DET);
990     QUADNODES(0,i)=Q1*((1./3.)*Q2*Q3-(1./2.)*(Q2+Q3)+1.);
991     QUADNODES(1,i)=Q2*((1./3.)*Q1*Q3-(1./2.)*(Q1+Q3)+1.);
992     QUADNODES(2,i)=Q3*((1./3.)*Q1*Q2-(1./2.)*(Q1+Q2)+1.);
993     }
994     }
995     #undef DIM
996     }
997    
998     /**************************************************************/
999    
1000     /* get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */
1001     /* as a X-product of a 1D scheme. */
1002    
1003     void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
1004 jgs 150 char error_msg[LenErrorMsg_MAX];
1005 jgs 82 int numQuadNodes1d,i,j,l;
1006 gross 1028 double *quadNodes1d=NULL,*quadWeights1d=NULL;
1007     bool_t set=FALSE;
1008 jgs 82 #define DIM 2
1009    
1010 gross 1028 quadNodes1d=TMPMEMALLOC(numQuadNodes,double);
1011     quadWeights1d=TMPMEMALLOC(numQuadNodes,double);
1012     if (! ( Finley_checkPtr(quadNodes1d) || Finley_checkPtr(quadWeights1d) ) ) {
1013     /* find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */
1014    
1015     for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
1016     if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
1017 jgs 82
1018 gross 1028 /* get 1D scheme: */
1019    
1020     Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
1021 jgs 82
1022 gross 1028 /* make 2D scheme: */
1023 jgs 82
1024 gross 1028 if (Finley_noError()) {
1025     l=0;
1026     for (i=0;i<numQuadNodes1d;i++) {
1027     for (j=0;j<numQuadNodes1d;j++) {
1028     QUADNODES(0,l)=quadNodes1d[i];
1029     QUADNODES(1,l)=quadNodes1d[j];
1030     QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];
1031     l++;
1032     }
1033     }
1034     set=TRUE;
1035     break;
1036     }
1037     }
1038     }
1039     if (!set) {
1040     sprintf(error_msg,"Finley_Quad_getNodesRec: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
1041     Finley_setError(VALUE_ERROR,error_msg);
1042     }
1043     TMPMEMFREE(quadNodes1d);
1044     TMPMEMFREE(quadWeights1d);
1045     }
1046     #undef DIM
1047 jgs 82 }
1048    
1049     /**************************************************************/
1050    
1051     /* get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
1052     /* as a X-product of a 1D scheme. */
1053    
1054     void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
1055 jgs 150 char error_msg[LenErrorMsg_MAX];
1056 jgs 82 int numQuadNodes1d,i,j,k,l;
1057 gross 1028 double *quadNodes1d=NULL,*quadWeights1d=NULL;
1058 artak 2034 bool_t set=FALSE;
1059 jgs 82 #define DIM 3
1060    
1061     /* find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
1062    
1063 gross 1028 quadNodes1d=TMPMEMALLOC(numQuadNodes,double);
1064     quadWeights1d=TMPMEMALLOC(numQuadNodes,double);
1065     if (! ( Finley_checkPtr(quadNodes1d) || Finley_checkPtr(quadWeights1d) ) ) {
1066     for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
1067     if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
1068 jgs 82
1069 gross 1028 /* get 1D scheme: */
1070 jgs 82
1071 gross 1028 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
1072 jgs 82
1073 gross 1028 /* make 3D scheme: */
1074 jgs 82
1075 gross 1028 if (Finley_noError()) {
1076     l=0;
1077     for (i=0;i<numQuadNodes1d;i++) {
1078     for (j=0;j<numQuadNodes1d;j++) {
1079     for (k=0;k<numQuadNodes1d;k++) {
1080     QUADNODES(0,l)=quadNodes1d[i];
1081     QUADNODES(1,l)=quadNodes1d[j];
1082     QUADNODES(2,l)=quadNodes1d[k];
1083     QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];
1084     l++;
1085     }
1086     }
1087     }
1088     set=TRUE;
1089     break;
1090     }
1091     }
1092     }
1093     if (!set) {
1094     sprintf(error_msg,"Finley_Quad_getNodesHex: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
1095     Finley_setError(VALUE_ERROR,error_msg);
1096     }
1097     TMPMEMFREE(quadNodes1d);
1098     TMPMEMFREE(quadWeights1d);
1099 jgs 82 }
1100     #undef DIM
1101     }
1102    
1103     /**************************************************************/
1104    
1105     /* get a quadrature scheme with numQuadNodes quadrature nodes for a point. As there */
1106     /* in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
1107     /* an error. */
1108    
1109     void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
1110 gross 2748 if (numQuadNodes==0) {
1111     return;
1112     } else {
1113     Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: Illegal number of quadrature nodes.");
1114 jgs 82 }
1115     }
1116    
1117     /**************************************************************/
1118    
1119     /* get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
1120     /* The nodes and weights are set from a table. */
1121    
1122     void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
1123     switch(numQuadNodes) {
1124     case 1:
1125     quadNodes[0]=0.5;
1126     quadWeights[0]=1.;
1127     break;
1128    
1129     case 2:
1130     quadNodes[0]=(1.-.577350269189626)/2.;
1131     quadNodes[1]=(1.+.577350269189626)/2.;
1132     quadWeights[0]=.5;
1133     quadWeights[1]=.5;
1134     break;
1135    
1136     case 3:
1137     quadNodes[0]=(1.-.774596669241483)/2.;
1138     quadNodes[1]=.5;
1139     quadNodes[2]=(1.+.774596669241483)/2.;
1140     quadWeights[0]=5./18.;
1141     quadWeights[1]=4./ 9.;
1142     quadWeights[2]=5./18.;
1143     break;
1144    
1145     case 4:
1146     quadNodes[0]=(1.-.861136311594053)/2.;
1147     quadNodes[1]=(1.-.339981043584856)/2.;
1148     quadNodes[2]=(1.+.339981043584856)/2.;
1149     quadNodes[3]=(1.+.861136311594053)/2.;
1150     quadWeights[0]=.347854845137454/2.;
1151     quadWeights[1]=.652145154862546/2.;
1152     quadWeights[2]=.652145154862546/2.;
1153     quadWeights[3]=.347854845137454/2.;
1154     break;
1155    
1156     case 5:
1157     quadNodes[0]=(1.-.906179845938664)/2.;
1158     quadNodes[1]=(1.-.538469310105683)/2.;
1159     quadNodes[2]= .5;
1160     quadNodes[3]=(1.+.538469310105683)/2.;
1161     quadNodes[4]=(1.+.906179845938664)/2.;
1162     quadWeights[0]=.236926885056189/2.;
1163     quadWeights[1]=.478628670499366/2.;
1164     quadWeights[2]=.568888888888889/2.;
1165     quadWeights[3]=.478628670499366/2.;
1166     quadWeights[4]=.236926885056189/2.;
1167     break;
1168    
1169     case 6:
1170     quadNodes[0]=(1.-.932469514203152)/2.;
1171     quadNodes[1]=(1.-.661209386466265)/2.;
1172     quadNodes[2]=(1.-.238619186083197)/2.;
1173     quadNodes[3]=(1.+.238619186083197)/2.;
1174     quadNodes[4]=(1.+.661209386466265)/2.;
1175     quadNodes[5]=(1.+.932469514203152)/2.;
1176     quadWeights[0]=.171324492379170/2.;
1177     quadWeights[1]=.360761573048139/2.;
1178     quadWeights[2]=.467913934572691/2.;
1179     quadWeights[3]=.467913934572691/2.;
1180     quadWeights[4]=.360761573048139/2.;
1181     quadWeights[5]=.171324492379170/2.;
1182     break;
1183    
1184     case 7:
1185     quadNodes[0]=(1.-.949107912342759)/2.;
1186     quadNodes[1]=(1.-.741531185599394)/2.;
1187     quadNodes[2]=(1.-.405845151377397)/2.;
1188     quadNodes[3]=0.5;
1189     quadNodes[4]=(1.+.405845151377397)/2.;
1190     quadNodes[5]=(1.+.741531185599394)/2.;
1191     quadNodes[6]=(1.+.949107912342759)/2.;
1192     quadWeights[0]= .129484966168870/2.;
1193     quadWeights[1]= .279705391489277/2.;
1194     quadWeights[2]= .381830050505119/2.;
1195     quadWeights[3]= .417959183673469/2.;
1196     quadWeights[4]= .381830050505119/2.;
1197     quadWeights[5]= .279705391489277/2.;
1198     quadWeights[6]= .129484966168870/2.;
1199     break;
1200    
1201     case 8:
1202     quadNodes[0]=(1.-.960289856497536)/2.;
1203     quadNodes[1]=(1.-.796666477413627)/2.;
1204     quadNodes[2]=(1.-.525532409916329)/2.;
1205     quadNodes[3]=(1.-.183434642495650)/2.;
1206     quadNodes[4]=(1.+.183434642495650)/2.;
1207     quadNodes[5]=(1.+.525532409916329)/2.;
1208     quadNodes[6]=(1.+.796666477413627)/2.;
1209     quadNodes[7]=(1.+.960289856497536)/2.;
1210     quadWeights[0]= .101228536290376/2.;
1211     quadWeights[1]= .222381034453374/2.;
1212     quadWeights[2]= .313706645877887/2.;
1213     quadWeights[3]= .362683783378362/2.;
1214     quadWeights[4]= .362683783378362/2.;
1215     quadWeights[5]= .313706645877887/2.;
1216     quadWeights[6]= .222381034453374/2.;
1217     quadWeights[7]= .101228536290376/2.;
1218     break;
1219    
1220     case 9:
1221     quadNodes[0]=(1.-.968160239507626)/2.;
1222     quadNodes[1]=(1.-.836031107326636)/2.;
1223     quadNodes[2]=(1.-.613371432700590)/2.;
1224     quadNodes[3]=(1.-.324253423403809)/2.;
1225     quadNodes[4]= .5;
1226     quadNodes[5]=(1.+.324253423403809)/2.;
1227     quadNodes[6]=(1.+.613371432700590)/2.;
1228     quadNodes[7]=(1.+.836031107326636)/2.;
1229     quadNodes[8]=(1.+.968160239507626)/2.;
1230     quadWeights[0]= .081274388361574/2.;
1231     quadWeights[1]= .180648160694857/2.;
1232     quadWeights[2]= .260610696402935/2.;
1233     quadWeights[3]= .312347077040003/2.;
1234     quadWeights[4]= .330239355001260/2.;
1235     quadWeights[5]= .312347077040003/2.;
1236     quadWeights[6]= .260610696402935/2.;
1237     quadWeights[7]= .180648160694857/2.;
1238     quadWeights[8]= .081274388361574/2.;
1239     break;
1240    
1241     case 10:
1242     quadNodes[0]=(1.-.973906528517172)/2.;
1243     quadNodes[1]=(1.-.865063366688985)/2.;
1244     quadNodes[2]=(1.-.679409568299024)/2.;
1245     quadNodes[3]=(1.-.433395394129247)/2.;
1246     quadNodes[4]=(1.-.148874338981631)/2.;
1247     quadNodes[5]=(1.+.148874338981631)/2.;
1248     quadNodes[6]=(1.+.433395394129247)/2.;
1249     quadNodes[7]=(1.+.679409568299024)/2.;
1250     quadNodes[8]=(1.+.865063366688985)/2.;
1251     quadNodes[9]=(1.+.973906528517172)/2.;
1252     quadWeights[0]= .066671344308688/2.;
1253     quadWeights[1]= .149451349150581/2.;
1254     quadWeights[2]= .219086362515982/2.;
1255     quadWeights[3]= .269266719309996/2.;
1256     quadWeights[4]= .295524224714753/2.;
1257     quadWeights[5]= .295524224714753/2.;
1258     quadWeights[6]= .269266719309996/2.;
1259     quadWeights[7]= .219086362515982/2.;
1260     quadWeights[8]= .149451349150581/2.;
1261     quadWeights[9]= .066671344308688/2.;
1262     break;
1263    
1264     default:
1265 gross 2748 Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesLine: Negative intergration order.");
1266 jgs 82 break;
1267     }
1268     }
1269    
1270    
1271     /**************************************************************/
1272    
1273     /* The following functions Finley_Quad_getNumNodes* return the nmber of quadrature points needed to */
1274     /* achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
1275     /* to consider the accuracy reduction through the construction process. */
1276    
1277    
1278     int Finley_Quad_getNumNodesPoint(int order) {
1279     return 0;
1280     }
1281    
1282     int Finley_Quad_getNumNodesLine(int order) {
1283 jgs 150 char error_msg[LenErrorMsg_MAX];
1284 jgs 82 if (order <0 ) {
1285 gross 964 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
1286 jgs 82 return -1;
1287     } else {
1288     if (order > 2*MAX_numQuadNodesLine-1) {
1289 gross 964 sprintf(error_msg,"Finley_Quad_getNumNodesPoint: requested integration order %d on line is too large (>%d).",
1290 jgs 82 order,2*MAX_numQuadNodesLine-1);
1291 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
1292 jgs 82 return -1;
1293     } else {
1294 jgs 150 Finley_resetError();
1295 jgs 82 return order/2+1;
1296     }
1297     }
1298     }
1299    
1300     int Finley_Quad_getNumNodesTri(int order) {
1301     int numQuadNodesLine;
1302 gross 1072 if (order<=1) {
1303 jgs 82 return 1;
1304 gross 1342 } else if (order<=2){
1305 btully 1170 return 3;
1306 gross 1342 } else if (order<=3){
1307 btully 1170 return 4;
1308 gross 1342 } else if (order<=4){
1309     return 6;
1310     } else if (order<=5){
1311     return 7;
1312     } else if (order<=6){
1313     return 12;
1314     } else if (order<=7){
1315     return 13;
1316     } else if (order<=8){
1317     return 16;
1318     } else if (order<=9){
1319     return 19;
1320 jgs 82 } else {
1321     numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
1322 jgs 150 if (Finley_noError()) {
1323 jgs 82 return numQuadNodesLine*numQuadNodesLine;
1324     } else {
1325     return -1;
1326     }
1327     }
1328     }
1329    
1330     int Finley_Quad_getNumNodesRec(int order) {
1331     int numQuadNodesLine;
1332     numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
1333 jgs 150 if (Finley_noError()) {
1334 jgs 82 return numQuadNodesLine*numQuadNodesLine;
1335     } else {
1336     return -1;
1337     }
1338     }
1339    
1340     int Finley_Quad_getNumNodesTet(int order) {
1341     int numQuadNodesLine;
1342 gross 1072 if (order<=1) {
1343     return 1;
1344 gross 1342 } else if (order<=2){
1345 btully 1170 return 4;
1346 gross 1342 } else if (order<=3){
1347 btully 1170 return 5;
1348 gross 1342 } else if (order<=4){
1349     return 11;
1350     } else if (order<=5){
1351     return 15;
1352     } else if (order<=6){
1353     return 24;
1354     } else if (order<=7){
1355     return 31;
1356     } else if (order<=8){
1357     return 45;
1358 jgs 82 } else {
1359 gross 1072 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
1360     if (Finley_noError()) {
1361     return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1362     } else {
1363     return -1;
1364     }
1365 jgs 82 }
1366     }
1367    
1368     int Finley_Quad_getNumNodesHex(int order) {
1369     int numQuadNodesLine;
1370     numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
1371 jgs 150 if (Finley_noError()) {
1372 jgs 82 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1373     } else {
1374     return -1;
1375     }
1376     }
1377 gross 2748 dim_t Finley_Quad_MacroPoint(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1378     dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1379     {
1380     return 0;
1381    
1382     }
1383     dim_t Finley_Quad_MacroLine(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1384     dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1385     {
1386     #define DIM 1
1387     dim_t s,q,i;
1388     register double x0, w;
1389     const double f=1./((double)numSubElements);
1390    
1391     if (new_len < numSubElements*numQuadNodes) {
1392     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroLine: array for new qurature scheme is too small");
1393     }
1394     for (q=0; q<numQuadNodes; ++q) {
1395    
1396     x0=quadNodes[INDEX2(0,q,DIM)];
1397     w=f*quadWeights[q];
1398    
1399     for (s=0; s<numSubElements; ++s) {
1400     new_quadWeights[INDEX2(q,s, numQuadNodes)] =w;
1401     new_quadNodes[INDEX3(0,q,s, DIM,numQuadNodes)] =(x0+s)*f;
1402     for (i=0;i<numF;i++) new_dFdv[INDEX4(i,0,q,s, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,numF, q,DIM)]*f;
1403     }
1404    
1405     }
1406     #undef DIM
1407     return numSubElements*numQuadNodes;
1408     }
1409     #define HALF 0.5
1410     #define TWO 2.
1411     dim_t Finley_Quad_MacroTri(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1412     dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1413     {
1414    
1415     #define DIM 2
1416     dim_t q,i;
1417     register double x0, x1, w, df0, df1;
1418    
1419     if (new_len < numSubElements*numQuadNodes) {
1420     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTri: array for new qurature scheme is too small");
1421     return -1;
1422     }
1423     if (numSubElements==1) {
1424    
1425     for (q=0; q<numQuadNodes; ++q) {
1426    
1427     x0=quadNodes[INDEX2(0,q,DIM)];
1428     x1=quadNodes[INDEX2(1,q,DIM)];
1429     w=quadWeights[q];
1430    
1431     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1432     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =x0;
1433     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =x1;
1434     for (i=0;i<numF;i++) {
1435     new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];
1436     new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];
1437     }
1438    
1439     }
1440    
1441     } else if (numSubElements==4) {
1442     const double f = 0.25;
1443     for (q=0; q<numQuadNodes; ++q) {
1444    
1445     x0=quadNodes[INDEX2(0,q,DIM)];
1446     x1=quadNodes[INDEX2(1,q,DIM)];
1447     w=f*quadWeights[q];
1448    
1449     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1450     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =HALF*x0;
1451     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =HALF*x1;
1452    
1453     new_quadWeights[INDEX2(q,1,numQuadNodes)] =w;
1454     new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)] =HALF*x0;
1455     new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)] =HALF*(x1+1);
1456    
1457     new_quadWeights[INDEX2(q,2,numQuadNodes)] =w;
1458     new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)] =HALF*(x0+1);
1459     new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)] =HALF*x1;
1460    
1461     new_quadWeights[INDEX2(q,3,numQuadNodes)] =w;
1462     new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)] =HALF*(1-x0);
1463     new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)] =HALF*(1-x1);
1464    
1465     for (i=0;i<numF;i++) {
1466     df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1467     df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1468    
1469     new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1470     new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1471    
1472     new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1473     new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1474    
1475     new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1476     new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1477    
1478     new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = -df0;
1479     new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = -df1;
1480     }
1481    
1482    
1483     }
1484     } else {
1485     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTri: unable to create quadrature scheme for macro element.");
1486     return -1;
1487     }
1488     #undef DIM
1489     return numSubElements*numQuadNodes;
1490     }
1491    
1492     dim_t Finley_Quad_MacroRec(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1493     dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1494     {
1495    
1496     #define DIM 2
1497     dim_t q,i;
1498     register double x0, x1, w, df0, df1;
1499    
1500     if (new_len < numSubElements*numQuadNodes) {
1501     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroRec: array for new qurature scheme is too small");
1502     return -1;
1503     }
1504     if (numSubElements==1) {
1505    
1506     for (q=0; q<numQuadNodes; ++q) {
1507    
1508     x0=quadNodes[INDEX2(0,q,DIM)];
1509     x1=quadNodes[INDEX2(1,q,DIM)];
1510     w=quadWeights[q];
1511    
1512     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1513     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =x0;
1514     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =x1;
1515     for (i=0;i<numF;i++) {
1516     new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0, q,numF, DIM)];
1517     new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1, q,numF, DIM)];
1518     }
1519    
1520     }
1521    
1522     } else if (numSubElements==4) {
1523     const double f = 0.25;
1524     for (q=0; q<numQuadNodes; ++q) {
1525    
1526     x0=quadNodes[INDEX2(0,q,DIM)];
1527     x1=quadNodes[INDEX2(1,q,DIM)];
1528     w=f*quadWeights[q];
1529    
1530     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1531     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =HALF*x0;
1532     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =HALF*x1;
1533    
1534     new_quadWeights[INDEX2(q,1,numQuadNodes)] =w;
1535     new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)] =HALF*x0;
1536     new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)] =HALF*(x1+1);
1537    
1538     new_quadWeights[INDEX2(q,2,numQuadNodes)] =w;
1539     new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)] =HALF*(x0+1);
1540     new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)] =HALF*x1;
1541    
1542     new_quadWeights[INDEX2(q,3,numQuadNodes)] =w;
1543     new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)] =HALF*(x0+1);
1544     new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)] =HALF*(x1+1);
1545    
1546     for (i=0;i<numF;i++) {
1547     df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1548     df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1549    
1550     new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1551     new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1552    
1553     new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1554     new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1555    
1556     new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1557     new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1558    
1559     new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = df0;
1560     new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1561     }
1562    
1563     }
1564     } else {
1565     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroRec: unable to create quadrature scheme for macro element.");
1566     return -1;
1567     }
1568     #undef DIM
1569     return numSubElements*numQuadNodes;
1570     }
1571    
1572    
1573     dim_t Finley_Quad_MacroTet(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1574     dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1575     {
1576     #define DIM 3
1577     dim_t q, i;
1578     register double x0, x1, x2, w, df0, df1, df2;
1579    
1580     if (new_len < numSubElements*numQuadNodes) {
1581     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTet: array for new qurature scheme is too small");
1582     return -1;
1583     }
1584     if (numSubElements==1) {
1585    
1586     for (q=0; q<numQuadNodes; ++q) {
1587    
1588     x0=quadNodes[INDEX2(0,q,DIM)];
1589     x1=quadNodes[INDEX2(1,q,DIM)];
1590     x2=quadNodes[INDEX2(2,q,DIM)];
1591     w=quadWeights[q];
1592    
1593     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1594     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =x0;
1595     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =x1;
1596     new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)] =x2;
1597    
1598     for (i=0;i<numF;i++) {
1599     new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];
1600     new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];
1601     new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q, numF, DIM)];
1602     }
1603     }
1604    
1605     } else if (numSubElements==8) {
1606     const double f = 0.125;
1607     for (q=0; q<numQuadNodes; ++q) {
1608    
1609     x0=quadNodes[INDEX2(0,q,DIM)];
1610     x1=quadNodes[INDEX2(1,q,DIM)];
1611     x2=quadNodes[INDEX2(2,q,DIM)];
1612     w=f*quadWeights[q];
1613    
1614     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1615     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =HALF*x0;
1616     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =HALF*x1;
1617     new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)] =HALF*x2;
1618    
1619     new_quadWeights[INDEX2(q,1,numQuadNodes)] =w;
1620     new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)] =HALF*(x0+1);
1621     new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)] =HALF*x1;
1622     new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)] =HALF*x2;
1623    
1624     new_quadWeights[INDEX2(q,2,numQuadNodes)] =w;
1625     new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)] =HALF*x0;
1626     new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)] =HALF*(x1+1);
1627     new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)] =HALF*x2;
1628    
1629     new_quadWeights[INDEX2(q,3,numQuadNodes)] =w;
1630     new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)] =HALF*x0;
1631     new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)] =HALF*x1;
1632     new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)] =HALF*(x2+1);
1633    
1634     new_quadWeights[INDEX2(q,4,numQuadNodes)] =w;
1635     new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)] =HALF*(1-x1);
1636     new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)] =HALF*(x0+x1);
1637     new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)] =HALF*x2;
1638    
1639     new_quadWeights[INDEX2(q,5,numQuadNodes)] =w;
1640     new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)] =HALF*(1-x0-x2);
1641     new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)] =HALF*(1-x1);
1642     new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)] =HALF*(x0+x1);
1643    
1644     new_quadWeights[INDEX2(q,6,numQuadNodes)] =w;
1645     new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)] =HALF*x2;
1646     new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)] =HALF*(1-x0-x2);
1647     new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)] =HALF*(1-x1);
1648    
1649     new_quadWeights[INDEX2(q,7,numQuadNodes)] =w;
1650     new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)] =HALF*(x0+x2);
1651     new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)] =HALF*x1;
1652     new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)] =HALF*(1-x0-x1);
1653    
1654     for (i=0;i<numF;i++) {
1655     df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1656     df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1657     df2=dFdv[INDEX3(i,2,q, numF, DIM)]*TWO;
1658    
1659     new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1660     new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1661     new_dFdv[INDEX4(i,2,q,0, numF,DIM,numQuadNodes)] = df2;
1662    
1663     new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1664     new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1665     new_dFdv[INDEX4(i,2,q,1, numF,DIM,numQuadNodes)] = df2;
1666    
1667     new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1668     new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1669     new_dFdv[INDEX4(i,2,q,2, numF,DIM,numQuadNodes)] = df2;
1670    
1671     new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = df0;
1672     new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1673     new_dFdv[INDEX4(i,2,q,3, numF,DIM,numQuadNodes)] = df2;
1674    
1675     new_dFdv[INDEX4(i,0,q,4, numF,DIM,numQuadNodes)] = df0-df1;
1676     new_dFdv[INDEX4(i,1,q,4, numF,DIM,numQuadNodes)] = df0;
1677     new_dFdv[INDEX4(i,2,q,4, numF,DIM,numQuadNodes)] = df2;
1678    
1679     new_dFdv[INDEX4(i,0,q,5, numF,DIM,numQuadNodes)] = -df2;
1680     new_dFdv[INDEX4(i,1,q,5, numF,DIM,numQuadNodes)] = df0-df2-df1;
1681     new_dFdv[INDEX4(i,2,q,5, numF,DIM,numQuadNodes)] = df0-df2;
1682    
1683     new_dFdv[INDEX4(i,0,q,6, numF,DIM,numQuadNodes)] = -df0+df2;
1684     new_dFdv[INDEX4(i,1,q,6, numF,DIM,numQuadNodes)] = -df0;
1685     new_dFdv[INDEX4(i,2,q,6, numF,DIM,numQuadNodes)] = -df1;
1686    
1687     new_dFdv[INDEX4(i,0,q,7, numF,DIM,numQuadNodes)] = df2;
1688     new_dFdv[INDEX4(i,1,q,7, numF,DIM,numQuadNodes)] = -df0+df1+df2;
1689     new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = -df0+df2;
1690     }
1691    
1692     }
1693     } else {
1694     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTet: unable to create quadrature scheme for macro element.");
1695     return -1;
1696     }
1697     #undef DIM
1698     return numSubElements*numQuadNodes;
1699     }
1700    
1701    
1702     dim_t Finley_Quad_MacroHex(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1703     dim_t new_len, double* new_quadNodes, double* new_quadWeights , double* new_dFdv)
1704     {
1705    
1706     #define DIM 3
1707     dim_t q, i;
1708     register double x0, x1, x2, w, df0, df1, df2;
1709    
1710     if (new_len < numSubElements*numQuadNodes) {
1711     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroHex: array for new qurature scheme is too small");
1712     return -1;
1713     }
1714     if (numSubElements==1) {
1715    
1716     for (q=0; q<numQuadNodes; ++q) {
1717    
1718     x0=quadNodes[INDEX2(0,q,DIM)];
1719     x1=quadNodes[INDEX2(1,q,DIM)];
1720     x2=quadNodes[INDEX2(2,q,DIM)];
1721     w=quadWeights[q];
1722    
1723     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1724     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =x0;
1725     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =x1;
1726     new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)] =x2;
1727    
1728     for (i=0;i<numF;i++) {
1729     new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q,numF, DIM)];
1730     new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q,numF, DIM)];
1731     new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q,numF, DIM)];
1732     }
1733     }
1734    
1735     } else if (numSubElements==8) {
1736     const double f = 0.125;
1737     for (q=0; q<numQuadNodes; ++q) {
1738    
1739     x0=quadNodes[INDEX2(0,q,DIM)];
1740     x1=quadNodes[INDEX2(1,q,DIM)];
1741     x2=quadNodes[INDEX2(2,q,DIM)];
1742     w=f*quadWeights[q];
1743    
1744     new_quadWeights[INDEX2(q,0,numQuadNodes)] =w;
1745     new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)] =HALF*x0;
1746     new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)] =HALF*x1;
1747     new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)] =HALF*x2;
1748    
1749     new_quadWeights[INDEX2(q,1,numQuadNodes)] =w;
1750     new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)] =HALF*(x0+1);
1751     new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)] =HALF*x1;
1752     new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)] =HALF*x2;
1753    
1754     new_quadWeights[INDEX2(q,2,numQuadNodes)] =w;
1755     new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)] =HALF*x0;
1756     new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)] =HALF*(x1+1);
1757     new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)] =HALF*x2;
1758    
1759     new_quadWeights[INDEX2(q,3,numQuadNodes)] =w;
1760     new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)] =HALF*(x0+1);
1761     new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)] =HALF*(x1+1);
1762     new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)] =HALF*x2;
1763    
1764     new_quadWeights[INDEX2(q,4,numQuadNodes)] =w;
1765     new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)] =HALF*x0;
1766     new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)] =HALF*x1;
1767     new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)] =HALF*(x2+1);
1768    
1769     new_quadWeights[INDEX2(q,5,numQuadNodes)] =w;
1770     new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)] =HALF*(x0+1);
1771     new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)] =HALF*x1;
1772     new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)] =HALF*(x2+1);
1773    
1774     new_quadWeights[INDEX2(q,6,numQuadNodes)] =w;
1775     new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)] =HALF*x0;
1776     new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)] =HALF*(x1+1);
1777     new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)] =HALF*(x2+1);
1778    
1779     new_quadWeights[INDEX2(q,7,numQuadNodes)] =w;
1780     new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)] =HALF*(x0+1);
1781     new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)] =HALF*(x1+1);
1782     new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)] =HALF*(x2+1);
1783    
1784     for (i=0;i<numF;i++) {
1785     df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1786     df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1787     df2=dFdv[INDEX3(i,2,q, numF, DIM)]*TWO;
1788    
1789     new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1790     new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1791     new_dFdv[INDEX4(i,2,q,0, numF,DIM,numQuadNodes)] = df2;
1792    
1793     new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1794     new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1795     new_dFdv[INDEX4(i,2,q,1, numF,DIM,numQuadNodes)] = df2;
1796    
1797     new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1798     new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1799     new_dFdv[INDEX4(i,2,q,2, numF,DIM,numQuadNodes)] = df2;
1800    
1801     new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = df0;
1802     new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1803     new_dFdv[INDEX4(i,2,q,3, numF,DIM,numQuadNodes)] = df2;
1804    
1805     new_dFdv[INDEX4(i,0,q,4, numF,DIM,numQuadNodes)] = df0;
1806     new_dFdv[INDEX4(i,1,q,4, numF,DIM,numQuadNodes)] = df1;
1807     new_dFdv[INDEX4(i,2,q,4, numF,DIM,numQuadNodes)] = df2;
1808    
1809     new_dFdv[INDEX4(i,0,q,5, numF,DIM,numQuadNodes)] = df0;
1810     new_dFdv[INDEX4(i,1,q,5, numF,DIM,numQuadNodes)] = df1;
1811     new_dFdv[INDEX4(i,2,q,5, numF,DIM,numQuadNodes)] = df2;
1812    
1813     new_dFdv[INDEX4(i,0,q,6, numF,DIM,numQuadNodes)] = df0;
1814     new_dFdv[INDEX4(i,1,q,6, numF,DIM,numQuadNodes)] = df1;
1815     new_dFdv[INDEX4(i,2,q,6, numF,DIM,numQuadNodes)] = df2;
1816    
1817     new_dFdv[INDEX4(i,0,q,7, numF,DIM,numQuadNodes)] = df0;
1818     new_dFdv[INDEX4(i,1,q,7, numF,DIM,numQuadNodes)] = df1;
1819     new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = df2;
1820     }
1821    
1822     }
1823     } else {
1824     Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroHex: unable to create quadrature scheme for macro element.");
1825     return -1;
1826     }
1827     #undef DIM
1828     return numSubElements*numQuadNodes;
1829     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26