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