LORENE
ope_ptens_rr_mat.C
1/*
2 * Builds the operator for Ope_pois_tens_rr
3 *
4 * (see file ope_elementary.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2004 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char ope_ptens_rr_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/ope_ptens_rr_mat.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $" ;
29
30/*
31 * $Id: ope_ptens_rr_mat.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
32 * $Log: ope_ptens_rr_mat.C,v $
33 * Revision 1.3 2014/10/13 08:53:34 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.2 2014/10/06 15:16:13 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.1 2004/12/23 16:30:15 j_novak
40 * New files and class for the solution of the rr component of the tensor Poisson
41 * equation.
42 *
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/ope_ptens_rr_mat.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
46 *
47 */
48
49//fichiers includes
50#include <cstdlib>
51
52#include "matrice.h"
53#include "type_parite.h"
54#include "proto.h"
55
56/*
57 * Routine caluclant l'operateur suivant :
58 *
59 * -R_CHEB : r^2d2sdx2+6*rdsdx+6-l*(l+1)Id
60 *
61 * -R_CHEBP et R_CHEBI : d2sdx2+6/r dsdx +(6-l(l+1))/r^2
62 *
63 * -R_CHEBU : d2sdx2+(6-l(l+1))/x^2
64 *
65 * Entree :
66 * -n nbre de points en r
67 * -l voire operateur. (doit etre >=2)
68 * -echelle utile uniquement pour R_CHEB : represente beta/alpha
69 * (cf mapping)
70 * - puis : exposant de multiplication dans la ZEC
71 * - base_r : base de developpement
72 * Sortie :
73 * La fonction renvoie la matrice.
74 */
75
76namespace Lorene {
77Matrice _ope_ptens_rr_mat_pas_prevu(int, int, double, int) ;
78Matrice _ope_ptens_rr_mat_r_chebp(int, int, double, int) ;
79Matrice _ope_ptens_rr_mat_r_chebi(int, int, double, int) ;
80Matrice _ope_ptens_rr_mat_r_chebu(int, int, double, int) ;
81Matrice _ope_ptens_rr_mat_r_cheb(int, int, double, int) ;
82
83 //-----------------------------------
84 // Routine pour les cas non prevus --
85 //-----------------------------------
86
87Matrice _ope_ptens_rr_mat_pas_prevu(int n, int l, double echelle, int puis) {
88 cout << "laplacien tens_rr pas prevu..." << endl ;
89 cout << "n : " << n << endl ;
90 cout << "l : " << l << endl ;
91 cout << "puissance : " << puis << endl ;
92 cout << "echelle : " << echelle << endl ;
93 abort() ;
94 exit(-1) ;
95 Matrice res(1, 1) ;
96 return res;
97}
98
99
100 //-------------------------
101 //---- CAS R_CHEBP ----
102 //-------------------------
103
104
105Matrice _ope_ptens_rr_mat_r_chebp (int n, int l, double, int) {
106
107 const int nmax = 100 ;// Nombre de Matrices stockees
108 static Matrice* tab[nmax] ; // les matrices calculees
109 static int nb_dejafait = 0 ; // nbre de matrices calculees
110 static int l_dejafait[nmax] ;
111 static int nr_dejafait[nmax] ;
112
113 int indice = -1 ;
114
115 // On determine si la matrice a deja ete calculee :
116 for (int conte=0 ; conte<nb_dejafait ; conte ++)
117 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
118 indice = conte ;
119
120 // Calcul a faire :
121 if (indice == -1) {
122 if (nb_dejafait >= nmax) {
123 cout << "_ope_ptens_rr_mat_r_chebp : trop de matrices" << endl ;
124 abort() ;
125 exit (-1) ;
126 }
127
128
129 l_dejafait[nb_dejafait] = l ;
130 nr_dejafait[nb_dejafait] = n ;
131
132 Matrice dd(n, n) ;
133 dd.set_etat_qcq() ;
134 Matrice xd(n, n) ;
135 xd.set_etat_qcq() ;
136 Matrice xx(n, n) ;
137 xx.set_etat_qcq() ;
138
139 double* vect = new double[n] ;
140
141 for (int i=0 ; i<n ; i++) {
142 for (int j=0 ; j<n ; j++)
143 vect[j] = 0 ;
144 vect[i] = 1 ;
145 d2sdx2_1d (n, &vect, R_CHEBP) ;
146
147 for (int j=0 ; j<n ; j++)
148 dd.set(j, i) = vect[j] ;
149 }
150
151 for (int i=0 ; i<n ; i++) {
152 for (int j=0 ; j<n ; j++)
153 vect[j] = 0 ;
154 vect[i] = 1 ;
155 sxdsdx_1d (n, &vect, R_CHEBP) ;
156 for (int j=0 ; j<n ; j++)
157 xd.set(j, i) = vect[j] ;
158
159 }
160
161 for (int i=0 ; i<n ; i++) {
162 for (int j=0 ; j<n ; j++)
163 vect[j] = 0 ;
164 vect[i] = 1 ;
165 sx2_1d (n, &vect, R_CHEBP) ;
166 for (int j=0 ; j<n ; j++)
167 xx.set(j, i) = vect[j] ;
168 }
169
170 delete [] vect ;
171
172 Matrice res(n, n) ;
173 res = dd+6*xd+(6-l*(l+1))*xx ;
174 tab[nb_dejafait] = new Matrice(res) ;
175 nb_dejafait ++ ;
176 return res ;
177 }
178
179 // Cas ou le calcul a deja ete effectue :
180 else
181 return *tab[indice] ;
182}
183
184
185
186 //------------------------
187 //-- CAS R_CHEBI ----
188 //------------------------
189
190
191Matrice _ope_ptens_rr_mat_r_chebi (int n, int l, double, int) {
192
193 const int nmax = 100 ;// Nombre de Matrices stockees
194 static Matrice* tab[nmax] ; // les matrices calculees
195 static int nb_dejafait = 0 ; // nbre de matrices calculees
196 static int l_dejafait[nmax] ;
197 static int nr_dejafait[nmax] ;
198
199 int indice = -1 ;
200
201 // On determine si la matrice a deja ete calculee :
202 for (int conte=0 ; conte<nb_dejafait ; conte ++)
203 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
204 indice = conte ;
205
206 // Calcul a faire :
207 if (indice == -1) {
208 if (nb_dejafait >= nmax) {
209 cout << "_ope_ptens_rr_mat_r_chebi : trop de matrices" << endl ;
210 abort() ;
211 exit (-1) ;
212 }
213
214
215 l_dejafait[nb_dejafait] = l ;
216 nr_dejafait[nb_dejafait] = n ;
217
218 Matrice dd(n, n) ;
219 dd.set_etat_qcq() ;
220 Matrice xd(n, n) ;
221 xd.set_etat_qcq() ;
222 Matrice xx(n, n) ;
223 xx.set_etat_qcq() ;
224
225 double* vect = new double[n] ;
226
227 for (int i=0 ; i<n ; i++) {
228 for (int j=0 ; j<n ; j++)
229 vect[j] = 0 ;
230 vect[i] = 1 ;
231 d2sdx2_1d (n, &vect, R_CHEBI) ; // appel dans le cas impair
232 for (int j=0 ; j<n ; j++)
233 dd.set(j, i) = vect[j] ;
234 }
235
236 for (int i=0 ; i<n ; i++) {
237 for (int j=0 ; j<n ; j++)
238 vect[j] = 0 ;
239 vect[i] = 1 ;
240 sxdsdx_1d (n, &vect, R_CHEBI) ;
241 for (int j=0 ; j<n ; j++)
242 xd.set(j, i) = vect[j] ;
243 }
244
245 for (int i=0 ; i<n ; i++) {
246 for (int j=0 ; j<n ; j++)
247 vect[j] = 0 ;
248 vect[i] = 1 ;
249 sx2_1d (n, &vect, R_CHEBI) ;
250 for (int j=0 ; j<n ; j++)
251 xx.set(j, i) = vect[j] ;
252 }
253
254 delete [] vect ;
255
256 Matrice res(n, n) ;
257 res = dd+6*xd+(6-l*(l+1))*xx ;
258 tab[nb_dejafait] = new Matrice(res) ;
259 nb_dejafait ++ ;
260 return res ;
261 }
262
263 // Cas ou le calcul a deja ete effectue :
264 else
265 return *tab[indice] ;
266}
267
268
269
270
271 //------------------------
272 //---- CAS R_CHEBU ----
273 //------------------------
274
275Matrice _ope_ptens_rr_mat_r_chebu( int n, int l, double, int puis) {
276
277 if (puis != 4) {
278 cout << "_ope_ptens_rr_mat_r_chebu : only the case dzpuis = 4 "
279 << '\n' << "is implemented! \n"
280 << "dzpuis = " << puis << endl ;
281 abort() ;
282 }
283
284 const int nmax = 200 ;// Nombre de Matrices stockees
285 static Matrice* tab[nmax] ; // les matrices calculees
286 static int nb_dejafait = 0 ; // nbre de matrices calculees
287 static int l_dejafait[nmax] ;
288 static int nr_dejafait[nmax] ;
289
290 int indice = -1 ;
291
292 // On determine si la matrice a deja ete calculee :
293 for (int conte=0 ; conte<nb_dejafait ; conte ++)
294 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
295 indice = conte ;
296
297 // Calcul a faire :
298 if (indice == -1) {
299 if (nb_dejafait >= nmax) {
300 cout << "_ope_ptens_rr_mat_r_chebu : trop de matrices" << endl ;
301 abort() ;
302 exit (-1) ;
303 }
304
305 l_dejafait[nb_dejafait] = l ;
306 nr_dejafait[nb_dejafait] = n ;
307
308 Matrice dd(n, n) ;
309 dd.set_etat_qcq() ;
310 Matrice xd(n, n) ;
311 xd.set_etat_qcq() ;
312 Matrice xx(n, n) ;
313 xx.set_etat_qcq() ;
314
315 double* vect = new double[n] ;
316
317 for (int i=0 ; i<n ; i++) {
318 for (int j=0 ; j<n ; j++)
319 vect[j] = 0 ;
320 vect[i] = 1 ;
321 d2sdx2_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
322 for (int j=0 ; j<n ; j++)
323 dd.set(j, i) = vect[j] ;
324 }
325
326 for (int i=0 ; i<n ; i++) {
327 for (int j=0 ; j<n ; j++)
328 vect[j] = 0 ;
329 vect[i] = 1 ;
330 dsdx_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
331 sxm1_1d_cheb (n, vect) ;
332 for (int j=0 ; j<n ; j++)
333 xd.set(j, i) = vect[j] ;
334 }
335
336 for (int i=0 ; i<n ; i++) {
337 for (int j=0 ; j<n ; j++)
338 vect[j] = 0 ;
339 vect[i] = 1 ;
340 sx2_1d (n, &vect, R_CHEBU) ;
341 for (int j=0 ; j<n ; j++)
342 xx.set(j, i) = vect[j] ;
343 }
344
345 delete [] vect ;
346
347 Matrice res(n, n) ;
348 res = dd - 4*xd + (6 -l*(l+1))*xx ;
349 tab[nb_dejafait] = new Matrice(res) ;
350 nb_dejafait ++ ;
351 return res ;
352 }
353
354 // Cas ou le calcul a deja ete effectue :
355 else
356 return *tab[indice] ;
357}
358
359
360 //-----------------------
361 //---- CAS R_CHEB ----
362 //-----------------------
363
364
365Matrice _ope_ptens_rr_mat_r_cheb (int n, int l, double echelle, int) {
366
367 const int nmax = 200 ;// Nombre de Matrices stockees
368 static Matrice* tab[nmax] ; // les matrices calculees
369 static int nb_dejafait = 0 ; // nbre de matrices calculees
370 static int l_dejafait[nmax] ;
371 static int nr_dejafait[nmax] ;
372 static double vieux_echelle = 0;
373
374 // Si on a change l'echelle : on detruit tout :
375 if (vieux_echelle != echelle) {
376 for (int i=0 ; i<nb_dejafait ; i++) {
377 l_dejafait[i] = -1 ;
378 nr_dejafait[i] = -1 ;
379 delete tab[i] ;
380 }
381
382 nb_dejafait = 0 ;
383 vieux_echelle = echelle ;
384 }
385
386 int indice = -1 ;
387
388 // On determine si la matrice a deja ete calculee :
389 for (int conte=0 ; conte<nb_dejafait ; conte ++)
390 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
391 indice = conte ;
392
393 // Calcul a faire :
394 if (indice == -1) {
395 if (nb_dejafait >= nmax) {
396 cout << "_ope_ptens_rr_mat_r_cheb : trop de matrices" << endl ;
397 abort() ;
398 exit (-1) ;
399 }
400
401
402 l_dejafait[nb_dejafait] = l ;
403 nr_dejafait[nb_dejafait] = n ;
404
405 Matrice dd(n, n) ;
406 dd.set_etat_qcq() ;
407 Matrice xd(n, n) ;
408 xd.set_etat_qcq() ;
409 Matrice xx(n, n) ;
410 xx.set_etat_qcq() ;
411
412 double* vect = new double[n] ;
413
414 for (int i=0 ; i<n ; i++) {
415 for (int j=0 ; j<n ; j++)
416 vect[j] = 0 ;
417 vect[i] = 1 ;
418 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
419 for (int j=0 ; j<n ; j++)
420 dd.set(j, i) = vect[j]*echelle*echelle ;
421 }
422
423 for (int i=0 ; i<n ; i++) {
424 for (int j=0 ; j<n ; j++)
425 vect[j] = 0 ;
426 vect[i] = 1 ;
427 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
428 multx_1d (n, &vect, R_CHEB) ;
429 for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
430 dd.set(j, i) += 2*echelle*vect[j] ;
431 }
432
433 for (int i=0 ; i<n ; i++) {
434 for (int j=0 ; j<n ; j++)
435 vect[j] = 0 ;
436 vect[i] = 1 ;
437 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
438 multx_1d (n, &vect, R_CHEB) ;
439 multx_1d (n, &vect, R_CHEB) ;
440 for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
441 dd.set(j, i) += vect[j] ;
442 }
443
444 for (int i=0 ; i<n ; i++) {
445 for (int j=0 ; j<n ; j++)
446 vect[j] = 0 ;
447 vect[i] = 1 ;
448 sxdsdx_1d (n, &vect, R_CHEB) ;
449 for (int j=0 ; j<n ; j++)
450 xd.set(j, i) = vect[j]*echelle ;
451 }
452
453 for (int i=0 ; i<n ; i++) {
454 for (int j=0 ; j<n ; j++)
455 vect[j] = 0 ;
456 vect[i] = 1 ;
457 sxdsdx_1d (n, &vect, R_CHEB) ;
458 multx_1d (n, &vect, R_CHEB) ;
459 for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
460 xd.set(j, i) += vect[j] ;
461 }
462
463 for (int i=0 ; i<n ; i++) {
464 for (int j=0 ; j<n ; j++)
465 vect[j] = 0 ;
466 vect[i] = 1 ;
467 sx2_1d (n, &vect, R_CHEB) ;
468 for (int j=0 ; j<n ; j++)
469 xx.set(j, i) = vect[j] ;
470 }
471
472 delete [] vect ;
473
474 Matrice res(n, n) ;
475 res = dd + 6*xd + (6 - l*(l+1))*xx ;
476 tab[nb_dejafait] = new Matrice(res) ;
477 nb_dejafait ++ ;
478 return res ;
479 }
480
481 // Cas ou le calcul a deja ete effectue :
482 else
483 return *tab[indice] ;
484}
485
486
487 //----------------------------
488 //--- La routine a appeler ---
489 //----------------------------
490
491Matrice ope_ptens_rr_mat(int n, int l, double echelle, int puis, int base_r)
492{
493
494 // Routines de derivation
495 static Matrice (*ope_ptens_rr_mat[MAX_BASE])(int, int, double, int) ;
496 static int nap = 0 ;
497
498 // Premier appel
499 if (nap==0) {
500 nap = 1 ;
501 for (int i=0 ; i<MAX_BASE ; i++) {
502 ope_ptens_rr_mat[i] = _ope_ptens_rr_mat_pas_prevu ;
503 }
504 // Les routines existantes
505 ope_ptens_rr_mat[R_CHEB >> TRA_R] = _ope_ptens_rr_mat_r_cheb ;
506 ope_ptens_rr_mat[R_CHEBU >> TRA_R] = _ope_ptens_rr_mat_r_chebu ;
507 ope_ptens_rr_mat[R_CHEBP >> TRA_R] = _ope_ptens_rr_mat_r_chebp ;
508 ope_ptens_rr_mat[R_CHEBI >> TRA_R] = _ope_ptens_rr_mat_r_chebi ;
509 }
510 assert (l>1) ;
511 Matrice res(ope_ptens_rr_mat[base_r](n, l, echelle, puis)) ;
512 return res ;
513}
514
515}
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#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