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 
21 
22 
23 /*
24  * $Id: ope_poisson_2d_non_dege.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_non_dege.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  //------------------------------------
36  // Routine pour les cas non prevus --
37  //-----------------------------------
38 
39 namespace Lorene {
40 Matrice _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 
53 Matrice _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 
123 Matrice _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 
194 Matrice _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  //-------------------
264 Matrice _poisson_2d_non_dege_r_chebu_quatre (const Matrice&, int) ;
265 Matrice _poisson_2d_non_dege_r_chebu_trois (const Matrice&, int) ;
266 Matrice _poisson_2d_non_dege_r_chebu_deux (const Matrice&, int) ;
267 
268 Matrice _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 ;
285 Matrice _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 ;
347 Matrice _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 ;
392 Matrice _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  }
478  non_dege = new Matrice(poisson_2d_non_dege[base_r](*ope_cl, l_quant,
479  beta/alpha, dzpuis)) ;
480 }
481 }
double alpha
Parameter of the associated mapping.
Matrice * ope_cl
Pointer on the banded-matrix of the operator.
virtual void do_non_dege() const
Computes the non-degenerated matrix of the operator.
double beta
Parameter of the associated mapping.
int dzpuis
the associated dzpuis, if in the compactified domain.
Lorene prototypes.
Definition: app_hor.h:67
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
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:263
virtual void do_ope_cl() const
Computes the banded-matrix of the operator.
#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
Matrice * non_dege
Pointer on the non-degenerated matrix of the operator.
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166