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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show 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 /* $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