LORENE
solh.C
1 /*
2  * Copyright (c) 1999-2001 Philippe Grandclement
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char 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 $" ;
24 
25 /*
26  * $Id: solh.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $
27  * $Log: solh.C,v $
28  * Revision 1.10 2014/10/13 08:53:30 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.9 2014/10/06 15:16:10 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.8 2008/02/18 13:53:43 j_novak
35  * Removal of special indentation instructions.
36  *
37  * Revision 1.7 2007/12/20 09:11:09 jl_cornou
38  * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
39  *
40  * Revision 1.6 2007/12/13 15:48:46 jl_cornou
41  * *** empty log message ***
42  *
43  * Revision 1.5 2007/12/12 12:30:49 jl_cornou
44  * *** empty log message ***
45  *
46  * Revision 1.4 2004/10/05 15:44:21 j_novak
47  * Minor speed enhancements.
48  *
49  * Revision 1.3 2003/01/31 08:49:58 e_gourgoulhon
50  * Increased the number nmax of stored matrices from 100 to 200.
51  *
52  * Revision 1.2 2002/10/16 14:37:12 j_novak
53  * Reorganization of #include instructions of standard C++, in order to
54  * use experimental version 3 of gcc.
55  *
56  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
57  * LORENE
58  *
59  * Revision 2.13 2000/01/18 14:15:50 phil
60  * enleve assert sur nobre de points min en r
61  *
62  * Revision 2.12 2000/01/04 18:59:39 phil
63  * Double nmax
64  *
65  * Revision 2.11 1999/10/11 14:29:37 phil
66  * &-> &&
67  *
68  * Revision 2.10 1999/09/30 09:23:25 phil
69  * remplacement des && en &
70  *
71  * Revision 2.9 1999/09/17 15:26:25 phil
72  * correction definition de NMAX
73  *
74  * Revision 2.8 1999/06/23 12:35:00 phil
75  * *** empty log message ***
76  *
77  * Revision 2.7 1999/04/28 10:49:19 phil
78  * augmentation de nmax a 50
79  *
80  * Revision 2.6 1999/04/27 13:12:34 phil
81  * *** empty log message ***
82  *
83  * Revision 2.5 1999/04/19 14:07:02 phil
84  * *** empty log message ***
85  *
86  * Revision 2.4 1999/04/16 13:19:07 phil
87  * *** empty log message ***
88  *
89  * Revision 2.3 1999/04/16 13:17:02 phil
90  * *** empty log message ***
91  *
92  * Revision 2.2 1999/04/14 13:57:05 phil
93  * Sauvegarde des solutions deja calculees
94  *
95  * Revision 2.1 1999/04/07 14:39:23 phil
96  * esthetique
97  *
98  * Revision 2.0 1999/04/07 14:11:18 phil
99  * *** empty log message ***
100  *
101  *
102  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solh.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $
103  *
104  */
105 
106 //fichiers includes
107 #include <cstdio>
108 #include <cstdlib>
109 #include <cmath>
110 
111 #include "proto.h"
112 #include "matrice.h"
113 #include "type_parite.h"
114 
115 
116 /*
117  *
118  * Renvoie une ou 2 solutions homogenes
119  * Si base_r = R_CHEB deux solutions (x+echelle)^l dans (0, *) et
120  * 1/(x+echelle)^(l+1) dans (1, *)
121  * Si base_r = R_CHEBU 1 solution (x-1)^l+1 dans (*)
122  * Si base_r = R_CHEBP ou R_CHEBI x^l dans (*)
123  *
124  * Entree :
125  * n : nbre de points en r
126  * l : nbre quantique associe
127  * echelle : cf ci-dessus, utile que dans le cas R_CHEB
128  * base_r : base de decomposition
129  *
130  * Sortie :
131  * Tbl contenant les coefficient de la ou des solution homogenes
132  *
133  */
134 
135  //------------------------------------
136  // Routine pour les cas non prevus --
137  //------------------------------------
138 namespace Lorene {
139 Tbl _solh_pas_prevu (int n, int l, double echelle) {
140 
141  cout << " Solution homogene pas prevue ..... : "<< endl ;
142  cout << " N : " << n << endl ;
143  cout << " l : " << l << endl ;
144  cout << " echelle : " << echelle << endl ;
145  abort() ;
146  exit(-1) ;
147  Tbl res(1) ;
148  return res;
149 }
150 
151 
152  //-------------------
153  //-- R_CHEB ------
154  //-------------------
155 
156 Tbl _solh_r_cheb (int n, int l, double echelle) {
157 
158  const int nmax = 200 ; // Nombre de Tbl stockes
159  static Tbl* tab[nmax] ; // les Tbl calcules
160  static int nb_dejafait = 0 ; // nbre de Tbl calcules
161  static int l_dejafait[nmax] ;
162  static int nr_dejafait[nmax] ;
163  static double vieux_echelle = 0;
164 
165  // Si on a change l'echelle : on detruit tout :
166  if (vieux_echelle != echelle) {
167  for (int i=0 ; i<nb_dejafait ; i++) {
168  l_dejafait[i] = -1 ;
169  nr_dejafait[i] = -1 ;
170  delete tab[i] ;
171  }
172  nb_dejafait = 0 ;
173  vieux_echelle = echelle ;
174  }
175 
176  int indice = -1 ;
177 
178  // On determine si la matrice a deja ete calculee :
179  for (int conte=0 ; conte<nb_dejafait ; conte ++)
180  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
181  indice = conte ;
182 
183  // Calcul a faire :
184  if (indice == -1) {
185  if (nb_dejafait >= nmax) {
186  cout << "_solh_r_cheb : trop de Tbl" << endl ;
187  abort() ;
188  exit (-1) ;
189  }
190 
191 
192  l_dejafait[nb_dejafait] = l ;
193  nr_dejafait[nb_dejafait] = n ;
194 
195  // assert (l < n) ;
196 
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] ;
201 
202  int * deg = new int[3] ;
203  deg[0] = 1 ;
204  deg[1] = 1 ;
205  deg[2] = n ;
206 
207  //Construction de la premiere solution homogene :
208  // cad celle polynomiale.
209 
210  if (l==0) {
211  pres->set(0, 0) = 1 ;
212  for (int i=1 ; i<n ; i++)
213  pres->set(0, i) = 0 ;
214  }
215  else {
216  for (int i=0 ; i<n ; i++)
217  coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), double(l)) ;
218 
219  cfrcheb(deg, deg, coloc, deg, coloc) ;
220  for (int i=0 ; i<n ;i++)
221  pres->set(0, i) = coloc[i] ;
222  }
223 
224 
225  // construction de la seconde solution homogene :
226  // cad celle fractionnelle.
227  for (int i=0 ; i<n ; i++)
228  coloc[i] = 1/pow(echelle-cos(M_PI*i/(n-1)), double(l+1)) ;
229 
230  cfrcheb(deg, deg, coloc, deg, coloc) ;
231  for (int i=0 ; i<n ;i++)
232  pres->set(1, i) = coloc[i] ;
233 
234 
235  delete [] coloc ;
236  delete [] deg ;
237  indice = nb_dejafait ;
238  nb_dejafait ++ ;
239  }
240 
241  return *tab[indice] ;
242 }
243 
244 
245  //-------------------
246  //-- R_JACO02 ------
247  //-------------------
248 
249 Tbl _solh_r_jaco02 (int n, int l, double echelle) {
250 
251  const int nmax = 200 ; // Nombre de Tbl stockes
252  static Tbl* tab[nmax] ; // les Tbl calcules
253  static int nb_dejafait = 0 ; // nbre de Tbl calcules
254  static int l_dejafait[nmax] ;
255  static int nr_dejafait[nmax] ;
256  static double vieux_echelle = 0;
257 
258  // Si on a change l'echelle : on detruit tout :
259  if (vieux_echelle != echelle) {
260  for (int i=0 ; i<nb_dejafait ; i++) {
261  l_dejafait[i] = -1 ;
262  nr_dejafait[i] = -1 ;
263  delete tab[i] ;
264  }
265  nb_dejafait = 0 ;
266  vieux_echelle = echelle ;
267  }
268 
269  int indice = -1 ;
270 
271  // On determine si la matrice a deja ete calculee :
272  for (int conte=0 ; conte<nb_dejafait ; conte ++)
273  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
274  indice = conte ;
275 
276  // Calcul a faire :
277  if (indice == -1) {
278  if (nb_dejafait >= nmax) {
279  cout << "_solh_r_jaco02 : trop de Tbl" << endl ;
280  abort() ;
281  exit (-1) ;
282  }
283 
284 
285  l_dejafait[nb_dejafait] = l ;
286  nr_dejafait[nb_dejafait] = n ;
287 
288  // assert (l < n) ;
289 
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] ;
294 
295  int * deg = new int[3] ;
296  deg[0] = 1 ;
297  deg[1] = 1 ;
298  deg[2] = n ;
299 
300 
301  double* zeta = pointsgausslobatto(n-1) ;
302  //Construction de la premiere solution homogene :
303  // cad celle polynomiale.
304 
305  if (l==0) {
306  pres->set(0, 0) = 1 ;
307  for (int i=1 ; i<n ; i++)
308  pres->set(0, i) = 0 ;
309  }
310  else {
311  for (int i=0 ; i<n ; i++)
312  coloc[i] = pow((echelle + zeta[i]), double(l)) ;
313 
314  cfrjaco02(deg, deg, coloc, deg, coloc) ;
315  for (int i=0 ; i<n ;i++)
316  pres->set(0, i) = coloc[i] ;
317  }
318 
319 
320  // construction de la seconde solution homogene :
321  // cad celle fractionnelle.
322  for (int i=0 ; i<n ; i++)
323  coloc[i] = 1/pow((echelle + zeta[i]), double(l+1)) ;
324 
325  cfrjaco02(deg, deg, coloc, deg, coloc) ;
326  for (int i=0 ; i<n ;i++)
327  pres->set(1, i) = coloc[i] ;
328 
329 
330  delete [] coloc ;
331  delete [] deg ;
332  indice = nb_dejafait ;
333  nb_dejafait ++ ;
334  }
335 
336  return *tab[indice] ;
337 }
338 
339 
340 
341  //-------------------
342  //-- R_CHEBP ------
343  //-------------------
344 
345 Tbl _solh_r_chebp (int n, int l, double) {
346 
347  const int nmax = 200 ; // Nombre de Tbl stockes
348  static Tbl* tab[nmax] ; // les Tbl calcules
349  static int nb_dejafait = 0 ; // nbre de Tbl calcules
350  static int l_dejafait[nmax] ;
351  static int nr_dejafait[nmax] ;
352 
353  int indice = -1 ;
354 
355  // On determine si la matrice a deja ete calculee :
356  for (int conte=0 ; conte<nb_dejafait ; conte ++)
357  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
358  indice = conte ;
359 
360  // Calcul a faire :
361  if (indice == -1) {
362  if (nb_dejafait >= nmax) {
363  cout << "_solh_r_chebp : trop de Tbl" << endl ;
364  abort() ;
365  exit (-1) ;
366  }
367 
368 
369  l_dejafait[nb_dejafait] = l ;
370  nr_dejafait[nb_dejafait] = n ;
371 
372  assert (div(l, 2).rem ==0) ;
373 
374  // assert (l < 2*n-1) ;
375 
376  tab[nb_dejafait] = new Tbl(n) ;
377  Tbl* pres = tab[nb_dejafait] ;
378  pres->set_etat_qcq() ;
379  double* coloc = new double[n] ;
380 
381  int * deg = new int[3] ;
382  deg[0] = 1 ;
383  deg[1] = 1 ;
384  deg[2] = n ;
385 
386  if (l==0) {
387  pres->set(0) = 1 ;
388  for (int i=1 ; i<n ; i++)
389  pres->set(i) = 0 ;
390  }
391  else {
392  for (int i=0 ; i<n ; i++)
393  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
394 
395  cfrchebp(deg, deg, coloc, deg, coloc) ;
396  for (int i=0 ; i<n ;i++)
397  pres->set(i) = coloc[i] ;
398  }
399 
400  delete [] coloc ;
401  delete [] deg ;
402  indice = nb_dejafait ;
403  nb_dejafait ++ ;
404  }
405 
406  return *tab[indice] ;
407 }
408 
409 
410  //-------------------
411  //-- R_CHEBI -----
412  //-------------------
413 
414 Tbl _solh_r_chebi (int n, int l, double) {
415 
416  const int nmax = 200 ; // Nombre de Tbl stockes
417  static Tbl* tab[nmax] ; // les Tbl calcules
418  static int nb_dejafait = 0 ; // nbre de Tbl calcules
419  static int l_dejafait[nmax] ;
420  static int nr_dejafait[nmax] ;
421 
422  int indice = -1 ;
423 
424  // On determine si la matrice a deja ete calculee :
425  for (int conte=0 ; conte<nb_dejafait ; conte ++)
426  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
427  indice = conte ;
428 
429  // Calcul a faire :
430  if (indice == -1) {
431  if (nb_dejafait >= nmax) {
432  cout << "_solh_r_chebi : trop de Tbl" << endl ;
433  abort() ;
434  exit (-1) ;
435  }
436 
437 
438  l_dejafait[nb_dejafait] = l ;
439  nr_dejafait[nb_dejafait] = n ;
440 
441 
442  assert (div(l, 2).rem == 1) ;
443  // assert (l < 2*n) ;
444 
445  tab[nb_dejafait] = new Tbl(n) ;
446  Tbl* pres = tab[nb_dejafait] ;
447  pres->set_etat_qcq() ;
448  double* coloc = new double[n] ;
449 
450  int * deg = new int[3] ;
451  deg[0] = 1 ;
452  deg[1] = 1 ;
453  deg[2] = n ;
454 
455  if (l==1) {
456  pres->set(0) = 1 ;
457  for (int i=1 ; i<n ; i++)
458  pres->set(i) = 0 ;
459  }
460  else {
461  for (int i=0 ; i<n ; i++)
462  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
463 
464  cfrchebi(deg, deg, coloc, deg, coloc) ;
465  for (int i=0 ; i<n ;i++)
466  pres->set(i) = coloc[i] ;
467  }
468 
469  delete [] coloc ;
470  delete [] deg ;
471  indice = nb_dejafait ;
472  nb_dejafait ++ ;
473  }
474 
475  return *tab[indice] ;
476 }
477 
478 
479 
480  //-------------------
481  //-- R_CHEBU -----
482  //-------------------
483 
484 Tbl _solh_r_chebu (int n, int l, double) {
485 
486  const int nmax = 200 ; // Nombre de Tbl stockes
487  static Tbl* tab[nmax] ; // les Tbl calcules
488  static int nb_dejafait = 0 ; // nbre de Tbl calcules
489  static int l_dejafait[nmax] ;
490  static int nr_dejafait[nmax] ;
491 
492  int indice = -1 ;
493 
494  // On determine si la matrice a deja ete calculee :
495  for (int conte=0 ; conte<nb_dejafait ; conte ++)
496  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
497  indice = conte ;
498 
499  // Calcul a faire :
500  if (indice == -1) {
501  if (nb_dejafait >= nmax) {
502  cout << "_solh_r_chebu : trop de Tbl" << endl ;
503  abort() ;
504  exit (-1) ;
505  }
506 
507 
508  l_dejafait[nb_dejafait] = l ;
509  nr_dejafait[nb_dejafait] = n ;
510 
511 
512  // assert (l < n-1) ;
513  tab[nb_dejafait] = new Tbl(n) ;
514  Tbl* pres = tab[nb_dejafait] ;
515  pres->set_etat_qcq() ;
516  double* coloc = new double[n] ;
517 
518  int * deg = new int[3] ;
519  deg[0] = 1 ;
520  deg[1] = 1 ;
521  deg[2] = n ;
522 
523  for (int i=0 ; i<n ; i++)
524  coloc[i] = pow(-1-cos(M_PI*i/(n-1)), double(l+1)) ;
525 
526  cfrcheb(deg, deg, coloc, deg, coloc) ;
527  for (int i=0 ; i<n ;i++)
528  pres->set(i) = coloc[i] ;
529 
530  delete [] coloc ;
531  delete [] deg ;
532  indice = nb_dejafait ;
533  nb_dejafait ++ ;
534  }
535 
536  return *tab[indice] ;
537 }
538 
539 
540 
541 
542  //-------------------
543  //-- Fonction ----
544  //-------------------
545 
546 
547 Tbl solh(int n, int l, double echelle, int base_r) {
548 
549  // Routines de derivation
550  static Tbl (*solh[MAX_BASE])(int, int, double) ;
551  static int nap = 0 ;
552 
553  // Premier appel
554  if (nap==0) {
555  nap = 1 ;
556  for (int i=0 ; i<MAX_BASE ; i++) {
557  solh[i] = _solh_pas_prevu ;
558  }
559  // Les routines existantes
560  solh[R_CHEB >> TRA_R] = _solh_r_cheb ;
561  solh[R_CHEBU >> TRA_R] = _solh_r_chebu ;
562  solh[R_CHEBP >> TRA_R] = _solh_r_chebp ;
563  solh[R_CHEBI >> TRA_R] = _solh_r_chebi ;
564  solh[R_JACO02 >> TRA_R] = _solh_r_jaco02 ;
565  }
566 
567  return solh[base_r](n, l, echelle) ;
568 }
569 }
Lorene prototypes.
Definition: app_hor.h:64
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166