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
23char 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
76namespace Lorene {
77Matrice _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
93Matrice _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
190Matrice _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 //----------------------------
303Matrice 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}
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64