LORENE
pointsgausslobatto.C
1 #include <math.h>
2 #include "tbl.h"
3 #include <cassert>
4 
5 namespace Lorene {
6 double* jacobi(int , double) ;
7 
8 double* 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:67