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