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 
21 
22 
23 /*
24  * $Id: ope_poisson_2d_mat.C,v 1.4 2016/12/05 16:18:12 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_mat.C,v 1.4 2016/12/05 16:18:12 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 
38 namespace Lorene {
39 Matrice _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 
53 Matrice _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 
139 Matrice _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 
223 Matrice _poisson_2d_mat_r_chebu_deux(int,int) ;
224 Matrice _poisson_2d_mat_r_chebu_trois(int,int) ;
225 Matrice _poisson_2d_mat_r_chebu_quatre(int,int) ;
226 
227 Matrice _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
248 Matrice _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
326 Matrice _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
410 Matrice _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 
490 Matrice _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
619  static Matrice (*poisson_2d_mat[MAX_BASE])(int, int, double, double, int);
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  }
634  ope_mat = new Matrice(poisson_2d_mat[base_r](nr, l_quant, alpha, beta, dzpuis)) ;
635 }
636 }
double alpha
Parameter of the associated mapping.
double beta
Parameter of the associated mapping.
int dzpuis
the associated dzpuis, if in the compactified domain.
Lorene prototypes.
Definition: app_hor.h:67
Matrice * ope_mat
Pointer on the matrix representation of the operator.
int l_quant
quantum number
int base_r
Radial basis of decomposition.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Matrix handling.
Definition: matrice.h:152
virtual void do_ope_mat() const
Computes the matrix of the operator.
int nr
Number of radial points.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166