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