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 
29 char 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 
130 namespace Lorene {
131 void c_est_pas_fait(char * ) ;
132 
133 void ipasprevu_r(const int*, const int*, double*, const int*, double*) ;
134 void ipasprevu_t(const int*, const int*, double*, const int*, double*) ;
135 void ipasprevu_p(const int* , const int* , const int* , double* , double* ) ;
136 void ibase_non_def_r(const int*, const int*, double*, const int*, double*) ;
137 void ibase_non_def_t(const int*, const int*, double*, const int*, double*) ;
138 void ibase_non_def_p(const int* , const int* , const int* , double* , double* ) ;
139 
140 void 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) ;
274  assert(base_t < MAX_BASE) ;
275  assert(base_p < MAX_BASE_2) ;
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
327  invcf_t[base_t]( deg, dimc, trav, dimc, trav ) ;
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 
341 void 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 
347 void 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 
353 void 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 
359 void 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 
364 void 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 
370 void 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
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
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
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
double * t
The array of double.
Definition: tbl.h:173
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
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)
Definition: type_parite.h:146
#define TRA_T
Translation en Theta, used for a bitwise shift (in hex)
Definition: type_parite.h:160
#define T_LEG_MP
fct. de Legendre associees avec m pair
Definition: type_parite.h:238
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define T_LEG_PI
fct. de Legendre associees paires avec m impair
Definition: type_parite.h:224
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define T_LEG
fct. de Legendre associees
Definition: type_parite.h:236
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
#define T_LEG_P
fct. de Legendre associees paires
Definition: type_parite.h:216
#define T_LEG_IP
fct. de Legendre associees impaires avec m pair
Definition: type_parite.h:222
#define NONDEF
base inconnue
Definition: type_parite.h:150
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define T_LEG_MI
fct. de Legendre associees avec m impair
Definition: type_parite.h:240
#define MSQ_T
Extraction de l'info sur Theta.
Definition: type_parite.h:154
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
#define T_LEG_II
fct. de Legendre associees impaires avec m impair
Definition: type_parite.h:226
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
#define T_LEG_I
fct. de Legendre associees impaires
Definition: type_parite.h:220
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
#define T_LEG_PP
fct. de Legendre associees paires avec m pair
Definition: type_parite.h:218
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
Definition: type_parite.h:162
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Lorene prototypes.
Definition: app_hor.h:64