LORENE
pointsgausslobatto.C
1#include <math.h>
2#include "tbl.h"
3#include <cassert>
4
5namespace Lorene {
6double* jacobi(int , double) ;
7
8double* pointsgausslobatto(int n) {
9
10 int nmax = 10 ;
11 int ndiv = (n+1)*(n+1)+5 ;
12 double pas = double(2)/double(ndiv-1) ;
13 double eps = 2.5E-12 ;
14
15 Tbl subdiv(ndiv) ;
16 subdiv.set_etat_qcq();
17 Tbl cs(n-1,2) ;
18 cs.set_etat_qcq() ;
19 double* xx = new double[nmax] ;
20 double* yy ;
21 double* zz ;
22 int i,k,l ;
23
24 for (i = 0 ; i < ndiv ; i++) {
25 subdiv.set(i) = double(-1) + pas * double(i) ;
26 }
27
28 double a = (2*double(n)+3)/double((n+1)*(n+1)) ;
29 double b = - 1 - a ;
30
31 int j=0;
32
33 for (i = 1 ; i < ndiv-3 ; i++) {
34 yy = jacobi(n+1 , subdiv(i)) ;
35 zz = jacobi(n+1 , subdiv(i+1)) ;
36 double omega1 = yy[n+1] + a * yy[n] + b * yy[n-1] ;
37 double omega2 = zz[n+1] + a * zz[n] + b * zz[n-1] ;
38 if (omega1*omega2 <= 0) {
39 cs.set(j,0) = subdiv(i) ;
40 cs.set(j,1) = subdiv(i+1) ;
41 j++;
42 }
43 delete [] yy ;
44 delete [] zz ;
45}
46
47
48
49 double* pointsgl = new double[j+2];
50 assert(j==n-1) ;
51 pointsgl[0] = -1;
52
53 for (l = 0 ; l < j ; l++) {
54 xx[0] = cs(l,0) ;
55 xx[1] = cs(l,1) ;
56 for (k = 2 ; k < nmax ; k++) {
57 yy = jacobi(n+1 , xx[k-2]) ;
58 zz = jacobi(n+1 , xx[k-1]) ;
59 double omega1 = yy[n+1] + a * yy[n] + b * yy[n-1] ;
60 double omega2 = zz[n+1] + a * zz[n] + b * zz[n-1] ;
61 if (( fabs(xx[k-2]-xx[k-1])>=eps )&&(fabs(omega2)>=eps )) {
62 xx[k]=(xx[k-2]*omega2 -xx[k-1]*omega1)/(omega2-omega1) ;
63 }
64 else {
65 xx[k]=xx[k-1];
66 }
67 delete [] yy ;
68 delete [] zz ;
69 }
70 pointsgl[l+1] = xx[nmax-1] ;
71 }
72 delete [] xx ;
73
74 pointsgl[n] = 1 ;
75
76 return pointsgl ;
77}
78}
Lorene prototypes.
Definition app_hor.h:64