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
30char 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 //-------------------------------------------------------------//
67namespace Lorene {
68 // version for an arbitrary point in (xi,theta',phi') //
69 //-------------------------------------------------------------//
70
71double Mtbl_cf::val_point_asymy(int l, double x, double theta, double phi)
72 const {
73
74// Routines de sommation
75static void (*som_r[MAX_BASE])
76 (double*, const int, const int, const int, const double, double*) ;
77static void (*som_tet[MAX_BASE])
78 (double*, const int, const int, const double, double*) ;
79static void (*som_phi[MAX_BASE_2])
80 (double*, const int, const double, double*) ;
81static 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
203double Mtbl_cf::val_point_jk_asymy(int l, double x, int j0, int k0) const {
204
205// Routines de sommation
206static void (*som_r[MAX_BASE])
207 (double*, const int, const int, const int, const double, double*) ;
208static 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}
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition base_val.h:331
const Tbl & phi_functions(int l, int np) const
Values of the phi basis functions at the phi collocation points.
const Tbl & theta_functions(int l, int nt) const
Values of the theta basis functions at the theta collocation points.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.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...
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition mtbl_cf.h:192
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition mtbl_cf.h:196
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
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 ...
Basic array class.
Definition tbl.h:161
int get_etat() const
Gives the logical state.
Definition tbl.h:394
double * t
The array of double.
Definition tbl.h:173
#define MAX_BASE_2
Smaller maximum bases used for phi (and higher dimensions for now)
#define TRA_T
Translation en Theta, used for a bitwise shift (in hex)
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define T_SIN_P
dev. sin seulement, harmoniques paires
#define MSQ_R
Extraction de l'info sur R.
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define T_COS_P
dev. cos seulement, harmoniques paires
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
#define P_COSSIN
dev. standart
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define MSQ_T
Extraction de l'info sur Theta.
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
#define R_CHEB
base de Chebychev ordinaire (fin)
#define T_COS
dev. cos seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define MSQ_P
Extraction de l'info sur Phi.
#define T_SIN
dev. sin seulement
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
Lorene prototypes.
Definition app_hor.h:64
void som_r_cheb_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition som_asymy.C:80
void som_tet_cossin_cp_asymy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CP ///.
Definition som_asymy.C:281
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
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_tet_cossin_ci_asymy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CI ///.
Definition som_asymy.C:325
void som_r_chebu_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition som_asymy.C:128