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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 201 - (hide annotations)
Wed Nov 23 04:10:21 2005 UTC (13 years, 8 months ago) by jgs
Original Path: trunk/finley/src/finley/Quadrature.c
File MIME type: text/plain
File size: 16822 byte(s)
copy finleyC and CPPAdapter to finley and finley/CPPAdapter to
facilitate scons builds

1 jgs 150 /*
2     ******************************************************************************
3     * *
4     * COPYRIGHT ACcESS 2003,2004,2005 - All Rights Reserved *
5     * *
6     * This software is the property of ACcESS. No part of this code *
7     * may be copied in any form or by any means without the expressed written *
8     * consent of ACcESS. Copying, use or modification of this software *
9     * by any unauthorised person is illegal unless that person has a software *
10     * license agreement with ACcESS. *
11     * *
12     ******************************************************************************
13     */
14 jgs 82
15     /**************************************************************/
16    
17     /* Finley: */
18    
19     /**************************************************************/
20    
21     /* Author: gross@access.edu.au */
22     /* Version: $Id$ */
23    
24     /**************************************************************/
25    
26     #include "Quadrature.h"
27    
28     #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
29     #define QUADWEIGHTS(_I_) quadWeights[_I_]
30    
31     /**************************************************************/
32    
33     /* get a quadrature scheme with numQuadNodes quadrature nodes for the tri */
34     /* as a queezed scheme on a quad [0,1]^2 */
35    
36     void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
37     int i;
38     double Q1,Q2;
39     #define DIM 2
40    
41     /* the easy case: */
42    
43     if (numQuadNodes==1) {
44     QUADNODES(0,0)=1./3.;
45     QUADNODES(1,0)=1./3.;
46     QUADWEIGHTS(0)= .5;
47     } else {
48    
49     /* get scheme on [0.1]^2 */
50    
51     Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
52 jgs 150 if (! Finley_noError()) return;
53 jgs 82
54     /* squeeze it: */
55    
56     for (i=0;i<numQuadNodes;i++) {
57     Q1=QUADNODES(0,i);
58     Q2=QUADNODES(1,i);
59     QUADWEIGHTS(i)=QUADWEIGHTS(i)*(1.-(1./2.)*(Q1+Q2));
60     QUADNODES(0,i)=Q1*(1.-(1./2.)*Q2);
61     QUADNODES(1,i)=Q2*(1.-(1./2.)*Q1);
62     }
63     }
64     #undef DIM
65     }
66    
67     /**************************************************************/
68    
69     /* get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
70     /* as a queezed scheme on a hex [0,1]^3 */
71    
72     void Finley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
73     int i;
74     double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;
75     #define DIM 3
76    
77     /* the easy case: */
78    
79     if (numQuadNodes==1) {
80     QUADNODES(0,0)= .25;
81     QUADNODES(1,0)= .25;
82     QUADNODES(2,0)= .25;
83     QUADWEIGHTS(0)=1./6.;
84     } else {
85    
86     /* get scheme on [0.1]^3 */
87    
88     Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
89 jgs 150 if (! Finley_noError()) return;
90 jgs 82
91     /* squeeze it: */
92    
93     for (i=0;i<numQuadNodes;i++) {
94     Q1=QUADNODES(0,i);
95     Q2=QUADNODES(1,i);
96     Q3=QUADNODES(2,i);
97    
98     JA11= (1./3.)*Q2*Q3-(1./2.)*(Q2+Q3) +1.;
99     JA12= (1./3.)*Q1*Q3-(1./2.)*Q1;
100     JA13= (1./3.)*Q1*Q2-(1./2.)*Q1;
101     JA21= (1./3.)*Q2*Q3-(1./2.)*Q2;
102     JA22= (1./3.)*Q1*Q3-(1./2.)*(Q1+Q3) +1.;
103     JA23= (1./3.)*Q1*Q2-(1./2.)*Q2;
104     JA31= (1./3.)*Q2*Q3-(1./2.)*Q3;
105     JA32= (1./3.)*Q1*Q3-(1./2.)*Q3;
106     JA33= (1./3.)*Q1*Q2-(1./2.)*(Q1+Q2) +1.;
107     DET=JA11*JA22*JA33+JA12*JA23*JA31+JA13*JA21*JA32-JA13*JA22*JA31-JA11*JA23*JA32-JA12*JA21*JA33;
108     quadWeights[i]=quadWeights[i]*ABS(DET);
109     QUADNODES(0,i)=Q1*((1./3.)*Q2*Q3-(1./2.)*(Q2+Q3)+1.);
110     QUADNODES(1,i)=Q2*((1./3.)*Q1*Q3-(1./2.)*(Q1+Q3)+1.);
111     QUADNODES(2,i)=Q3*((1./3.)*Q1*Q2-(1./2.)*(Q1+Q2)+1.);
112     }
113     }
114     #undef DIM
115     }
116    
117     /**************************************************************/
118    
119     /* get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */
120     /* as a X-product of a 1D scheme. */
121    
122     void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
123 jgs 150 char error_msg[LenErrorMsg_MAX];
124 jgs 82 int numQuadNodes1d,i,j,l;
125     double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];
126     #define DIM 2
127    
128     /* find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */
129    
130     for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
131     if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
132    
133     /* get 1D scheme: */
134    
135     Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
136    
137     /* make 2D scheme: */
138    
139 jgs 150 if (Finley_noError()) {
140 jgs 82 l=0;
141     for (i=0;i<numQuadNodes1d;i++) {
142     for (j=0;j<numQuadNodes1d;j++) {
143     QUADNODES(0,l)=quadNodes1d[i];
144     QUADNODES(1,l)=quadNodes1d[j];
145     QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];
146     l++;
147     }
148     }
149     }
150     return;
151     }
152     }
153 jgs 150 sprintf(error_msg,"__FILE__: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
154     Finley_setError(VALUE_ERROR,error_msg);
155 jgs 82 #undef DIM
156     }
157    
158     /**************************************************************/
159    
160     /* get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
161     /* as a X-product of a 1D scheme. */
162    
163     void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
164 jgs 150 char error_msg[LenErrorMsg_MAX];
165 jgs 82 int numQuadNodes1d,i,j,k,l;
166     double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];
167     #define DIM 3
168    
169     /* find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
170    
171     for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
172     if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
173    
174     /* get 1D scheme: */
175    
176     Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
177    
178     /* make 3D scheme: */
179    
180 jgs 150 if (Finley_noError()) {
181 jgs 82 l=0;
182     for (i=0;i<numQuadNodes1d;i++) {
183     for (j=0;j<numQuadNodes1d;j++) {
184     for (k=0;k<numQuadNodes1d;k++) {
185     QUADNODES(0,l)=quadNodes1d[i];
186     QUADNODES(1,l)=quadNodes1d[j];
187     QUADNODES(2,l)=quadNodes1d[k];
188     QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];
189     l++;
190     }
191     }
192     }
193     }
194    
195     return;
196     }
197     }
198 jgs 150 sprintf(error_msg,"__FILE__: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
199     Finley_setError(VALUE_ERROR,error_msg);
200 jgs 82 #undef DIM
201     }
202    
203     /**************************************************************/
204    
205     /* get a quadrature scheme with numQuadNodes quadrature nodes for a point. As there */
206     /* in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
207     /* an error. */
208    
209     void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
210     if (numQuadNodes!=0) {
211 jgs 150 Finley_setError(VALUE_ERROR,"__FILE__: There is no quadrature scheme on points.");
212 jgs 82 }
213     }
214    
215     /**************************************************************/
216    
217     /* get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
218     /* The nodes and weights are set from a table. */
219    
220     void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
221     switch(numQuadNodes) {
222     case 1:
223     quadNodes[0]=0.5;
224     quadWeights[0]=1.;
225     break;
226    
227     case 2:
228     quadNodes[0]=(1.-.577350269189626)/2.;
229     quadNodes[1]=(1.+.577350269189626)/2.;
230     quadWeights[0]=.5;
231     quadWeights[1]=.5;
232     break;
233    
234     case 3:
235     quadNodes[0]=(1.-.774596669241483)/2.;
236     quadNodes[1]=.5;
237     quadNodes[2]=(1.+.774596669241483)/2.;
238     quadWeights[0]=5./18.;
239     quadWeights[1]=4./ 9.;
240     quadWeights[2]=5./18.;
241     break;
242    
243     case 4:
244     quadNodes[0]=(1.-.861136311594053)/2.;
245     quadNodes[1]=(1.-.339981043584856)/2.;
246     quadNodes[2]=(1.+.339981043584856)/2.;
247     quadNodes[3]=(1.+.861136311594053)/2.;
248     quadWeights[0]=.347854845137454/2.;
249     quadWeights[1]=.652145154862546/2.;
250     quadWeights[2]=.652145154862546/2.;
251     quadWeights[3]=.347854845137454/2.;
252     break;
253    
254     case 5:
255     quadNodes[0]=(1.-.906179845938664)/2.;
256     quadNodes[1]=(1.-.538469310105683)/2.;
257     quadNodes[2]= .5;
258     quadNodes[3]=(1.+.538469310105683)/2.;
259     quadNodes[4]=(1.+.906179845938664)/2.;
260     quadWeights[0]=.236926885056189/2.;
261     quadWeights[1]=.478628670499366/2.;
262     quadWeights[2]=.568888888888889/2.;
263     quadWeights[3]=.478628670499366/2.;
264     quadWeights[4]=.236926885056189/2.;
265     break;
266    
267     case 6:
268     quadNodes[0]=(1.-.932469514203152)/2.;
269     quadNodes[1]=(1.-.661209386466265)/2.;
270     quadNodes[2]=(1.-.238619186083197)/2.;
271     quadNodes[3]=(1.+.238619186083197)/2.;
272     quadNodes[4]=(1.+.661209386466265)/2.;
273     quadNodes[5]=(1.+.932469514203152)/2.;
274     quadWeights[0]=.171324492379170/2.;
275     quadWeights[1]=.360761573048139/2.;
276     quadWeights[2]=.467913934572691/2.;
277     quadWeights[3]=.467913934572691/2.;
278     quadWeights[4]=.360761573048139/2.;
279     quadWeights[5]=.171324492379170/2.;
280     break;
281    
282     case 7:
283     quadNodes[0]=(1.-.949107912342759)/2.;
284     quadNodes[1]=(1.-.741531185599394)/2.;
285     quadNodes[2]=(1.-.405845151377397)/2.;
286     quadNodes[3]=0.5;
287     quadNodes[4]=(1.+.405845151377397)/2.;
288     quadNodes[5]=(1.+.741531185599394)/2.;
289     quadNodes[6]=(1.+.949107912342759)/2.;
290     quadWeights[0]= .129484966168870/2.;
291     quadWeights[1]= .279705391489277/2.;
292     quadWeights[2]= .381830050505119/2.;
293     quadWeights[3]= .417959183673469/2.;
294     quadWeights[4]= .381830050505119/2.;
295     quadWeights[5]= .279705391489277/2.;
296     quadWeights[6]= .129484966168870/2.;
297     break;
298    
299     case 8:
300     quadNodes[0]=(1.-.960289856497536)/2.;
301     quadNodes[1]=(1.-.796666477413627)/2.;
302     quadNodes[2]=(1.-.525532409916329)/2.;
303     quadNodes[3]=(1.-.183434642495650)/2.;
304     quadNodes[4]=(1.+.183434642495650)/2.;
305     quadNodes[5]=(1.+.525532409916329)/2.;
306     quadNodes[6]=(1.+.796666477413627)/2.;
307     quadNodes[7]=(1.+.960289856497536)/2.;
308     quadWeights[0]= .101228536290376/2.;
309     quadWeights[1]= .222381034453374/2.;
310     quadWeights[2]= .313706645877887/2.;
311     quadWeights[3]= .362683783378362/2.;
312     quadWeights[4]= .362683783378362/2.;
313     quadWeights[5]= .313706645877887/2.;
314     quadWeights[6]= .222381034453374/2.;
315     quadWeights[7]= .101228536290376/2.;
316     break;
317    
318     case 9:
319     quadNodes[0]=(1.-.968160239507626)/2.;
320     quadNodes[1]=(1.-.836031107326636)/2.;
321     quadNodes[2]=(1.-.613371432700590)/2.;
322     quadNodes[3]=(1.-.324253423403809)/2.;
323     quadNodes[4]= .5;
324     quadNodes[5]=(1.+.324253423403809)/2.;
325     quadNodes[6]=(1.+.613371432700590)/2.;
326     quadNodes[7]=(1.+.836031107326636)/2.;
327     quadNodes[8]=(1.+.968160239507626)/2.;
328     quadWeights[0]= .081274388361574/2.;
329     quadWeights[1]= .180648160694857/2.;
330     quadWeights[2]= .260610696402935/2.;
331     quadWeights[3]= .312347077040003/2.;
332     quadWeights[4]= .330239355001260/2.;
333     quadWeights[5]= .312347077040003/2.;
334     quadWeights[6]= .260610696402935/2.;
335     quadWeights[7]= .180648160694857/2.;
336     quadWeights[8]= .081274388361574/2.;
337     break;
338    
339     case 10:
340     quadNodes[0]=(1.-.973906528517172)/2.;
341     quadNodes[1]=(1.-.865063366688985)/2.;
342     quadNodes[2]=(1.-.679409568299024)/2.;
343     quadNodes[3]=(1.-.433395394129247)/2.;
344     quadNodes[4]=(1.-.148874338981631)/2.;
345     quadNodes[5]=(1.+.148874338981631)/2.;
346     quadNodes[6]=(1.+.433395394129247)/2.;
347     quadNodes[7]=(1.+.679409568299024)/2.;
348     quadNodes[8]=(1.+.865063366688985)/2.;
349     quadNodes[9]=(1.+.973906528517172)/2.;
350     quadWeights[0]= .066671344308688/2.;
351     quadWeights[1]= .149451349150581/2.;
352     quadWeights[2]= .219086362515982/2.;
353     quadWeights[3]= .269266719309996/2.;
354     quadWeights[4]= .295524224714753/2.;
355     quadWeights[5]= .295524224714753/2.;
356     quadWeights[6]= .269266719309996/2.;
357     quadWeights[7]= .219086362515982/2.;
358     quadWeights[8]= .149451349150581/2.;
359     quadWeights[9]= .066671344308688/2.;
360     break;
361    
362     default:
363 jgs 150 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
364 jgs 82 break;
365     }
366     }
367     /**************************************************************/
368     /* */
369     /* the following function are used define the meshes on the surface in the xy-plane */
370    
371     /* triangle surface on a tetrahedron */
372     void Finley_Quad_getNodesTriOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
373     Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesTri);
374     }
375     /* rectangular surface on a hexahedron */
376     void Finley_Quad_getNodesRecOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
377     Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesRec);
378     }
379     /* line surface on a triangle or rectangle */
380     void Finley_Quad_getNodesLineOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
381     Finley_Quad_makeNodesOnFace(2,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesLine);
382     }
383     /* point surface on a line */
384     void Finley_Quad_getNodesPointOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
385     Finley_Quad_makeNodesOnFace(1,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesPoint);
386     }
387    
388     void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {
389     int q,i;
390     double quadNodesOnFace[numQuadNodes*(dim-1)];
391    
392     #define DIM dim
393     getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);
394 jgs 150 if (! Finley_noError()) return;
395 jgs 82
396     for (q=0;q<numQuadNodes;q++) {
397     for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];
398     QUADNODES(dim-1,q)=0;
399     }
400     #undef DIM
401     }
402    
403     /**************************************************************/
404    
405     /* The following functions Finley_Quad_getNumNodes* return the nmber of quadrature points needed to */
406     /* achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
407     /* to consider the accuracy reduction through the construction process. */
408    
409    
410     int Finley_Quad_getNumNodesPoint(int order) {
411     if (order <0 ) {
412 jgs 150 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
413 jgs 82 return -1;
414     } else {
415     return 0;
416     }
417     }
418    
419     int Finley_Quad_getNumNodesLine(int order) {
420 jgs 150 char error_msg[LenErrorMsg_MAX];
421 jgs 82 if (order <0 ) {
422 jgs 150 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
423 jgs 82 return -1;
424     } else {
425     if (order > 2*MAX_numQuadNodesLine-1) {
426 jgs 150 sprintf(error_msg,"__FILE__: requested integration order %d on line is too large (>%d).",
427 jgs 82 order,2*MAX_numQuadNodesLine-1);
428 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
429 jgs 82 return -1;
430     } else {
431 jgs 150 Finley_resetError();
432 jgs 82 return order/2+1;
433     }
434     }
435     }
436    
437     int Finley_Quad_getNumNodesTri(int order) {
438     int numQuadNodesLine;
439     if (order==1) {
440     return 1;
441     } else {
442     numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
443 jgs 150 if (Finley_noError()) {
444 jgs 82 return numQuadNodesLine*numQuadNodesLine;
445     } else {
446     return -1;
447     }
448     }
449     }
450    
451     int Finley_Quad_getNumNodesRec(int order) {
452     int numQuadNodesLine;
453     numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
454 jgs 150 if (Finley_noError()) {
455 jgs 82 return numQuadNodesLine*numQuadNodesLine;
456     } else {
457     return -1;
458     }
459     }
460    
461     int Finley_Quad_getNumNodesTet(int order) {
462     int numQuadNodesLine;
463     numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
464 jgs 150 if (Finley_noError()) {
465 jgs 82 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
466     } else {
467     return -1;
468     }
469     }
470    
471     int Finley_Quad_getNumNodesHex(int order) {
472     int numQuadNodesLine;
473     numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
474 jgs 150 if (Finley_noError()) {
475 jgs 82 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
476     } else {
477     return -1;
478     }
479     }
480     /*
481     * $Log$
482 jgs 150 * Revision 1.2 2005/09/15 03:44:23 jgs
483     * Merge of development branch dev-02 back to main trunk on 2005-09-15
484 jgs 82 *
485 jgs 150 * Revision 1.1.1.1.6.1 2005/09/07 06:26:20 gross
486     * the solver from finley are put into the standalone package paso now
487     *
488     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
489     * initial import of project esys2
490     *
491 jgs 82 * Revision 1.2 2004/08/03 04:49:06 gross
492     * bug in Quadrature.c fixed
493     *
494     * Revision 1.1.1.1 2004/06/24 04:00:40 johng
495     * Initial version of eys using boost-python.
496     *
497     *
498     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26