LORENE
valeur_val_propre_1d.C
1/*
2 * Copyright (c) 2004 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 valeur_val_propre_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_val_propre_1d.C,v 1.3 2014/10/13 08:53:51 j_novak Exp $" ;
24
25
26
27/*
28 * $Id: valeur_val_propre_1d.C,v 1.3 2014/10/13 08:53:51 j_novak Exp $
29 * $Log: valeur_val_propre_1d.C,v $
30 * Revision 1.3 2014/10/13 08:53:51 j_novak
31 * Lorene classes and functions now belong to the namespace Lorene.
32 *
33 * Revision 1.2 2014/10/06 15:13:24 j_novak
34 * Modified #include directives to use c++ syntax.
35 *
36 * Revision 1.1 2004/08/24 09:14:52 p_grandclement
37 * Addition of some new operators, like Poisson in 2d... It now requieres the
38 * GSL library to work.
39 *
40 * Also, the way a variable change is stored by a Param_elliptic is changed and
41 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
42 * will requiere some modification. (It should concern only the ones about monopoles)
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_val_propre_1d.C,v 1.3 2014/10/13 08:53:51 j_novak Exp $
46 *
47 */
48
49// headers du C
50#include <cassert>
51#include <cstdlib>
52
53// headers Lorene
54
55#include "type_parite.h"
56#include "valeur.h"
57#include "matrice.h"
58#include "proto.h"
59
60namespace Lorene {
61
62//****************************************************************
63// CAS PAIR
64//****************************************************************
65
66void rotate_propre_pair (Valeur& so, bool sens) {
67
68 so.coef() ;
69 so.set_etat_cf_qcq() ;
70
71 static int nt_courant = 0 ;
72 static Matrice* passage = 0x0 ;
73 static Matrice* inv = 0x0 ;
74
75 int nt = so.mg->get_nt(0) ;
76
77
78 if (nt != nt_courant) {
79 // On doit calculer les matrices... Pas de la tarte...
80 if (passage != 0x0)
81 delete passage ;
82 if (inv != 0x0)
83 delete inv ;
84
85 nt_courant = nt ;
86
87 Matrice ope (nt, nt) ;
88 ope.set_etat_qcq() ;
89 for (int i=0 ; i<nt ; i++)
90 for (int j=0 ; j<nt ; j++)
91 ope.set(i,j) = 0 ;
92
93 double c_courant = 0 ;
94 for (int j=0 ; j<nt ; j++) {
95 ope.set(0, j) = 2*j ;
96 for (int i=1 ; i<j ; i++)
97 ope.set(i,j) = 4*j ;
98 ope.set(j,j) = c_courant ;
99 c_courant -= 8*j + 2 ;
100 }
101
102 passage = new Matrice (ope.vect_propre()) ;
103 passage->set_band(nt-1,0) ;
104 passage->set_lu() ;
105 // Un peu nul pour calculer l'inverse mais bon...
106 inv = new Matrice (nt, nt) ;
107 inv->set_etat_qcq() ;
108
109 Tbl auxi (nt) ;
110 auxi.set_etat_qcq() ;
111
112 for (int i=0 ; i<nt ; i++) {
113 for (int j=0 ; j<nt ; j++)
114 auxi.set(j) = 0 ;
115 auxi.set(i) = 1 ;
116 Tbl sortie (passage->inverse(auxi)) ;
117 for (int j=0 ; j<nt ; j++)
118 inv->set(j,i) = sortie(j) ;
119 }
120 }
121
122 // Fin du calcul des matrices...
123 for (int l=0 ; l<so.mg->get_nzone() ; l++)
124 if (so.c_cf->t[l]->get_etat() != ETATZERO)
125 for (int k=0 ; k<so.mg->get_np(l) ; k++)
126 for (int i=0 ; i<so.mg->get_nr(l) ; i++) {
127
128 Tbl coefs(nt) ;
129 coefs.set_etat_qcq() ;
130 for (int j=0 ; j<nt ; j++)
131 coefs.set(j) = (*so.c_cf)(l,k,j,i) ;
132 Tbl prod (nt) ;
133 prod.set_etat_qcq() ;
134 for (int j=0 ; j<nt ; j++)
135 prod.set(j) = 0 ;
136
137 if (sens) {
138 //calcul direct
139 for (int j=0 ; j<nt ; j++)
140 for (int jb=0 ; jb<nt ; jb++)
141 prod.set(j) += (*inv)(j, jb) * coefs(jb) ;
142 }
143 else {
144 //calcul inverse
145 for (int j=0 ; j<nt ; j++)
146 for (int jb=0 ; jb<nt ; jb++)
147 prod.set(j) += (*passage)(j, jb) * coefs(jb) ;
148 }
149
150 for (int j=0 ; j<nt ; j++)
151 so.c_cf->set(l,k,j,i) = prod(j) ;
152
153 }
154
155 // On met les bonnes bases :
156 int base_tet = so.base.get_base_t(0) ;
157 Base_val new_base (so.base) ;
158
159 if (sens) {
160 switch (base_tet) {
161 case T_COS_P:
162 new_base.set_base_t(T_CL_COS_P) ;
163 break ;
164 case T_SIN_P:
165 new_base.set_base_t(T_CL_SIN_P) ;
166 break ;
167 default:
168 cout << "Problem in rotate_propre_pair" << endl ;
169 abort() ;
170 break ;
171 }
172 }
173 else {
174 switch (base_tet) {
175 case T_CL_COS_P:
176 new_base.set_base_t(T_COS_P) ;
177 break ;
178 case T_CL_SIN_P:
179 new_base.set_base_t(T_SIN_P) ;
180 break ;
181 default:
182 cout << "Problem in rotate_propre_pair" << endl ;
183 abort() ;
184 break ;
185 }
186 }
187
188 so.c_cf->base = new_base ;
189 so.base = new_base ;
190}
191
192//****************************************************************
193// CAS IMPAIR
194//****************************************************************
195
196void rotate_propre_impair (Valeur& so, bool sens) {
197 so.coef() ;
198 so.set_etat_cf_qcq() ;
199
200 static int nt_courant = 0 ;
201 static Matrice* passage = 0x0 ;
202 static Matrice* inv = 0x0 ;
203
204 int nt = so.mg->get_nt(0) ;
205
206
207 if (nt != nt_courant) {
208 // On doit calculer les matrices... Pas de la tarte...
209 if (passage != 0x0)
210 delete passage ;
211 if (inv != 0x0)
212 delete inv ;
213
214 nt_courant = nt ;
215
216 Matrice ope (nt, nt) ;
217 ope.set_etat_qcq() ;
218 for (int i=0 ; i<nt ; i++)
219 for (int j=0 ; j<nt ; j++)
220 ope.set(i,j) = 0 ;
221
222 double c_courant = 0 ;
223 for (int j=0 ; j<nt ; j++) {
224 for (int i=0 ; i<j ; i++)
225 ope.set(i,j) = 2+4*j ;
226 ope.set(j,j) = c_courant ;
227 c_courant -= 8*j + 6 ;
228 }
229
230 passage = new Matrice (ope.vect_propre()) ;
231 passage->set_band(nt-1,0) ;
232 passage->set_lu() ;
233 // Un peu nul pour calculer l'inverse mais bon...
234 inv = new Matrice (nt, nt) ;
235 inv->set_etat_qcq() ;
236
237 Tbl auxi (nt) ;
238 auxi.set_etat_qcq() ;
239
240 for (int i=0 ; i<nt ; i++) {
241 for (int j=0 ; j<nt ; j++)
242 auxi.set(j) = 0 ;
243 auxi.set(i) = 1 ;
244 Tbl sortie (passage->inverse(auxi)) ;
245 for (int j=0 ; j<nt ; j++)
246 inv->set(j,i) = sortie(j) ;
247 }
248 }
249
250 // Fin du calcul des matrices...
251
252 for (int l=0 ; l<so.mg->get_nzone() ; l++)
253 if (so.c_cf->t[l]->get_etat() != ETATZERO)
254 for (int k=0 ; k<so.mg->get_np(l) ; k++)
255 for (int i=0 ; i<so.mg->get_nr(l) ; i++) {
256
257 Tbl coefs(nt) ;
258 coefs.set_etat_qcq() ;
259 for (int j=0 ; j<nt ; j++)
260 coefs.set(j) = (*so.c_cf)(l,k,j,i) ;
261 Tbl prod (nt) ;
262 prod.set_etat_qcq() ;
263 for (int j=0 ; j<nt ; j++)
264 prod.set(j) = 0 ;
265
266 if (sens) {
267 //calcul direct
268 for (int j=0 ; j<nt ; j++)
269 for (int jb=0 ; jb<nt ; jb++)
270 prod.set(j) += (*inv)(j, jb) * coefs(jb) ;
271 }
272 else {
273 //calcul inverse
274 for (int j=0 ; j<nt ; j++)
275 for (int jb=0 ; jb<nt ; jb++)
276 prod.set(j) += (*passage)(j, jb) * coefs(jb) ;
277 }
278
279 for (int j=0 ; j<nt ; j++)
280 so.c_cf->set(l,k,j,i) = prod(j) ;
281
282 }
283
284 // On met les bonnes bases :
285 int base_tet = so.base.get_base_t(0) ;
286 Base_val new_base (so.base) ;
287
288 if (sens) {
289 switch (base_tet) {
290 case T_COS_I:
291 new_base.set_base_t(T_CL_COS_I) ;
292 break ;
293 case T_SIN_I:
294 new_base.set_base_t(T_CL_SIN_I) ;
295 break ;
296 default:
297 cout << "Problem in rotate_propre_impair" << endl ;
298 abort() ;
299 break ;
300 }
301 }
302 else {
303 switch (base_tet) {
304 case T_CL_COS_I:
305 new_base.set_base_t(T_COS_I) ;
306 break ;
307 case T_CL_SIN_I:
308 new_base.set_base_t(T_SIN_I) ;
309 break ;
310 default:
311 cout << "Problem in rotate_propre_impair" << endl ;
312 abort() ;
313 break ;
314 }
315 }
316
317 so.c_cf->base = new_base ;
318 so.base = new_base ;
319}
320
321
323
324 // On recupere la base en tet :
325 int nz = mg->get_nzone() ;
326 int base_tet = base.get_base_t(0) ;
327 // On verifie que c'est bien la meme partout...
328 for (int i=1 ; i<nz ; i++)
330
331 switch (base_tet) {
332 case T_COS_P:
333 rotate_propre_pair (*this, true) ;
334 break ;
335 case T_SIN_P:
336 rotate_propre_pair (*this, true) ;
337 break ;
338 case T_COS_I:
339 rotate_propre_impair (*this, true) ;
340 break ;
341 case T_SIN_I:
342 rotate_propre_impair (*this, true) ;
343 break ;
344 case T_CL_COS_P:
345 break ;
346 case T_CL_SIN_P:
347 break ;
348 case T_CL_COS_I:
349 break ;
350 case T_CL_SIN_I:
351 break ;
352 default:
353 cout << "Unknown basis in Valeur::val_propre_1d" << endl ;
354 abort() ;
355 break ;
356 }
357}
359
360 // On recupere la base en tet :
361 int nz = mg->get_nzone() ;
362 int base_tet = base.get_base_t(0) ;
363 // On verifie que c'est bien la meme partout...
364 for (int i=1 ; i<nz ; i++)
366
367 switch (base_tet) {
368 case T_CL_COS_P:
369 rotate_propre_pair (*this, false) ;
370 break ;
371 case T_CL_SIN_P:
372 rotate_propre_pair (*this, false) ;
373 break ;
374 case T_CL_COS_I:
375 rotate_propre_impair (*this, false) ;
376 break ;
377 case T_CL_SIN_I:
378 rotate_propre_impair (*this, false) ;
379 break ;
380 case T_COS_P:
381 break ;
382 case T_SIN_P:
383 break ;
384 case T_COS_I:
385 break ;
386 case T_SIN_I:
387 break ;
388 default:
389 cout << "Unknown basis in Valeur::val_propre_1d_i" << endl ;
390 abort() ;
391 break ;
392 }
393}
394
395
396}
Bases of the spectral expansions.
Definition base_val.h:322
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition base_val.h:411
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
Basic array class.
Definition tbl.h:161
Values and coefficients of a (real-value) function.
Definition valeur.h:287
void val_propre_1d_i()
Inverse transformation of val_propre_1d.
friend void rotate_propre_impair(Valeur &, bool)
Friend fonction.
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
friend void rotate_propre_pair(Valeur &, bool)
Friend fonction.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
void val_propre_1d()
Set the basis to the eigenvalues of .
#define T_CL_COS_I
CL of odd cosines.
#define T_CL_SIN_I
CL of odd sines.
#define T_SIN_P
dev. sin seulement, harmoniques paires
#define T_CL_SIN_P
CL of even sines.
#define T_COS_P
dev. cos seulement, harmoniques paires
#define T_CL_COS_P
CL of even cosines.
#define T_SIN_I
dev. sin seulement, harmoniques impaires
#define T_COS_I
dev. cos seulement, harmoniques impaires
Lorene prototypes.
Definition app_hor.h:64