23char solh_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solh.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $" ;
113#include "type_parite.h"
139Tbl _solh_pas_prevu (
int n,
int l,
double echelle) {
141 cout <<
" Solution homogene pas prevue ..... : "<< endl ;
142 cout <<
" N : " << n << endl ;
143 cout <<
" l : " << l << endl ;
144 cout <<
" echelle : " << echelle << endl ;
156Tbl _solh_r_cheb (
int n,
int l,
double echelle) {
158 const int nmax = 200 ;
159 static Tbl* tab[nmax] ;
160 static int nb_dejafait = 0 ;
161 static int l_dejafait[nmax] ;
162 static int nr_dejafait[nmax] ;
163 static double vieux_echelle = 0;
166 if (vieux_echelle != echelle) {
167 for (
int i=0 ; i<nb_dejafait ; i++) {
169 nr_dejafait[i] = -1 ;
173 vieux_echelle = echelle ;
179 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
180 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
185 if (nb_dejafait >= nmax) {
186 cout <<
"_solh_r_cheb : trop de Tbl" << endl ;
192 l_dejafait[nb_dejafait] = l ;
193 nr_dejafait[nb_dejafait] = n ;
197 tab[nb_dejafait] =
new Tbl(2, n) ;
198 Tbl* pres = tab[nb_dejafait] ;
199 pres->set_etat_qcq() ;
200 double* coloc =
new double[n] ;
202 int * deg =
new int[3] ;
211 pres->set(0, 0) = 1 ;
212 for (
int i=1 ; i<n ; i++)
213 pres->set(0, i) = 0 ;
216 for (
int i=0 ; i<n ; i++)
217 coloc[i] =
pow(echelle-
cos(M_PI*i/(n-1)),
double(l)) ;
219 cfrcheb(deg, deg, coloc, deg, coloc) ;
220 for (
int i=0 ; i<n ;i++)
221 pres->set(0, i) = coloc[i] ;
227 for (
int i=0 ; i<n ; i++)
228 coloc[i] = 1/
pow(echelle-
cos(M_PI*i/(n-1)),
double(l+1)) ;
230 cfrcheb(deg, deg, coloc, deg, coloc) ;
231 for (
int i=0 ; i<n ;i++)
232 pres->set(1, i) = coloc[i] ;
237 indice = nb_dejafait ;
241 return *tab[indice] ;
249Tbl _solh_r_jaco02 (
int n,
int l,
double echelle) {
251 const int nmax = 200 ;
252 static Tbl* tab[nmax] ;
253 static int nb_dejafait = 0 ;
254 static int l_dejafait[nmax] ;
255 static int nr_dejafait[nmax] ;
256 static double vieux_echelle = 0;
259 if (vieux_echelle != echelle) {
260 for (
int i=0 ; i<nb_dejafait ; i++) {
262 nr_dejafait[i] = -1 ;
266 vieux_echelle = echelle ;
272 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
273 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
278 if (nb_dejafait >= nmax) {
279 cout <<
"_solh_r_jaco02 : trop de Tbl" << endl ;
285 l_dejafait[nb_dejafait] = l ;
286 nr_dejafait[nb_dejafait] = n ;
290 tab[nb_dejafait] =
new Tbl(2, n) ;
291 Tbl* pres = tab[nb_dejafait] ;
292 pres->set_etat_qcq() ;
293 double* coloc =
new double[n] ;
295 int * deg =
new int[3] ;
301 double* zeta = pointsgausslobatto(n-1) ;
306 pres->set(0, 0) = 1 ;
307 for (
int i=1 ; i<n ; i++)
308 pres->set(0, i) = 0 ;
311 for (
int i=0 ; i<n ; i++)
312 coloc[i] =
pow((echelle + zeta[i]),
double(l)) ;
314 cfrjaco02(deg, deg, coloc, deg, coloc) ;
315 for (
int i=0 ; i<n ;i++)
316 pres->set(0, i) = coloc[i] ;
322 for (
int i=0 ; i<n ; i++)
323 coloc[i] = 1/
pow((echelle + zeta[i]),
double(l+1)) ;
325 cfrjaco02(deg, deg, coloc, deg, coloc) ;
326 for (
int i=0 ; i<n ;i++)
327 pres->set(1, i) = coloc[i] ;
332 indice = nb_dejafait ;
336 return *tab[indice] ;
345Tbl _solh_r_chebp (
int n,
int l,
double) {
347 const int nmax = 200 ;
348 static Tbl* tab[nmax] ;
349 static int nb_dejafait = 0 ;
350 static int l_dejafait[nmax] ;
351 static int nr_dejafait[nmax] ;
356 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
357 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
362 if (nb_dejafait >= nmax) {
363 cout <<
"_solh_r_chebp : trop de Tbl" << endl ;
369 l_dejafait[nb_dejafait] = l ;
370 nr_dejafait[nb_dejafait] = n ;
372 assert (div(l, 2).rem ==0) ;
376 tab[nb_dejafait] =
new Tbl(n) ;
377 Tbl* pres = tab[nb_dejafait] ;
378 pres->set_etat_qcq() ;
379 double* coloc =
new double[n] ;
381 int * deg =
new int[3] ;
388 for (
int i=1 ; i<n ; i++)
392 for (
int i=0 ; i<n ; i++)
393 coloc[i] =
pow(
sin(M_PI*i/2/(n-1)),
double(l)) ;
395 cfrchebp(deg, deg, coloc, deg, coloc) ;
396 for (
int i=0 ; i<n ;i++)
397 pres->set(i) = coloc[i] ;
402 indice = nb_dejafait ;
406 return *tab[indice] ;
414Tbl _solh_r_chebi (
int n,
int l,
double) {
416 const int nmax = 200 ;
417 static Tbl* tab[nmax] ;
418 static int nb_dejafait = 0 ;
419 static int l_dejafait[nmax] ;
420 static int nr_dejafait[nmax] ;
425 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
426 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
431 if (nb_dejafait >= nmax) {
432 cout <<
"_solh_r_chebi : trop de Tbl" << endl ;
438 l_dejafait[nb_dejafait] = l ;
439 nr_dejafait[nb_dejafait] = n ;
442 assert (div(l, 2).rem == 1) ;
445 tab[nb_dejafait] =
new Tbl(n) ;
446 Tbl* pres = tab[nb_dejafait] ;
447 pres->set_etat_qcq() ;
448 double* coloc =
new double[n] ;
450 int * deg =
new int[3] ;
457 for (
int i=1 ; i<n ; i++)
461 for (
int i=0 ; i<n ; i++)
462 coloc[i] =
pow(
sin(M_PI*i/2/(n-1)),
double(l)) ;
464 cfrchebi(deg, deg, coloc, deg, coloc) ;
465 for (
int i=0 ; i<n ;i++)
466 pres->set(i) = coloc[i] ;
471 indice = nb_dejafait ;
475 return *tab[indice] ;
484Tbl _solh_r_chebu (
int n,
int l,
double) {
486 const int nmax = 200 ;
487 static Tbl* tab[nmax] ;
488 static int nb_dejafait = 0 ;
489 static int l_dejafait[nmax] ;
490 static int nr_dejafait[nmax] ;
495 for (
int conte=0 ; conte<nb_dejafait ; conte ++)
496 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
501 if (nb_dejafait >= nmax) {
502 cout <<
"_solh_r_chebu : trop de Tbl" << endl ;
508 l_dejafait[nb_dejafait] = l ;
509 nr_dejafait[nb_dejafait] = n ;
513 tab[nb_dejafait] =
new Tbl(n) ;
514 Tbl* pres = tab[nb_dejafait] ;
515 pres->set_etat_qcq() ;
516 double* coloc =
new double[n] ;
518 int * deg =
new int[3] ;
523 for (
int i=0 ; i<n ; i++)
524 coloc[i] =
pow(-1-
cos(M_PI*i/(n-1)),
double(l+1)) ;
526 cfrcheb(deg, deg, coloc, deg, coloc) ;
527 for (
int i=0 ; i<n ;i++)
528 pres->set(i) = coloc[i] ;
532 indice = nb_dejafait ;
536 return *tab[indice] ;
547Tbl solh(
int n,
int l,
double echelle,
int base_r) {
550 static Tbl (*solh[
MAX_BASE])(int, int, double) ;
557 solh[i] = _solh_pas_prevu ;
567 return solh[base_r](n, l, echelle) ;
Cmp sin(const Cmp &)
Sine.
Cmp pow(const Cmp &, int)
Power .
Cmp cos(const Cmp &)
Cosine.
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement