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 
28 char 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 
76 namespace Lorene {
77 Matrice _ope_ptens_rr_mat_pas_prevu(int, int, double, int) ;
78 Matrice _ope_ptens_rr_mat_r_chebp(int, int, double, int) ;
79 Matrice _ope_ptens_rr_mat_r_chebi(int, int, double, int) ;
80 Matrice _ope_ptens_rr_mat_r_chebu(int, int, double, int) ;
81 Matrice _ope_ptens_rr_mat_r_cheb(int, int, double, int) ;
82 
83  //-----------------------------------
84  // Routine pour les cas non prevus --
85  //-----------------------------------
86 
87 Matrice _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 
105 Matrice _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 
191 Matrice _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 
275 Matrice _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 
365 Matrice _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 
491 Matrice 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 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180