LORENE
laplacien_mat.C
1 /*
2  * Copyright (c) 1999-2001 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 as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 
24 
25 /*
26  * $Id: laplacien_mat.C,v 1.12 2016/12/05 16:18:09 j_novak Exp $
27  * $Log: laplacien_mat.C,v $
28  * Revision 1.12 2016/12/05 16:18:09 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.11 2014/10/13 08:53:29 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.10 2014/10/06 15:16:08 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.9 2007/12/11 15:28:22 jl_cornou
38  * Jacobi(0,2) polynomials partially implemented
39  *
40  * Revision 1.8 2007/06/06 07:43:28 p_grandclement
41  * nmax increased to 200
42  *
43  * Revision 1.7 2005/01/27 10:19:43 j_novak
44  * Now using Diff operators to build the matrices.
45  *
46  * Revision 1.6 2004/10/05 15:44:21 j_novak
47  * Minor speed enhancements.
48  *
49  * Revision 1.5 2004/02/06 10:53:54 j_novak
50  * New dzpuis = 5 -> dzpuis = 3 case (not ready yet).
51  *
52  * Revision 1.4 2003/01/31 08:49:58 e_gourgoulhon
53  * Increased the number nmax of stored matrices from 100 to 200.
54  *
55  * Revision 1.3 2002/10/16 14:37:11 j_novak
56  * Reorganization of #include instructions of standard C++, in order to
57  * use experimental version 3 of gcc.
58  *
59  * Revision 1.2 2002/10/09 12:47:32 j_novak
60  * Execution speed improved
61  *
62  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
63  * LORENE
64  *
65  * Revision 2.15 2000/05/22 13:36:33 phil
66  * ajout du cas dzpuis == 3
67  *
68  * Revision 2.14 2000/01/04 19:00:09 phil
69  * Double nmax
70  *
71  * Revision 2.13 1999/10/11 14:27:26 phil
72  * & -> &&
73  *
74  * Revision 2.12 1999/09/30 09:17:11 phil
75  * remplacement && en & et initialisation des variables statiques.
76  *
77  * Revision 2.11 1999/09/17 15:19:30 phil
78  * correction de definition de nmax
79  *
80  * Revision 2.10 1999/09/03 14:07:15 phil
81  * pas de modif
82  *
83  * Revision 2.9 1999/07/08 09:54:20 phil
84  * *** empty log message ***
85  *
86  * Revision 2.8 1999/07/07 10:02:39 phil
87  * Passage aux operateurs 1d
88  *
89  * Revision 2.7 1999/06/23 12:34:07 phil
90  * ajout de dzpuis = 2
91  *
92  * Revision 2.6 1999/04/28 10:45:54 phil
93  * augmentation de NMAX a 50
94  *
95  * Revision 2.5 1999/04/19 14:03:42 phil
96  * *** empty log message ***
97  *
98  * Revision 2.4 1999/04/16 13:15:52 phil
99  * *** empty log message ***
100  *
101  * Revision 2.3 1999/04/14 13:57:26 phil
102  * Sauvegarde des Matrices deja calculees
103  *
104  * Revision 2.2 1999/04/13 13:58:30 phil
105  * ajout proto.h
106  *
107  * Revision 2.1 1999/04/07 14:22:17 phil
108  * *** empty log message ***
109  *
110  * Revision 2.0 1999/04/07 14:09:41 phil
111  * *** empty log message ***
112  *
113  *
114  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/laplacien_mat.C,v 1.12 2016/12/05 16:18:09 j_novak Exp $
115  *
116  */
117 
118 //fichiers includes
119 #include <cstdio>
120 #include <cstdlib>
121 
122 #include "diff.h"
123 #include "proto.h"
124 
125 /*
126  * Routine calculant l'operateur suivant :
127  *
128  * -R_CHEB : r^2d2sdx2+2*rdsdx-l*(l+1)Id
129  *
130  * -R_CHEBP et R_CHEBI : d2sdx2+2/r dsdx-l(l+1)/r^2
131  *
132  * -R_CHEBU : d2sdx2-l(l+1)/x^2
133  *
134  * -R_JACO02 : d2sdx2 + 2/r dsdx - l(l+1)/r^2
135  *
136  * Entree :
137  * -n nbre de points en r
138  -l voire operateur.
139  -echelle utile uniquement pour R_CHEB : represente beta/alpha
140  (cf mapping)
141 
142  - puis : exposant de multiplication dans la ZEC
143  - base_r : base de developpement
144 
145  Sortie :
146  La fonction renvoie la matrice.
147 
148  */
149  //-----------------------------------
150  // Routine pour les cas non prevus --
151  //-----------------------------------
152 
153 namespace Lorene {
154 Matrice _laplacien_mat_pas_prevu(int n, int l, double echelle, int puis) {
155  cout << "laplacien pas prevu..." << endl ;
156  cout << "n : " << n << endl ;
157  cout << "l : " << l << endl ;
158  cout << "puissance : " << puis << endl ;
159  cout << "echelle : " << echelle << endl ;
160  abort() ;
161  exit(-1) ;
162  Matrice res(1, 1) ;
163  return res;
164 }
165 
166 
167  //-------------------------
168  //-- CAS R_JACO02 -----
169  //-------------------------
170 
171 
172 Matrice _laplacien_mat_r_jaco02 (int n, int l, double, int) {
173 
174  const int nmax = 200 ;// Nombre de Matrices stockees
175  static Matrice* tab[nmax] ; // les matrices calculees
176  static int nb_dejafait = 0 ; // nbre de matrices calculees
177  static int l_dejafait[nmax] ;
178  static int nr_dejafait[nmax] ;
179 
180  int indice = -1 ;
181 
182  // On determine si la matrice a deja ete calculee :
183  for (int conte=0 ; conte<nb_dejafait ; conte ++)
184  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
185  indice = conte ;
186 
187  // Calcul a faire :
188  if (indice == -1) {
189  if (nb_dejafait >= nmax) {
190  cout << "_laplacien_mat_r_jaco02 : trop de matrices" << endl ;
191  abort() ;
192  exit (-1) ;
193  }
194 
195 
196  l_dejafait[nb_dejafait] = l ;
197  nr_dejafait[nb_dejafait] = n ;
198 
199  Diff_dsdx2 d2(R_JACO02, n) ;
200  Diff_sxdsdx sxd(R_JACO02, n) ;
201  Diff_sx2 sx2(R_JACO02, n) ;
202 
203  tab[nb_dejafait] = new Matrice(d2 + 2.*sxd -(l*(l+1))*sx2) ;
204 
205  indice = nb_dejafait ;
206  nb_dejafait ++ ;
207  }
208 
209  return *tab[indice] ;
210 }
211 
212 
213  //-------------------------
214  //-- CAS R_CHEBP -----
215  //--------------------------
216 
217 
218 Matrice _laplacien_mat_r_chebp (int n, int l, double, int) {
219 
220  const int nmax = 200 ;// Nombre de Matrices stockees
221  static Matrice* tab[nmax] ; // les matrices calculees
222  static int nb_dejafait = 0 ; // nbre de matrices calculees
223  static int l_dejafait[nmax] ;
224  static int nr_dejafait[nmax] ;
225 
226  int indice = -1 ;
227 
228  // On determine si la matrice a deja ete calculee :
229  for (int conte=0 ; conte<nb_dejafait ; conte ++)
230  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
231  indice = conte ;
232 
233  // Calcul a faire :
234  if (indice == -1) {
235  if (nb_dejafait >= nmax) {
236  cout << "_laplacien_mat_r_chebp : trop de matrices" << endl ;
237  abort() ;
238  exit (-1) ;
239  }
240 
241 
242  l_dejafait[nb_dejafait] = l ;
243  nr_dejafait[nb_dejafait] = n ;
244 
245  Diff_dsdx2 d2(R_CHEBP, n) ;
246  Diff_sxdsdx sxd(R_CHEBP, n) ;
247  Diff_sx2 sx2(R_CHEBP, n) ;
248 
249  tab[nb_dejafait] = new Matrice(d2 + 2.*sxd -(l*(l+1))*sx2) ;
250 
251  indice = nb_dejafait ;
252  nb_dejafait ++ ;
253  }
254 
255  return *tab[indice] ;
256 }
257 
258 
259 
260  //------------------------
261  //-- CAS R_CHEBI ----
262  //------------------------
263 
264 
265 Matrice _laplacien_mat_r_chebi (int n, int l, double, int) {
266 
267  const int nmax = 200 ;// Nombre de Matrices stockees
268  static Matrice* tab[nmax] ; // les matrices calculees
269  static int nb_dejafait = 0 ; // nbre de matrices calculees
270  static int l_dejafait[nmax] ;
271  static int nr_dejafait[nmax] ;
272 
273  int indice = -1 ;
274 
275  // On determine si la matrice a deja ete calculee :
276  for (int conte=0 ; conte<nb_dejafait ; conte ++)
277  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
278  indice = conte ;
279 
280  // Calcul a faire :
281  if (indice == -1) {
282  if (nb_dejafait >= nmax) {
283  cout << "_laplacien_mat_r_chebi : trop de matrices" << endl ;
284  abort() ;
285  exit (-1) ;
286  }
287 
288 
289  l_dejafait[nb_dejafait] = l ;
290  nr_dejafait[nb_dejafait] = n ;
291 
292  Diff_dsdx2 d2(R_CHEBI, n) ;
293  Diff_sxdsdx sxd(R_CHEBI, n) ;
294  Diff_sx2 sx2(R_CHEBI, n) ;
295 
296  tab[nb_dejafait] = new Matrice(d2 + 2.*sxd - (l*(l+1))*sx2) ;
297  indice = nb_dejafait ;
298  nb_dejafait ++ ;
299  }
300 
301  return *tab[indice] ;
302 }
303 
304 
305 
306 
307  //-------------------------
308  //-- CAS R_CHEBU -----
309  //-------------------------
310 
311 Matrice _laplacien_mat_r_chebu( int n, int l, double, int puis) {
312  Matrice res(n, n) ;
313  res.set_etat_qcq() ;
314  switch (puis) {
315  case 4 :
316  res = _laplacien_mat_r_chebu_quatre (n, l) ;
317  break ;
318  case 3 :
319  res = _laplacien_mat_r_chebu_trois (n, l) ;
320  break ;
321  case 2 :
322  res = _laplacien_mat_r_chebu_deux (n, l) ;
323  break ;
324  case 5 :
325  res = _laplacien_mat_r_chebu_cinq (n, l) ;
326  break ;
327  default :
328  abort() ;
329  exit(-1) ;
330  }
331  return res ;
332 }
333 
334  // Cas ou dzpuis = 4
335 Matrice _laplacien_mat_r_chebu_quatre (int n, int l) {
336 
337  const int nmax = 200 ;// Nombre de Matrices stockees
338  static Matrice* tab[nmax] ; // les matrices calculees
339  static int nb_dejafait = 0 ; // nbre de matrices calculees
340  static int l_dejafait[nmax] ;
341  static int nr_dejafait[nmax] ;
342 
343  int indice = -1 ;
344 
345  // On determine si la matrice a deja ete calculee :
346  for (int conte=0 ; conte<nb_dejafait ; conte ++)
347  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
348  indice = conte ;
349 
350  // Calcul a faire :
351  if (indice == -1) {
352  if (nb_dejafait >= nmax) {
353  cout << "_laplacien_mat_r_chebu_quatre : trop de matrices" << endl ;
354  abort() ;
355  exit (-1) ;
356  }
357 
358 
359  l_dejafait[nb_dejafait] = l ;
360  nr_dejafait[nb_dejafait] = n ;
361 
362  Diff_dsdx2 dd(R_CHEBU, n) ;
363  Diff_sx2 xx(R_CHEBU, n) ;
364 
365  tab[nb_dejafait] = new Matrice(dd-(l*(l+1))*xx) ;
366  indice = nb_dejafait ;
367  nb_dejafait ++ ;
368  }
369 
370  return *tab[indice] ;
371 }
372 
373 // Cas ou dzpuis =3
374 Matrice _laplacien_mat_r_chebu_trois (int n, int l) {
375 
376  const int nmax = 200 ;// Nombre de Matrices stockees
377  static Matrice* tab[nmax] ; // les matrices calculees
378  static int nb_dejafait = 0 ; // nbre de matrices calculees
379  static int l_dejafait[nmax] ;
380  static int nr_dejafait[nmax] ;
381 
382  int indice = -1 ;
383 
384  // On determine si la matrice a deja ete calculee :
385  for (int conte=0 ; conte<nb_dejafait ; conte ++)
386  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
387  indice = conte ;
388 
389  // Calcul a faire :
390  if (indice == -1) {
391  if (nb_dejafait >= nmax) {
392  cout << "_laplacien_mat_r_chebu_trois : trop de matrices" << endl ;
393  abort() ;
394  exit (-1) ;
395  }
396 
397 
398  l_dejafait[nb_dejafait] = l ;
399  nr_dejafait[nb_dejafait] = n ;
400 
401  Diff_xdsdx2 xd2(R_CHEBU, n) ;
402  Diff_sx sx(R_CHEBU, n) ;
403 
404  tab[nb_dejafait] = new Matrice(xd2 -(l*(l+1))*sx) ;
405 
406  indice = nb_dejafait ;
407  nb_dejafait ++ ;
408  }
409 
410  return *tab[indice] ;
411 }
412 
413 
414  //Cas ou dzpuis = 2
415 Matrice _laplacien_mat_r_chebu_deux (int n, int l) {
416 
417  const int nmax = 200 ;// Nombre de Matrices stockees
418  static Matrice* tab[nmax] ; // les matrices calculees
419  static int nb_dejafait = 0 ; // nbre de matrices calculees
420  static int l_dejafait[nmax] ;
421  static int nr_dejafait[nmax] ;
422 
423  int indice = -1 ;
424 
425  // On determine si la matrice a deja ete calculee :
426  for (int conte=0 ; conte<nb_dejafait ; conte ++)
427  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
428  indice = conte ;
429 
430  // Calcul a faire :
431  if (indice == -1) {
432  if (nb_dejafait >= nmax) {
433  cout << "_laplacien_mat_r_chebu_deux : trop de matrices" << endl ;
434  abort() ;
435  exit (-1) ;
436  }
437 
438 
439  l_dejafait[nb_dejafait] = l ;
440  nr_dejafait[nb_dejafait] = n ;
441 
442  Diff_x2dsdx2 x2dd(R_CHEBU, n) ;
443  Diff_id id(R_CHEBU, n) ;
444 
445  tab[nb_dejafait] = new Matrice(x2dd - (l*(l+1))*id) ;
446 
447  indice = nb_dejafait ;
448  nb_dejafait ++ ;
449  }
450 
451  return *tab[indice] ;
452 }
453 
454  //Cas ou dzpuis = 5
455 Matrice _laplacien_mat_r_chebu_cinq (int n, int l) {
456 
457  const int nmax = 200 ;// Nombre de Matrices stockees
458  static Matrice* tab[nmax] ; // les matrices calculees
459  static int nb_dejafait = 0 ; // nbre de matrices calculees
460  static int l_dejafait[nmax] ;
461  static int nr_dejafait[nmax] ;
462 
463  int indice = -1 ;
464 
465  // On determine si la matrice a deja ete calculee :
466  for (int conte=0 ; conte<nb_dejafait ; conte ++)
467  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
468  indice = conte ;
469 
470  // Calcul a faire :
471  if (indice == -1) {
472  if (nb_dejafait >= nmax) {
473  cout << "_laplacien_mat_r_chebu_cinq : trop de matrices" << endl ;
474  abort() ;
475  exit (-1) ;
476  }
477 
478 
479  l_dejafait[nb_dejafait] = l ;
480  nr_dejafait[nb_dejafait] = n ;
481 
482  Diff_x2dsdx2 x2dd(R_CHEBU, n) ;
483  Diff_xdsdx xd1(R_CHEBU, n) ;
484  Diff_id id(R_CHEBU, n) ;
485 
486  tab[nb_dejafait] = new Matrice( x2dd + 6.*xd1 + (6-l*(l+1))*id ) ;
487 
488  indice = nb_dejafait ;
489  nb_dejafait ++ ;
490  }
491 
492  return *tab[indice] ;
493 }
494 
495  //-------------------------
496  //-- CAS R_CHEB -----
497  //-----------------------
498 
499 
500 Matrice _laplacien_mat_r_cheb (int n, int l, double echelle, int) {
501 
502  const int nmax = 200 ;// Nombre de Matrices stockees
503  static Matrice* tab[nmax] ; // les matrices calculees
504  static int nb_dejafait = 0 ; // nbre de matrices calculees
505  static int l_dejafait[nmax] ;
506  static int nr_dejafait[nmax] ;
507  static double vieux_echelle = 0;
508 
509  // Si on a change l'echelle : on detruit tout :
510  if (vieux_echelle != echelle) {
511  for (int i=0 ; i<nb_dejafait ; i++) {
512  l_dejafait[i] = -1 ;
513  nr_dejafait[i] = -1 ;
514  delete tab[i] ;
515  }
516 
517  nb_dejafait = 0 ;
518  vieux_echelle = echelle ;
519  }
520 
521  int indice = -1 ;
522 
523  // On determine si la matrice a deja ete calculee :
524  for (int conte=0 ; conte<nb_dejafait ; conte ++)
525  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
526  indice = conte ;
527 
528  // Calcul a faire :
529  if (indice == -1) {
530  if (nb_dejafait >= nmax) {
531  cout << "_laplacien_mat_r_cheb : trop de matrices" << endl ;
532  abort() ;
533  exit (-1) ;
534  }
535 
536 
537  l_dejafait[nb_dejafait] = l ;
538  nr_dejafait[nb_dejafait] = n ;
539 
540  Diff_dsdx2 d2(R_CHEB, n) ;
541  Diff_xdsdx2 xd2(R_CHEB, n) ;
542  Diff_x2dsdx2 x2d2(R_CHEB, n) ;
543  Diff_dsdx d1(R_CHEB, n) ;
544  Diff_xdsdx xd1(R_CHEB, n) ;
545  Diff_id id(R_CHEB, n) ;
546 
547  tab[nb_dejafait] =
548  new Matrice(x2d2 + (2*echelle)*xd2 + (echelle*echelle)*d2
549  + 2*xd1 + (2*echelle)*d1 - (l*(l+1))*id) ;
550  indice = nb_dejafait ;
551  nb_dejafait ++ ;
552  }
553 
554  return *tab[indice] ;
555 }
556 
557 
558  //--------------------------
559  //- La routine a appeler ---
560  //----------------------------
561 Matrice laplacien_mat(int n, int l, double echelle, int puis, int base_r)
562 {
563 
564  // Routines de derivation
565  static Matrice (*laplacien_mat[MAX_BASE])(int, int, double, int) ;
566  static int nap = 0 ;
567 
568  // Premier appel
569  if (nap==0) {
570  nap = 1 ;
571  for (int i=0 ; i<MAX_BASE ; i++) {
572  laplacien_mat[i] = _laplacien_mat_pas_prevu ;
573  }
574  // Les routines existantes
575  laplacien_mat[R_CHEB >> TRA_R] = _laplacien_mat_r_cheb ;
576  laplacien_mat[R_CHEBU >> TRA_R] = _laplacien_mat_r_chebu ;
577  laplacien_mat[R_CHEBP >> TRA_R] = _laplacien_mat_r_chebp ;
578  laplacien_mat[R_CHEBI >> TRA_R] = _laplacien_mat_r_chebi ;
579  laplacien_mat[R_JACO02 >> TRA_R] = _laplacien_mat_r_jaco02 ;
580  }
581 
582  return laplacien_mat[base_r](n, l, echelle, puis) ;
583 }
584 
585 }
Lorene prototypes.
Definition: app_hor.h:67
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#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
#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