LORENE
ope_poisson_2d_cl.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_cl_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_cl.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $" ;
22
23/*
24 * $Id: ope_poisson_2d_cl.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_cl.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
26 *
27 */
28#include <cmath>
29#include <cstdlib>
30
31#include "proto.h"
32#include "ope_elementary.h"
33
34// Version Matrice --> Matrice
35namespace Lorene {
36Matrice _cl_poisson_2d_pas_prevu (const Matrice & source, int, double,int) {
37 cout << "Combinaison lineaire pas prevu..." << endl ;
38 abort() ;
39 exit(-1) ;
40 return source;
41}
42
43
44 //-------------------
45 //-- R_CHEB ------
46 //-------------------
47
48Matrice _cl_poisson_2d_r_cheb (const Matrice &source, int l, double echelle, int) {
49 int n = source.get_dim(0) ;assert (n == source.get_dim(1)) ;
50
51 const int nmax = 100 ; // Nombre de Matrices stockees
52 static Matrice* tab[nmax] ; // les matrices calculees
53 static int nb_dejafait = 0 ; // nbre de matrices calculees
54 static int l_dejafait[nmax] ;
55 static int nr_dejafait[nmax] ;
56 static double vieux_echelle = 0 ;
57
58 // Si on a change l'echelle : on detruit tout :
59 if (vieux_echelle != echelle) {
60 for (int i=0 ; i<nb_dejafait ; i++) {
61 l_dejafait[i] = -1 ;
62 nr_dejafait[i] = -1 ;
63 delete tab[i] ;
64 }
65 nb_dejafait = 0 ;
66 vieux_echelle = echelle ;
67 }
68
69 int indice = -1 ;
70
71 // On determine si la matrice a deja ete calculee :
72 for (int conte=0 ; conte<nb_dejafait ; conte ++)
73 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
74 indice = conte ;
75
76 // Calcul a faire :
77 if (indice == -1) {
78 if (nb_dejafait >= nmax) {
79 cout << "_cl_poisson_r_cheb : trop de matrices" << endl ;
80 abort() ;
81 exit (-1) ;
82 }
83
84 l_dejafait[nb_dejafait] = l ;
85 nr_dejafait[nb_dejafait] = n ;
86
87 Matrice barre(source) ;
88 int dirac = 1 ;
89 for (int i=0 ; i<n-2 ; i++) {
90 for (int j=0 ; j<n ; j++)
91 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
92 /(i+1) ;
93 if (i==0) dirac = 0 ;
94 }
95
96 Matrice res(barre) ;
97 for (int i=0 ; i<n-4 ; i++)
98 for (int j=0 ; j<n ; j++)
99 res.set(i, j) = barre(i, j)-barre(i+2, j) ;
100 tab[nb_dejafait] = new Matrice(res) ;
101 nb_dejafait ++ ;
102 return res ;
103 }
104
105 // Cas ou le calcul a deja ete effectue :
106 else
107 return *tab[indice] ;
108}
109
110 //-------------------
111 //-- R_CHEBP -----
112 //-------------------
113
114
115Matrice _cl_poisson_2d_r_chebp (const Matrice &source, int l, double, int) {
116
117 int n = source.get_dim(0) ;
118 assert (n == source.get_dim(1)) ;
119
120 const int nmax = 100 ; // Nombre de Matrices stockees
121 static Matrice* tab[nmax] ; // les matrices calculees
122 static int nb_dejafait = 0 ; // nbre de matrices calculees
123 static int l_dejafait[nmax] ;
124 static int nr_dejafait[nmax] ;
125
126 int indice = -1 ;
127
128 // On determine si la matrice a deja ete calculee :
129 for (int conte=0 ; conte<nb_dejafait ; conte ++)
130 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
131 indice = conte ;
132
133 // Calcul a faire :
134 if (indice == -1) {
135 if (nb_dejafait >= nmax) {
136 cout << "_cl_poisson_2d_r_chebp : trop de matrices" << endl ;
137 abort() ;
138 exit (-1) ;
139 }
140
141 l_dejafait[nb_dejafait] = l ;
142 nr_dejafait[nb_dejafait] = n ;
143
144 Matrice barre(source) ;
145
146 int dirac = 1 ;
147 for (int i=0 ; i<n-2 ; i++) {
148 for (int j=0 ; j<n ; j++)
149 barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
150 if (i==0) dirac = 0 ;
151 }
152
153 Matrice tilde(barre) ;
154 for (int i=0 ; i<n-4 ; i++)
155 for (int j=0 ; j<n ; j++)
156 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
157
158 Matrice res(tilde) ;
159 for (int i=0 ; i<n-4 ; i++)
160 for (int j=0 ; j<n ; j++)
161 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
162 tab[nb_dejafait] = new Matrice(res) ;
163 nb_dejafait ++ ;
164 return res ;
165 }
166
167 // Cas ou le calcul a deja ete effectue :
168 else
169 return *tab[indice] ;
170}
171
172 //-------------------
173 //-- R_CHEBI -----
174 //-------------------
175
176
177Matrice _cl_poisson_2d_r_chebi (const Matrice &source, int l, double, int) {
178 int n = source.get_dim(0) ;
179 assert (n == source.get_dim(1)) ;
180
181
182 const int nmax = 100 ; // Nombre de Matrices stockees
183 static Matrice* tab[nmax] ; // les matrices calculees
184 static int nb_dejafait = 0 ; // nbre de matrices calculees
185 static int l_dejafait[nmax] ;
186 static int nr_dejafait[nmax] ;
187
188 int indice = -1 ;
189
190 // On determine si la matrice a deja ete calculee :
191 for (int conte=0 ; conte<nb_dejafait ; conte ++)
192 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
193 indice = conte ;
194
195 // Calcul a faire :
196 if (indice == -1) {
197 if (nb_dejafait >= nmax) {
198 cout << "_cl_poisson_2d_r_chebi : trop de matrices" << endl ;
199 abort() ;
200 exit (-1) ;
201 }
202
203 l_dejafait[nb_dejafait] = l ;
204 nr_dejafait[nb_dejafait] = n ;
205
206 Matrice barre(source) ;
207
208 for (int i=0 ; i<n-2 ; i++)
209 for (int j=0 ; j<n ; j++)
210 barre.set(i, j) = source(i, j)-source(i+2, j) ;
211
212 Matrice tilde(barre) ;
213 for (int i=0 ; i<n-4 ; i++)
214 for (int j=0 ; j<n ; j++)
215 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
216
217 Matrice res(tilde) ;
218 for (int i=0 ; i<n-4 ; i++)
219 for (int j=0 ; j<n ; j++)
220 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
221 tab[nb_dejafait] = new Matrice(res) ;
222 nb_dejafait ++ ;
223 return res ;
224 }
225
226 // Cas ou le calcul a deja ete effectue :
227 else
228 return *tab[indice] ;
229}
230 //-------------------
231 //-- R_CHEBU -----
232 //-------------------
233
234Matrice _cl_poisson_2d_r_chebu_quatre (const Matrice&, int) ;
235Matrice _cl_poisson_2d_r_chebu_trois (const Matrice&, int) ;
236Matrice _cl_poisson_2d_r_chebu_deux (const Matrice&, int) ;
237
238
239Matrice _cl_poisson_2d_r_chebu (const Matrice &source, int l, double, int puis) {
240 int n = source.get_dim(0) ;
241 assert (n == source.get_dim(1)) ;
242
243 Matrice res(n, n) ;
244 res.set_etat_qcq() ;
245
246 switch (puis) {
247 case 4 :
248 res = _cl_poisson_2d_r_chebu_quatre(source, l) ;
249 break ;
250 case 3 :
251 res = _cl_poisson_2d_r_chebu_trois (source, l) ;
252 break ;
253 case 2 :
254 res = _cl_poisson_2d_r_chebu_deux(source, l) ;
255 break ;
256 default :
257 abort() ;
258 exit(-1) ;
259 }
260
261 return res ;
262}
263
264
265// Cas dzpuis = 4
266Matrice _cl_poisson_2d_r_chebu_quatre (const Matrice &source, int l) {
267 int n = source.get_dim(0) ;
268 assert (n == source.get_dim(1)) ;
269
270
271 const int nmax = 200 ; // Nombre de Matrices stockees
272 static Matrice* tab[nmax] ; // les matrices calculees
273 static int nb_dejafait = 0 ; // nbre de matrices calculees
274 static int l_dejafait[nmax] ;
275 static int nr_dejafait[nmax] ;
276
277 int indice = -1 ;
278
279 // On determine si la matrice a deja ete calculee :
280 for (int conte=0 ; conte<nb_dejafait ; conte ++)
281 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
282 indice = conte ;
283
284 // Calcul a faire :
285 if (indice == -1) {
286 if (nb_dejafait >= nmax) {
287 cout << "_cl_poisson_2d_r_chebu_quatre : trop de matrices" << endl ;
288 abort() ;
289 exit (-1) ;
290 }
291
292 l_dejafait[nb_dejafait] = l ;
293 nr_dejafait[nb_dejafait] = n ;
294
295 Matrice barre(source) ;
296
297 int dirac = 1 ;
298 for (int i=0 ; i<n-2 ; i++) {
299 for (int j=0 ; j<n ; j++)
300 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
301 if (i==0) dirac = 0 ;
302 }
303
304 Matrice tilde(barre) ;
305 for (int i=0 ; i<n-4 ; i++)
306 for (int j=0 ; j<n ; j++)
307 tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
308
309 Matrice prime(tilde) ;
310 for (int i=0 ; i<n-4 ; i++)
311 for (int j=0 ; j<n ; j++)
312 prime.set(i, j) = (tilde(i, j)-tilde(i+1, j)) ;
313
314 Matrice res(prime) ;
315 for (int i=0 ; i<n-4 ; i++)
316 for (int j=0 ; j<n ; j++)
317 res.set(i, j) = (prime(i, j)-prime(i+2, j)) ;
318 tab[nb_dejafait] = new Matrice(res) ;
319 nb_dejafait ++ ;
320 return res ;
321 }
322
323 // Cas ou le calcul a deja ete effectue :
324 else
325 return *tab[indice] ;
326}
327
328// Cas dzpuis == 3
329Matrice _cl_poisson_2d_r_chebu_trois (const Matrice &source, int l) {
330 int n = source.get_dim(0) ;
331 assert (n == source.get_dim(1)) ;
332
333
334 const int nmax = 200 ; // Nombre de Matrices stockees
335 static Matrice* tab[nmax] ; // les matrices calculees
336 static int nb_dejafait = 0 ; // nbre de matrices calculees
337 static int l_dejafait[nmax] ;
338 static int nr_dejafait[nmax] ;
339
340 int indice = -1 ;
341
342 // On determine si la matrice a deja ete calculee :
343 for (int conte=0 ; conte<nb_dejafait ; conte ++)
344 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
345 indice = conte ;
346
347 // Calcul a faire :
348 if (indice == -1) {
349 if (nb_dejafait >= nmax) {
350 cout << "_cl_poisson_2d_r_chebu_trois : trop de matrices" << endl ;
351 abort() ;
352 exit (-1) ;
353 }
354
355 l_dejafait[nb_dejafait] = l ;
356 nr_dejafait[nb_dejafait] = n ;
357
358 Matrice barre(source) ;
359
360 int dirac = 1 ;
361 for (int i=0 ; i<n-2 ; i++) {
362 for (int j=0 ; j<n ; j++)
363 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
364 if (i==0) dirac = 0 ;
365 }
366
367 Matrice tilde(barre) ;
368 for (int i=0 ; i<n-2 ; i++)
369 for (int j=0 ; j<n ; j++)
370 tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
371
372 Matrice res(tilde) ;
373 for (int i=0 ; i<n-2 ; i++)
374 for (int j=0 ; j<n ; j++)
375 res.set(i, j) = (tilde(i, j)+tilde(i+1, j)) ;
376
377 tab[nb_dejafait] = new Matrice(res) ;
378 nb_dejafait ++ ;
379 return res ;
380 }
381
382 // Cas ou le calcul a deja ete effectue :
383 else
384 return *tab[indice] ;
385}
386
387
388//Cas dzpuis == 2
389Matrice _cl_poisson_2d_r_chebu_deux (const Matrice &source, int l) {
390 int n = source.get_dim(0) ;
391 assert (n == source.get_dim(1)) ;
392
393
394 const int nmax = 200 ; // Nombre de Matrices stockees
395 static Matrice* tab[nmax] ; // les matrices calculees
396 static int nb_dejafait = 0 ; // nbre de matrices calculees
397 static int l_dejafait[nmax] ;
398 static int nr_dejafait[nmax] ;
399
400 int indice = -1 ;
401
402 // On determine si la matrice a deja ete calculee :
403 for (int conte=0 ; conte<nb_dejafait ; conte ++)
404 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
405 indice = conte ;
406
407 // Calcul a faire :
408 if (indice == -1) {
409 if (nb_dejafait >= nmax) {
410 cout << "_cl_poisson_2d_r_chebu_deux : trop de matrices" << endl ;
411 abort() ;
412 exit (-1) ;
413 }
414
415 l_dejafait[nb_dejafait] = l ;
416 nr_dejafait[nb_dejafait] = n ;
417
418 Matrice barre(source) ;
419
420 int dirac = 1 ;
421 for (int i=0 ; i<n-2 ; i++) {
422 for (int j=0 ; j<n ; j++)
423 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
424 if (i==0) dirac = 0 ;
425 }
426
427 Matrice tilde(barre) ;
428 for (int i=0 ; i<n-4 ; i++)
429 for (int j=0 ; j<n ; j++)
430 tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
431
432 Matrice res(tilde) ;
433 for (int i=0 ; i<n-4 ; i++)
434 for (int j=0 ; j<n ; j++)
435 res.set(i, j) = (tilde(i, j)+tilde(i+1, j)) ;
436
437 return res ;
438 }
439
440 // Cas ou le calcul a deja ete effectue :
441 else
442 return *tab[indice] ;
443}
444
446 if (ope_mat == 0x0)
447 do_ope_mat() ;
448
449 if (ope_cl != 0x0)
450 delete ope_cl ;
451
452 // Routines de derivation
453 static Matrice (*cl_poisson_2d[MAX_BASE])(const Matrice&, int, double, int);
454 static int nap = 0 ;
455
456 // Premier appel
457 if (nap==0) {
458 nap = 1 ;
459 for (int i=0 ; i<MAX_BASE ; i++) {
460 cl_poisson_2d[i] = _cl_poisson_2d_pas_prevu ;
461 }
462 // Les routines existantes
463 cl_poisson_2d[R_CHEBP >> TRA_R] = _cl_poisson_2d_r_chebp ;
464 cl_poisson_2d[R_CHEBI >> TRA_R] = _cl_poisson_2d_r_chebi ;
465 cl_poisson_2d[R_CHEB >> TRA_R] = _cl_poisson_2d_r_cheb ;
466 cl_poisson_2d[R_CHEBU >> TRA_R] = _cl_poisson_2d_r_chebu ;
467 }
468 ope_cl = new Matrice(cl_poisson_2d[base_r](*ope_mat, l_quant, beta/alpha,
469 dzpuis)) ;
470}
471
472
473}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:260
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double beta
Parameter of the associated mapping.
double alpha
Parameter of the associated mapping.
int base_r
Radial basis of decomposition.
Matrice * ope_cl
Pointer on the banded-matrix of the operator.
int dzpuis
the associated dzpuis, if in the compactified domain.
virtual void do_ope_cl() const
Computes the banded-matrix of the operator.
int l_quant
quantum number
virtual void do_ope_mat() const
Computes the matrix of the operator.
#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