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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1072 - (hide annotations)
Thu Mar 29 06:44:30 2007 UTC (12 years, 4 months ago) by gross
File MIME type: text/plain
File size: 17572 byte(s)
PDE assemblage for reduced integration order + tests added.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26