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