ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/outil/src/ot_quadrature_gauss.cpp
Revision: 1019
Committed: Tue Jun 4 21:16:50 2019 UTC (6 years ago) by francois
File size: 19485 byte(s)
Log Message: 
restructuration de magic
outil est sorti de lib pour pouvoir etre utiliser en dehors de lib
template est merge avec outil
poly_occ et un sous projet de magic qui utilise le nouveau outil

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 "ot_quadrature_gauss.h"
10     //---------------------------------------------------------------------------
11     #pragma package(smart_init)
12    
13     void OT_QUADRATURE_GAUSS::gauss_legendre_points(double x1, double x2, double x[], double w[], int n)
14     {
15     int m,j,i;
16     double z1,z,xm,xl,pp,p3,p2,p1;
17     m=(n+1)/2;
18     xm=0.5*(x2+x1);
19     xl=0.5*(x2-x1);
20     for (i=1;i<=m;i++) {
21     z=cos(M_PI*(i-0.25)/(n+0.5));
22     do {
23     p1=1;
24     p2=0;
25     for (j=1;j<=n;j++) {
26     p3=p2;
27     p2=p1;
28     p1=((2*j-1)*z*p2-(j-1)*p3)/j;
29     }
30     pp=n*(z*p1-p2)/(z*z-1);
31     z1=z;
32     z=z1-p1/pp;
33     } while (fabs(z-z1) > EPS);
34     x[i]=xm-xl*z;
35     x[n+1-i]=xm+xl*z;
36     w[i]=2*xl/((1-z*z)*pp*pp);
37     w[n+1-i]=w[i];
38     }
39     }
40 francois 754
41    
42    
43    
44    
45     int OT_POINTS_GAUSS::get_nb_point_seg(int degre)
46     {
47     if (degre<2) return 1;
48     if (degre<4) return 2;
49     if (degre<6) return 3;
50     if (degre<8) return 4;
51     if (degre<10) return 5;
52     if (degre<12) return 6;
53     if (degre<14) return 7;
54     return 0;
55     }
56    
57 francois 757
58    
59     int OT_POINTS_GAUSS::get_degre_gauss_seg(int num)
60     {
61     if (num<2) return 1;
62     if (num<3) return 3;
63     if (num<4) return 7;
64     if (num<5) return 9;
65     if (num<6) return 11;
66     if (num<7) return 13;
67     if (num<8) return 15;
68     return 0;
69     }
70    
71    
72    
73 francois 754 int OT_POINTS_GAUSS::get_nb_point_tri(int degre)
74     {
75     if (degre<2) return 1;
76     if (degre<3) return 3;
77     if (degre<4) return 4;
78     if (degre<5) return 6;
79     if (degre<6) return 7;
80 francois 757 return 0;
81     }
82 francois 754
83 francois 757 int OT_POINTS_GAUSS::get_degre_gauss_tri(int num)
84     {
85     if (num<2) return 1;
86     if (num<4) return 2;
87     if (num<6) return 3;
88     if (num<7) return 4;
89     if (num<8) return 5;
90     return 0;
91     }
92    
93 francois 754 int OT_POINTS_GAUSS::get_nb_point_quad(int degre)
94     {
95     if (degre<3) return 3;
96     if (degre<4) return 4;
97     if (degre<6) return 7;
98     return 0;
99     }
100 francois 757
101     int OT_POINTS_GAUSS::get_degre_gauss_quad(int num)
102     {
103     if (num<4) return 2;
104     if (num<5) return 3;
105     if (num<8) return 5;
106     return 0;
107     }
108    
109    
110 francois 754 int OT_POINTS_GAUSS::get_nb_point_quad_prod(int degre)
111     {
112     if (degre<2) return 1;
113     if (degre<4) return 4;
114     if (degre<6) return 9;
115     return 0;
116     }
117 francois 757
118     int OT_POINTS_GAUSS::get_degre_gauss_quad_prod(int num)
119     {
120     if (num<2) return 1;
121     if (num<5) return 3;
122     if (num<10) return 5;
123     return 0;
124     }
125    
126 francois 754 int OT_POINTS_GAUSS::get_nb_point_tetra(int degre)
127     {
128     if (degre<2) return 1;
129     if (degre<3) return 4;
130     if (degre<4) return 5;
131     if (degre<6) return 15;
132 francois 757 return 0;
133     }
134 francois 754
135 francois 757 int OT_POINTS_GAUSS::get_degre_gauss_tetra(int num)
136     {
137     if (num<2) return 1;
138     if (num<5) return 2;
139     if (num<6) return 3;
140     if (num<16) return 5;
141     return 0;
142     }
143    
144 francois 754 int OT_POINTS_GAUSS::get_nb_point_hexa(int degre)
145     {
146     if (degre<3) return 4;
147     if (degre<4) return 6;
148     if (degre<6) return 14;
149     return 0;
150     }
151    
152 francois 757 int OT_POINTS_GAUSS::get_degre_gauss_hexa(int num)
153     {
154     if (num<5) return 2;
155     if (num<7) return 3;
156     if (num<15) return 5;
157     return 0;
158     }
159    
160 francois 754 int OT_POINTS_GAUSS::get_nb_point_hexa_prod(int degre)
161     {
162     if (degre<2) return 1;
163     if (degre<4) return 8;
164     if (degre<6) return 27;
165     return 0;
166     }
167    
168 francois 757 int OT_POINTS_GAUSS::get_degre_gauss_hexa_prod(int num)
169     {
170     if (num<2) return 1;
171     if (num<9) return 3;
172     if (num<28) return 5;
173     return 0;
174     }
175    
176 francois 876 int OT_POINTS_GAUSS::get_nb_point_penta(int degre)
177     {
178     if (degre<3) return 6;
179     if (degre<4) return 8;
180     if (degre<6) return 21;
181     return 0;
182     }
183 francois 757
184 francois 876 int OT_POINTS_GAUSS::get_degre_gauss_penta(int num)
185     {
186     if (num<7) return 2;
187     if (num<9) return 3;
188     if (num<22) return 5;
189     return 0;
190     }
191 francois 757
192    
193 francois 876
194 francois 754 void OT_POINTS_GAUSS::get_pt_gauss_seg(int degre,int num,double &w,double *u)
195     {
196     if (degre<2)
197     {
198     if (num==0) {w=2;u[0]=0;}
199     return;
200     }
201     if (degre<4)
202     {
203     if (num==0) {w=1;u[0]=0.577350269189626;}
204     if (num==1) {w=1;u[0]=-0.577350269189626;}
205     return;
206     }
207     if (degre<6)
208     {
209     if (num==0) {w=0.555555555555556;u[0]=-0.774596669241483;}
210     if (num==1) {w=0.888888888888889;u[0]=0.;}
211     if (num==2) {w=0.555555555555556;u[0]=0.774596669241483;}
212     return;
213     }
214     if (degre<8)
215     {
216     if (num==0) {w=0.652145154862546;u[0]=0.339981043584856;}
217     if (num==1) {w=0.652145154862546;u[0]=-0.339981043584856;}
218     if (num==2) {w=0.347854845137454;u[0]=0.861136311594053;}
219     if (num==3) {w=0.347854845137454;u[0]=-0.861136311594053;}
220     return;
221     }
222     // limite conformite code aster
223     if (degre<10)
224     {
225     if (num==0) {w=0.568888888888889;u[0]=0.;}
226     if (num==1) {w=0.478628670499366;u[0]=-0.538469310105683;}
227     if (num==2) {w=0.478628670499366;u[0]=0.538469310105683;}
228     if (num==3) {w=0.236926885056189;u[0]=-0.906179845938664;}
229     if (num==4) {w=0.236926885056189;u[0]=0.906179845938664;}
230     return;
231     }
232     if (degre<12)
233     {
234     if (num==0) {w=0.467913934572691;u[0]=-0.238619186083197;}
235     if (num==1) {w=0.467913934572691;u[0]=0.238619186083197;}
236     if (num==2) {w=0.360761573048139;u[0]=-0.661209386466265;}
237     if (num==3) {w=0.360761573048139;u[0]=0.661209386466265;}
238     if (num==4) {w=0.171324492379170;u[0]=-0.932469514203152;}
239     if (num==5) {w=0.171324492379170;u[0]=0.932469514203152;}
240     return;
241     }
242     if (degre<14)
243     {
244     if (num==0) {w=0.417959183673469;u[0]=0;}
245     if (num==1) {w=0.381830050505119;u[0]=-0.405845151377397;}
246     if (num==2) {w=0.381830050505119;u[0]=0.405845151377397;}
247     if (num==3) {w=0.279705391489277;u[0]=-0.741531185599394;}
248     if (num==4) {w=0.279705391489277;u[0]=0.741531185599394;}
249     if (num==5) {w=0.129484966168870;u[0]=-0.949107912342759;}
250     if (num==6) {w=0.129484966168870;u[0]=0.949107912342759;}
251     return;
252     }
253     }
254    
255    
256     void OT_POINTS_GAUSS::get_pt_gauss_tri(int degre, int num, double& w, double* uv)
257     {
258     if (degre<2)
259     {
260     if (num==0) {w=0.5;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
261     return;
262     }
263     if (degre<3)
264     {
265     if (num==0) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.166666666666667;}
266     if (num==1) {w=0.166666666666667;uv[0]=0.666666666666667;uv[1]=0.166666666666667;}
267     if (num==2) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.666666666666667;}
268     return;
269     }
270     if (degre<4)
271     {
272     if (num==0) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.2;}
273     if (num==1) {w=0.260416666666667;uv[0]=0.6;uv[1]=0.2;}
274     if (num==2) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.6;}
275     if (num==3) {w=-0.28125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
276     return;
277     }
278     if (degre<5)
279     {
280     if (num==0) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.091576213509771;}
281     if (num==1) {w=0.054975871827661;uv[0]=0.816847572980458;uv[1]=0.091576213509771;}
282     if (num==2) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.816847572980458;}
283     if (num==3) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.445948490915965;}
284     if (num==4) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.108103018168070;}
285     if (num==5) {w=0.111690794839005;uv[0]=0.108103018168070;uv[1]=0.445948490915965;}
286     return;
287     }
288     if (degre<6)
289     {
290     if (num==0) {w=0.1125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
291     if (num==1) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.470142064105115;}
292     if (num==2) {w=0.066197076394253;uv[0]=0.059715871789770;uv[1]=0.470142064105115;}
293     if (num==3) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.059715871789770;}
294     if (num==4) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.101286507323456;}
295     if (num==5) {w=0.062969590272413;uv[0]=0.797426985353088;uv[1]=0.101286507323456;}
296     if (num==6) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.797426985353088;}
297     return;
298     }
299     }
300    
301 francois 757
302    
303    
304    
305 francois 754 void OT_POINTS_GAUSS::get_pt_gauss_tet(int degre, int num, double& w, double* uvw)
306     {
307     if (degre<2)
308     {
309     if (num==0) {w=0.166666666666667;uvw[0]=0.25;uvw[1]=0.25;uvw[2]=0.25;}
310     return;
311     }
312     if (degre<3)
313     {
314     if (num==0) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.138196601125011;uvw[2]=0.138196601125011;}
315     if (num==1) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.138196601125011;uvw[2]=0.585410196624968;}
316     if (num==2) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.585410196624968;uvw[2]=0.138196601125011;}
317     if (num==3) {w=0.041666666666667;uvw[0]=0.585410196624968;uvw[1]=0.138196601125011;uvw[2]=0.138196601125011;}
318     return;
319     }
320     if (degre<4)
321     {
322     if (num==0) {w=-0.133333333333333;uvw[0]=0.25;uvw[1]=0.25;uvw[2]=0.25;}
323     if (num==1) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.166666666666667;uvw[2]=0.166666666666667;}
324     if (num==2) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.166666666666667;uvw[2]=0.5;}
325     if (num==3) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.5;uvw[2]=0.166666666666667;}
326     if (num==4) {w=0.075;uvw[0]=0.5;uvw[1]=0.166666666666667;uvw[2]=0.166666666666667;}
327     return;
328     }
329     if (degre<6)
330     {
331     double a=0.25;
332     double b1=0.3197936278;
333     double b2=0.0919710781;
334     double c2=0.7240867658;
335     double c1=0.0406191165;
336     double d=0.0563508327;
337     double e=0.4436491673;
338     double w1=0.0197530864;
339     double w3=0.011989514;
340     double w2=0.0115113679;
341     double w4=0.0088183422;
342     if (num==0) {w=w1;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
343     if (num==1) {w=w2;uvw[0]=b1;uvw[1]=b1;uvw[2]=b1;}
344     if (num==2) {w=w2;uvw[0]=b1;uvw[1]=b1;uvw[2]=c1;}
345     if (num==3) {w=w2;uvw[0]=b1;uvw[1]=c1;uvw[2]=b1;}
346     if (num==4) {w=w2;uvw[0]=c1;uvw[1]=b1;uvw[2]=b1;}
347     if (num==5) {w=w3;uvw[0]=b2;uvw[1]=b2;uvw[2]=b2;}
348     if (num==6) {w=w3;uvw[0]=b2;uvw[1]=b2;uvw[2]=c2;}
349     if (num==7) {w=w3;uvw[0]=b2;uvw[1]=c2;uvw[2]=b2;}
350     if (num==8) {w=w3;uvw[0]=c2;uvw[1]=b2;uvw[2]=b2;}
351     if (num==9) {w=w4;uvw[0]=d;uvw[1]=d;uvw[2]=e;}
352     if (num==10) {w=w4;uvw[0]=d;uvw[1]=e;uvw[2]=d;}
353     if (num==11) {w=w4;uvw[0]=e;uvw[1]=d;uvw[2]=d;}
354     if (num==12) {w=w4;uvw[0]=d;uvw[1]=e;uvw[2]=e;}
355     if (num==13) {w=w4;uvw[0]=e;uvw[1]=d;uvw[2]=e;}
356     if (num==14) {w=w4;uvw[0]=e;uvw[1]=e;uvw[2]=d;}
357     return;
358     }
359     }
360    
361     void OT_POINTS_GAUSS::get_pt_gauss_qua(int degre, int num, double& w, double* uv)
362     {
363     if (degre<3)
364     {
365     if (num==0) {w=1.333333333333333;uv[0]=0.816496580927726;uv[1]=0.;}
366     if (num==1) {w=1.333333333333333;uv[0]=-0.408248290463863;uv[1]=-0.707106781186548;}
367     if (num==2) {w=1.333333333333333;uv[0]=-0.408248290463863;uv[1]=0.707106781186548;}
368     return;
369     }
370     if (degre<4)
371     {
372     if (num==0) {w=1.;uv[0]=0.816496580927726;uv[1]=0.;}
373     if (num==1) {w=1.;uv[0]=-0.816496580927726;uv[1]=0.;}
374     if (num==2) {w=1.;uv[0]=0.;uv[1]=0.816496580927726;}
375     if (num==3) {w=1.;uv[0]=0.;uv[1]=-0.816496580927726;}
376     return;
377     }
378     if (degre<6)
379     {
380     if (num==0) {w=1.142857142857140;uv[0]=0.;uv[1]=0.;}
381     if (num==1) {w=0.317460317460317;uv[0]=0.;uv[1]=0.966091783079296;}
382     if (num==2) {w=0.317460317460317;uv[0]=0.;uv[1]=-0.966091783079296;}
383     if (num==3) {w=0.555555555555556;uv[0]=0.774596669241483;uv[1]=-0.577350269189626;}
384     if (num==4) {w=0.555555555555556;uv[0]=-0.774596669241483;uv[1]=-0.577350269189626;}
385     if (num==5) {w=0.555555555555556;uv[0]=-0.774596669241483;uv[1]=0.577350269189626;}
386     if (num==6) {w=0.555555555555556;uv[0]=0.774596669241483;uv[1]=0.577350269189626;}
387     return;
388     }
389     }
390    
391     void OT_POINTS_GAUSS::get_pt_gauss_qua_prod(int degre, int num, double& w, double* uv)
392     {
393     if (degre<2)
394     {
395     if (num==0) {w=4.0;uv[0]=0.;uv[1]=0.;}
396     return;
397     }
398     if (degre<4)
399     {
400     if (num==0) {w=1.;uv[0]=-0.577350269189626;uv[1]=-0.577350269189626;}
401 francois 876 if (num==1) {w=1.;uv[0]=0.577350269189626;uv[1]=-0.577350269189626;}
402     if (num==2) {w=1.;uv[0]=0.577350269189626;uv[1]=0.577350269189626;}
403     if (num==3) {w=1.;uv[0]=-0.577350269189626;uv[1]=0.577350269189626;}
404 francois 754 return;
405     }
406     if (degre<6)
407     {
408     double a=0.774596669241483;
409     if (num==0) {w=25./81.;uv[0]=-a;uv[1]=-a;}
410     if (num==1) {w=25./81.;uv[0]=a;uv[1]=-a;}
411     if (num==2) {w=25./81.;uv[0]=a;uv[1]=a;}
412     if (num==3) {w=25./81.;uv[0]=-a;uv[1]=a;}
413     if (num==4) {w=40./81.;uv[0]=0.;uv[1]=-a;}
414     if (num==5) {w=40./81.;uv[0]=a;uv[1]=0.;}
415     if (num==6) {w=40./81.;uv[0]=0.;uv[1]=a;}
416     if (num==7) {w=40./81.;uv[0]=-a;uv[1]=0.;}
417     if (num==8) {w=64./81.;uv[0]=0.;uv[1]=0.;}
418     return;
419     }
420     }
421    
422    
423     void OT_POINTS_GAUSS::get_pt_gauss_hex_prod(int degre, int num, double& w, double* uvw)
424     {
425     if (degre<2)
426     {
427     if (num==0) {w=8.0;uvw[0]=0.;uvw[1]=0.;uvw[2]=0.;}
428     return;
429     }
430     if (degre<4)
431     {
432     double a=0.577350269189626;
433     if (num==0) {w=1.0;uvw[0]=-a;uvw[1]=-a;uvw[2]=-a;}
434     if (num==1) {w=1.0;uvw[0]=-a;uvw[1]=-a;uvw[2]=a;}
435     if (num==2) {w=1.0;uvw[0]=-a;uvw[1]=a;uvw[2]=-a;}
436     if (num==3) {w=1.0;uvw[0]=-a;uvw[1]=a;uvw[2]=a;}
437     if (num==4) {w=1.0;uvw[0]=a;uvw[1]=-a;uvw[2]=-a;}
438     if (num==5) {w=1.0;uvw[0]=a;uvw[1]=-a;uvw[2]=a;}
439     if (num==6) {w=1.0;uvw[0]=a;uvw[1]=a;uvw[2]=-a;}
440     if (num==7) {w=1.0;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
441     return;
442     }
443     if (degre<6)
444     {
445     double a=0.774596669241483;
446     double c1=0.555555555555556;
447     double c2=0.888888888888889;
448     if (num==0) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=-a;uvw[2]=-a;}
449     if (num==1) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=-a;uvw[2]=0.;}
450     if (num==2) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=-a;uvw[2]=a;}
451     if (num==3) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=0.;uvw[2]=-a;}
452     if (num==4) {w=c1*c2*c2;uvw[0]=-a;uvw[1]=0.;uvw[2]=0.;}
453     if (num==5) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=0.;uvw[2]=a;}
454     if (num==6) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=a;uvw[2]=-a;}
455     if (num==7) {w=c1*c1*c2;uvw[0]=-a;uvw[1]=a;uvw[2]=0.;}
456     if (num==8) {w=c1*c1*c1;uvw[0]=-a;uvw[1]=a;uvw[2]=a;}
457     if (num==9) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=-a;uvw[2]=-a;}
458     if (num==10) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=-a;uvw[2]=0.;}
459     if (num==11) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=-a;uvw[2]=a;}
460     if (num==12) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=0.;uvw[2]=-a;}
461     if (num==13) {w=c2*c2*c2;uvw[0]=0.;uvw[1]=0.;uvw[2]=0.;}
462     if (num==14) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=0.;uvw[2]=a;}
463     if (num==15) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=a;uvw[2]=-a;}
464     if (num==16) {w=c1*c2*c2;uvw[0]=0.;uvw[1]=a;uvw[2]=0.;}
465     if (num==17) {w=c1*c1*c2;uvw[0]=0.;uvw[1]=a;uvw[2]=a;}
466     if (num==18) {w=c1*c1*c1;uvw[0]=a;uvw[1]=-a;uvw[2]=-a;}
467     if (num==19) {w=c1*c1*c2;uvw[0]=a;uvw[1]=-a;uvw[2]=0.;}
468     if (num==20) {w=c1*c1*c1;uvw[0]=a;uvw[1]=-a;uvw[2]=a;}
469     if (num==21) {w=c1*c1*c2;uvw[0]=a;uvw[1]=0.;uvw[2]=-a;}
470     if (num==22) {w=c1*c2*c2;uvw[0]=a;uvw[1]=0.;uvw[2]=0.;}
471     if (num==23) {w=c1*c1*c2;uvw[0]=a;uvw[1]=0.;uvw[2]=a;}
472     if (num==24) {w=c1*c1*c1;uvw[0]=a;uvw[1]=a;uvw[2]=-a;}
473     if (num==25) {w=c1*c1*c2;uvw[0]=a;uvw[1]=a;uvw[2]=0.;}
474     if (num==26) {w=c1*c1*c1;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
475     return;
476     }
477     }
478    
479    
480    
481     void OT_POINTS_GAUSS::get_pt_gauss_hex(int degre, int num, double& w, double* uvw)
482     {
483     if (degre<3)
484     {
485     if (num==0) {w=2.;uvw[0]=0.;uvw[1]=0.816496580927726;uvw[2]=-0,577350269189626;}
486     if (num==1) {w=2.;uvw[0]=0.;uvw[1]=-0.816496580927726;uvw[2]=-0,577350269189626;}
487     if (num==2) {w=2.;uvw[0]=-0.816496580927726;uvw[1]=0.;uvw[2]=0,577350269189626;}
488 francois 876 if (num==3) {w=2.;uvw[0]=0.816496580927726;uvw[1]=0.;uvw[2]=0,577350269189626;}
489 francois 754 return;
490     }
491     if (degre<4)
492     {
493     if (num==0) {w=1.333333333333333;uvw[0]=1.;uvw[1]=0.;uvw[2]=0.;}
494     if (num==1) {w=1.333333333333333;uvw[0]=-1.;uvw[1]=0.;uvw[2]=0.;}
495     if (num==2) {w=1.333333333333333;uvw[0]=0.;uvw[1]=1.;uvw[2]=0.;}
496     if (num==3) {w=1.333333333333333;uvw[0]=0.;uvw[1]=-1.;uvw[2]=0.;}
497     if (num==4) {w=1.333333333333333;uvw[0]=0.;uvw[1]=0.;uvw[2]=1.;}
498     if (num==5) {w=1.333333333333333;uvw[0]=0.;uvw[1]=0.;uvw[2]=-1.;}
499     return;
500     }
501     if (degre<6)
502     {
503     double a=0.795822425754221;
504     double b=0.758786910639328;
505     double w1=0.886426592797784;
506     double w2=0.335180055401662;
507     if (num==0) {w=w1;uvw[0]=a;uvw[1]=0.;uvw[2]=0.;}
508     if (num==1) {w=w1;uvw[0]=-a;uvw[1]=0.;uvw[2]=0.;}
509     if (num==2) {w=w1;uvw[0]=0.;uvw[1]=a;uvw[2]=0.;}
510     if (num==3) {w=w1;uvw[0]=0.;uvw[1]=-a;uvw[2]=0.;}
511     if (num==4) {w=w1;uvw[0]=0.;uvw[1]=0.;uvw[2]=a;}
512     if (num==5) {w=w1;uvw[0]=0.;uvw[1]=0.;uvw[2]=-a;}
513     if (num==6) {w=w2;uvw[0]=b;uvw[1]=b;uvw[2]=b;}
514     if (num==7) {w=w2;uvw[0]=b;uvw[1]=-b;uvw[2]=-b;}
515     if (num==8) {w=w2;uvw[0]=b;uvw[1]=b;uvw[2]=-b;}
516     if (num==9) {w=w2;uvw[0]=b;uvw[1]=-b;uvw[2]=b;}
517     if (num==10) {w=w2;uvw[0]=-b;uvw[1]=b;uvw[2]=b;}
518     if (num==11) {w=w2;uvw[0]=-b;uvw[1]=-b;uvw[2]=-b;}
519     if (num==12) {w=w2;uvw[0]=-b;uvw[1]=b;uvw[2]=-b;}
520     if (num==13) {w=w2;uvw[0]=-b;uvw[1]=-b;uvw[2]=b;}
521     return;
522     }
523    
524 francois 876 }
525    
526     void OT_POINTS_GAUSS::get_pt_gauss_penta(int degre, int num, double& w, double* uvw)
527     {
528     if (degre<3)
529     {
530 francois 885 double a=1./sqrt(3.);
531     if (num==0) {w=1./6.;uvw[0]=-a;uvw[1]=0.5;uvw[2]=0.5;}
532     if (num==1) {w=1./6.;uvw[0]=-a;uvw[1]=0.;uvw[2]=0.5;}
533     if (num==2) {w=1./6.;uvw[0]=-a;uvw[1]=0.5;uvw[2]=0.;}
534     if (num==3) {w=1./6.;uvw[0]=a;uvw[1]=0.5;uvw[2]=0.5;}
535     if (num==4) {w=1./6.;uvw[0]=a;uvw[1]=0.;uvw[2]=0.5;}
536     if (num==5) {w=1./6.;uvw[0]=a;uvw[1]=0.5;uvw[2]=0.;}
537 francois 876 return;
538     }
539     if (degre<4)
540     {
541     double a=1./sqrt(3);
542     if (num==0) {w=-27./96.;uvw[0]=-a;uvw[1]=1./3.;uvw[2]=1./3.;}
543     if (num==1) {w=25./96.;uvw[0]=-a;uvw[1]=0.6;uvw[2]=0.2;}
544     if (num==2) {w=25./96.;uvw[0]=-a;uvw[1]=0.2;uvw[2]=0.6;}
545     if (num==3) {w=25./96.;uvw[0]=-a;uvw[1]=0.2;uvw[2]=0.2;}
546     if (num==4) {w=-27./96.;uvw[0]=a;uvw[1]=1./3.;uvw[2]=1./3.;}
547     if (num==5) {w=25./96.;uvw[0]=a;uvw[1]=0.6;uvw[2]=0.2;}
548     if (num==6) {w=25./96.;uvw[0]=a;uvw[1]=0.2;uvw[2]=0.6;}
549     if (num==7) {w=25./96.;uvw[0]=a;uvw[1]=0.2;uvw[2]=0.2;}
550     return;
551     }
552     if (degre<6)
553     {
554     double alpha=sqrt(3./5.);
555     double c1=5./9.;
556     double c2=8./9.;
557     double a=(6+sqrt(15))/21;
558     double b=(6-sqrt(15))/21;
559     if (num==0) {w=c1*9./80.;uvw[0]=-alpha;uvw[1]=1./3.;uvw[2]=1./3.;}
560     if (num==1) {w=c1*(155.+sqrt(15))/2400.;uvw[0]=-alpha;uvw[1]=a;uvw[2]=a;}
561     if (num==2) {w=c1*(155.+sqrt(15))/2400.;uvw[0]=-alpha;uvw[1]=1.-2.*a;uvw[2]=a;}
562     if (num==3) {w=c1*(155.+sqrt(15))/2400.;uvw[0]=-alpha;uvw[1]=a;uvw[2]=1.-2.*a;}
563     if (num==4) {w=c1*(155.-sqrt(15))/2400.;uvw[0]=-alpha;uvw[1]=b;uvw[2]=b;}
564     if (num==5) {w=c1*(155.-sqrt(15))/2400.;uvw[0]=-alpha;uvw[1]=1.-2.*b;uvw[2]=b;}
565     if (num==6) {w=c1*(155.-sqrt(15))/2400.;uvw[0]=-alpha;uvw[1]=b;uvw[2]=1.-2.*b;}
566     if (num==7) {w=c2*9./80.;uvw[0]=0.;uvw[1]=1./3.;uvw[2]=1./3.;}
567     if (num==8) {w=c2*(155.+sqrt(15))/2400.;uvw[0]=0.;uvw[1]=a;uvw[2]=a;}
568     if (num==9) {w=c2*(155.+sqrt(15))/2400.;uvw[0]=0.;uvw[1]=1.-2.*a;uvw[2]=a;}
569     if (num==10) {w=c2*(155.+sqrt(15))/2400.;uvw[0]=0.;uvw[1]=a;uvw[2]=1.-2.*a;}
570     if (num==11) {w=c2*(155.-sqrt(15))/2400.;uvw[0]=0.;uvw[1]=b;uvw[2]=b;}
571     if (num==12) {w=c2*(155.-sqrt(15))/2400.;uvw[0]=0.;uvw[1]=1.-2.*b;uvw[2]=b;}
572     if (num==13) {w=c2*(155.-sqrt(15))/2400.;uvw[0]=0.;uvw[1]=b;uvw[2]=1.-2.*b;}
573     if (num==14) {w=c1*9./80.;uvw[0]=alpha;uvw[1]=1./3.;uvw[2]=1./3.;}
574     if (num==15) {w=c1*(155.+sqrt(15))/2400.;uvw[0]=alpha;uvw[1]=a;uvw[2]=a;}
575     if (num==16) {w=c1*(155.+sqrt(15))/2400.;uvw[0]=alpha;uvw[1]=1.-2.*a;uvw[2]=a;}
576     if (num==17) {w=c1*(155.+sqrt(15))/2400.;uvw[0]=alpha;uvw[1]=a;uvw[2]=1.-2.*a;}
577     if (num==18) {w=c1*(155.-sqrt(15))/2400.;uvw[0]=alpha;uvw[1]=b;uvw[2]=b;}
578     if (num==19) {w=c1*(155.-sqrt(15))/2400.;uvw[0]=alpha;uvw[1]=1.-2.*b;uvw[2]=b;}
579     if (num==20) {w=c1*(155.-sqrt(15))/2400.;uvw[0]=alpha;uvw[1]=b;uvw[2]=1.-2.*b;}
580    
581    
582     return;
583     }
584    
585 francois 754 }