LORENE
valeur_coef.C
1 /*
2  * Computation of the spectral coefficients.
3  */
4 
5 /*
6  * Copyright (c) 1999-2001 Eric Gourgoulhon
7  *
8  * This file is part of LORENE.
9  *
10  * LORENE is free software; you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation; either version 2 of the License, or
13  * (at your option) any later version.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 
27 
28 
29 
30 /*
31  * $Id: valeur_coef.C,v 1.19 2016/12/05 16:18:20 j_novak Exp $
32  * $Log: valeur_coef.C,v $
33  * Revision 1.19 2016/12/05 16:18:20 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.18 2014/10/13 08:53:49 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.17 2013/06/07 14:44:34 j_novak
40  * Coefficient computation for even Legendre basis.
41  *
42  * Revision 1.16 2013/06/05 15:06:11 j_novak
43  * Legendre bases are treated as standard bases, when the multi-grid
44  * (Mg3d) is built with BASE_LEG.
45  *
46  * Revision 1.15 2012/01/17 17:51:16 j_penner
47  * *** empty log message ***
48  *
49  * Revision 1.14 2012/01/17 15:07:57 j_penner
50  * using MAX_BASE_2 for the phi coordinate
51  *
52  * Revision 1.13 2009/10/23 12:56:29 j_novak
53  * New base T_LEG_MI
54  *
55  * Revision 1.12 2009/10/13 13:49:58 j_novak
56  * New base T_LEG_MP.
57  *
58  * Revision 1.11 2008/10/07 15:01:58 j_novak
59  * The case nt=1 is now treated separately.
60  *
61  * Revision 1.10 2008/05/24 15:09:02 j_novak
62  * Getting back to previous version, the new one was an error.
63  *
64  * Revision 1.9 2008/05/24 15:05:22 j_novak
65  * New method Scalar::match_tau to match the output of an explicit time-marching scheme with the tau method.
66  *
67  * Revision 1.8 2008/02/18 13:53:51 j_novak
68  * Removal of special indentation instructions.
69  *
70  * Revision 1.7 2007/12/11 15:28:25 jl_cornou
71  * Jacobi(0,2) polynomials partially implemented
72  *
73  * Revision 1.6 2004/11/23 15:17:19 m_forot
74  * Added the bases for the cases without any equatorial symmetry
75  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
76  *
77  * Revision 1.5 2003/09/17 12:30:22 j_novak
78  * New checks for changing to T_LEG* bases.
79  *
80  * Revision 1.4 2003/09/16 13:07:41 j_novak
81  * New files for coefficient trnasformation to/from the T_LEG_II base.
82  *
83  * Revision 1.3 2002/11/12 17:44:35 j_novak
84  * Added transformation function for T_COS basis.
85  *
86  * Revision 1.2 2002/10/16 14:37:15 j_novak
87  * Reorganization of #include instructions of standard C++, in order to
88  * use experimental version 3 of gcc.
89  *
90  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
91  * LORENE
92  *
93  * Revision 2.10 2000/10/04 14:41:26 eric
94  * Ajout des bases T_LEG_IP et T_LEG_PI
95  *
96  * Revision 2.9 2000/09/29 16:09:25 eric
97  * Mise a zero des coefficients k=1 et k=2 dans le cas np=1.
98  *
99  * Revision 2.8 2000/09/07 15:14:30 eric
100  * Ajout de la base P_COSSIN_I
101  *
102  * Revision 2.7 2000/08/16 10:32:53 eric
103  * Suppression de Mtbl::dzpuis.
104  * >> .
105  *
106  * Revision 2.6 1999/11/30 12:42:29 eric
107  * Valeur::base est desormais du type Base_val et non plus Base_val*.
108  *
109  * Revision 2.5 1999/11/24 16:06:13 eric
110  * Ajout du test de l'admissibilite FFT des nombres de degres de liberte.
111  *
112  * Revision 2.4 1999/10/28 07:43:14 eric
113  * Modif commentaires.
114  *
115  * Revision 2.3 1999/10/13 15:51:07 eric
116  * Ajout de la base dans l'appel au constructeur de Mtbl_cf.
117  *
118  * Revision 2.2 1999/06/22 14:21:23 phil
119  * Ajout de dzpuis
120  *
121  * Revision 2.1 1999/03/01 14:55:12 eric
122  * *** empty log message ***
123  *
124  * Revision 2.0 1999/02/22 15:40:37 hyc
125  * *** empty log message ***
126  *
127  *
128  * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_coef.C,v 1.19 2016/12/05 16:18:20 j_novak Exp $
129  *
130  */
131 #include<cmath>
132 
133 // Header Lorene
134 #include "mtbl.h"
135 #include "mtbl_cf.h"
136 #include "valeur.h"
137 #include "proto.h"
138 
139 // Prototypage local
140 namespace Lorene {
141 void pasprevu_r(const int*, const int*, double*, const int*, double*) ;
142 void pasprevu_t(const int*, const int*, double*, const int*, double*) ;
143 void pasprevu_p(const int* ,const int* , double* ) ;
144 
145 void base_non_def_r(const int*, const int*, double*, const int*, double*) ;
146 void base_non_def_t(const int*, const int*, double*, const int*, double*) ;
147 void base_non_def_p(const int* ,const int* , double* ) ;
148 
149 bool admissible_fft(int ) ;
150 
151 void Valeur::coef() const {
152 
153  // Variables statiques
154  static void (*coef_r[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
155  static void (*coef_t[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
156  static void (*coef_p[MAX_BASE_2])(const int* ,const int* , double* ) ;
157  static int premier_appel = 1 ;
158 
159  // Premier appel
160  if (premier_appel) {
161  premier_appel = 0 ;
162 
163  for (int i=0; i<MAX_BASE; i++) {
164  coef_r[i] = pasprevu_r ;
165  coef_t[i] = pasprevu_t ;
166  if(i%2 == 0){
167  coef_p[i/2] = pasprevu_p ;
168  }
169  }
170 
171  coef_r[NONDEF] = base_non_def_r ;
172  coef_r[R_CHEB >> TRA_R] = cfrcheb ;
173  coef_r[R_CHEBU >> TRA_R] = cfrcheb ;
174  coef_r[R_CHEBP >> TRA_R] = cfrchebp ;
175  coef_r[R_CHEBI >> TRA_R] = cfrchebi ;
176  coef_r[R_CHEBPIM_P >> TRA_R] = cfrchebpimp ;
177  coef_r[R_CHEBPIM_I >> TRA_R] = cfrchebpimi ;
178  coef_r[R_CHEBPI_P >> TRA_R] = cfrchebpip ;
179  coef_r[R_CHEBPI_I >> TRA_R] = cfrchebpii ;
180  coef_r[R_LEG >> TRA_R] = cfrleg ;
181  coef_r[R_LEGP >> TRA_R] = cfrlegp ;
182  coef_r[R_LEGI >> TRA_R] = cfrlegi ;
183  coef_r[R_JACO02 >> TRA_R] = cfrjaco02 ;
184 
185  coef_t[NONDEF] = base_non_def_t ;
186  coef_t[T_COS >> TRA_T] = cftcos ;
187  coef_t[T_SIN >> TRA_T] = cftsin ;
188  coef_t[T_COS_P >> TRA_T] = cftcosp ;
189  coef_t[T_COS_I >> TRA_T] = cftcosi ;
190  coef_t[T_SIN_P >> TRA_T] = cftsinp ;
191  coef_t[T_SIN_I >> TRA_T] = cftsini ;
192  coef_t[T_COSSIN_CP >> TRA_T] = cftcossincp ;
193  coef_t[T_COSSIN_SI >> TRA_T] = cftcossinsi ;
194  coef_t[T_COSSIN_SP >> TRA_T] = cftcossinsp ;
195  coef_t[T_COSSIN_CI >> TRA_T] = cftcossinci ;
196  coef_t[T_COSSIN_S >> TRA_T] = cftcossins ;
197  coef_t[T_COSSIN_C >> TRA_T] = cftcossinc ;
198  coef_t[T_LEG_P >> TRA_T] = cftlegp ;
199  coef_t[T_LEG_PP >> TRA_T] = cftlegpp ;
200  coef_t[T_LEG_I >> TRA_T] = cftlegi ;
201  coef_t[T_LEG_IP >> TRA_T] = cftlegip ;
202  coef_t[T_LEG_PI >> TRA_T] = cftlegpi ;
203  coef_t[T_LEG_II >> TRA_T] = cftlegii ;
204  coef_t[T_LEG_MP >> TRA_T] = cftlegmp ;
205  coef_t[T_LEG_MI >> TRA_T] = cftlegmi ;
206  coef_t[T_LEG >> TRA_T] = cftleg ;
207 
208  coef_p[NONDEF] = base_non_def_p ;
209  coef_p[P_COSSIN >> TRA_P] = cfpcossin ;
210  coef_p[P_COSSIN_P >> TRA_P] = cfpcossin ;
211  coef_p[P_COSSIN_I >> TRA_P] = cfpcossini ;
212 
213  } // fin des operation de premier appel
214 
215 
216  //-------------------//
217  // DEBUT DU CALCUL //
218  //-------------------//
219 
220  // Tout null ?
221  if (etat == ETATZERO) {
222  return ;
223  }
224 
225  // Protection
226  assert(etat != ETATNONDEF) ;
227 
228  // Peut-etre rien a faire
229  if (c_cf != 0x0) {
230  return ;
231  }
232 
233  // Il faut bosser
234  assert(c != 0x0) ; // ..si on peut
235  c_cf = new Mtbl_cf(mg, base) ;
236  c_cf->set_etat_qcq() ;
237 
238  // Boucles sur les zones
239  int nz = mg->get_nzone() ;
240  for (int l=0; l<nz; l++) {
241 
242  // Initialisation des valeurs de this->c_cf avec celle de this->c :
243  const Tbl* f = (c->t)[l] ;
244  Tbl* cf = (c_cf->t)[l] ;
245 
246  if (f->get_etat() == ETATZERO) {
247  cf->set_etat_zero() ;
248  continue ; // on ne fait rien si le tbl = 0
249  }
250 
251  cf->set_etat_qcq() ;
252 
253  int np = f->get_dim(2) ;
254  int nt = f->get_dim(1) ;
255  int nr = f->get_dim(0) ;
256 
257  int np_c = cf->get_dim(2) ;
258  int nt_c = cf->get_dim(1) ;
259  int nr_c = cf->get_dim(0) ;
260 
261  // Attention a ce qui suit... (deg et dim)
262  int deg[3] ;
263  deg[0] = np ;
264  deg[1] = nt ;
265  deg[2] = nr ;
266 
267  int dim[3] ;
268  dim[0] = np_c ;
269  dim[1] = nt_c ;
270  dim[2] = nr_c ;
271 
272  int nrnt = nr * nt ;
273  int nrnt_c = nr_c * nt_c ;
274 
275  for (int i=0; i<np ; i++) {
276  for (int j=0; j<nt ; j++) {
277  for (int k=0; k<nr ; k++) {
278  int index = nrnt * i + nr * j + k ;
279  int index_c = nrnt_c * i + nr_c * j + k ;
280  (cf->t)[index_c] = (f->t)[index] ;
281  }
282  }
283  }
284 
285  // On recupere les bases en r, theta et phi :
286  int base_r = ( base.b[l] & MSQ_R ) >> TRA_R ;
287  int base_t = ( base.b[l] & MSQ_T ) >> TRA_T ;
288  int base_p = ( base.b[l] & MSQ_P ) >> TRA_P ;
289  int vbase_t = base.b[l] & MSQ_T ;
290  int vbase_p = base.b[l] & MSQ_P ;
291 
292  assert(base_r < MAX_BASE) ;
293  assert(base_t < MAX_BASE) ;
294  assert(base_p < MAX_BASE_2) ;
295 
296  // Transformation en phi
297  // ---------------------
298  if ( np > 1 ) {
299  assert( admissible_fft(np) ) ;
300 
301  coef_p[base_p]( deg, dim, (cf->t) ) ;
302  }
303  else{ // Cas np=1 : mise a zero des coefficients k=1 et k=2 :
304  for (int i=nrnt; i<3*nrnt; i++) {
305  cf->t[i] = 0 ;
306  }
307  }
308 
309  // Transformation en theta:
310  // ------------------------
311 
312  if ( nt > 1 ) {
313  assert( admissible_fft(nt-1) ) ;
314  bool pair = ( (vbase_t == T_LEG_PP) || (vbase_t == T_LEG_IP)
315  || (vbase_t == T_LEG_MP) ) ;
316  bool impair = ( (vbase_t == T_LEG_PI) || (vbase_t == T_LEG_II)
317  || (vbase_t == T_LEG_MI) ) ;
318 
319  if ((pair && (vbase_p == P_COSSIN_I)) ||
320  (impair && (vbase_p == P_COSSIN_P)) )
321  pasprevu_t(deg, dim, (cf->t), dim, (cf->t) ) ;
322  else
323  coef_t[base_t](deg, dim, (cf->t), dim, (cf->t)) ;
324  }
325  else {
326  if ((vbase_t == T_LEG_PP) || (vbase_t == T_LEG_PI) ||
327  (vbase_t == T_LEG_IP) || (vbase_t == T_LEG_II) ||
328  (vbase_t == T_LEG_P) || (vbase_t == T_LEG_I) ||
329  (vbase_t == T_LEG) || (vbase_t == T_LEG_MP)
330  || (vbase_t == T_LEG_MI) ) {
331 
332  *c_cf->t[l] *=sqrt(2.) ;
333  }
334  }
335 
336 
337  // Transformation en r:
338  // --------------------
339  if ( nr > 1 ) {
340  assert( admissible_fft(nr-1) || (mg->get_colloc_r(l) != BASE_CHEB) ) ;
341  coef_r[base_r](deg, dim, (cf->t), dim, (cf->t)) ;
342  }
343 
344  } // fin de la boucle sur les differentes zones
345 
346 }
347 
348  //------------------------//
349  // Les machins pas prevus //
350  //------------------------//
351 
352 void pasprevu_r(const int*, const int*, double*, const int*, double*) {
353  cout << "Valeur::coef: the required expansion basis in r " << endl ;
354  cout << " is not implemented !" << endl ;
355  abort() ;
356 }
357 
358 void pasprevu_t(const int*, const int*, double*, const int*, double*) {
359  cout << "Valeur::coef: the required expansion basis in theta " << endl ;
360  cout << " is not implemented !" << endl ;
361  abort() ;
362 }
363 
364 void pasprevu_p(const int*, const int*, double*) {
365  cout << "Valeur::coef: the required expansion basis in phi " << endl ;
366  cout << " is not implemented !" << endl ;
367  abort() ;
368 }
369 
370 void base_non_def_r(const int*, const int*, double*, const int*, double*) {
371  cout << "Valeur::coef: the expansion basis in r is undefined !" << endl ;
372  abort() ;
373 }
374 
375 void base_non_def_t(const int*, const int*, double*, const int*, double*) {
376  cout << "Valeur::coef: the expansion basis in theta is undefined !" << endl ;
377  abort() ;
378 }
379 
380 void base_non_def_p(const int*, const int*, double*) {
381  cout << "Valeur::coef: the expansion basis in phi is undefined !" << endl ;
382  abort() ;
383 }
384 
385 }
#define T_LEG
fct. de Legendre associees
Definition: type_parite.h:236
#define MAX_BASE_2
Smaller maximum bases used for phi (and higher dimensions for now)
Definition: type_parite.h:146
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
#define T_LEG_MP
fct. de Legendre associees avec m pair
Definition: type_parite.h:238
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
#define T_LEG_PI
fct. de Legendre associees paires avec m impair
Definition: type_parite.h:224
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
Definition: type_parite.h:162
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
#define T_LEG_MI
fct. de Legendre associees avec m impair
Definition: type_parite.h:240
Lorene prototypes.
Definition: app_hor.h:67
#define TRA_T
Translation en Theta, used for a bitwise shift (in hex)
Definition: type_parite.h:160
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
int get_colloc_r(int l) const
Returns the type of collocation points used in domain no.
Definition: grilles.h:528
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition: valeur.h:302
#define T_LEG_I
fct. de Legendre associees impaires
Definition: type_parite.h:220
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
double * t
The array of double.
Definition: tbl.h:176
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:334
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
#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 set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl_cf.C:303
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
#define T_LEG_IP
fct. de Legendre associees impaires avec m pair
Definition: type_parite.h:222
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:350
#define T_LEG_P
fct. de Legendre associees paires
Definition: type_parite.h:216
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
#define NONDEF
base inconnue
Definition: type_parite.h:150
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
bool admissible_fft(int)
Checks whether or not a given number of degrees of freedom is compatible with the FFT transform...
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: valeur.h:305
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
Basic array class.
Definition: tbl.h:164
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
#define T_LEG_II
fct. de Legendre associees impaires avec m impair
Definition: type_parite.h:226
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
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
#define T_LEG_PP
fct. de Legendre associees paires avec m pair
Definition: type_parite.h:218
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166