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 
30 char 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  //-------------------------------------------------------------//
111 namespace Lorene {
112  // version for an arbitrary point in (xi,theta',phi') //
113  //-------------------------------------------------------------//
114 
115 double Mtbl_cf::val_point(int l, double x, double theta, double phi) const {
116 
117 // Routines de sommation
118 static void (*som_r[MAX_BASE])
119  (double*, const int, const int, const int, const double, double*) ;
120 static void (*som_tet[MAX_BASE])
121  (double*, const int, const int, const double, double*) ;
122 static void (*som_phi[MAX_BASE_2])
123  (double*, const int, const double, double*) ;
124 static 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 
258 double Mtbl_cf::val_point_jk(int l, double x, int j0, int k0) const {
259 
260 // Routines de sommation
261 static void (*som_r[MAX_BASE])
262  (double*, const int, const int, const int, const double, double*) ;
263 static 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 
428 double 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 
452 double 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 }
#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:462
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex)
Definition: type_parite.h:162
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene prototypes.
Definition: app_hor.h:64
const Tbl & operator()(int l) const
Read-only of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:306
#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:192
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
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...
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
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...
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
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...
#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:196
#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
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
double * t
The array of double.
Definition: tbl.h:173
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
#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 T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
#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
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 ...
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:452
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
#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:200
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition: type_parite.h:249
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
#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
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166