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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Quadrature.c
File MIME type: text/plain
File size: 15806 byte(s)
Initial revision

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26