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