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