LORENE
valeur_coef_i.C
1/*
2 * Computations of the values at the collocation points from the spectral
3 * coefficients
4 *
5 */
6
7/*
8 * Copyright (c) 1999-2003 Eric Gourgoulhon
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28
29char valeur_coef_i_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_coef_i.C,v 1.17 2014/10/13 08:53:49 j_novak Exp $" ;
30
31/*
32 * $Id: valeur_coef_i.C,v 1.17 2014/10/13 08:53:49 j_novak Exp $
33 * $Log: valeur_coef_i.C,v $
34 * Revision 1.17 2014/10/13 08:53:49 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.16 2014/10/06 15:13:22 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.15 2013/06/06 15:31:33 j_novak
41 * Functions to compute Legendre coefficients (not fully tested yet).
42 *
43 * Revision 1.14 2012/01/17 15:08:02 j_penner
44 * using MAX_BASE_2 for the phi coordinate
45 *
46 * Revision 1.13 2009/10/23 12:56:29 j_novak
47 * New base T_LEG_MI
48 *
49 * Revision 1.12 2009/10/13 13:49:58 j_novak
50 * New base T_LEG_MP.
51 *
52 * Revision 1.11 2009/10/08 16:23:14 j_novak
53 * Addition of new bases T_COS and T_SIN.
54 *
55 * Revision 1.10 2008/10/07 15:01:58 j_novak
56 * The case nt=1 is now treated separately.
57 *
58 * Revision 1.9 2007/12/11 15:28:25 jl_cornou
59 * Jacobi(0,2) polynomials partially implemented
60 *
61 * Revision 1.8 2004/11/23 15:17:19 m_forot
62 * Added the bases for the cases without any equatorial symmetry
63 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
64 *
65 * Revision 1.7 2004/08/24 09:14:52 p_grandclement
66 * Addition of some new operators, like Poisson in 2d... It now requieres the
67 * GSL library to work.
68 *
69 * Also, the way a variable change is stored by a Param_elliptic is changed and
70 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
71 * will requiere some modification. (It should concern only the ones about monopoles)
72 *
73 * Revision 1.6 2003/10/13 20:51:25 e_gourgoulhon
74 * Replaced malloc by new
75 *
76 * Revision 1.5 2003/09/17 12:30:22 j_novak
77 * New checks for changing to T_LEG* bases.
78 *
79 * Revision 1.4 2003/09/16 13:07:41 j_novak
80 * New files for coefficient trnasformation to/from the T_LEG_II base.
81 *
82 * Revision 1.3 2002/11/12 17:44:35 j_novak
83 * Added transformation function for T_COS basis.
84 *
85 * Revision 1.2 2002/10/16 14:37:15 j_novak
86 * Reorganization of #include instructions of standard C++, in order to
87 * use experimental version 3 of gcc.
88 *
89 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
90 * LORENE
91 *
92 * Revision 2.9 2000/10/04 14:41:40 eric
93 * Ajout des bases T_LEG_IP et T_LEG_PI
94 *
95 * Revision 2.8 2000/09/07 15:14:44 eric
96 * Ajout de la base P_COSSIN_I
97 *
98 * Revision 2.7 2000/08/16 10:33:04 eric
99 * Suppression de Mtbl::dzpuis.
100 *
101 * Revision 2.6 1999/12/16 16:41:48 phil
102 * *** empty log message ***
103 *
104 * Revision 2.5 1999/11/30 12:43:03 eric
105 * Valeur::base est desormais du type Base_val et non plus Base_val*.
106 *
107 * Revision 2.4 1999/10/28 07:44:20 eric
108 * Modif commentaires.
109 *
110 * Revision 2.3 1999/10/13 15:51:33 eric
111 * Anglisation.
112 *
113 * Revision 2.2 1999/06/22 15:10:02 phil
114 * ajout de dzpuis
115 *
116 * Revision 2.1 1999/02/22 15:40:25 hyc
117 * *** empty log message ***
118 *
119 *
120 * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_coef_i.C,v 1.17 2014/10/13 08:53:49 j_novak Exp $
121 *
122 */
123
124#include <cmath>
125
126// Header Lorene
127#include "valeur.h"
128#include "proto.h"
129
130namespace Lorene {
131void c_est_pas_fait(char * ) ;
132
133void ipasprevu_r(const int*, const int*, double*, const int*, double*) ;
134void ipasprevu_t(const int*, const int*, double*, const int*, double*) ;
135void ipasprevu_p(const int* , const int* , const int* , double* , double* ) ;
136void ibase_non_def_r(const int*, const int*, double*, const int*, double*) ;
137void ibase_non_def_t(const int*, const int*, double*, const int*, double*) ;
138void ibase_non_def_p(const int* , const int* , const int* , double* , double* ) ;
139
140void Valeur::coef_i() const {
141
142 // Variables statiques
143 static void (*invcf_r[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
144 static void (*invcf_t[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
145 static void (*invcf_p[MAX_BASE_2])(const int* , const int* , const int*, double* , double* ) ;
146 static int premier_appel = 1 ;
147
148 // Premier appel
149 if (premier_appel) {
150 premier_appel = 0 ;
151
152 for (int i=0; i<MAX_BASE; i++) {
153 invcf_r[i] = ipasprevu_r ;
154 invcf_t[i] = ipasprevu_t ;
155 if(i%2==0){
156 invcf_p[i/2] = ipasprevu_p ; // saves a loop
157 }
158 }
159
160 invcf_r[NONDEF] = ibase_non_def_r ;
161 invcf_r[R_CHEB >> TRA_R] = circheb ;
162 invcf_r[R_CHEBU >> TRA_R] = circheb ;
163 invcf_r[R_CHEBP >> TRA_R] = circhebp ;
164 invcf_r[R_CHEBI >> TRA_R] = circhebi ;
165 invcf_r[R_CHEBPIM_P >> TRA_R] = circhebpimp ;
166 invcf_r[R_CHEBPIM_I >> TRA_R] = circhebpimi ;
167 invcf_r[R_CHEBPI_P >> TRA_R] = circhebpip ;
168 invcf_r[R_CHEBPI_I >> TRA_R] = circhebpii ;
169 invcf_r[R_LEG >> TRA_R] = cirleg ;
170 invcf_r[R_LEGP >> TRA_R] = cirlegp ;
171 invcf_r[R_LEGI >> TRA_R] = cirlegi ;
172 invcf_r[R_JACO02 >> TRA_R] = cirjaco02 ;
173
174 invcf_t[NONDEF] = ibase_non_def_t ;
175 invcf_t[T_COS >> TRA_T] = citcos ;
176 invcf_t[T_SIN >> TRA_T] = citsin ;
177 invcf_t[T_COS_P >> TRA_T] = citcosp ;
178 invcf_t[T_COS_I >> TRA_T] = citcosi ;
179 invcf_t[T_SIN_P >> TRA_T] = citsinp ;
180 invcf_t[T_SIN_I >> TRA_T] = citsini ;
181 invcf_t[T_COSSIN_CP >> TRA_T] = citcossincp ;
182 invcf_t[T_COSSIN_SI >> TRA_T] = citcossinsi ;
183 invcf_t[T_COSSIN_SP >> TRA_T] = citcossinsp ;
184 invcf_t[T_COSSIN_CI >> TRA_T] = citcossinci ;
185 invcf_t[T_COSSIN_S >> TRA_T] = citcossins ;
186 invcf_t[T_COSSIN_C >> TRA_T] = citcossinc ;
187 invcf_t[T_LEG_P >> TRA_T] = citlegp ;
188 invcf_t[T_LEG_PP >> TRA_T] = citlegpp ;
189 invcf_t[T_LEG_I >> TRA_T] = citlegi ;
190 invcf_t[T_LEG_IP >> TRA_T] = citlegip ;
191 invcf_t[T_LEG_PI >> TRA_T] = citlegpi ;
192 invcf_t[T_LEG_II >> TRA_T] = citlegii ;
193 invcf_t[T_LEG_MP >> TRA_T] = citlegmp ;
194 invcf_t[T_LEG_MI >> TRA_T] = citlegmi ;
195 invcf_t[T_LEG >> TRA_T] = citleg ;
196
197 invcf_p[NONDEF] = ibase_non_def_p ;
198 invcf_p[P_COSSIN >> TRA_P] = cipcossin ;
199 invcf_p[P_COSSIN_P >> TRA_P] = cipcossin ;
200 invcf_p[P_COSSIN_I >> TRA_P] = cipcossini ;
201
202 } // fin des operation de premier appel
203
204 //------------------//
205 // DEBUT DU CALCUL //
206 //------------------//
207
208 // Tout null ?
209 if (etat == ETATZERO) {
210 return ;
211 }
212
213 // Protection
214 assert(etat != ETATNONDEF) ;
215
216 // Peut-etre rien a faire
217 if (c != 0x0) {
218 return ;
219 }
220
221 // Il faut bosser
222 assert(c_cf != 0x0) ; // ..si on peut
223 assert(c_cf->base == base) ; // Consistence des bases
224
225 c = new Mtbl(mg) ;
226 c->set_etat_qcq() ;
227
228 // Boucles sur les zones
229 int nz = mg->get_nzone() ;
230 for (int l=0; l<nz; l++) {
231
232 // Initialisation des valeurs de this->c_cf avec celle de this->c :
233 Tbl* f = (c->t)[l] ;
234 const Tbl* cf = (c_cf->t)[l] ;
235
236 if (cf->get_etat() == ETATZERO) {
237 f->set_etat_zero() ;
238 continue ; // on ne fait rien si le tbl(cf) = 0
239 }
240
241 f->set_etat_qcq() ;
242
243 int np = f->get_dim(2) ;
244 int nt = f->get_dim(1) ;
245 int nr = f->get_dim(0) ;
246
247 int np_c = cf->get_dim(2) ;
248 int nt_c = cf->get_dim(1) ;
249 int nr_c = cf->get_dim(0) ;
250
251 // Attention a ce qui suit... (deg et dim)
252 int deg[3] ;
253 deg[0] = np ;
254 deg[1] = nt ;
255 deg[2] = nr ;
256
257 int dimc[3] ;
258 dimc[0] = np_c ;
259 dimc[1] = nt_c ;
260 dimc[2] = nr_c ;
261
262 // Allocation de l'espace memoire pour le tableau de travail trav
263 int ntot = cf->get_taille() ;
264 double* trav = new double[ntot] ;
265
266 // On recupere les bases en r, theta et phi :
267 int base_r = ( base.b[l] & MSQ_R ) >> TRA_R ;
268 int base_t = ( base.b[l] & MSQ_T ) >> TRA_T ;
269 int base_p = ( base.b[l] & MSQ_P ) >> TRA_P ;
270 int vbase_t = base.b[l] & MSQ_T ;
271 int vbase_p = base.b[l] & MSQ_P ;
272
273 assert(base_r < MAX_BASE) ;
276
277 // Transformation inverse en r:
278 if ( nr == 1 ) {
279 for (int i=0; i<ntot; i++) {
280 trav[i] = cf->t[i] ; // simple recopie cf --> trav
281 }
282 }
283 else {
284 invcf_r[base_r]( deg, dimc, (cf->t), dimc, trav ) ;
285 }
286
287 // Partie angulaire
288 if ( np == 1) {
289 if (nt==1) {
290 for (int i=0 ; i<f->get_taille() ; i++)
291 f->t[i] = trav[i] ;
292 if ((vbase_t == T_LEG_PP) || (vbase_t == T_LEG_PI) ||
293 (vbase_t == T_LEG_IP) || (vbase_t == T_LEG_II) ||
294 (vbase_t == T_LEG_P) || (vbase_t == T_LEG_I) ||
295 (vbase_t == T_LEG_MP) || (vbase_t == T_LEG_MI) ||
296 (vbase_t == T_LEG) ) {
297
298 *f /=sqrt(2.) ;
299 }
300 }
301
302 else {
303 bool pair = ( (vbase_t == T_LEG_PP) || (vbase_t == T_LEG_IP)
304 || (vbase_t == T_LEG_MP) ) ;
305 bool impair = ( (vbase_t == T_LEG_PI) || (vbase_t == T_LEG_II)
306 || (vbase_t == T_LEG_MI)) ;
307
308 if ((pair && (vbase_p == P_COSSIN_I)) ||
309 (impair && (vbase_p == P_COSSIN_P)) )
310 ipasprevu_t(deg, dimc, trav, deg, (f->t) ) ;
311 else
312 invcf_t[base_t]( deg, dimc, trav, deg, (f->t) ) ;
313 }
314 }
315 else {
316 // Cas 3-D
317 // ...... Transformation inverse en theta:
318 bool pair = ( (vbase_t == T_LEG_PP) || (vbase_t == T_LEG_IP)
319 || (vbase_t == T_LEG_MP) ) ;
320 bool impair = ( (vbase_t == T_LEG_PI) || (vbase_t == T_LEG_II)
321 || (vbase_t == T_LEG_MI) ) ;
322
323 if ((pair && (vbase_p == P_COSSIN_I)) ||
324 (impair && (vbase_p == P_COSSIN_P)) )
325 ipasprevu_t(deg, dimc, trav, dimc, trav ) ;
326 else
328 // ...... Transformation inverse en phi:
329 invcf_p[base_p]( deg, dimc, deg, trav, (f->t) ) ;
330 }
331 // Menage
332 delete [] trav ;
333 } // fin de la boucle sur les differentes zones
334
335}
336
337 //------------------------//
338 // Les machins pas prevus //
339 //------------------------//
340
341void ipasprevu_r(const int*, const int*, double*, const int*, double*) {
342 cout << "Valeur::coef_i: the required expansion basis in r " << endl ;
343 cout << " is not implemented !" << endl ;
344 abort() ;
345}
346
347void ipasprevu_t(const int*, const int*, double*, const int*, double* ) {
348 cout << "Valeur::coef_i: the required expansion basis in theta " << endl ;
349 cout << " is not implemented !" << endl ;
350 abort() ;
351}
352
353void ipasprevu_p(const int*, const int*, const int*, double*, double* ) {
354 cout << "Valeur::coef_i: the required expansion basis in phi " << endl ;
355 cout << " is not implemented !" << endl ;
356 abort() ;
357}
358
359void ibase_non_def_r(const int*, const int*, double*, const int*, double*) {
360 cout << "Valeur::coef_i: the expansion basis in r is undefined !" << endl ;
361 abort() ;
362}
363
364void ibase_non_def_t(const int*, const int*, double*, const int*, double*) {
365 cout << "Valeur::coef_i: the expansion basis in theta is undefined !"
366 << endl ;
367 abort() ;
368}
369
370void ibase_non_def_p(const int*, const int*, const int*, double*, double*) {
371 cout << "Valeur::coef_i: the expansion basis in phi is undefined !" << endl ;
372 abort() ;
373}
374
375}
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
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Multi-domain array.
Definition mtbl.h:118
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
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
void coef_i() const
Computes the physical value of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
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.
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Lorene prototypes.
Definition app_hor.h:64