ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/outil/src/ot_quadrature_gauss.cpp
Revision: 754
Committed: Thu Oct 29 21:49:15 2015 UTC (9 years, 8 months ago) by francois
Original Path: magic/lib/outil/src/ot_quadrature_gauss.cpp
File size: 15153 byte(s)
Log Message:
Centralisation et unification de la définition des points de gauss des elements finis et arrangement dans le meme ordre que Aster

File Contents

# User Rev Content
1 francois 283 //---------------------------------------------------------------------------
2    
3     #include <math.h>
4     #include <stdio.h>
5     #define EPS 3.0e-11
6    
7     #pragma hdrstop
8    
9     #include "gestionversion.h"
10     #include "ot_quadrature_gauss.h"
11     //---------------------------------------------------------------------------
12     #pragma package(smart_init)
13    
14     void OT_QUADRATURE_GAUSS::gauss_legendre_points(double x1, double x2, double x[], double w[], int n)
15     {
16     int m,j,i;
17     double z1,z,xm,xl,pp,p3,p2,p1;
18     m=(n+1)/2;
19     xm=0.5*(x2+x1);
20     xl=0.5*(x2-x1);
21     for (i=1;i<=m;i++) {
22     z=cos(M_PI*(i-0.25)/(n+0.5));
23     do {
24     p1=1;
25     p2=0;
26     for (j=1;j<=n;j++) {
27     p3=p2;
28     p2=p1;
29     p1=((2*j-1)*z*p2-(j-1)*p3)/j;
30     }
31     pp=n*(z*p1-p2)/(z*z-1);
32     z1=z;
33     z=z1-p1/pp;
34     } while (fabs(z-z1) > EPS);
35     x[i]=xm-xl*z;
36     x[n+1-i]=xm+xl*z;
37     w[i]=2*xl/((1-z*z)*pp*pp);
38     w[n+1-i]=w[i];
39     }
40     }
41 francois 754
42    
43    
44    
45    
46     int OT_POINTS_GAUSS::get_nb_point_seg(int degre)
47     {
48     if (degre<2) return 1;
49     if (degre<4) return 2;
50     if (degre<6) return 3;
51     if (degre<8) return 4;
52     if (degre<10) return 5;
53     if (degre<12) return 6;
54     if (degre<14) return 7;
55     return 0;
56     }
57    
58     int OT_POINTS_GAUSS::get_nb_point_tri(int degre)
59     {
60     if (degre<2) return 1;
61     if (degre<3) return 3;
62     if (degre<4) return 4;
63     if (degre<5) return 6;
64     if (degre<6) return 7;
65     return 0; }
66    
67     int OT_POINTS_GAUSS::get_nb_point_quad(int degre)
68     {
69     if (degre<3) return 3;
70     if (degre<4) return 4;
71     if (degre<6) return 7;
72     return 0;
73     }
74     int OT_POINTS_GAUSS::get_nb_point_quad_prod(int degre)
75     {
76     if (degre<2) return 1;
77     if (degre<4) return 4;
78     if (degre<6) return 9;
79     return 0;
80     }
81     int OT_POINTS_GAUSS::get_nb_point_tetra(int degre)
82     {
83     if (degre<2) return 1;
84     if (degre<3) return 4;
85     if (degre<4) return 5;
86     if (degre<6) return 15;
87     return 0; }
88    
89     int OT_POINTS_GAUSS::get_nb_point_hexa(int degre)
90     {
91     if (degre<3) return 4;
92     if (degre<4) return 6;
93     if (degre<6) return 14;
94     return 0;
95     }
96    
97     int OT_POINTS_GAUSS::get_nb_point_hexa_prod(int degre)
98     {
99     if (degre<2) return 1;
100     if (degre<4) return 8;
101     if (degre<6) return 27;
102     return 0;
103     }
104    
105     void OT_POINTS_GAUSS::get_pt_gauss_seg(int degre,int num,double &w,double *u)
106     {
107     if (degre<2)
108     {
109     if (num==0) {w=2;u[0]=0;}
110     return;
111     }
112     if (degre<4)
113     {
114     if (num==0) {w=1;u[0]=0.577350269189626;}
115     if (num==1) {w=1;u[0]=-0.577350269189626;}
116     return;
117     }
118     if (degre<6)
119     {
120     if (num==0) {w=0.555555555555556;u[0]=-0.774596669241483;}
121     if (num==1) {w=0.888888888888889;u[0]=0.;}
122     if (num==2) {w=0.555555555555556;u[0]=0.774596669241483;}
123     return;
124     }
125     if (degre<8)
126     {
127     if (num==0) {w=0.652145154862546;u[0]=0.339981043584856;}
128     if (num==1) {w=0.652145154862546;u[0]=-0.339981043584856;}
129     if (num==2) {w=0.347854845137454;u[0]=0.861136311594053;}
130     if (num==3) {w=0.347854845137454;u[0]=-0.861136311594053;}
131     return;
132     }
133     // limite conformite code aster
134     if (degre<10)
135     {
136     if (num==0) {w=0.568888888888889;u[0]=0.;}
137     if (num==1) {w=0.478628670499366;u[0]=-0.538469310105683;}
138     if (num==2) {w=0.478628670499366;u[0]=0.538469310105683;}
139     if (num==3) {w=0.236926885056189;u[0]=-0.906179845938664;}
140     if (num==4) {w=0.236926885056189;u[0]=0.906179845938664;}
141     return;
142     }
143     if (degre<12)
144     {
145     if (num==0) {w=0.467913934572691;u[0]=-0.238619186083197;}
146     if (num==1) {w=0.467913934572691;u[0]=0.238619186083197;}
147     if (num==2) {w=0.360761573048139;u[0]=-0.661209386466265;}
148     if (num==3) {w=0.360761573048139;u[0]=0.661209386466265;}
149     if (num==4) {w=0.171324492379170;u[0]=-0.932469514203152;}
150     if (num==5) {w=0.171324492379170;u[0]=0.932469514203152;}
151     return;
152     }
153     if (degre<14)
154     {
155     if (num==0) {w=0.417959183673469;u[0]=0;}
156     if (num==1) {w=0.381830050505119;u[0]=-0.405845151377397;}
157     if (num==2) {w=0.381830050505119;u[0]=0.405845151377397;}
158     if (num==3) {w=0.279705391489277;u[0]=-0.741531185599394;}
159     if (num==4) {w=0.279705391489277;u[0]=0.741531185599394;}
160     if (num==5) {w=0.129484966168870;u[0]=-0.949107912342759;}
161     if (num==6) {w=0.129484966168870;u[0]=0.949107912342759;}
162     return;
163     }
164     }
165    
166    
167     void OT_POINTS_GAUSS::get_pt_gauss_tri(int degre, int num, double& w, double* uv)
168     {
169     if (degre<2)
170     {
171     if (num==0) {w=0.5;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
172     return;
173     }
174     if (degre<3)
175     {
176     if (num==0) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.166666666666667;}
177     if (num==1) {w=0.166666666666667;uv[0]=0.666666666666667;uv[1]=0.166666666666667;}
178     if (num==2) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.666666666666667;}
179     return;
180     }
181     if (degre<4)
182     {
183     if (num==0) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.2;}
184     if (num==1) {w=0.260416666666667;uv[0]=0.6;uv[1]=0.2;}
185     if (num==2) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.6;}
186     if (num==3) {w=-0.28125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
187     return;
188     }
189     if (degre<5)
190     {
191     if (num==0) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.091576213509771;}
192     if (num==1) {w=0.054975871827661;uv[0]=0.816847572980458;uv[1]=0.091576213509771;}
193     if (num==2) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.816847572980458;}
194     if (num==3) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.445948490915965;}
195     if (num==4) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.108103018168070;}
196     if (num==5) {w=0.111690794839005;uv[0]=0.108103018168070;uv[1]=0.445948490915965;}
197     return;
198     }
199     if (degre<6)
200     {
201     if (num==0) {w=0.1125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
202     if (num==1) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.470142064105115;}
203     if (num==2) {w=0.066197076394253;uv[0]=0.059715871789770;uv[1]=0.470142064105115;}
204     if (num==3) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.059715871789770;}
205     if (num==4) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.101286507323456;}
206     if (num==5) {w=0.062969590272413;uv[0]=0.797426985353088;uv[1]=0.101286507323456;}
207     if (num==6) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.797426985353088;}
208     return;
209     }
210     }
211    
212     void OT_POINTS_GAUSS::get_pt_gauss_tet(int degre, int num, double& w, double* uvw)
213     {
214     if (degre<2)
215     {
216     if (num==0) {w=0.166666666666667;uvw[0]=0.25;uvw[1]=0.25;uvw[2]=0.25;}
217     return;
218     }
219     if (degre<3)
220     {
221     if (num==0) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.138196601125011;uvw[2]=0.138196601125011;}
222     if (num==1) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.138196601125011;uvw[2]=0.585410196624968;}
223     if (num==2) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.585410196624968;uvw[2]=0.138196601125011;}
224     if (num==3) {w=0.041666666666667;uvw[0]=0.585410196624968;uvw[1]=0.138196601125011;uvw[2]=0.138196601125011;}
225     return;
226     }
227     if (degre<4)
228     {
229     if (num==0) {w=-0.133333333333333;uvw[0]=0.25;uvw[1]=0.25;uvw[2]=0.25;}
230     if (num==1) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.166666666666667;uvw[2]=0.166666666666667;}
231     if (num==2) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.166666666666667;uvw[2]=0.5;}
232     if (num==3) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.5;uvw[2]=0.166666666666667;}
233     if (num==4) {w=0.075;uvw[0]=0.5;uvw[1]=0.166666666666667;uvw[2]=0.166666666666667;}
234     return;
235     }
236     if (degre<6)
237     {
238     double a=0.25;
239     double b1=0.3197936278;
240     double b2=0.0919710781;
241     double c2=0.7240867658;
242     double c1=0.0406191165;
243     double d=0.0563508327;
244     double e=0.4436491673;
245     double w1=0.0197530864;
246     double w3=0.011989514;
247     double w2=0.0115113679;
248     double w4=0.0088183422;
249     if (num==0) {w=w1;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
250     if (num==1) {w=w2;uvw[0]=b1;uvw[1]=b1;uvw[2]=b1;}
251     if (num==2) {w=w2;uvw[0]=b1;uvw[1]=b1;uvw[2]=c1;}
252     if (num==3) {w=w2;uvw[0]=b1;uvw[1]=c1;uvw[2]=b1;}
253     if (num==4) {w=w2;uvw[0]=c1;uvw[1]=b1;uvw[2]=b1;}
254     if (num==5) {w=w3;uvw[0]=b2;uvw[1]=b2;uvw[2]=b2;}
255     if (num==6) {w=w3;uvw[0]=b2;uvw[1]=b2;uvw[2]=c2;}
256     if (num==7) {w=w3;uvw[0]=b2;uvw[1]=c2;uvw[2]=b2;}
257     if (num==8) {w=w3;uvw[0]=c2;uvw[1]=b2;uvw[2]=b2;}
258     if (num==9) {w=w4;uvw[0]=d;uvw[1]=d;uvw[2]=e;}
259     if (num==10) {w=w4;uvw[0]=d;uvw[1]=e;uvw[2]=d;}
260     if (num==11) {w=w4;uvw[0]=e;uvw[1]=d;uvw[2]=d;}
261     if (num==12) {w=w4;uvw[0]=d;uvw[1]=e;uvw[2]=e;}
262     if (num==13) {w=w4;uvw[0]=e;uvw[1]=d;uvw[2]=e;}
263     if (num==14) {w=w4;uvw[0]=e;uvw[1]=e;uvw[2]=d;}
264     return;
265     }
266     }
267    
268     void OT_POINTS_GAUSS::get_pt_gauss_qua(int degre, int num, double& w, double* uv)
269     {
270     if (degre<3)
271     {
272     if (num==0) {w=1.333333333333333;uv[0]=0.816496580927726;uv[1]=0.;}
273     if (num==1) {w=1.333333333333333;uv[0]=-0.408248290463863;uv[1]=-0.707106781186548;}
274     if (num==2) {w=1.333333333333333;uv[0]=-0.408248290463863;uv[1]=0.707106781186548;}
275     return;
276     }
277     if (degre<4)
278     {
279     if (num==0) {w=1.;uv[0]=0.816496580927726;uv[1]=0.;}
280     if (num==1) {w=1.;uv[0]=-0.816496580927726;uv[1]=0.;}
281     if (num==2) {w=1.;uv[0]=0.;uv[1]=0.816496580927726;}
282     if (num==3) {w=1.;uv[0]=0.;uv[1]=-0.816496580927726;}
283     return;
284     }
285     if (degre<6)
286     {
287     if (num==0) {w=1.142857142857140;uv[0]=0.;uv[1]=0.;}
288     if (num==1) {w=0.317460317460317;uv[0]=0.;uv[1]=0.966091783079296;}
289     if (num==2) {w=0.317460317460317;uv[0]=0.;uv[1]=-0.966091783079296;}
290     if (num==3) {w=0.555555555555556;uv[0]=0.774596669241483;uv[1]=-0.577350269189626;}
291     if (num==4) {w=0.555555555555556;uv[0]=-0.774596669241483;uv[1]=-0.577350269189626;}
292     if (num==5) {w=0.555555555555556;uv[0]=-0.774596669241483;uv[1]=0.577350269189626;}
293     if (num==6) {w=0.555555555555556;uv[0]=0.774596669241483;uv[1]=0.577350269189626;}
294     return;
295     }
296     }
297    
298     void OT_POINTS_GAUSS::get_pt_gauss_qua_prod(int degre, int num, double& w, double* uv)
299     {
300     if (degre<2)
301     {
302     if (num==0) {w=4.0;uv[0]=0.;uv[1]=0.;}
303     return;
304     }
305     if (degre<4)
306     {
307     if (num==0) {w=1.;uv[0]=-0.577350269189626;uv[1]=-0.577350269189626;}
308     if (num==0) {w=1.;uv[0]=0.577350269189626;uv[1]=-0.577350269189626;}
309     if (num==0) {w=1.;uv[0]=0.577350269189626;uv[1]=0.577350269189626;}
310     if (num==0) {w=1.;uv[0]=-0.577350269189626;uv[1]=0.577350269189626;}
311     return;
312     }
313     if (degre<6)
314     {
315     double a=0.774596669241483;
316     if (num==0) {w=25./81.;uv[0]=-a;uv[1]=-a;}
317     if (num==1) {w=25./81.;uv[0]=a;uv[1]=-a;}
318     if (num==2) {w=25./81.;uv[0]=a;uv[1]=a;}
319     if (num==3) {w=25./81.;uv[0]=-a;uv[1]=a;}
320     if (num==4) {w=40./81.;uv[0]=0.;uv[1]=-a;}
321     if (num==5) {w=40./81.;uv[0]=a;uv[1]=0.;}
322     if (num==6) {w=40./81.;uv[0]=0.;uv[1]=a;}
323     if (num==7) {w=40./81.;uv[0]=-a;uv[1]=0.;}
324     if (num==8) {w=64./81.;uv[0]=0.;uv[1]=0.;}
325     return;
326     }
327     }
328    
329    
330     void OT_POINTS_GAUSS::get_pt_gauss_hex_prod(int degre, int num, double& w, double* uvw)
331     {
332     if (degre<2)
333     {
334     if (num==0) {w=8.0;uvw[0]=0.;uvw[1]=0.;uvw[2]=0.;}
335     return;
336     }
337     if (degre<4)
338     {
339     double a=0.577350269189626;
340     if (num==0) {w=1.0;uvw[0]=-a;uvw[1]=-a;uvw[2]=-a;}
341     if (num==1) {w=1.0;uvw[0]=-a;uvw[1]=-a;uvw[2]=a;}
342     if (num==2) {w=1.0;uvw[0]=-a;uvw[1]=a;uvw[2]=-a;}
343     if (num==3) {w=1.0;uvw[0]=-a;uvw[1]=a;uvw[2]=a;}
344     if (num==4) {w=1.0;uvw[0]=a;uvw[1]=-a;uvw[2]=-a;}
345     if (num==5) {w=1.0;uvw[0]=a;uvw[1]=-a;uvw[2]=a;}
346     if (num==6) {w=1.0;uvw[0]=a;uvw[1]=a;uvw[2]=-a;}
347     if (num==7) {w=1.0;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
348     return;
349     }
350     if (degre<6)
351     {
352     double a=0.774596669241483;
353     double c1=0.555555555555556;
354     double c2=0.888888888888889;
355     if (num==0) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=-a;uvw[2]=-a;}
356     if (num==1) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=-a;uvw[2]=0.;}
357     if (num==2) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=-a;uvw[2]=a;}
358     if (num==3) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=0.;uvw[2]=-a;}
359     if (num==4) {w=c1*c2*c2;uvw[0]=-a;uvw[1]=0.;uvw[2]=0.;}
360     if (num==5) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=0.;uvw[2]=a;}
361     if (num==6) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=a;uvw[2]=-a;}
362     if (num==7) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=a;uvw[2]=0.;}
363     if (num==8) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=a;uvw[2]=a;}
364     if (num==9) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=-a;uvw[2]=-a;}
365     if (num==10) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=-a;uvw[2]=0.;}
366     if (num==11) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=-a;uvw[2]=a;}
367     if (num==12) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=0.;uvw[2]=-a;}
368     if (num==13) {w=c2*c2*c2;uvw[0]=0.;uvw[1]=0.;uvw[2]=0.;}
369     if (num==14) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=0.;uvw[2]=a;}
370     if (num==15) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=a;uvw[2]=-a;}
371     if (num==16) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=a;uvw[2]=0.;}
372     if (num==17) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=a;uvw[2]=a;}
373     if (num==18) {w=c1*c1*c1;uvw[0]=a;uvw[1]=-a;uvw[2]=-a;}
374     if (num==19) {w=c1*c1*c2;uvw[0]=a;uvw[1]=-a;uvw[2]=0.;}
375     if (num==20) {w=c1*c1*c1;uvw[0]=a;uvw[1]=-a;uvw[2]=a;}
376     if (num==21) {w=c1*c1*c2;uvw[0]=a;uvw[1]=0.;uvw[2]=-a;}
377     if (num==22) {w=c1*c2*c2;uvw[0]=a;uvw[1]=0.;uvw[2]=0.;}
378     if (num==23) {w=c1*c1*c2;uvw[0]=a;uvw[1]=0.;uvw[2]=a;}
379     if (num==24) {w=c1*c1*c1;uvw[0]=a;uvw[1]=a;uvw[2]=-a;}
380     if (num==25) {w=c1*c1*c2;uvw[0]=a;uvw[1]=a;uvw[2]=0.;}
381     if (num==26) {w=c1*c1*c1;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
382     return;
383     }
384     }
385    
386    
387    
388     void OT_POINTS_GAUSS::get_pt_gauss_hex(int degre, int num, double& w, double* uvw)
389     {
390     if (degre<3)
391     {
392     if (num==0) {w=2.;uvw[0]=0.;uvw[1]=0.816496580927726;uvw[2]=-0,577350269189626;}
393     if (num==1) {w=2.;uvw[0]=0.;uvw[1]=-0.816496580927726;uvw[2]=-0,577350269189626;}
394     if (num==2) {w=2.;uvw[0]=-0.816496580927726;uvw[1]=0.;uvw[2]=0,577350269189626;}
395     if (num==2) {w=2.;uvw[0]=0.816496580927726;uvw[1]=0.;uvw[2]=0,577350269189626;}
396     return;
397     }
398     if (degre<4)
399     {
400     if (num==0) {w=1.333333333333333;uvw[0]=1.;uvw[1]=0.;uvw[2]=0.;}
401     if (num==1) {w=1.333333333333333;uvw[0]=-1.;uvw[1]=0.;uvw[2]=0.;}
402     if (num==2) {w=1.333333333333333;uvw[0]=0.;uvw[1]=1.;uvw[2]=0.;}
403     if (num==3) {w=1.333333333333333;uvw[0]=0.;uvw[1]=-1.;uvw[2]=0.;}
404     if (num==4) {w=1.333333333333333;uvw[0]=0.;uvw[1]=0.;uvw[2]=1.;}
405     if (num==5) {w=1.333333333333333;uvw[0]=0.;uvw[1]=0.;uvw[2]=-1.;}
406     return;
407     }
408     if (degre<6)
409     {
410     double a=0.795822425754221;
411     double b=0.758786910639328;
412     double w1=0.886426592797784;
413     double w2=0.335180055401662;
414     if (num==0) {w=w1;uvw[0]=a;uvw[1]=0.;uvw[2]=0.;}
415     if (num==1) {w=w1;uvw[0]=-a;uvw[1]=0.;uvw[2]=0.;}
416     if (num==2) {w=w1;uvw[0]=0.;uvw[1]=a;uvw[2]=0.;}
417     if (num==3) {w=w1;uvw[0]=0.;uvw[1]=-a;uvw[2]=0.;}
418     if (num==4) {w=w1;uvw[0]=0.;uvw[1]=0.;uvw[2]=a;}
419     if (num==5) {w=w1;uvw[0]=0.;uvw[1]=0.;uvw[2]=-a;}
420     if (num==6) {w=w2;uvw[0]=b;uvw[1]=b;uvw[2]=b;}
421     if (num==7) {w=w2;uvw[0]=b;uvw[1]=-b;uvw[2]=-b;}
422     if (num==8) {w=w2;uvw[0]=b;uvw[1]=b;uvw[2]=-b;}
423     if (num==9) {w=w2;uvw[0]=b;uvw[1]=-b;uvw[2]=b;}
424     if (num==10) {w=w2;uvw[0]=-b;uvw[1]=b;uvw[2]=b;}
425     if (num==11) {w=w2;uvw[0]=-b;uvw[1]=-b;uvw[2]=-b;}
426     if (num==12) {w=w2;uvw[0]=-b;uvw[1]=b;uvw[2]=-b;}
427     if (num==13) {w=w2;uvw[0]=-b;uvw[1]=-b;uvw[2]=b;}
428     return;
429     }
430    
431     }