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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 69473 byte(s)
Updated the copyright header.


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Quadrature.cpp:5567-5588 /branches/diaplayground/finley/src/Quadrature.cpp:4940-5147 /branches/lapack2681/finley/src/Quadrature.cpp:2682-2741 /branches/pasowrap/finley/src/Quadrature.cpp:3661-3674 /branches/py3_attempt2/finley/src/Quadrature.cpp:3871-3891 /branches/restext/finley/src/Quadrature.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Quadrature.cpp:3669-3791 /branches/stage3.0/finley/src/Quadrature.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Quadrature.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Quadrature.cpp:5898-6118 /release/3.0/finley/src/Quadrature.cpp:2591-2601 /release/4.0/finley/src/Quadrature.cpp:5380-5406 /trunk/finley/src/Quadrature.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26