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
31char 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 //-------------------------------------------------------------//
65namespace Lorene {
66 // version for an arbitrary point in (xi,theta',phi') //
67 //-------------------------------------------------------------//
68
69double Mtbl_cf::val_point_symy(int l, double x, double theta, double phi)
70 const {
71
72// Routines de sommation
73static void (*som_r[MAX_BASE])
74 (double*, const int, const int, const int, const double, double*) ;
75static void (*som_tet[MAX_BASE])
76 (double*, const int, const int, const double, double*) ;
77static void (*som_phi[MAX_BASE_2])
78 (double*, const int, const double, double*) ;
79static 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
201double Mtbl_cf::val_point_jk_symy(int l, double x, int j0, int k0) const {
202
203// Routines de sommation
204static void (*som_r[MAX_BASE])
205 (double*, const int, const int, const int, const double, double*) ;
206static 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}
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_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 ...
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...
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
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_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition som_symy.C:82
void som_tet_cossin_ci_symy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CI ///.
Definition som_symy.C:393
void som_r_chebu_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition som_symy.C:145
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
void som_tet_cossin_cp_symy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CP ///.
Definition som_symy.C:338
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