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