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 char 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 
60 namespace Lorene {
61 
62 //****************************************************************
63 // CAS PAIR
64 //****************************************************************
65 
66 void 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 
196 void 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++)
329  assert (base_tet == base.get_base_t(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++)
365  assert (base_tet == base.get_base_t(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 }
friend void rotate_propre_pair(Valeur &, bool)
Friend fonction.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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:712
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
Lorene prototypes.
Definition: app_hor.h:64
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
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
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
void set_base_t(int base_t)
Sets the expansion basis for functions in all domains.
Definition: base_val.C:195
#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:361
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition: valeur.h:292
#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:305
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:448
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:364
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
void val_propre_1d_i()
Inverse transformation of val_propre_1d.
Bases of the spectral expansions.
Definition: base_val.h:322
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
#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:507
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205