LORENE
mtbl_cf_vp_symy.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 symmetric 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 
31 char mtbl_cf_vp_symy_C[] = "$Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_vp_symy.C,v 1.3 2014/10/13 08:53:09 j_novak Exp $" ;
32 
33 /*
34  * $Id: mtbl_cf_vp_symy.C,v 1.3 2014/10/13 08:53:09 j_novak Exp $
35  * $Log: mtbl_cf_vp_symy.C,v $
36  * Revision 1.3 2014/10/13 08:53:09 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.2 2012/01/17 15:09:22 j_penner
40  * using MAX_BASE_2 for the phi coordinate
41  *
42  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
43  * LORENE
44  *
45  * Revision 2.2 2000/09/08 16:07:36 eric
46  * Ajout de la base P_COSSIN_I
47  *
48  * Revision 2.1 2000/03/06 15:57:34 eric
49  * *** empty log message ***
50  *
51  * Revision 2.0 2000/03/06 10:27:02 eric
52  * *** empty log message ***
53  *
54  *
55  * $Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_vp_symy.C,v 1.3 2014/10/13 08:53:09 j_novak Exp $
56  *
57  */
58 
59 
60 // Headers Lorene
61 #include "mtbl_cf.h"
62 #include "proto.h"
63 
64  //-------------------------------------------------------------//
65 namespace Lorene {
66  // version for an arbitrary point in (xi,theta',phi') //
67  //-------------------------------------------------------------//
68 
69 double Mtbl_cf::val_point_symy(int l, double x, double theta, double phi)
70  const {
71 
72 // Routines de sommation
73 static void (*som_r[MAX_BASE])
74  (double*, const int, const int, const int, const double, double*) ;
75 static void (*som_tet[MAX_BASE])
76  (double*, const int, const int, const double, double*) ;
77 static void (*som_phi[MAX_BASE_2])
78  (double*, const int, const double, double*) ;
79 static int premier_appel = 1 ;
80 
81 // Initialisations au premier appel
82 // --------------------------------
83  if (premier_appel == 1) {
84 
85  premier_appel = 0 ;
86 
87  for (int i=0 ; i<MAX_BASE ; i++) {
88  if(i%2==0){
89  som_phi[i/2] = som_phi_pas_prevu ;
90  }
91  som_tet[i] = som_tet_pas_prevu ;
92  som_r[i] = som_r_pas_prevu ;
93  }
94 
95  som_r[R_CHEB >> TRA_R] = som_r_cheb_symy ;
96  som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
97  som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
98  som_r[R_CHEBU >> TRA_R] = som_r_chebu_symy ;
101 
102  som_tet[T_COS >> TRA_T] = som_tet_cos ;
103  som_tet[T_SIN >> TRA_T] = som_tet_sin ;
104  som_tet[T_COS_P >> TRA_T] = som_tet_cos_p ;
105  som_tet[T_SIN_P >> TRA_T] = som_tet_sin_p ;
108 
109  som_phi[P_COSSIN >> TRA_P] = som_phi_cossin_symy ;
110  som_phi[P_COSSIN_P >> TRA_P] = som_phi_cossin_p ;
111  som_phi[P_COSSIN_I >> TRA_P] = som_phi_cossin_i ;
112 
113  } // fin des operations de premier appel
114 
115 
116  assert (etat != ETATNONDEF) ;
117 
118  double resu ; // valeur de retour
119 
120 // Cas ou tous les coefficients sont nuls :
121  if (etat == ETATZERO ) {
122  resu = 0 ;
123  return resu ;
124  }
125 
126 // Nombre de points en phi, theta et r :
127  int np = mg->get_np(l) ;
128  int nt = mg->get_nt(l) ;
129  int nr = mg->get_nr(l) ;
130 
131 // Bases de developpement :
132  int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
133  int base_t = (base.b[l] & MSQ_T) >> TRA_T ;
134  int base_p = (base.b[l] & MSQ_P) >> TRA_P ;
135 
136 //--------------------------------------
137 // Calcul de la valeur au point demande
138 //--------------------------------------
139 
140 // Pointeur sur le tableau contenant les coefficients:
141 
142  assert(etat == ETATQCQ) ;
143  Tbl* tbcf = t[l] ;
144 
145  if (tbcf->get_etat() == ETATZERO ) {
146  resu = 0 ;
147  return resu ;
148  }
149 
150 
151  assert(tbcf->get_etat() == ETATQCQ) ;
152 
153  double* cf = tbcf->t ;
154 
155  // Tableaux de travail
156  double* trp = new double [np+2] ;
157  double* trtp = new double [(np+2)*nt] ;
158 
159  if (nr == 1) {
160 
161 // Cas particulier nr = 1 (Fonction purement angulaire) :
162 // ----------------------
163 
164  som_tet[base_t](cf, nt, np, theta, trp) ; // sommation sur theta
165  }
166  else {
167 
168 // Cas general
169 // -----------
170 
171  som_r[base_r](cf, nr, nt, np, x, trtp) ; // sommation sur r
172  som_tet[base_t](trtp, nt, np, theta, trp) ; // sommation sur theta
173  }
174 
175 // Sommation sur phi
176 // -----------------
177 
178  if (np == 1) {
179  resu = trp[0] ; // cas axisymetrique
180  }
181  else {
182  som_phi[base_p](trp, np, phi, &resu) ; // sommation sur phi
183  }
184 
185  // Menage
186  delete [] trp ;
187  delete [] trtp ;
188 
189  // Termine
190  return resu ;
191 
192 }
193 
194 
195 
196  //-------------------------------------------------------------//
197  // version for an arbitrary point in xi //
198  // but collocation point in (theta',phi') //
199  //-------------------------------------------------------------//
200 
201 double Mtbl_cf::val_point_jk_symy(int l, double x, int j0, int k0) const {
202 
203 // Routines de sommation
204 static void (*som_r[MAX_BASE])
205  (double*, const int, const int, const int, const double, double*) ;
206 static int premier_appel = 1 ;
207 
208 // Initialisations au premier appel
209 // --------------------------------
210  if (premier_appel == 1) {
211 
212  premier_appel = 0 ;
213 
214  for (int i=0 ; i<MAX_BASE ; i++) {
215  som_r[i] = som_r_pas_prevu ;
216  }
217 
218  som_r[R_CHEB >> TRA_R] = som_r_cheb_symy ;
219  som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
220  som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
221  som_r[R_CHEBU >> TRA_R] = som_r_chebu_symy ;
224 
225  } // fin des operations de premier appel
226 
227  assert (etat != ETATNONDEF) ;
228 
229  double resu ; // valeur de retour
230 
231 // Cas ou tous les coefficients sont nuls :
232  if (etat == ETATZERO ) {
233  resu = 0 ;
234  return resu ;
235  }
236 
237 // Nombre de points en phi, theta et r :
238  int np = mg->get_np(l) ;
239  int nt = mg->get_nt(l) ;
240  int nr = mg->get_nr(l) ;
241 
242 // Bases de developpement :
243  int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
244 
245 //------------------------------------------------------------------------
246 // Valeurs des fonctions de base en phi aux points de collocation en phi
247 // et des fonctions de base en theta aux points de collocation en theta
248 //------------------------------------------------------------------------
249 
250  Tbl tab_phi = base.phi_functions(l, np) ;
251  Tbl tab_theta = base.theta_functions(l, nt) ;
252 
253 
254 //--------------------------------------
255 // Calcul de la valeur au point demande
256 //--------------------------------------
257 
258 // Pointeur sur le tableau contenant les coefficients:
259 
260  assert(etat == ETATQCQ) ;
261  Tbl* tbcf = t[l] ;
262 
263  if (tbcf->get_etat() == ETATZERO ) {
264  resu = 0 ;
265  return resu ;
266  }
267 
268 
269  assert(tbcf->get_etat() == ETATQCQ) ;
270 
271  double* cf = tbcf->t ;
272 
273  // Tableau de travail
274  double* coef_tp = new double [(np+2)*nt] ;
275 
276 
277 // 1/ Sommation sur r
278 // ------------------
279 
280  som_r[base_r](cf, nr, nt, np, x, coef_tp) ;
281 
282 
283 // 2/ Sommation sur theta et phi
284 // -----------------------------
285  double* pi = coef_tp ; // pointeur courant sur les coef en theta et phi
286 
287 // Sommation sur le premier phi, k=0
288 
289  double somt = 0 ;
290  for (int j=0 ; j<nt ; j++) {
291  somt += (*pi) * tab_theta(0, j, j0) ;
292  pi++ ; // theta suivant
293  }
294  resu = somt * tab_phi(0, k0) ;
295 
296  if (np > 1) { // sommation sur phi
297 
298  // On saute le phi suivant (sin(0)), k=1
299  pi += nt ;
300 
301  // Sommation sur le reste des phi (pour k=2,...,np)
302 
303  int base_t = base.b[l] & MSQ_T ;
304 
305  switch (base_t) {
306 
307  case T_COSSIN_CP : {
308 
309  for (int k=2 ; k<np+1 ; k+=2) { // k+=2 : on saute les sin(m phi)
310  int m_par = (k/2)%2 ; // 0 pour m pair, 1 pour m impair
311  somt = 0 ;
312  for (int j=0 ; j<nt ; j++) {
313  somt += (*pi) * tab_theta(m_par, j, j0) ;
314  pi++ ; // theta suivant
315  }
316  resu += somt * tab_phi(k, k0) ;
317 
318  // On saute le sin(k*phi) :
319  pi += nt ;
320  }
321  break ;
322  }
323 
324  default: {
325  cout << "Mtbl_cf::val_point_jk_symy: unknown theta basis ! "
326  << endl ;
327  abort () ;
328  }
329 
330  } // fin des cas sur base_t
331 
332  } // fin du cas np > 1
333 
334 
335  // Menage
336  delete [] coef_tp ;
337 
338  // Termine
339  return resu ;
340 
341 }
342 
343 
344 }
#define MAX_BASE_2
Smaller maximum bases used for phi (and higher dimensions for now)
Definition: type_parite.h:146
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
void som_tet_cossin_cp_symy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CP ///.
Definition: som_symy.C:338
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
void som_r_chebu_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition: som_symy.C:145
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
Definition: type_parite.h:162
void som_r_cheb_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition: som_symy.C:82
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
void som_r_chebpim_i_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_I ///.
Definition: som_symy.C:269
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
void som_r_chebpim_p_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_P ///.
Definition: som_symy.C:206
#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
#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
#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
double val_point_symy(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 ...
#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
void som_tet_cossin_ci_symy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CI ///.
Definition: som_symy.C:393
#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
double val_point_jk_symy(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 MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
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