LORENE
sh_ptens_rr.C
1 /*
2  * Methods for computing the homogeneous solutions for Ope_pois_tens_rr.
3  *
4  * (see file Ope_elementary for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char sh_ptens_rr_C[] = "$Heade$" ;
29 
30 /*
31  * $Id: sh_ptens_rr.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
32  * $Log: sh_ptens_rr.C,v $
33  * Revision 1.3 2014/10/13 08:53:34 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.2 2014/10/06 15:16:13 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.1 2004/12/23 16:30:15 j_novak
40  * New files and class for the solution of the rr component of the tensor Poisson
41  * equation.
42  *
43  *
44  *
45  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/sh_ptens_rr.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
46  *
47  */
48 
49 //fichiers includes
50 #include <cstdlib>
51 #include <cmath>
52 
53 #include "matrice.h"
54 #include "type_parite.h"
55 #include "proto.h"
56 
57 /*
58  *
59  * Renvoie une ou 2 solution homogene
60  * Si base_r = R_CHEB deux solutions (x+echelle)^(l-2) dans (0, *) et
61  * 1/(x+echelle)^(l+3) dans (1, *)
62  * Si base_r = R_CHEBU 1 solution (x-1)^(l+3) dans (*)
63  * Si base_r = R_CHEBP ou R_CHEBI x^(l-2) dans (*)
64  * l est necessairement > 2...
65  *
66  * Entree :
67  * n : nbre de points en r
68  * l : nbre quantique associe
69  * echelle : cf ci-dessus, utile que dans le cas R_CHEB
70  * base_r : base de decomposition
71  *
72  * Sortie :
73  * Tbl contenant les coefficient de la ou des solution homogenes
74  *
75  */
76 
77 namespace Lorene {
78 Tbl _sh_ptens_rr_pas_prevu (int, int, double) ;
79 Tbl _sh_ptens_rr_cheb (int, int, double) ;
80 Tbl _sh_ptens_rr_chebp (int, int, double) ;
81 Tbl _sh_ptens_rr_chebi (int, int, double) ;
82 Tbl _sh_ptens_rr_chebu (int, int, double) ;
83 
84  //------------------------------------
85  // Routine pour les cas non prevus --
86  //------------------------------------
87 Tbl _sh_ptens_rr_pas_prevu (int n, int l, double echelle) {
88 
89  cout << " Solution homogene pas prevue ..... : "<< endl ;
90  cout << " N : " << n << endl ;
91  cout << " l : " << l << endl ;
92  cout << " echelle : " << echelle << endl ;
93  abort() ;
94  exit(-1) ;
95  Tbl res(1) ;
96  return res;
97 }
98 
99 
100  //-------------------
101  //-- R_CHEB ------
102  //-------------------
103 
104 Tbl _sh_ptens_rr_cheb (int n, int l, double echelle) {
105 
106  const int nmax = 200 ; // Nombre de Tbl stockes
107  static Tbl* tab[nmax] ; // les Tbl calcules
108  static int nb_dejafait = 0 ; // nbre de Tbl calcules
109  static int l_dejafait[nmax] ;
110  static int nr_dejafait[nmax] ;
111  static double vieux_echelle = 0;
112 
113  // Si on a change l'echelle : on detruit tout :
114  if (vieux_echelle != echelle) {
115  for (int i=0 ; i<nb_dejafait ; i++) {
116  l_dejafait[i] = -1 ;
117  nr_dejafait[i] = -1 ;
118  delete tab[i] ;
119  }
120  nb_dejafait = 0 ;
121  vieux_echelle = echelle ;
122  }
123 
124  int indice = -1 ;
125 
126  // On determine si la matrice a deja ete calculee :
127  for (int conte=0 ; conte<nb_dejafait ; conte ++)
128  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
129  indice = conte ;
130 
131  // Calcul a faire :
132  if (indice == -1) {
133  if (nb_dejafait >= nmax) {
134  cout << "_sh_ptens_rr_cheb : trop de Tbl" << endl ;
135  abort() ;
136  exit (-1) ;
137  }
138 
139 
140  l_dejafait[nb_dejafait] = l ;
141  nr_dejafait[nb_dejafait] = n ;
142 
143  Tbl res(2, n) ;
144  res.set_etat_qcq() ;
145  double* coloc = new double[n] ;
146 
147  int * deg = new int[3] ;
148  deg[0] = 1 ;
149  deg[1] = 1 ;
150  deg[2] = n ;
151 
152  //Construction de la premiere solution homogene :
153  // cad celle polynomiale.
154 
155  if (l==2) {
156  res.set(0, 0) = 1 ;
157  for (int i=1 ; i<n ; i++)
158  res.set(0, i) = 0 ;
159  }
160  else {
161  for (int i=0 ; i<n ; i++)
162  coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), double(l-2)) ;
163 
164  cfrcheb(deg, deg, coloc, deg, coloc) ;
165  for (int i=0 ; i<n ;i++)
166  res.set(0, i) = coloc[i] ;
167  }
168 
169 
170  // construction de la seconde solution homogene :
171  // cad celle fractionnelle.
172  for (int i=0 ; i<n ; i++)
173  coloc[i] = 1/pow(echelle-cos(M_PI*i/(n-1)), double(l+3)) ;
174 
175  cfrcheb(deg, deg, coloc, deg, coloc) ;
176  for (int i=0 ; i<n ;i++)
177  res.set(1, i) = coloc[i] ;
178 
179  delete [] coloc ;
180  delete [] deg ;
181  tab[nb_dejafait] = new Tbl(res) ;
182  nb_dejafait ++ ;
183  return res ;
184  }
185 
186  else return *tab[indice] ;
187 }
188 
189  //-------------------
190  //-- R_CHEBP ------
191  //-------------------
192 
193 Tbl _sh_ptens_rr_chebp (int n, int l, double) {
194 
195  const int nmax = 200 ; // Nombre de Tbl stockes
196  static Tbl* tab[nmax] ; // les Tbl calcules
197  static int nb_dejafait = 0 ; // nbre de Tbl calcules
198  static int l_dejafait[nmax] ;
199  static int nr_dejafait[nmax] ;
200 
201  int indice = -1 ;
202 
203  // On determine si la matrice a deja ete calculee :
204  for (int conte=0 ; conte<nb_dejafait ; conte ++)
205  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
206  indice = conte ;
207 
208  // Calcul a faire :
209  if (indice == -1) {
210  if (nb_dejafait >= nmax) {
211  cout << "_sh_ptens_rr_chebp : trop de Tbl" << endl ;
212  abort() ;
213  exit (-1) ;
214  }
215 
216 
217  l_dejafait[nb_dejafait] = l ;
218  nr_dejafait[nb_dejafait] = n ;
219 
220  Tbl res(n) ;
221  res.set_etat_qcq() ;
222  double* coloc = new double[n] ;
223 
224  int * deg = new int[3] ;
225  deg[0] = 1 ;
226  deg[1] = 1 ;
227  deg[2] = n ;
228 
229  assert (l % 2 == 0) ;
230  if (l==0) {
231  res.set(0) = 1 ;
232  for (int i=1 ; i<n ; i++)
233  res.set(i) = 0 ;
234  }
235  else {
236  for (int i=0 ; i<n ; i++)
237  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l-2)) ;
238 
239  cfrchebp(deg, deg, coloc, deg, coloc) ;
240  for (int i=0 ; i<n ;i++)
241  res.set(i) = coloc[i] ;
242  }
243 
244  delete [] coloc ;
245  delete [] deg ;
246  tab[nb_dejafait] = new Tbl(res) ;
247  nb_dejafait ++ ;
248  return res ;
249  }
250 
251  else return *tab[indice] ;
252 }
253 
254 
255  //-------------------
256  //-- R_CHEBI -----
257  //-------------------
258 
259 Tbl _sh_ptens_rr_chebi (int n, int l, double) {
260 
261  const int nmax = 200 ; // Nombre de Tbl stockes
262  static Tbl* tab[nmax] ; // les Tbl calcules
263  static int nb_dejafait = 0 ; // nbre de Tbl calcules
264  static int l_dejafait[nmax] ;
265  static int nr_dejafait[nmax] ;
266 
267  int indice = -1 ;
268 
269  // On determine si la matrice a deja ete calculee :
270  for (int conte=0 ; conte<nb_dejafait ; conte ++)
271  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
272  indice = conte ;
273 
274  // Calcul a faire :
275  if (indice == -1) {
276  if (nb_dejafait >= nmax) {
277  cout << "_sh_ptens_rr_chebi : trop de Tbl" << endl ;
278  abort() ;
279  exit (-1) ;
280  }
281 
282 
283  l_dejafait[nb_dejafait] = l ;
284  nr_dejafait[nb_dejafait] = n ;
285 
286 
287  assert (l%2 == 1) ;
288 
289  Tbl res(n) ;
290  res.set_etat_qcq() ;
291  double* coloc = new double[n] ;
292 
293  int * deg = new int[3] ;
294  deg[0] = 1 ;
295  deg[1] = 1 ;
296  deg[2] = n ;
297 
298  for (int i=0 ; i<n ; i++)
299  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l-2)) ;
300 
301  cfrchebi(deg, deg, coloc, deg, coloc) ;
302  for (int i=0 ; i<n ;i++)
303  res.set(i) = coloc[i] ;
304 
305 
306  delete [] coloc ;
307  delete [] deg ;
308  tab[nb_dejafait] = new Tbl(res) ;
309  nb_dejafait ++ ;
310  return res ;
311  }
312 
313  else return *tab[indice] ;
314 }
315 
316 
317 
318  //-------------------
319  //-- R_CHEBU -----
320  //-------------------
321 
322 Tbl _sh_ptens_rr_chebu (int n, int l, double) {
323 
324  const int nmax = 200 ; // Nombre de Tbl stockes
325  static Tbl* tab[nmax] ; // les Tbl calcules
326  static int nb_dejafait = 0 ; // nbre de Tbl calcules
327  static int l_dejafait[nmax] ;
328  static int nr_dejafait[nmax] ;
329 
330  int indice = -1 ;
331 
332  // On determine si la matrice a deja ete calculee :
333  for (int conte=0 ; conte<nb_dejafait ; conte ++)
334  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
335  indice = conte ;
336 
337  // Calcul a faire :
338  if (indice == -1) {
339  if (nb_dejafait >= nmax) {
340  cout << "_sh_ptens_rr_chebu : trop de Tbl" << endl ;
341  abort() ;
342  exit (-1) ;
343  }
344 
345 
346  l_dejafait[nb_dejafait] = l ;
347  nr_dejafait[nb_dejafait] = n ;
348 
349  Tbl res(n) ;
350  res.set_etat_qcq() ;
351  double* coloc = new double[n] ;
352 
353  int * deg = new int[3] ;
354  deg[0] = 1 ;
355  deg[1] = 1 ;
356  deg[2] = n ;
357 
358  for (int i=0 ; i<n ; i++)
359  coloc[i] = pow(-1-cos(M_PI*i/(n-1)), double(l+3)) ;
360 
361  cfrcheb(deg, deg, coloc, deg, coloc) ;
362  for (int i=0 ; i<n ;i++)
363  res.set(i) = coloc[i] ;
364 
365  delete [] coloc ;
366  delete [] deg ;
367  tab[nb_dejafait] = new Tbl(res) ;
368  nb_dejafait ++ ;
369  return res ;
370  }
371 
372  else return *tab[indice] ;
373 }
374 
375 
376 
377 
378  //-------------------
379  //-- Fonction ----
380  //-------------------
381 
382 
383 Tbl sh_ptens_rr(int n, int l, double echelle, int base_r) {
384 
385  // Routines de derivation
386  static Tbl (*sh_ptens_rr[MAX_BASE])(int, int, double) ;
387  static int nap = 0 ;
388 
389  // Premier appel
390  if (nap==0) {
391  nap = 1 ;
392  for (int i=0 ; i<MAX_BASE ; i++) {
393  sh_ptens_rr[i] = _sh_ptens_rr_pas_prevu ;
394  }
395  // Les routines existantes
396  sh_ptens_rr[R_CHEB >> TRA_R] = _sh_ptens_rr_cheb ;
397  sh_ptens_rr[R_CHEBU >> TRA_R] = _sh_ptens_rr_chebu ;
398  sh_ptens_rr[R_CHEBP >> TRA_R] = _sh_ptens_rr_chebp ;
399  sh_ptens_rr[R_CHEBI >> TRA_R] = _sh_ptens_rr_chebi ;
400  }
401 
402  assert (l > 1) ;
403 
404  Tbl res(sh_ptens_rr[base_r](n, l, echelle)) ;
405  return res ;
406 }
407 }
Lorene prototypes.
Definition: app_hor.h:64
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