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
28char 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
77namespace Lorene {
78Tbl _sh_ptens_rr_pas_prevu (int, int, double) ;
79Tbl _sh_ptens_rr_cheb (int, int, double) ;
80Tbl _sh_ptens_rr_chebp (int, int, double) ;
81Tbl _sh_ptens_rr_chebi (int, int, double) ;
82Tbl _sh_ptens_rr_chebu (int, int, double) ;
83
84 //------------------------------------
85 // Routine pour les cas non prevus --
86 //------------------------------------
87Tbl _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
104Tbl _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
193Tbl _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
259Tbl _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
322Tbl _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
383Tbl 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}
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#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
Lorene prototypes.
Definition app_hor.h:64