LORENE
ope_poisson_2d_mat.C
1/*
2 * Copyright (c) 2003 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_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_mat.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $" ;
22
23/*
24 * $Id: ope_poisson_2d_mat.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_mat.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 // Routine pour les cas non prevus --
36 //-----------------------------------
37
38namespace Lorene {
39Matrice _poisson_2d_mat_pas_prevu(int, int, double, double, int) {
40 cout << "laplacien pas prevu..." << endl ;
41 abort() ;
42 exit(-1) ;
43 Matrice res(1, 1) ;
44 return res;
45}
46
47
48 //-------------------------
49 //-- CAS R_CHEBP -----
50 //--------------------------
51
52
53Matrice _poisson_2d_mat_r_chebp (int n, int l, double, double, int) {
54
55 const int nmax = 100 ;// Nombre de Matrices stockees
56 static Matrice* tab[nmax] ; // les matrices calculees
57 static int nb_dejafait = 0 ; // nbre de matrices calculees
58 static int l_dejafait[nmax] ;
59 static int nr_dejafait[nmax] ;
60
61 int indice = -1 ;
62
63 // On determine si la matrice a deja ete calculee :
64 for (int conte=0 ; conte<nb_dejafait ; conte ++)
65 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
66 indice = conte ;
67
68 // Calcul a faire :
69 if (indice == -1) {
70 if (nb_dejafait >= nmax) {
71 cout << "_poisson_2d_mat_r_chebp : trop de matrices" << endl ;
72 abort() ;
73 exit (-1) ;
74 }
75
76
77 l_dejafait[nb_dejafait] = l ;
78 nr_dejafait[nb_dejafait] = n ;
79
80 Matrice dd(n, n) ;
81 dd.set_etat_qcq() ;
82 Matrice xd(n, n) ;
83 xd.set_etat_qcq() ;
84 Matrice xx(n, n) ;
85 xx.set_etat_qcq() ;
86
87 double* vect = new double[n] ;
88
89 for (int i=0 ; i<n ; i++) {
90 for (int j=0 ; j<n ; j++)
91 vect[j] = 0 ;
92 vect[i] = 1 ;
93 d2sdx2_1d (n, &vect, R_CHEBP) ;
94
95 for (int j=0 ; j<n ; j++)
96 dd.set(j, i) = vect[j] ;
97 }
98
99 for (int i=0 ; i<n ; i++) {
100 for (int j=0 ; j<n ; j++)
101 vect[j] = 0 ;
102 vect[i] = 1 ;
103 sxdsdx_1d (n, &vect, R_CHEBP) ;
104 for (int j=0 ; j<n ; j++)
105 xd.set(j, i) = vect[j] ;
106
107 }
108
109 for (int i=0 ; i<n ; i++) {
110 for (int j=0 ; j<n ; j++)
111 vect[j] = 0 ;
112 vect[i] = 1 ;
113 sx2_1d (n, &vect, R_CHEBP) ;
114 for (int j=0 ; j<n ; j++)
115 xx.set(j, i) = vect[j] ;
116 }
117
118 delete [] vect ;
119
120 Matrice res(n, n) ;
121 res = dd+xd-l*l*xx ;
122 tab[nb_dejafait] = new Matrice(res) ;
123 nb_dejafait ++ ;
124 return res ;
125 }
126
127 // Cas ou le calcul a deja ete effectue :
128 else
129 return *tab[indice] ;
130}
131
132
133
134 //------------------------
135 //-- CAS R_CHEBI ----
136 //------------------------
137
138
139Matrice _poisson_2d_mat_r_chebi (int n, int l, double, double, int) {
140
141 const int nmax = 100 ;// Nombre de Matrices stockees
142 static Matrice* tab[nmax] ; // les matrices calculees
143 static int nb_dejafait = 0 ; // nbre de matrices calculees
144 static int l_dejafait[nmax] ;
145 static int nr_dejafait[nmax] ;
146
147 int indice = -1 ;
148
149 // On determine si la matrice a deja ete calculee :
150 for (int conte=0 ; conte<nb_dejafait ; conte ++)
151 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
152 indice = conte ;
153
154 // Calcul a faire :
155 if (indice == -1) {
156 if (nb_dejafait >= nmax) {
157 cout << "_poisson_2d_mat_r_chebi : trop de matrices" << endl ;
158 abort() ;
159 exit (-1) ;
160 }
161
162
163 l_dejafait[nb_dejafait] = l ;
164 nr_dejafait[nb_dejafait] = n ;
165
166 Matrice dd(n, n) ;
167 dd.set_etat_qcq() ;
168 Matrice xd(n, n) ;
169 xd.set_etat_qcq() ;
170 Matrice xx(n, n) ;
171 xx.set_etat_qcq() ;
172
173 double* vect = new double[n] ;
174
175 for (int i=0 ; i<n ; i++) {
176 for (int j=0 ; j<n ; j++)
177 vect[j] = 0 ;
178 vect[i] = 1 ;
179 d2sdx2_1d (n, &vect, R_CHEBI) ; // appel dans le cas impair
180 for (int j=0 ; j<n ; j++)
181 dd.set(j, i) = vect[j] ;
182 }
183
184 for (int i=0 ; i<n ; i++) {
185 for (int j=0 ; j<n ; j++)
186 vect[j] = 0 ;
187 vect[i] = 1 ;
188 sxdsdx_1d (n, &vect, R_CHEBI) ;
189 for (int j=0 ; j<n ; j++)
190 xd.set(j, i) = vect[j] ;
191 }
192
193 for (int i=0 ; i<n ; i++) {
194 for (int j=0 ; j<n ; j++)
195 vect[j] = 0 ;
196 vect[i] = 1 ;
197 sx2_1d (n, &vect, R_CHEBI) ;
198 for (int j=0 ; j<n ; j++)
199 xx.set(j, i) = vect[j] ;
200 }
201
202 delete [] vect ;
203
204 Matrice res(n, n) ;
205 res = dd+xd-l*l*xx ;
206 tab[nb_dejafait] = new Matrice(res) ;
207 nb_dejafait ++ ;
208 return res ;
209 }
210
211 // Cas ou le calcul a deja ete effectue :
212 else
213 return *tab[indice] ;
214}
215
216
217
218
219 //-------------------------
220 //-- CAS R_CHEBU -----
221 //-------------------------
222
223Matrice _poisson_2d_mat_r_chebu_deux(int,int) ;
224Matrice _poisson_2d_mat_r_chebu_trois(int,int) ;
225Matrice _poisson_2d_mat_r_chebu_quatre(int,int) ;
226
227Matrice _poisson_2d_mat_r_chebu( int n, int l, double, double, int puis) {
228 Matrice res(n, n) ;
229 res.set_etat_qcq() ;
230 switch (puis) {
231 case 4 :
232 res = _poisson_2d_mat_r_chebu_quatre (n, l) ;
233 break ;
234 case 3 :
235 res = _poisson_2d_mat_r_chebu_trois (n, l) ;
236 break ;
237 case 2 :
238 res = _poisson_2d_mat_r_chebu_deux (n, l) ;
239 break ;
240 default :
241 abort() ;
242 exit(-1) ;
243 }
244 return res ;
245}
246
247 // Cas ou dzpuis = 4
248Matrice _poisson_2d_mat_r_chebu_quatre (int n, int l) {
249
250 const int nmax = 200 ;// Nombre de Matrices stockees
251 static Matrice* tab[nmax] ; // les matrices calculees
252 static int nb_dejafait = 0 ; // nbre de matrices calculees
253 static int l_dejafait[nmax] ;
254 static int nr_dejafait[nmax] ;
255
256 int indice = -1 ;
257
258 // On determine si la matrice a deja ete calculee :
259 for (int conte=0 ; conte<nb_dejafait ; conte ++)
260 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
261 indice = conte ;
262
263 // Calcul a faire :
264 if (indice == -1) {
265 if (nb_dejafait >= nmax) {
266 cout << "_poisson_2d_mat_r_chebu_quatre : trop de matrices" << endl ;
267 abort() ;
268 exit (-1) ;
269 }
270
271
272 l_dejafait[nb_dejafait] = l ;
273 nr_dejafait[nb_dejafait] = n ;
274
275 Matrice dd(n, n) ;
276 dd.set_etat_qcq() ;
277 Matrice dx (n,n) ;
278 dx.set_etat_qcq() ;
279 Matrice xx(n, n) ;
280 xx.set_etat_qcq() ;
281
282 double* vect = new double[n] ;
283
284 for (int i=0 ; i<n ; i++) {
285 for (int j=0 ; j<n ; j++)
286 vect[j] = 0 ;
287 vect[i] = 1 ;
288 d2sdx2_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
289 for (int j=0 ; j<n ; j++)
290 dd.set(j, i) = vect[j] ;
291 }
292
293 for (int i=0 ; i<n ; i++) {
294 for (int j=0 ; j<n ; j++)
295 vect[j] = 0 ;
296 vect[i] = 1 ;
297 sxdsdx_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
298 for (int j=0 ; j<n ; j++)
299 dx.set(j, i) = vect[j] ;
300 }
301
302 for (int i=0 ; i<n ; i++) {
303 for (int j=0 ; j<n ; j++)
304 vect[j] = 0 ;
305 vect[i] = 1 ;
306 sx2_1d (n, &vect, R_CHEBU) ;
307 for (int j=0 ; j<n ; j++)
308 xx.set(j, i) = vect[j] ;
309 }
310
311 delete [] vect ;
312
313 Matrice res(n, n) ;
314 res = dd+dx-l*l*xx ;
315 tab[nb_dejafait] = new Matrice(res) ;
316 nb_dejafait ++ ;
317 return res ;
318 }
319
320 // Cas ou le calcul a deja ete effectue :
321 else
322 return *tab[indice] ;
323}
324
325// Cas ou dzpuis =3
326Matrice _poisson_2d_mat_r_chebu_trois (int n, int l) {
327
328 const int nmax = 200 ;// Nombre de Matrices stockees
329 static Matrice* tab[nmax] ; // les matrices calculees
330 static int nb_dejafait = 0 ; // nbre de matrices calculees
331 static int l_dejafait[nmax] ;
332 static int nr_dejafait[nmax] ;
333
334 int indice = -1 ;
335
336 // On determine si la matrice a deja ete calculee :
337 for (int conte=0 ; conte<nb_dejafait ; conte ++)
338 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
339 indice = conte ;
340
341 // Calcul a faire :
342 if (indice == -1) {
343 if (nb_dejafait >= nmax) {
344 cout << "_poisson_2d_mat_r_chebu_trois : trop de matrices" << endl ;
345 abort() ;
346 exit (-1) ;
347 }
348
349
350 l_dejafait[nb_dejafait] = l ;
351 nr_dejafait[nb_dejafait] = n ;
352
353 Matrice dd(n, n) ;
354 dd.set_etat_qcq() ;
355 Matrice dx(n, n) ;
356 dx.set_etat_qcq() ;
357 Matrice xx(n, n) ;
358 xx.set_etat_qcq() ;
359
360 double* vect = new double[n] ;
361 double* auxi = new double[n] ;
362
363 for (int i=0 ; i<n ; i++) {
364
365 for (int j=0 ; j<n ; j++)
366 vect[j] = 0 ;
367 vect[i] = 1 ;
368 d2sdx2_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
369 mult_xm1_1d_cheb (n, vect, auxi) ;
370 for (int j=0 ; j<n ; j++)
371 dd.set(j, i) = auxi[j] ;
372 }
373
374 for (int i=0 ; i<n ; i++) {
375 for (int j=0 ; j<n ; j++)
376 vect[j] = 0 ;
377 vect[i] = 1 ;
378 dsdx_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
379 for (int j=0 ; j<n ; j++)
380 dx.set(j, i) = vect[j] ;
381 }
382
383
384 for (int i=0 ; i<n ; i++) {
385 for (int j=0 ; j<n ; j++)
386 vect[j] = 0 ;
387 vect[i] = 1 ;
388 sxm1_1d_cheb (n, vect) ;
389 for (int j=0 ; j<n ; j++)
390 xx.set(j, i) = vect[j] ;
391 }
392
393 delete [] vect ;
394 delete [] auxi ;
395
396 Matrice res(n, n) ;
397 res = dd+dx-l*l*xx ;
398 tab[nb_dejafait] = new Matrice(res) ;
399 nb_dejafait ++ ;
400 return res ;
401 }
402
403 // Cas ou le calcul a deja ete effectue :
404 else
405 return *tab[indice] ;
406}
407
408
409 //Cas ou dzpuis = 2
410Matrice _poisson_2d_mat_r_chebu_deux (int n, int l) {
411
412 const int nmax = 200 ;// Nombre de Matrices stockees
413 static Matrice* tab[nmax] ; // les matrices calculees
414 static int nb_dejafait = 0 ; // nbre de matrices calculees
415 static int l_dejafait[nmax] ;
416 static int nr_dejafait[nmax] ;
417
418 int indice = -1 ;
419
420 // On determine si la matrice a deja ete calculee :
421 for (int conte=0 ; conte<nb_dejafait ; conte ++)
422 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
423 indice = conte ;
424
425 // Calcul a faire :
426 if (indice == -1) {
427 if (nb_dejafait >= nmax) {
428 cout << "_poisson_2d_mat_r_chebu_deux : trop de matrices" << endl ;
429 abort() ;
430 exit (-1) ;
431 }
432
433
434 l_dejafait[nb_dejafait] = l ;
435 nr_dejafait[nb_dejafait] = n ;
436
437 Matrice dx (n,n) ;
438 dx.set_etat_qcq() ;
439 Matrice res(n, n) ;
440 res.set_etat_qcq() ;
441
442
443 double* vect = new double[n] ;
444
445 double* x2vect = new double[n] ;
446
447 for (int i=0 ; i<n ; i++) {
448 for (int j=0 ; j<n ; j++)
449 vect[j] = 0 ;
450 vect[i] = 1 ;
451 d2sdx2_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
452 mult2_xm1_1d_cheb (n, vect, x2vect) ; // multiplication par (x-1)^2
453 for (int j=0 ; j<n ; j++)
454 res.set(j, i) = x2vect[j] ;
455 }
456
457 for (int i=0 ; i<n ; i++) {
458 for (int j=0 ; j<n ; j++)
459 vect[j] = 0 ;
460 vect[i] = 1 ;
461 dsdx_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
462 mult_xm1_1d_cheb (n, vect, x2vect) ; // multiplication par (x-1)
463 for (int j=0 ; j<n ; j++)
464 dx.set(j, i) = x2vect[j] ;
465 }
466
467 delete [] vect ;
468 delete [] x2vect ;
469
470 res = res + dx ;
471
472 for (int i=0 ; i<n ; i++)
473 res.set(i, i) -= l*l ;
474
475 tab[nb_dejafait] = new Matrice(res) ;
476 nb_dejafait ++ ;
477 return res ;
478 }
479
480 // Cas ou le calcul a deja ete effectue :
481 else
482 return *tab[indice] ;
483}
484
485 //-------------------------
486 //-- CAS R_CHEB -----
487 //-----------------------
488
489
490Matrice _poisson_2d_mat_r_cheb (int n, int l, double alf, double bet, int) {
491
492 double echelle = bet / alf ;
493
494 const int nmax = 200 ;// Nombre de Matrices stockees
495 static Matrice* tab[nmax] ; // les matrices calculees
496 static int nb_dejafait = 0 ; // nbre de matrices calculees
497 static int l_dejafait[nmax] ;
498 static int nr_dejafait[nmax] ;
499 static double vieux_echelle = 0;
500
501 // Si on a change l'echelle : on detruit tout :
502 if (vieux_echelle != echelle) {
503 for (int i=0 ; i<nb_dejafait ; i++) {
504 l_dejafait[i] = -1 ;
505 nr_dejafait[i] = -1 ;
506 delete tab[i] ;
507 }
508
509 nb_dejafait = 0 ;
510 vieux_echelle = echelle ;
511 }
512
513 int indice = -1 ;
514
515 // On determine si la matrice a deja ete calculee :
516 for (int conte=0 ; conte<nb_dejafait ; conte ++)
517 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
518 indice = conte ;
519
520 // Calcul a faire :
521 if (indice == -1) {
522 if (nb_dejafait >= nmax) {
523 cout << "_poisson_2d_mat_r_cheb : trop de matrices" << endl ;
524 abort() ;
525 exit (-1) ;
526 }
527
528
529 l_dejafait[nb_dejafait] = l ;
530 nr_dejafait[nb_dejafait] = n ;
531
532 Matrice dd(n, n) ;
533 dd.set_etat_qcq() ;
534 Matrice xd(n, n) ;
535 xd.set_etat_qcq() ;
536 Matrice xx(n, n) ;
537 xx.set_etat_qcq() ;
538
539 double* vect = new double[n] ;
540
541 for (int i=0 ; i<n ; i++) {
542 for (int j=0 ; j<n ; j++)
543 vect[j] = 0 ;
544 vect[i] = 1 ;
545 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
546 for (int j=0 ; j<n ; j++)
547 dd.set(j, i) = vect[j]*echelle*echelle ;
548 }
549
550 for (int i=0 ; i<n ; i++) {
551 for (int j=0 ; j<n ; j++)
552 vect[j] = 0 ;
553 vect[i] = 1 ;
554 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
555 multx_1d (n, &vect, R_CHEB) ;
556 for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
557 dd.set(j, i) += 2*echelle*vect[j] ;
558 }
559
560 for (int i=0 ; i<n ; i++) {
561 for (int j=0 ; j<n ; j++)
562 vect[j] = 0 ;
563 vect[i] = 1 ;
564 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
565 multx_1d (n, &vect, R_CHEB) ;
566 multx_1d (n, &vect, R_CHEB) ;
567 for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
568 dd.set(j, i) += vect[j] ;
569 }
570
571 for (int i=0 ; i<n ; i++) {
572 for (int j=0 ; j<n ; j++)
573 vect[j] = 0 ;
574 vect[i] = 1 ;
575 sxdsdx_1d (n, &vect, R_CHEB) ;
576 for (int j=0 ; j<n ; j++)
577 xd.set(j, i) = vect[j]*echelle ;
578 }
579
580 for (int i=0 ; i<n ; i++) {
581 for (int j=0 ; j<n ; j++)
582 vect[j] = 0 ;
583 vect[i] = 1 ;
584 sxdsdx_1d (n, &vect, R_CHEB) ;
585 multx_1d (n, &vect, R_CHEB) ;
586 for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
587 xd.set(j, i) += vect[j] ;
588 }
589
590 for (int i=0 ; i<n ; i++) {
591 for (int j=0 ; j<n ; j++)
592 vect[j] = 0 ;
593 vect[i] = 1 ;
594 sx2_1d (n, &vect, R_CHEB) ;
595 for (int j=0 ; j<n ; j++)
596 xx.set(j, i) = vect[j] ;
597 }
598
599 delete [] vect ;
600
601 Matrice res(n, n) ;
602 res = dd+xd-l*l*xx ;
603 tab[nb_dejafait] = new Matrice(res) ;
604 nb_dejafait ++ ;
605 return res ;
606 }
607
608 // Cas ou le calcul a deja ete effectue :
609 else
610 return *tab[indice] ;
611}
612
613
615 if (ope_mat != 0x0)
616 delete ope_mat ;
617
618 // Routines de derivation
620 static int nap = 0 ;
621
622 // Premier appel
623 if (nap==0) {
624 nap = 1 ;
625 for (int i=0 ; i<MAX_BASE ; i++) {
626 poisson_2d_mat[i] = _poisson_2d_mat_pas_prevu ;
627 }
628 // Les routines existantes
629 poisson_2d_mat[R_CHEB >> TRA_R] = _poisson_2d_mat_r_cheb ;
630 poisson_2d_mat[R_CHEBP >> TRA_R] = _poisson_2d_mat_r_chebp ;
631 poisson_2d_mat[R_CHEBI >> TRA_R] = _poisson_2d_mat_r_chebi ;
632 poisson_2d_mat[R_CHEBU >> TRA_R] = _poisson_2d_mat_r_chebu ;
633 }
635}
636}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
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.
int nr
Number of radial points.
int dzpuis
the associated dzpuis, if in the compactified domain.
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