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
23char 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
150namespace Lorene {
151Matrice _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
169Matrice _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
215Matrice _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
262Matrice _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
308Matrice _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
332Matrice _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
371Matrice _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
412Matrice _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
452Matrice _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
497Matrice _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 //----------------------------
558Matrice 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}
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#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