LORENE
ope_poisson_2d_solh.C
1/*
2 * Copyright (c) 2004 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 version 2
8 * as published by the Free Software Foundation.
9 *
10 * LORENE is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with LORENE; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 */
20
21char ope_poisson_2d_solh_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_solh.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $" ;
22
23/*
24 * $Id: ope_poisson_2d_solh.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_solh.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $
26 *
27 */
28#include <cmath>
29#include <cstdlib>
30
31#include "proto.h"
32#include "ope_elementary.h"
33
34
35 //------------------------------------
36 // Routine pour les cas non prevus --
37 //------------------------------------
38namespace Lorene {
39Tbl _solh_poisson_2d_pas_prevu (int, int, double, double, Tbl&) {
40
41 cout << " Solution homogene pas prevue ..... : "<< endl ;
42 exit(-1) ;
43 Tbl res(1) ;
44 return res;
45}
46
47
48 //-------------------
49 //-- R_CHEB ------
50 //-------------------
51
52Tbl _solh_poisson_2d_r_cheb (int n, int l, double alpha, double beta,
53 Tbl& val_lim) {
54
55
56 double echelle = beta / alpha ;
57
58 // CASE l != 0
59 if (l > 0) {
60 val_lim.set(0,0) = pow(echelle-1, double(l)) ;
61 val_lim.set(0,1) = double(l) * pow(echelle-1, double(l-1))/alpha ;
62 val_lim.set(0,2) = pow(echelle+1, double(l)) ;
63 val_lim.set(0,3) = double(l) * pow(echelle+1, double(l-1))/alpha ;
64
65
66 val_lim.set(1,0) = pow(echelle-1, -double(l)) ;
67 val_lim.set(1,1) = -double(l) * pow(echelle-1, -double(l+1))/alpha ;
68 val_lim.set(1,2) = pow(echelle+1, -double(l)) ;
69 val_lim.set(1,3) = -double(l) * pow(echelle+1, -double(l+1))/alpha ;
70 }
71 // CASE l =0
72 else {
73 val_lim.set(0,0) = 1. ;
74 val_lim.set(0,1) = 0. ;
75 val_lim.set(0,2) = 1. ;
76 val_lim.set(0,3) = 0. ;
77
78 val_lim.set(1,0) = log(echelle-1) ;
79 val_lim.set(1,1) = 1./(echelle-1)/alpha ;
80 val_lim.set(1,2) = log(echelle+1) ;
81 val_lim.set(1,3) = 1./(echelle+1)/alpha ;
82 }
83
84 const int nmax = 200 ; // Nombre de Tbl stockes
85 static Tbl* tab[nmax] ; // les Tbl calcules
86 static int nb_dejafait = 0 ; // nbre de Tbl calcules
87 static int l_dejafait[nmax] ;
88 static int nr_dejafait[nmax] ;
89 static double vieux_echelle = 0;
90
91 // Si on a change l'echelle : on detruit tout :
92 if (vieux_echelle != echelle) {
93 for (int i=0 ; i<nb_dejafait ; i++) {
94 l_dejafait[i] = -1 ;
95 nr_dejafait[i] = -1 ;
96 delete tab[i] ;
97 }
98 nb_dejafait = 0 ;
99 vieux_echelle = echelle ;
100 }
101
102 int indice = -1 ;
103
104 // On determine si la matrice a deja ete calculee :
105 for (int conte=0 ; conte<nb_dejafait ; conte ++)
106 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
107 indice = conte ;
108
109 // Calcul a faire :
110 if (indice == -1) {
111 if (nb_dejafait >= nmax) {
112 cout << "_solh_poisson_2d_r_cheb : trop de Tbl" << endl ;
113 abort() ;
114 exit (-1) ;
115 }
116
117
118 l_dejafait[nb_dejafait] = l ;
119 nr_dejafait[nb_dejafait] = n ;
120
121 Tbl res(2, n) ;
122 res.set_etat_qcq() ;
123 double* coloc = new double[n] ;
124
125 int * deg = new int[3] ;
126 deg[0] = 1 ;
127 deg[1] = 1 ;
128 deg[2] = n ;
129
130 //Construction de la premiere solution homogene :
131 // cad celle polynomiale.
132
133 for (int i=0 ; i<n ; i++)
134 if (l!=0)
135 coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), double(l)) ;
136 else
137 coloc[i] = 1 ;
138
139 cfrcheb(deg, deg, coloc, deg, coloc) ;
140 for (int i=0 ; i<n ;i++)
141 res.set(0, i) = coloc[i] ;
142
143 // construction de la seconde solution homogene :
144 // cad celle fractionnelle (ou log dans le cas l==0)
145 for (int i=0 ; i<n ; i++)
146 if (l != 0)
147 coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), -double(l)) ;
148 else
149 coloc[i] = log(echelle-cos(M_PI*i/(n-1))) ;
150
151 cfrcheb(deg, deg, coloc, deg, coloc) ;
152 for (int i=0 ; i<n ;i++)
153 res.set(1, i) = coloc[i] ;
154
155
156 delete [] coloc ;
157 delete [] deg ;
158 tab[nb_dejafait] = new Tbl(res) ;
159 nb_dejafait ++ ;
160 return res ;
161 }
162
163 else return *tab[indice] ;
164}
165
166 //-------------------
167 //-- R_CHEBP ------
168 //-------------------
169
170Tbl _solh_poisson_2d_r_chebp (int n, int l, double alpha,
171 double, Tbl& val_lim) {
172
173 val_lim.set(0,0) = (l!=0) ? 1 : 0 ;
174 val_lim.set(0,1) = (l!=1) ? 0 : 1 ;
175 val_lim.set(0,2) = 1. ;
176 val_lim.set(0,3) = double(l)/alpha ;
177
178 const int nmax = 200 ; // Nombre de Tbl stockes
179 static Tbl* tab[nmax] ; // les Tbl calcules
180 static int nb_dejafait = 0 ; // nbre de Tbl calcules
181 static int l_dejafait[nmax] ;
182 static int nr_dejafait[nmax] ;
183
184 int indice = -1 ;
185
186 // On determine si la matrice a deja ete calculee :
187 for (int conte=0 ; conte<nb_dejafait ; conte ++)
188 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
189 indice = conte ;
190
191 // Calcul a faire :
192 if (indice == -1) {
193 if (nb_dejafait >= nmax) {
194 cout << "_solh_poisson_2d_r_chebp : trop de Tbl" << endl ;
195 abort() ;
196 exit (-1) ;
197 }
198
199
200 l_dejafait[nb_dejafait] = l ;
201 nr_dejafait[nb_dejafait] = n ;
202
203 assert (div(l, 2).rem ==0) ;
204
205 Tbl res(n) ;
206 res.set_etat_qcq() ;
207 double* coloc = new double[n] ;
208
209 int * deg = new int[3] ;
210 deg[0] = 1 ;
211 deg[1] = 1 ;
212 deg[2] = n ;
213
214 for (int i=0 ; i<n ; i++)
215 if (l!=0)
216 coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
217 else
218 coloc[i] = 1 ;
219
220 cfrchebp(deg, deg, coloc, deg, coloc) ;
221 for (int i=0 ; i<n ;i++)
222 res.set(i) = coloc[i] ;
223
224 delete [] coloc ;
225 delete [] deg ;
226 tab[nb_dejafait] = new Tbl(res) ;
227 nb_dejafait ++ ;
228 return res ;
229 }
230
231 else return *tab[indice] ;
232}
233
234
235 //-------------------
236 //-- R_CHEBI -----
237 //-------------------
238
239Tbl _solh_poisson_2d_r_chebi (int n, int l,
240 double alpha, double, Tbl& val_lim) {
241
242 val_lim.set(0,0) = 0 ;
243 val_lim.set(0,1) = (l!=1) ? 0 : 1 ;
244 val_lim.set(0,2) = 1. ;
245 val_lim.set(0,3) = double(l)/alpha ;
246
247 const int nmax = 200 ; // Nombre de Tbl stockes
248 static Tbl* tab[nmax] ; // les Tbl calcules
249 static int nb_dejafait = 0 ; // nbre de Tbl calcules
250 static int l_dejafait[nmax] ;
251 static int nr_dejafait[nmax] ;
252
253 int indice = -1 ;
254
255 // On determine si la matrice a deja ete calculee :
256 for (int conte=0 ; conte<nb_dejafait ; conte ++)
257 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
258 indice = conte ;
259
260 // Calcul a faire :
261 if (indice == -1) {
262 if (nb_dejafait >= nmax) {
263 cout << "_solh_r_chebi : trop de Tbl" << endl ;
264 abort() ;
265 exit (-1) ;
266 }
267
268
269 l_dejafait[nb_dejafait] = l ;
270 nr_dejafait[nb_dejafait] = n ;
271
272
273 assert (div(l, 2).rem == 1) ;
274
275 Tbl res(n) ;
276 res.set_etat_qcq() ;
277 double* coloc = new double[n] ;
278
279 int * deg = new int[3] ;
280 deg[0] = 1 ;
281 deg[1] = 1 ;
282 deg[2] = n ;
283
284 for (int i=0 ; i<n ; i++)
285 coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
286
287 cfrchebi(deg, deg, coloc, deg, coloc) ;
288 for (int i=0 ; i<n ;i++)
289 res.set(i) = coloc[i] ;
290
291 delete [] coloc ;
292 delete [] deg ;
293 tab[nb_dejafait] = new Tbl(res) ;
294 nb_dejafait ++ ;
295 return res ;
296 }
297
298 else return *tab[indice] ;
299}
300
301
302
303 //-------------------
304 //-- R_CHEBU -----
305 //-------------------
306
307Tbl _solh_poisson_2d_r_chebu (int n, int l, double alpha,
308 double, Tbl& val_lim) {
309
310
311 if (l==0) {
312 cout << "Case l=0 in 2D Poisson not defined in the external compactified domain..." << endl ;
313 abort() ;
314 }
315
316 val_lim.set(0,0) = pow(-2., double(l)) ;
317 val_lim.set(0,1) = -double(l)*alpha*pow(-2, double(l+1.)) ;
318 val_lim.set(0,2) = 0 ;
319 val_lim.set(0,3) = 0 ;
320
321 const int nmax = 200 ; // Nombre de Tbl stockes
322 static Tbl* tab[nmax] ; // les Tbl calcules
323 static int nb_dejafait = 0 ; // nbre de Tbl calcules
324 static int l_dejafait[nmax] ;
325 static int nr_dejafait[nmax] ;
326
327 int indice = -1 ;
328
329 // On determine si la matrice a deja ete calculee :
330 for (int conte=0 ; conte<nb_dejafait ; conte ++)
331 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
332 indice = conte ;
333
334 // Calcul a faire :
335 if (indice == -1) {
336 if (nb_dejafait >= nmax) {
337 cout << "_solh_poisson_2d_r_chebu : trop de Tbl" << endl ;
338 abort() ;
339 exit (-1) ;
340 }
341
342
343 l_dejafait[nb_dejafait] = l ;
344 nr_dejafait[nb_dejafait] = n ;
345
346
347 // assert (l < n-1) ;
348 Tbl res(n) ;
349 res.set_etat_qcq() ;
350 double* coloc = new double[n] ;
351
352 int * deg = new int[3] ;
353 deg[0] = 1 ;
354 deg[1] = 1 ;
355 deg[2] = n ;
356
357 for (int i=0 ; i<n ; i++)
358 coloc[i] = pow(-1-cos(M_PI*i/(n-1)), double(l)) ;
359
360 cfrcheb(deg, deg, coloc, deg, coloc) ;
361 for (int i=0 ; i<n ;i++)
362 res.set(i) = coloc[i] ;
363
364 delete [] coloc ;
365 delete [] deg ;
366 tab[nb_dejafait] = new Tbl(res) ;
367 nb_dejafait ++ ;
368 return res ;
369 }
370
371 else return *tab[indice] ;
372}
373
374
376
377 // Routines de derivation
378 static Tbl (*solh_poisson_2d[MAX_BASE]) (int, int, double, double, Tbl&) ;
379 static int nap = 0 ;
380
381 // Premier appel
382 if (nap==0) {
383 nap = 1 ;
384 for (int i=0 ; i<MAX_BASE ; i++) {
385 solh_poisson_2d[i] = _solh_poisson_2d_pas_prevu ;
386 }
387 // Les routines existantes
388 solh_poisson_2d[R_CHEB >> TRA_R] = _solh_poisson_2d_r_cheb ;
389 solh_poisson_2d[R_CHEBP >> TRA_R] = _solh_poisson_2d_r_chebp ;
390 solh_poisson_2d[R_CHEBI >> TRA_R] = _solh_poisson_2d_r_chebi ;
391 solh_poisson_2d[R_CHEBU >> TRA_R] = _solh_poisson_2d_r_chebu ;
392 }
393
394 Tbl val_lim (2,4) ;
395 val_lim.set_etat_qcq() ;
397
398 s_one_minus = val_lim(0,0) ;
399 ds_one_minus = val_lim(0,1) ;
400 s_one_plus = val_lim(0,2) ;
401 ds_one_plus = val_lim(0,3) ;
402
403 s_two_minus = val_lim(1,0) ;
404 ds_two_minus = val_lim(1,1) ;
405 s_two_plus = val_lim(1,2) ;
406 ds_two_plus = val_lim(1,3) ;
407
408 return res ;
409}
410}
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
double ds_two_minus
Value of the derivative of the second homogeneous solution at the inner boundary.
double s_two_plus
Value of the second homogeneous solution at the outer boundary.
double s_one_minus
Value of the first homogeneous solution at the inner boundary.
double beta
Parameter of the associated mapping.
double ds_one_plus
Value of the derivative of the first homogeneous solution at the outer boundary.
double ds_one_minus
Value of the derivative of the first homogeneous solution at the inner boundary.
double alpha
Parameter of the associated mapping.
double s_two_minus
Value of the second homogeneous solution at the inner boundary.
int base_r
Radial basis of decomposition.
double s_one_plus
Value of the first homogeneous solution at the outer boundary.
double ds_two_plus
Value of the derivative of the second homogeneous solution at the outer boundary.
int nr
Number of radial points.
int l_quant
quantum number
virtual Tbl get_solh() const
Computes the homogeneous solutions(s).
Basic array class.
Definition tbl.h:161
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
#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