23char lap_cpt_mat_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/lap_cpt_mat.C,v 1.5 2014/10/13 08:53:29 j_novak Exp $" ;
60#include "type_parite.h"
77Matrice _lap_cpt_mat_pas_prevu(
int n,
int l) {
78 cout <<
"laplacien * (1-r^2/R_0^2) pas prevu..." << endl ;
79 cout <<
"n : " << n << endl ;
80 cout <<
"l : " << l << endl ;
93Matrice _lap_cpt_mat_r_chebp (
int n,
int l) {
95 const int nmax = 200 ;
96 static Matrice* tab[nmax] ;
97 static int nb_dejafait = 0 ;
98 static int l_dejafait[nmax] ;
99 static int nr_dejafait[nmax] ;
104 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
105 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
110 if (nb_dejafait >= nmax) {
111 cout <<
"_laplacien_nul_mat_r_chebp : trop de matrices" << endl ;
117 l_dejafait[nb_dejafait] = l ;
118 nr_dejafait[nb_dejafait] = n ;
120 Matrice res(n-1, n-1) ;
124 double* xdsdx =
new double[n] ;
125 double* x2d2sdx2 =
new double[n] ;
126 double* d2sdx2 =
new double[n] ;
127 double* sxdsdx =
new double[n] ;
128 double* sx2 =
new double [n] ;
130 for (
int i=0 ; i< n-1 ; i++) {
131 for (
int j=0 ; j<n ; j++)
135 xdsdx_1d (n, &xdsdx,
R_CHEBP) ;
138 for (
int j=0 ; j<n ; j++)
143 d2sdx2_1d(n, &x2d2sdx2,
R_CHEBP) ;
144 for (
int j=0 ; j<n ; j++)
145 d2sdx2[j] = x2d2sdx2[j] ;
146 multx2_1d(n, &x2d2sdx2,
R_CHEBP) ;
149 for (
int j=0 ; j<n ; j++)
153 sxdsdx_1d(n, &sxdsdx,
R_CHEBP) ;
156 for (
int j=0 ; j<n ; j++)
162 for (
int j=0 ; j<n-1 ; j++)
163 res.set(j, i) = (d2sdx2[j] + 2*sxdsdx[j] - l*(l+1)*sx2[j])
164 - (x2d2sdx2[j]+2*xdsdx[j]) ;
165 res.set(i, i) += l*(l+1) ;
167 res.set(i+1, i) += l*(l+1) ;
175 tab[nb_dejafait] =
new Matrice(res) ;
180 return *tab[indice] ;
190Matrice _lap_cpt_mat_r_chebi (
int n,
int l) {
192 const int nmax = 200 ;
193 static Matrice* tab[nmax] ;
194 static int nb_dejafait = 0 ;
195 static int l_dejafait[nmax] ;
196 static int nr_dejafait[nmax] ;
201 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
202 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
207 if (nb_dejafait >= nmax) {
208 cout <<
"_laplacien_nul_mat_r_chebp : trop de matrices" << endl ;
214 l_dejafait[nb_dejafait] = l ;
215 nr_dejafait[nb_dejafait] = n ;
218 int taille = (l == 1) ? n : n-1 ;
219 Matrice res(taille, taille) ;
223 double* xdsdx =
new double[n] ;
224 double* x2d2sdx2 =
new double[n] ;
225 double* d2sdx2 =
new double[n] ;
226 double* sxdsdx =
new double[n] ;
227 double* sx2 =
new double [n] ;
231 for (
int i=0 ; i<taille ; i++) {
243 for (
int j=0 ; j<n ; j++)
247 xdsdx[i+1] = f_deux ;
248 xdsdx_1d (n, &xdsdx,
R_CHEBI) ;
251 for (
int j=0 ; j<n ; j++)
255 x2d2sdx2[i+1] = f_deux ;
257 d2sdx2_1d(n, &x2d2sdx2,
R_CHEBI) ;
258 for (
int j=0 ; j<n ; j++)
259 d2sdx2[j] = x2d2sdx2[j] ;
260 multx2_1d(n, &x2d2sdx2,
R_CHEBI) ;
263 for (
int j=0 ; j<n ; j++)
267 sxdsdx[i+1] = f_deux ;
268 sxdsdx_1d(n, &sxdsdx,
R_CHEBI) ;
271 for (
int j=0 ; j<n ; j++)
278 for (
int j=0 ; j<taille ; j++)
279 res.set(j, i) = (d2sdx2[j] + 2*sxdsdx[j] - l*(l+1)*sx2[j])
280 - (x2d2sdx2[j]+2*xdsdx[j]) ;
281 res.set(i, i) += l*(l+1)*f_un ;
283 res.set(i+1, i) += l*(l+1)*f_deux ;
291 tab[nb_dejafait] =
new Matrice(res) ;
296 return *tab[indice] ;
303Matrice lap_cpt_mat(
int n,
int l,
int base_r)
307 static Matrice (*lap_cpt_mat[
MAX_BASE])(int, int) ;
314 lap_cpt_mat[i] = _lap_cpt_mat_pas_prevu ;
321 Matrice res(lap_cpt_mat[base_r](n, l)) ;
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEBP
base de Cheb. paire (rare) seulement