LORENE
mtbl_cf_val_point.C
1/*
2 * Member functions of the Mtbl_cf class for computing the value of a field
3 * at an arbitrary point
4 *
5 * (see file mtbl_cf.h for the documentation).
6 */
7
8/*
9 * Copyright (c) 1999-2002 Eric Gourgoulhon
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char mtbl_cf_val_point_C[] = "$Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_val_point.C,v 1.15 2014/10/13 08:53:09 j_novak Exp $" ;
31
32/*
33 * $Id: mtbl_cf_val_point.C,v 1.15 2014/10/13 08:53:09 j_novak Exp $
34 * $Log: mtbl_cf_val_point.C,v $
35 * Revision 1.15 2014/10/13 08:53:09 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.14 2013/06/13 14:18:18 j_novak
39 * Inclusion of new bases R_LEG, R_LEGP and R_LEGI.
40 *
41 * Revision 1.13 2012/01/17 15:09:05 j_penner
42 * using MAX_BASE_2 for the phi coordinate
43 *
44 * Revision 1.12 2009/10/08 16:21:16 j_novak
45 * Addition of new bases T_COS and T_SIN.
46 *
47 * Revision 1.11 2007/12/20 09:11:08 jl_cornou
48 * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
49 *
50 * Revision 1.10 2007/12/11 15:28:16 jl_cornou
51 * Jacobi(0,2) polynomials partially implemented
52 *
53 * Revision 1.9 2007/10/23 17:15:13 j_novak
54 * Added the bases T_COSSIN_C and T_COSSIN_S
55 *
56 * Revision 1.8 2006/06/06 14:57:01 j_novak
57 * Summation functions for angular coefficients at xi=+/-1.
58 *
59 * Revision 1.7 2006/05/30 08:30:15 n_vasset
60 * Implementation of sine-like bases (T_SIN_P, T_SIN_I, T_COSSIN_SI, etc...).
61 *
62 * Revision 1.6 2005/05/27 14:55:00 j_novak
63 * Added new bases T_COSSIN_CI and T_COS_I
64 *
65 * Revision 1.5 2005/02/16 15:10:39 m_forot
66 * Correct the case T_COSSIN_C
67 *
68 * Revision 1.4 2004/12/17 13:35:03 m_forot
69 * Add the case T_LEG
70 *
71 * Revision 1.3 2002/05/11 12:37:31 e_gourgoulhon
72 * Added basis T_COSSIN_SI.
73 *
74 * Revision 1.2 2002/05/05 16:22:33 e_gourgoulhon
75 * Added the case of the theta basis T_COSSIN_SP.
76 *
77 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
78 * LORENE
79 *
80 * Revision 1.7 2000/09/08 16:26:43 eric
81 * Ajout de la base T_SIN_I.
82 *
83 * Revision 1.6 2000/09/08 16:07:16 eric
84 * Ajout de la base P_COSSIN_I
85 *
86 * Revision 1.5 2000/09/06 14:00:19 eric
87 * Ajout de la base T_COS_I.
88 *
89 * Revision 1.4 1999/12/29 13:11:42 eric
90 * Ajout de la fonction val_point_jk.
91 *
92 * Revision 1.3 1999/12/07 15:10:45 eric
93 * *** empty log message ***
94 *
95 * Revision 1.2 1999/12/07 14:52:34 eric
96 * Changement ordre des arguments (phi,theta,xi) --> (xi,theta,phi)
97 *
98 * Revision 1.1 1999/12/06 16:47:39 eric
99 * Initial revision
100 *
101 *
102 * $Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_val_point.C,v 1.15 2014/10/13 08:53:09 j_novak Exp $
103 *
104 */
105
106// Headers Lorene
107#include "mtbl_cf.h"
108#include "proto.h"
109
110 //-------------------------------------------------------------//
111namespace Lorene {
112 // version for an arbitrary point in (xi,theta',phi') //
113 //-------------------------------------------------------------//
114
115double Mtbl_cf::val_point(int l, double x, double theta, double phi) const {
116
117// Routines de sommation
118static void (*som_r[MAX_BASE])
119 (double*, const int, const int, const int, const double, double*) ;
120static void (*som_tet[MAX_BASE])
121 (double*, const int, const int, const double, double*) ;
122static void (*som_phi[MAX_BASE_2])
123 (double*, const int, const double, double*) ;
124static int premier_appel = 1 ;
125
126// Initialisations au premier appel
127// --------------------------------
128 if (premier_appel == 1) {
129
130 premier_appel = 0 ;
131
132 for (int i=0 ; i<MAX_BASE ; i++) {
133 if(i%2 == 0){
134 som_phi[i/2] = som_phi_pas_prevu ;
135 }
136 som_tet[i] = som_tet_pas_prevu ;
137 som_r[i] = som_r_pas_prevu ;
138 }
139
140 som_r[R_CHEB >> TRA_R] = som_r_cheb ;
141 som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
142 som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
143 som_r[R_CHEBU >> TRA_R] = som_r_chebu ;
144 som_r[R_CHEBPIM_P >> TRA_R] = som_r_chebpim_p ;
145 som_r[R_CHEBPIM_I >> TRA_R] = som_r_chebpim_i ;
146 som_r[R_CHEBPI_P >> TRA_R] = som_r_chebpi_p ;
147 som_r[R_CHEBPI_I >> TRA_R] = som_r_chebpi_i ;
148 som_r[R_LEG >> TRA_R] = som_r_leg ;
149 som_r[R_LEGP >> TRA_R] = som_r_legp ;
150 som_r[R_LEGI >> TRA_R] = som_r_legi ;
151 som_r[R_JACO02 >> TRA_R] = som_r_jaco02 ;
152
153 som_tet[T_COS >> TRA_T] = som_tet_cos ;
154 som_tet[T_SIN >> TRA_T] = som_tet_sin ;
155 som_tet[T_COS_P >> TRA_T] = som_tet_cos_p ;
156 som_tet[T_COS_I >> TRA_T] = som_tet_cos_i ;
157 som_tet[T_SIN_P >> TRA_T] = som_tet_sin_p ;
158 som_tet[T_SIN_I >> TRA_T] = som_tet_sin_i ;
159 som_tet[T_COSSIN_CP >> TRA_T] = som_tet_cossin_cp ;
160 som_tet[T_COSSIN_CI >> TRA_T] = som_tet_cossin_ci ;
161 som_tet[T_COSSIN_SP >> TRA_T] = som_tet_cossin_sp ;
162 som_tet[T_COSSIN_SI >> TRA_T] = som_tet_cossin_si ;
163 som_tet[T_COSSIN_C >> TRA_T] = som_tet_cossin_c ;
164 som_tet[T_COSSIN_S >> TRA_T] = som_tet_cossin_s ;
165
166 som_phi[P_COSSIN >> TRA_P] = som_phi_cossin ;
167 som_phi[P_COSSIN_P >> TRA_P] = som_phi_cossin_p ;
168 som_phi[P_COSSIN_I >> TRA_P] = som_phi_cossin_i ;
169
170 } // fin des operations de premier appel
171
172
173 assert (etat != ETATNONDEF) ;
174
175 double resu ; // valeur de retour
176
177// Cas ou tous les coefficients sont nuls :
178 if (etat == ETATZERO ) {
179 resu = 0 ;
180 return resu ;
181 }
182
183// Nombre de points en phi, theta et r :
184 int np = mg->get_np(l) ;
185 int nt = mg->get_nt(l) ;
186 int nr = mg->get_nr(l) ;
187
188// Bases de developpement :
189 int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
190 int base_t = (base.b[l] & MSQ_T) >> TRA_T ;
191 int base_p = (base.b[l] & MSQ_P) >> TRA_P ;
192
193//--------------------------------------
194// Calcul de la valeur au point demande
195//--------------------------------------
196
197// Pointeur sur le tableau contenant les coefficients:
198
199 assert(etat == ETATQCQ) ;
200 Tbl* tbcf = t[l] ;
201
202 if (tbcf->get_etat() == ETATZERO ) {
203 resu = 0 ;
204 return resu ;
205 }
206
207
208 assert(tbcf->get_etat() == ETATQCQ) ;
209
210 double* cf = tbcf->t ;
211
212 // Tableaux de travail
213 double* trp = new double [np+2] ;
214 double* trtp = new double [(np+2)*nt] ;
215
216 if (nr == 1) {
217
218// Cas particulier nr = 1 (Fonction purement angulaire) :
219// ----------------------
220
221 som_tet[base_t](cf, nt, np, theta, trp) ; // sommation sur theta
222 }
223 else {
224
225// Cas general
226// -----------
227
228 som_r[base_r](cf, nr, nt, np, x, trtp) ; // sommation sur r
229 som_tet[base_t](trtp, nt, np, theta, trp) ; // sommation sur theta
230 }
231
232// Sommation sur phi
233// -----------------
234
235 if (np == 1) {
236 resu = trp[0] ; // cas axisymetrique
237 }
238 else {
239 som_phi[base_p](trp, np, phi, &resu) ; // sommation sur phi
240 }
241
242 // Menage
243 delete [] trp ;
244 delete [] trtp ;
245
246 // Termine
247 return resu ;
248
249}
250
251
252
253 //-------------------------------------------------------------//
254 // version for an arbitrary point in xi //
255 // but collocation point in (theta',phi') //
256 //-------------------------------------------------------------//
257
258double Mtbl_cf::val_point_jk(int l, double x, int j0, int k0) const {
259
260// Routines de sommation
261static void (*som_r[MAX_BASE])
262 (double*, const int, const int, const int, const double, double*) ;
263static int premier_appel = 1 ;
264
265// Initialisations au premier appel
266// --------------------------------
267 if (premier_appel == 1) {
268
269 premier_appel = 0 ;
270
271 for (int i=0 ; i<MAX_BASE ; i++) {
272 som_r[i] = som_r_pas_prevu ;
273 }
274
275 som_r[R_CHEB >> TRA_R] = som_r_cheb ;
276 som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
277 som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
278 som_r[R_CHEBU >> TRA_R] = som_r_chebu ;
279 som_r[R_CHEBPIM_P >> TRA_R] = som_r_chebpim_p ;
280 som_r[R_CHEBPIM_I >> TRA_R] = som_r_chebpim_i ;
281 som_r[R_CHEBPI_P >> TRA_R] = som_r_chebpi_p ;
282 som_r[R_CHEBPI_I >> TRA_R] = som_r_chebpi_i ;
283 som_r[R_JACO02 >> TRA_R] = som_r_jaco02 ;
284
285 } // fin des operations de premier appel
286
287 assert (etat != ETATNONDEF) ;
288
289 double resu ; // valeur de retour
290
291// Cas ou tous les coefficients sont nuls :
292 if (etat == ETATZERO ) {
293 resu = 0 ;
294 return resu ;
295 }
296
297// Nombre de points en phi, theta et r :
298 int np = mg->get_np(l) ;
299 int nt = mg->get_nt(l) ;
300 int nr = mg->get_nr(l) ;
301
302// Bases de developpement :
303 int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
304
305//------------------------------------------------------------------------
306// Valeurs des fonctions de base en phi aux points de collocation en phi
307// et des fonctions de base en theta aux points de collocation en theta
308//------------------------------------------------------------------------
309
310 Tbl tab_phi = base.phi_functions(l, np) ;
311 Tbl tab_theta = base.theta_functions(l, nt) ;
312
313
314//--------------------------------------
315// Calcul de la valeur au point demande
316//--------------------------------------
317
318// Pointeur sur le tableau contenant les coefficients:
319
320 assert(etat == ETATQCQ) ;
321 Tbl* tbcf = t[l] ;
322
323 if (tbcf->get_etat() == ETATZERO ) {
324 resu = 0 ;
325 return resu ;
326 }
327
328
329 assert(tbcf->get_etat() == ETATQCQ) ;
330
331 double* cf = tbcf->t ;
332
333 // Tableau de travail
334 double* coef_tp = new double [(np+2)*nt] ;
335
336
337// 1/ Sommation sur r
338// ------------------
339
340 som_r[base_r](cf, nr, nt, np, x, coef_tp) ;
341
342
343// 2/ Sommation sur theta et phi
344// -----------------------------
345 double* pi = coef_tp ; // pointeur courant sur les coef en theta et phi
346
347// Sommation sur le premier phi, k=0
348
349 double somt = 0 ;
350 for (int j=0 ; j<nt ; j++) {
351 somt += (*pi) * tab_theta(0, j, j0) ;
352 pi++ ; // theta suivant
353 }
354 resu = somt * tab_phi(0, k0) ;
355
356 if (np > 1) { // sommation sur phi
357
358 // On saute le phi suivant (sin(0)), k=1
359 pi += nt ;
360
361 // Sommation sur le reste des phi (pour k=2,...,np)
362
363 int base_t = base.b[l] & MSQ_T ;
364
365 switch (base_t) {
366
367 case T_COS :
368 case T_SIN :
369 case T_SIN_P :
370 case T_SIN_I :
371 case T_COS_I :
372 case T_COS_P : {
373
374 for (int k=2 ; k<np+1 ; k++) {
375 somt = 0 ;
376 for (int j=0 ; j<nt ; j++) {
377 somt += (*pi) * tab_theta(0, j, j0) ;
378 pi++ ; // theta suivant
379 }
380 resu += somt * tab_phi(k, k0) ;
381 }
382 break ;
383 }
384
385 case T_COSSIN_C :
386 case T_COSSIN_S :
387 case T_COSSIN_SP :
388 case T_COSSIN_SI :
389 case T_COSSIN_CI :
390 case T_COSSIN_CP : {
391
392 for (int k=2 ; k<np+1 ; k++) {
393 int m_par = (k/2)%2 ; // 0 pour m pair, 1 pour m impair
394 somt = 0 ;
395 for (int j=0 ; j<nt ; j++) {
396 somt += (*pi) * tab_theta(m_par, j, j0) ;
397 pi++ ; // theta suivant
398 }
399 resu += somt * tab_phi(k, k0) ;
400 }
401 break ;
402 }
403
404 default: {
405 cout << "Mtbl_cf::val_point_jk: unknown theta basis ! " << endl ;
406 abort () ;
407 }
408
409 } // fin des cas sur base_t
410
411 } // fin du cas np > 1
412
413
414 // Menage
415 delete [] coef_tp ;
416
417 // Termine
418 return resu ;
419
420}
421
422
423 //-------------------------------------------------------------//
424 // version for xi = 1 //
425 // and collocation point in (theta',phi') //
426 //-------------------------------------------------------------//
427
428double Mtbl_cf::val_out_bound_jk(int l, int j0, int k0) const {
429
430#ifndef NDEBUG
431// Bases de developpement :
432 int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
433 assert((base_r == R_CHEB) || (base_r == R_CHEBU) || (base_r == R_CHEBP)
434 || (base_r == R_CHEBI) || (base_r == R_CHEBPIM_P) || (base_r == R_CHEBPIM_I)
435 || (base_r == R_CHEBPI_P) || (base_r == R_CHEBPI_I) || (base_r == R_JACO02)) ;
436#endif
437
438 int nr = mg->get_nr(l) ;
439 double resu = 0 ;
440 for (int i=0; i<nr; i++)
441 resu += operator()(l, k0, j0, i) ;
442
443 return resu ;
444}
445
446
447 //-------------------------------------------------------------//
448 // version for xi = -1 //
449 // and collocation point in (theta',phi') //
450 //-------------------------------------------------------------//
451
452double Mtbl_cf::val_in_bound_jk(int l, int j0, int k0) const {
453
454#ifndef NDEBUG
455// Bases de developpement :
456 int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
457 assert((base_r == R_CHEB) || (base_r == R_CHEBU)) ;
458#endif
459
460 int nr = mg->get_nr(l) ;
461 double resu = 0 ;
462 int pari = 1 ;
463 for (int i=0; i<nr; i++) {
464 resu += pari*operator()(l, k0, j0, i) ;
465 pari = - pari ;
466 }
467
468 return resu ;
469}
470
471}
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
double val_point_jk(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...
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
const Tbl & operator()(int l) const
Read-only of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:306
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(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_out_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
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 R_LEGP
base de Legendre paire (rare) seulement
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
#define R_LEGI
base de Legendre impaire (rare) seulement
#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 T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
#define R_LEG
base de Legendre ordinaire (fin)
#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_SIN_I
dev. sin seulement, harmoniques impaires
#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 T_COS_I
dev. cos seulement, harmoniques impaires
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Lorene prototypes.
Definition app_hor.h:64