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