LORENE
lap_cpt_mat.C
1 /*
2  * Copyright (c) 2000-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 lap_cpt_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/lap_cpt_mat.C,v 1.5 2014/10/13 08:53:29 j_novak Exp $" ;
24 
25 /*
26  * $Id: lap_cpt_mat.C,v 1.5 2014/10/13 08:53:29 j_novak Exp $
27  * $Log: lap_cpt_mat.C,v $
28  * Revision 1.5 2014/10/13 08:53:29 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.4 2014/10/06 15:16:08 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.3 2007/06/21 20:06:31 k_taniguchi
35  * nmax increased to 200
36  *
37  * Revision 1.2 2002/10/16 14:37:11 j_novak
38  * Reorganization of #include instructions of standard C++, in order to
39  * use experimental version 3 of gcc.
40  *
41  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
42  * LORENE
43  *
44  * Revision 2.0 2000/03/16 16:23:08 phil
45  * *** empty log message ***
46  *
47  *
48  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/lap_cpt_mat.C,v 1.5 2014/10/13 08:53:29 j_novak Exp $
49  *
50  */
51 
52 
53 
54 
55 //fichiers includes
56 #include <cstdio>
57 #include <cstdlib>
58 
59 #include "matrice.h"
60 #include "type_parite.h"
61 #include "proto.h"
62 
63 /*
64  * Routine renvoyant la matrice de l'operateur (1-x^2)*Laplacien f = s
65  * Pour l != 1 le resultat est donne en s est donne en Chebyshev et
66  * f en Gelerkin (T_i + T_{i+1} pour l pair et (2*i+3)T_i + (2*i+1)T_{i+1} pour
67  * l impair.
68  * Pour l=1 pas de probleme de singularite on reste donc en Chebyshev.
69  */
70 
71 
72  //-----------------------------------
73  // Routine pour les cas non prevus --
74  //-----------------------------------
75 
76 namespace Lorene {
77 Matrice _lap_cpt_mat_pas_prevu(int n, int l) {
78  cout << "laplacien * (1-r^2/R_0^2) pas prevu..." << endl ;
79  cout << "n : " << n << endl ;
80  cout << "l : " << l << endl ;
81  abort() ;
82  exit(-1) ;
83  Matrice res(1, 1) ; // On ne passe jamais ici de toute facon !
84  return res;
85 }
86 
87 
88  //-------------------------
89  //-- CAS R_CHEBP -----
90  //--------------------------
91 
92 
93 Matrice _lap_cpt_mat_r_chebp (int n, int l) {
94 
95  const int nmax = 200 ;// Nombre de Matrices stockees
96  static Matrice* tab[nmax] ; // les matrices calculees
97  static int nb_dejafait = 0 ; // nbre de matrices calculees
98  static int l_dejafait[nmax] ;
99  static int nr_dejafait[nmax] ;
100 
101  int indice = -1 ;
102 
103  // On determine si la matrice a deja ete calculee :
104  for (int conte=0 ; conte<nb_dejafait ; conte ++)
105  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
106  indice = conte ;
107 
108  // Calcul a faire :
109  if (indice == -1) {
110  if (nb_dejafait >= nmax) {
111  cout << "_laplacien_nul_mat_r_chebp : trop de matrices" << endl ;
112  abort() ;
113  exit (-1) ;
114  }
115 
116 
117  l_dejafait[nb_dejafait] = l ;
118  nr_dejafait[nb_dejafait] = n ;
119 
120  Matrice res(n-1, n-1) ;
121  res.set_etat_qcq() ;
122 
123 
124  double* xdsdx = new double[n] ;
125  double* x2d2sdx2 = new double[n] ;
126  double* d2sdx2 = new double[n] ;
127  double* sxdsdx = new double[n] ;
128  double* sx2 = new double [n] ;
129 
130  for (int i=0 ; i< n-1 ; i++) {
131  for (int j=0 ; j<n ; j++)
132  xdsdx[j] = 0 ;
133  xdsdx[i] = 1 ;
134  xdsdx[i+1] = 1 ;
135  xdsdx_1d (n, &xdsdx, R_CHEBP) ;
136 
137 
138  for (int j=0 ; j<n ; j++)
139  x2d2sdx2[j] = 0 ;
140  x2d2sdx2[i] = 1 ;
141  x2d2sdx2[i+1] = 1 ;
142 
143  d2sdx2_1d(n, &x2d2sdx2, R_CHEBP) ;
144  for (int j=0 ; j<n ; j++)
145  d2sdx2[j] = x2d2sdx2[j] ;
146  multx2_1d(n, &x2d2sdx2, R_CHEBP) ;
147 
148 
149  for (int j=0 ; j<n ; j++)
150  sxdsdx[j] = 0 ;
151  sxdsdx[i] = 1 ;
152  sxdsdx[i+1] = 1 ;
153  sxdsdx_1d(n, &sxdsdx, R_CHEBP) ;
154 
155 
156  for (int j=0 ; j<n ; j++)
157  sx2[j] = 0 ;
158  sx2[i] = 1 ;
159  sx2[i+1] = 1 ;
160  sx2_1d(n, &sx2, R_CHEBP) ;
161 
162  for (int j=0 ; j<n-1 ; j++)
163  res.set(j, i) = (d2sdx2[j] + 2*sxdsdx[j] - l*(l+1)*sx2[j])
164  - (x2d2sdx2[j]+2*xdsdx[j]) ;
165  res.set(i, i) += l*(l+1) ;
166  if (i < n-2)
167  res.set(i+1, i) += l*(l+1) ;
168  }
169  delete [] d2sdx2 ;
170  delete [] x2d2sdx2 ;
171  delete [] sxdsdx ;
172  delete [] xdsdx ;
173  delete [] sx2 ;
174 
175  tab[nb_dejafait] = new Matrice(res) ;
176  nb_dejafait ++ ;
177  return res ;
178  }
179  else
180  return *tab[indice] ;
181 }
182 
183 
184 
185  //------------------------
186  //-- CAS R_CHEBI ----
187  //------------------------
188 
189 
190 Matrice _lap_cpt_mat_r_chebi (int n, int l) {
191 
192  const int nmax = 200 ;// Nombre de Matrices stockees
193  static Matrice* tab[nmax] ; // les matrices calculees
194  static int nb_dejafait = 0 ; // nbre de matrices calculees
195  static int l_dejafait[nmax] ;
196  static int nr_dejafait[nmax] ;
197 
198  int indice = -1 ;
199 
200  // On determine si la matrice a deja ete calculee :
201  for (int conte=0 ; conte<nb_dejafait ; conte ++)
202  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
203  indice = conte ;
204 
205  // Calcul a faire :
206  if (indice == -1) {
207  if (nb_dejafait >= nmax) {
208  cout << "_laplacien_nul_mat_r_chebp : trop de matrices" << endl ;
209  abort() ;
210  exit (-1) ;
211  }
212 
213 
214  l_dejafait[nb_dejafait] = l ;
215  nr_dejafait[nb_dejafait] = n ;
216 
217  // Non degenere si l = 1
218  int taille = (l == 1) ? n : n-1 ;
219  Matrice res(taille, taille) ;
220  res.set_etat_qcq() ;
221 
222 
223  double* xdsdx = new double[n] ;
224  double* x2d2sdx2 = new double[n] ;
225  double* d2sdx2 = new double[n] ;
226  double* sxdsdx = new double[n] ;
227  double* sx2 = new double [n] ;
228 
229  int f_un, f_deux ;
230 
231  for (int i=0 ; i<taille ; i++) {
232 
233  // Gelerkin ou Chebyshev ????
234  if (taille == n) {
235  f_un = 1 ;
236  f_deux = 0 ;
237  }
238  else {
239  f_un = 2*i+3 ;
240  f_deux = 2*i+1 ;
241  }
242 
243  for (int j=0 ; j<n ; j++)
244  xdsdx[j] = 0 ;
245  xdsdx[i] = f_un ;
246  if (i+1 < n)
247  xdsdx[i+1] = f_deux ;
248  xdsdx_1d (n, &xdsdx, R_CHEBI) ;
249 
250 
251  for (int j=0 ; j<n ; j++)
252  x2d2sdx2[j] = 0 ;
253  x2d2sdx2[i] = f_un ;
254  if (i+1 < n)
255  x2d2sdx2[i+1] = f_deux ;
256 
257  d2sdx2_1d(n, &x2d2sdx2, R_CHEBI) ;
258  for (int j=0 ; j<n ; j++)
259  d2sdx2[j] = x2d2sdx2[j] ;
260  multx2_1d(n, &x2d2sdx2, R_CHEBI) ;
261 
262 
263  for (int j=0 ; j<n ; j++)
264  sxdsdx[j] = 0 ;
265  sxdsdx[i] = f_un ;
266  if (i+1 < n)
267  sxdsdx[i+1] = f_deux ;
268  sxdsdx_1d(n, &sxdsdx, R_CHEBI) ;
269 
270 
271  for (int j=0 ; j<n ; j++)
272  sx2[j] = 0 ;
273  sx2[i] = f_un ;
274  if (i+1 < n)
275  sx2[i+1] = f_deux ;
276  sx2_1d(n, &sx2, R_CHEBI) ;
277 
278  for (int j=0 ; j<taille ; j++)
279  res.set(j, i) = (d2sdx2[j] + 2*sxdsdx[j] - l*(l+1)*sx2[j])
280  - (x2d2sdx2[j]+2*xdsdx[j]) ;
281  res.set(i, i) += l*(l+1)*f_un ;
282  if (i < taille-1)
283  res.set(i+1, i) += l*(l+1)*f_deux ;
284  }
285  delete [] d2sdx2 ;
286  delete [] x2d2sdx2 ;
287  delete [] sxdsdx ;
288  delete [] xdsdx ;
289  delete [] sx2 ;
290 
291  tab[nb_dejafait] = new Matrice(res) ;
292  nb_dejafait ++ ;
293  return res ;
294  }
295  else
296  return *tab[indice] ;
297 
298 }
299 
300  //--------------------------
301  //- La routine a appeler ---
302  //----------------------------
303 Matrice lap_cpt_mat(int n, int l, int base_r)
304 {
305 
306  // Routines de derivation
307  static Matrice (*lap_cpt_mat[MAX_BASE])(int, int) ;
308  static int nap = 0 ;
309 
310  // Premier appel
311  if (nap==0) {
312  nap = 1 ;
313  for (int i=0 ; i<MAX_BASE ; i++) {
314  lap_cpt_mat[i] = _lap_cpt_mat_pas_prevu ;
315  }
316  // Les routines existantes
317  lap_cpt_mat[R_CHEBP >> TRA_R] = _lap_cpt_mat_r_chebp ;
318  lap_cpt_mat[R_CHEBI >> TRA_R] = _lap_cpt_mat_r_chebi ;
319  }
320 
321  Matrice res(lap_cpt_mat[base_r](n, l)) ;
322  return res ;
323 }
324 
325 }
Lorene prototypes.
Definition: app_hor.h:64
#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 MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144