/[escript]/branches/RW_WIN32/finley/src/finleyC/Quadrature.c
ViewVC logotype

Contents of /branches/RW_WIN32/finley/src/finleyC/Quadrature.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 186 - (show annotations)
Sun Nov 20 10:51:16 2005 UTC (14 years ago) by robwdcock
File MIME type: text/plain
File size: 17544 byte(s)
PARTIAL WIN32 PORT
+ Paso, ext_math and mmio libraries build using scons in top level directory
+ Scons build system currently hard coded against RW's third-party library paths.
+ Still lots TODO to complete port and build system 
1 /*
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
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 if (! Finley_noError()) return;
53
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 if (! Finley_noError()) return;
90
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 char error_msg[LenErrorMsg_MAX];
124 int numQuadNodes1d,i,j,l;
125 #define DIM 2
126 /* Win32 Refactor */
127 double *quadNodes1d = (numQuadNodes>0) ? TMPMEMALLOC(numQuadNodes,double) : (double*)NULL;
128 double *quadWeights1d = (numQuadNodes>0) ? TMPMEMALLOC(numQuadNodes,double) : (double*)NULL;
129
130 /* find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */
131
132 for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
133 if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
134
135 /* get 1D scheme: */
136
137 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
138
139 /* make 2D scheme: */
140
141 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 }
152 /* Win32 Refactor */
153 TMPMEMFREE(quadNodes1d);
154 TMPMEMFREE(quadWeights1d);
155 return;
156 }
157 }
158 sprintf(error_msg,"__FILE__: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
159 Finley_setError(VALUE_ERROR,error_msg);
160
161 /* Win32 Refactor */
162 TMPMEMFREE(quadNodes1d);
163 TMPMEMFREE(quadWeights1d);
164 #undef DIM
165 }
166
167 /**************************************************************/
168
169 /* get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
170 /* as a X-product of a 1D scheme. */
171
172 void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
173 char error_msg[LenErrorMsg_MAX];
174 int numQuadNodes1d,i,j,k,l;
175 #define DIM 3
176 /* Win32 Refactor */
177 double *quadNodes1d = (numQuadNodes>0) ? TMPMEMALLOC(numQuadNodes,double) : (double*)NULL;
178 double *quadWeights1d = (numQuadNodes>0) ? TMPMEMALLOC(numQuadNodes,double) : (double*)NULL;
179
180 /* find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
181
182 for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
183 if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
184
185 /* get 1D scheme: */
186
187 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
188
189 /* make 3D scheme: */
190
191 if (Finley_noError()) {
192 l=0;
193 for (i=0;i<numQuadNodes1d;i++) {
194 for (j=0;j<numQuadNodes1d;j++) {
195 for (k=0;k<numQuadNodes1d;k++) {
196 QUADNODES(0,l)=quadNodes1d[i];
197 QUADNODES(1,l)=quadNodes1d[j];
198 QUADNODES(2,l)=quadNodes1d[k];
199 QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];
200 l++;
201 }
202 }
203 }
204 }
205 /* Win32 Refactor */
206 TMPMEMFREE(quadNodes1d);
207 TMPMEMFREE(quadWeights1d);
208
209 return;
210 }
211 }
212 sprintf(error_msg,"__FILE__: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
213 Finley_setError(VALUE_ERROR,error_msg);
214 /* Win32 Refactor */
215 TMPMEMFREE(quadNodes1d);
216 TMPMEMFREE(quadWeights1d);
217
218 #undef DIM
219 }
220
221 /**************************************************************/
222
223 /* get a quadrature scheme with numQuadNodes quadrature nodes for a point. As there */
224 /* in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
225 /* an error. */
226
227 void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
228 if (numQuadNodes!=0) {
229 Finley_setError(VALUE_ERROR,"__FILE__: There is no quadrature scheme on points.");
230 }
231 }
232
233 /**************************************************************/
234
235 /* get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
236 /* The nodes and weights are set from a table. */
237
238 void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
239 switch(numQuadNodes) {
240 case 1:
241 quadNodes[0]=0.5;
242 quadWeights[0]=1.;
243 break;
244
245 case 2:
246 quadNodes[0]=(1.-.577350269189626)/2.;
247 quadNodes[1]=(1.+.577350269189626)/2.;
248 quadWeights[0]=.5;
249 quadWeights[1]=.5;
250 break;
251
252 case 3:
253 quadNodes[0]=(1.-.774596669241483)/2.;
254 quadNodes[1]=.5;
255 quadNodes[2]=(1.+.774596669241483)/2.;
256 quadWeights[0]=5./18.;
257 quadWeights[1]=4./ 9.;
258 quadWeights[2]=5./18.;
259 break;
260
261 case 4:
262 quadNodes[0]=(1.-.861136311594053)/2.;
263 quadNodes[1]=(1.-.339981043584856)/2.;
264 quadNodes[2]=(1.+.339981043584856)/2.;
265 quadNodes[3]=(1.+.861136311594053)/2.;
266 quadWeights[0]=.347854845137454/2.;
267 quadWeights[1]=.652145154862546/2.;
268 quadWeights[2]=.652145154862546/2.;
269 quadWeights[3]=.347854845137454/2.;
270 break;
271
272 case 5:
273 quadNodes[0]=(1.-.906179845938664)/2.;
274 quadNodes[1]=(1.-.538469310105683)/2.;
275 quadNodes[2]= .5;
276 quadNodes[3]=(1.+.538469310105683)/2.;
277 quadNodes[4]=(1.+.906179845938664)/2.;
278 quadWeights[0]=.236926885056189/2.;
279 quadWeights[1]=.478628670499366/2.;
280 quadWeights[2]=.568888888888889/2.;
281 quadWeights[3]=.478628670499366/2.;
282 quadWeights[4]=.236926885056189/2.;
283 break;
284
285 case 6:
286 quadNodes[0]=(1.-.932469514203152)/2.;
287 quadNodes[1]=(1.-.661209386466265)/2.;
288 quadNodes[2]=(1.-.238619186083197)/2.;
289 quadNodes[3]=(1.+.238619186083197)/2.;
290 quadNodes[4]=(1.+.661209386466265)/2.;
291 quadNodes[5]=(1.+.932469514203152)/2.;
292 quadWeights[0]=.171324492379170/2.;
293 quadWeights[1]=.360761573048139/2.;
294 quadWeights[2]=.467913934572691/2.;
295 quadWeights[3]=.467913934572691/2.;
296 quadWeights[4]=.360761573048139/2.;
297 quadWeights[5]=.171324492379170/2.;
298 break;
299
300 case 7:
301 quadNodes[0]=(1.-.949107912342759)/2.;
302 quadNodes[1]=(1.-.741531185599394)/2.;
303 quadNodes[2]=(1.-.405845151377397)/2.;
304 quadNodes[3]=0.5;
305 quadNodes[4]=(1.+.405845151377397)/2.;
306 quadNodes[5]=(1.+.741531185599394)/2.;
307 quadNodes[6]=(1.+.949107912342759)/2.;
308 quadWeights[0]= .129484966168870/2.;
309 quadWeights[1]= .279705391489277/2.;
310 quadWeights[2]= .381830050505119/2.;
311 quadWeights[3]= .417959183673469/2.;
312 quadWeights[4]= .381830050505119/2.;
313 quadWeights[5]= .279705391489277/2.;
314 quadWeights[6]= .129484966168870/2.;
315 break;
316
317 case 8:
318 quadNodes[0]=(1.-.960289856497536)/2.;
319 quadNodes[1]=(1.-.796666477413627)/2.;
320 quadNodes[2]=(1.-.525532409916329)/2.;
321 quadNodes[3]=(1.-.183434642495650)/2.;
322 quadNodes[4]=(1.+.183434642495650)/2.;
323 quadNodes[5]=(1.+.525532409916329)/2.;
324 quadNodes[6]=(1.+.796666477413627)/2.;
325 quadNodes[7]=(1.+.960289856497536)/2.;
326 quadWeights[0]= .101228536290376/2.;
327 quadWeights[1]= .222381034453374/2.;
328 quadWeights[2]= .313706645877887/2.;
329 quadWeights[3]= .362683783378362/2.;
330 quadWeights[4]= .362683783378362/2.;
331 quadWeights[5]= .313706645877887/2.;
332 quadWeights[6]= .222381034453374/2.;
333 quadWeights[7]= .101228536290376/2.;
334 break;
335
336 case 9:
337 quadNodes[0]=(1.-.968160239507626)/2.;
338 quadNodes[1]=(1.-.836031107326636)/2.;
339 quadNodes[2]=(1.-.613371432700590)/2.;
340 quadNodes[3]=(1.-.324253423403809)/2.;
341 quadNodes[4]= .5;
342 quadNodes[5]=(1.+.324253423403809)/2.;
343 quadNodes[6]=(1.+.613371432700590)/2.;
344 quadNodes[7]=(1.+.836031107326636)/2.;
345 quadNodes[8]=(1.+.968160239507626)/2.;
346 quadWeights[0]= .081274388361574/2.;
347 quadWeights[1]= .180648160694857/2.;
348 quadWeights[2]= .260610696402935/2.;
349 quadWeights[3]= .312347077040003/2.;
350 quadWeights[4]= .330239355001260/2.;
351 quadWeights[5]= .312347077040003/2.;
352 quadWeights[6]= .260610696402935/2.;
353 quadWeights[7]= .180648160694857/2.;
354 quadWeights[8]= .081274388361574/2.;
355 break;
356
357 case 10:
358 quadNodes[0]=(1.-.973906528517172)/2.;
359 quadNodes[1]=(1.-.865063366688985)/2.;
360 quadNodes[2]=(1.-.679409568299024)/2.;
361 quadNodes[3]=(1.-.433395394129247)/2.;
362 quadNodes[4]=(1.-.148874338981631)/2.;
363 quadNodes[5]=(1.+.148874338981631)/2.;
364 quadNodes[6]=(1.+.433395394129247)/2.;
365 quadNodes[7]=(1.+.679409568299024)/2.;
366 quadNodes[8]=(1.+.865063366688985)/2.;
367 quadNodes[9]=(1.+.973906528517172)/2.;
368 quadWeights[0]= .066671344308688/2.;
369 quadWeights[1]= .149451349150581/2.;
370 quadWeights[2]= .219086362515982/2.;
371 quadWeights[3]= .269266719309996/2.;
372 quadWeights[4]= .295524224714753/2.;
373 quadWeights[5]= .295524224714753/2.;
374 quadWeights[6]= .269266719309996/2.;
375 quadWeights[7]= .219086362515982/2.;
376 quadWeights[8]= .149451349150581/2.;
377 quadWeights[9]= .066671344308688/2.;
378 break;
379
380 default:
381 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
382 break;
383 }
384 }
385 /**************************************************************/
386 /* */
387 /* the following function are used define the meshes on the surface in the xy-plane */
388
389 /* triangle surface on a tetrahedron */
390 void Finley_Quad_getNodesTriOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
391 Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesTri);
392 }
393 /* rectangular surface on a hexahedron */
394 void Finley_Quad_getNodesRecOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
395 Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesRec);
396 }
397 /* line surface on a triangle or rectangle */
398 void Finley_Quad_getNodesLineOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
399 Finley_Quad_makeNodesOnFace(2,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesLine);
400 }
401 /* point surface on a line */
402 void Finley_Quad_getNodesPointOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
403 Finley_Quad_makeNodesOnFace(1,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesPoint);
404 }
405
406 void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {
407 int q,i;
408 /* Win32 Refactor */
409 double *quadNodesOnFace = (numQuadNodes*(dim-1)>0) ? TMPMEMALLOC(numQuadNodes*(dim-1),double) : (double*)NULL;
410
411 #define DIM dim
412 getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);
413 if (! Finley_noError()) return;
414
415 for (q=0;q<numQuadNodes;q++) {
416 for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];
417 QUADNODES(dim-1,q)=0;
418 }
419 #undef DIM
420 TMPMEMFREE(quadNodesOnFace);
421 }
422
423 /**************************************************************/
424
425 /* The following functions Finley_Quad_getNumNodes* return the nmber of quadrature points needed to */
426 /* achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
427 /* to consider the accuracy reduction through the construction process. */
428
429
430 int Finley_Quad_getNumNodesPoint(int order) {
431 if (order <0 ) {
432 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
433 return -1;
434 } else {
435 return 0;
436 }
437 }
438
439 int Finley_Quad_getNumNodesLine(int order) {
440 char error_msg[LenErrorMsg_MAX];
441 if (order <0 ) {
442 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
443 return -1;
444 } else {
445 if (order > 2*MAX_numQuadNodesLine-1) {
446 sprintf(error_msg,"__FILE__: requested integration order %d on line is too large (>%d).",
447 order,2*MAX_numQuadNodesLine-1);
448 Finley_setError(VALUE_ERROR,error_msg);
449 return -1;
450 } else {
451 Finley_resetError();
452 return order/2+1;
453 }
454 }
455 }
456
457 int Finley_Quad_getNumNodesTri(int order) {
458 int numQuadNodesLine;
459 if (order==1) {
460 return 1;
461 } else {
462 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
463 if (Finley_noError()) {
464 return numQuadNodesLine*numQuadNodesLine;
465 } else {
466 return -1;
467 }
468 }
469 }
470
471 int Finley_Quad_getNumNodesRec(int order) {
472 int numQuadNodesLine;
473 numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
474 if (Finley_noError()) {
475 return numQuadNodesLine*numQuadNodesLine;
476 } else {
477 return -1;
478 }
479 }
480
481 int Finley_Quad_getNumNodesTet(int order) {
482 int numQuadNodesLine;
483 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
484 if (Finley_noError()) {
485 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
486 } else {
487 return -1;
488 }
489 }
490
491 int Finley_Quad_getNumNodesHex(int order) {
492 int numQuadNodesLine;
493 numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
494 if (Finley_noError()) {
495 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
496 } else {
497 return -1;
498 }
499 }
500 /*
501 * $Log$
502 * Revision 1.2 2005/09/15 03:44:23 jgs
503 * Merge of development branch dev-02 back to main trunk on 2005-09-15
504 *
505 * Revision 1.1.1.1.6.1 2005/09/07 06:26:20 gross
506 * the solver from finley are put into the standalone package paso now
507 *
508 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
509 * initial import of project esys2
510 *
511 * Revision 1.2 2004/08/03 04:49:06 gross
512 * bug in Quadrature.c fixed
513 *
514 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
515 * Initial version of eys using boost-python.
516 *
517 *
518 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26