LORENE
mtbl_cf_vp_asymy.C
1 /*
2  * Member functions of the Mtbl_cf class for computing the value of a field
3  * at an arbitrary point, when the field is antisymmetric with respect to the
4  * y=0 plane.
5  *
6  * (see file mtbl_cf.h for the documentation).
7  */
8 
9 /*
10  * Copyright (c) 1999-2001 Eric Gourgoulhon
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 char mtbl_cf_vp_asymy_C[] = "$Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_vp_asymy.C,v 1.3 2014/10/13 08:53:09 j_novak Exp $" ;
31 
32 /*
33  * $Id: mtbl_cf_vp_asymy.C,v 1.3 2014/10/13 08:53:09 j_novak Exp $
34  * $Log: mtbl_cf_vp_asymy.C,v $
35  * Revision 1.3 2014/10/13 08:53:09 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.2 2012/01/17 15:09:14 j_penner
39  * using MAX_BASE_2 for the phi coordinate
40  *
41  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
42  * LORENE
43  *
44  * Revision 2.3 2000/09/08 16:07:26 eric
45  * Ajout de la base P_COSSIN_I
46  *
47  * Revision 2.2 2000/03/06 15:59:03 eric
48  * *** empty log message ***
49  *
50  * Revision 2.1 2000/03/06 15:57:29 eric
51  * *** empty log message ***
52  *
53  * Revision 2.0 2000/03/06 10:26:54 eric
54  * *** empty log message ***
55  *
56  *
57  * $Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_vp_asymy.C,v 1.3 2014/10/13 08:53:09 j_novak Exp $
58  *
59  */
60 
61 
62 // Headers Lorene
63 #include "mtbl_cf.h"
64 #include "proto.h"
65 
66  //-------------------------------------------------------------//
67 namespace Lorene {
68  // version for an arbitrary point in (xi,theta',phi') //
69  //-------------------------------------------------------------//
70 
71 double Mtbl_cf::val_point_asymy(int l, double x, double theta, double phi)
72  const {
73 
74 // Routines de sommation
75 static void (*som_r[MAX_BASE])
76  (double*, const int, const int, const int, const double, double*) ;
77 static void (*som_tet[MAX_BASE])
78  (double*, const int, const int, const double, double*) ;
79 static void (*som_phi[MAX_BASE_2])
80  (double*, const int, const double, double*) ;
81 static int premier_appel = 1 ;
82 
83 // Initialisations au premier appel
84 // --------------------------------
85  if (premier_appel == 1) {
86 
87  premier_appel = 0 ;
88 
89  for (int i=0 ; i<MAX_BASE ; i++) {
90  if(i%2==0){
91  som_phi[i/2] = som_phi_pas_prevu ;
92  }
93  som_tet[i] = som_tet_pas_prevu ;
94  som_r[i] = som_r_pas_prevu ;
95  }
96 
97  som_r[R_CHEB >> TRA_R] = som_r_cheb_asymy ;
98  som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
99  som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
100  som_r[R_CHEBU >> TRA_R] = som_r_chebu_asymy ;
103 
104  som_tet[T_COS >> TRA_T] = som_tet_cos ;
105  som_tet[T_SIN >> TRA_T] = som_tet_sin ;
106  som_tet[T_COS_P >> TRA_T] = som_tet_cos_p ;
107  som_tet[T_SIN_P >> TRA_T] = som_tet_sin_p ;
110 
111  som_phi[P_COSSIN >> TRA_P] = som_phi_cossin_asymy ;
112  som_phi[P_COSSIN_P >> TRA_P] = som_phi_cossin_p ;
113  som_phi[P_COSSIN_I >> TRA_P] = som_phi_cossin_i ;
114 
115  } // fin des operations de premier appel
116 
117 
118  assert (etat != ETATNONDEF) ;
119 
120  double resu ; // valeur de retour
121 
122 // Cas ou tous les coefficients sont nuls :
123  if (etat == ETATZERO ) {
124  resu = 0 ;
125  return resu ;
126  }
127 
128 // Nombre de points en phi, theta et r :
129  int np = mg->get_np(l) ;
130  int nt = mg->get_nt(l) ;
131  int nr = mg->get_nr(l) ;
132 
133 // Bases de developpement :
134  int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
135  int base_t = (base.b[l] & MSQ_T) >> TRA_T ;
136  int base_p = (base.b[l] & MSQ_P) >> TRA_P ;
137 
138 //--------------------------------------
139 // Calcul de la valeur au point demande
140 //--------------------------------------
141 
142 // Pointeur sur le tableau contenant les coefficients:
143 
144  assert(etat == ETATQCQ) ;
145  Tbl* tbcf = t[l] ;
146 
147  if (tbcf->get_etat() == ETATZERO ) {
148  resu = 0 ;
149  return resu ;
150  }
151 
152 
153  assert(tbcf->get_etat() == ETATQCQ) ;
154 
155  double* cf = tbcf->t ;
156 
157  // Tableaux de travail
158  double* trp = new double [np+2] ;
159  double* trtp = new double [(np+2)*nt] ;
160 
161  if (nr == 1) {
162 
163 // Cas particulier nr = 1 (Fonction purement angulaire) :
164 // ----------------------
165 
166  som_tet[base_t](cf, nt, np, theta, trp) ; // sommation sur theta
167  }
168  else {
169 
170 // Cas general
171 // -----------
172 
173  som_r[base_r](cf, nr, nt, np, x, trtp) ; // sommation sur r
174  som_tet[base_t](trtp, nt, np, theta, trp) ; // sommation sur theta
175  }
176 
177 // Sommation sur phi
178 // -----------------
179 
180  if (np == 1) {
181  resu = trp[0] ; // cas axisymetrique
182  }
183  else {
184  som_phi[base_p](trp, np, phi, &resu) ; // sommation sur phi
185  }
186 
187  // Menage
188  delete [] trp ;
189  delete [] trtp ;
190 
191  // Termine
192  return resu ;
193 
194 }
195 
196 
197 
198  //-------------------------------------------------------------//
199  // version for an arbitrary point in xi //
200  // but collocation point in (theta',phi') //
201  //-------------------------------------------------------------//
202 
203 double Mtbl_cf::val_point_jk_asymy(int l, double x, int j0, int k0) const {
204 
205 // Routines de sommation
206 static void (*som_r[MAX_BASE])
207  (double*, const int, const int, const int, const double, double*) ;
208 static int premier_appel = 1 ;
209 
210 // Initialisations au premier appel
211 // --------------------------------
212  if (premier_appel == 1) {
213 
214  premier_appel = 0 ;
215 
216  for (int i=0 ; i<MAX_BASE ; i++) {
217  som_r[i] = som_r_pas_prevu ;
218  }
219 
220  som_r[R_CHEB >> TRA_R] = som_r_cheb_asymy ;
221  som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
222  som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
223  som_r[R_CHEBU >> TRA_R] = som_r_chebu_asymy ;
226 
227  } // fin des operations de premier appel
228 
229  assert (etat != ETATNONDEF) ;
230 
231  double resu ; // valeur de retour
232 
233 // Cas ou tous les coefficients sont nuls :
234  if (etat == ETATZERO ) {
235  resu = 0 ;
236  return resu ;
237  }
238 
239 // Nombre de points en phi, theta et r :
240  int np = mg->get_np(l) ;
241  int nt = mg->get_nt(l) ;
242  int nr = mg->get_nr(l) ;
243 
244 // Bases de developpement :
245  int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
246 
247 //------------------------------------------------------------------------
248 // Valeurs des fonctions de base en phi aux points de collocation en phi
249 // et des fonctions de base en theta aux points de collocation en theta
250 //------------------------------------------------------------------------
251 
252  Tbl tab_phi = base.phi_functions(l, np) ;
253  Tbl tab_theta = base.theta_functions(l, nt) ;
254 
255 
256 //--------------------------------------
257 // Calcul de la valeur au point demande
258 //--------------------------------------
259 
260 // Pointeur sur le tableau contenant les coefficients:
261 
262  assert(etat == ETATQCQ) ;
263  Tbl* tbcf = t[l] ;
264 
265  if (tbcf->get_etat() == ETATZERO ) {
266  resu = 0 ;
267  return resu ;
268  }
269 
270 
271  assert(tbcf->get_etat() == ETATQCQ) ;
272 
273  double* cf = tbcf->t ;
274 
275  // Tableau de travail
276  double* coef_tp = new double [(np+2)*nt] ;
277 
278 
279 // 1/ Sommation sur r
280 // ------------------
281 
282  som_r[base_r](cf, nr, nt, np, x, coef_tp) ;
283 
284 
285 // 2/ Sommation sur theta et phi
286 // -----------------------------
287  double* pi = coef_tp ; // pointeur courant sur les coef en theta et phi
288 
289  resu = 0 ;
290 
291  if (np > 1) { // sommation sur phi
292 
293  // On saute les coef k=0 (cos(0 phi), k=1 (sin(0 phi) et k=2
294  pi += 3*nt ;
295 
296  // Sommation sur le reste des phi (pour k=3,...,np)
297 
298  int base_t = base.b[l] & MSQ_T ;
299 
300  switch (base_t) {
301 
302  case T_COSSIN_CP : {
303 
304  for (int k=3 ; k<np+1 ; k+=2) { // k+=2 : on saute les cos(m phi)
305  int m_par = (k/2)%2 ; // 0 pour m pair, 1 pour m impair
306  double somt = 0 ;
307  for (int j=0 ; j<nt ; j++) {
308  somt += (*pi) * tab_theta(m_par, j, j0) ;
309  pi++ ; // theta suivant
310  }
311  resu += somt * tab_phi(k, k0) ;
312 
313  // On saute le cos(k*phi) :
314  pi += nt ;
315  }
316  break ;
317  }
318 
319  default: {
320  cout << "Mtbl_cf::val_point_jk_asymy: unknown theta basis ! "
321  << endl ;
322  abort () ;
323  }
324 
325  } // fin des cas sur base_t
326 
327  } // fin du cas np > 1
328 
329 
330  // Menage
331  delete [] coef_tp ;
332 
333  // Termine
334  return resu ;
335 
336 }
337 
338 
339 }
void som_tet_cossin_cp_asymy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CP ///.
Definition: som_asymy.C:281
#define MAX_BASE_2
Smaller maximum bases used for phi (and higher dimensions for now)
Definition: type_parite.h:146
void som_r_chebpim_p_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_P ///.
Definition: som_asymy.C:177
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
Definition: type_parite.h:162
Lorene prototypes.
Definition: app_hor.h:64
#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
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: mtbl_cf.h:192
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
void som_r_chebu_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition: som_asymy.C:128
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: mtbl_cf.h:196
#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
double val_point_jk_asymy(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
double * t
The array of double.
Definition: tbl.h:173
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
#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
void som_tet_cossin_ci_asymy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CI ///.
Definition: som_asymy.C:325
void som_r_chebpim_i_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_I ///.
Definition: som_asymy.C:227
void som_r_cheb_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition: som_asymy.C:80
#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
const Tbl & theta_functions(int l, int nt) const
Values of the theta basis functions at the theta collocation points.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
const Tbl & phi_functions(int l, int np) const
Values of the phi basis functions at the phi collocation points.
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
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 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_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
double val_point_asymy(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
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
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166