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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


1 /*
2 ************************************************************
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 */
12
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 if (! Finley_noError()) return;
51
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 if (! Finley_noError()) return;
88
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 char error_msg[LenErrorMsg_MAX];
122 int numQuadNodes1d,i,j,l;
123 double *quadNodes1d=NULL,*quadWeights1d=NULL;
124 bool_t set=FALSE;
125 #define DIM 2
126
127 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
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 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 }
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 char error_msg[LenErrorMsg_MAX];
173 int numQuadNodes1d,i,j,k,l;
174 double *quadNodes1d=NULL,*quadWeights1d=NULL;
175 bool_t set;
176 #define DIM 3
177
178 /* find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
179
180 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
186 /* get 1D scheme: */
187
188 Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
189
190 /* make 3D scheme: */
191
192 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 }
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 Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: There is no quadrature scheme on points.");
229 }
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 Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
381 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 double *quadNodesOnFace=NULL;
408 #define DIM dim
409 quadNodesOnFace=TMPMEMALLOC(numQuadNodes*(dim-1),double);
410
411 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 }
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 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
434 return -1;
435 } else {
436 return 0;
437 }
438 }
439
440 int Finley_Quad_getNumNodesLine(int order) {
441 char error_msg[LenErrorMsg_MAX];
442 if (order <0 ) {
443 Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");
444 return -1;
445 } else {
446 if (order > 2*MAX_numQuadNodesLine-1) {
447 sprintf(error_msg,"Finley_Quad_getNumNodesPoint: requested integration order %d on line is too large (>%d).",
448 order,2*MAX_numQuadNodesLine-1);
449 Finley_setError(VALUE_ERROR,error_msg);
450 return -1;
451 } else {
452 Finley_resetError();
453 return order/2+1;
454 }
455 }
456 }
457
458 int Finley_Quad_getNumNodesTri(int order) {
459 int numQuadNodesLine;
460 if (order<=1) {
461 return 1;
462 } else {
463 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
464 if (Finley_noError()) {
465 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 if (Finley_noError()) {
476 return numQuadNodesLine*numQuadNodesLine;
477 } else {
478 return -1;
479 }
480 }
481
482 int Finley_Quad_getNumNodesTet(int order) {
483 int numQuadNodesLine;
484 if (order<=1) {
485 return 1;
486 } else {
487 numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
488 if (Finley_noError()) {
489 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
490 } else {
491 return -1;
492 }
493 }
494 }
495
496 int Finley_Quad_getNumNodesHex(int order) {
497 int numQuadNodesLine;
498 numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
499 if (Finley_noError()) {
500 return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
501 } else {
502 return -1;
503 }
504 }
505 /*
506 * $Log$
507 * 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 *
510 * 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 * 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