mtbl_cf_val_point.C
/*
* Member functions of the Mtbl_cf class for computing the value of a field
* at an arbitrary point
*
* (see file mtbl_cf.h for the documentation).
*/
/*
* Copyright (c) 1999-2002 Eric Gourgoulhon
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: mtbl_cf_val_point.C,v 1.15 2014/10/13 08:53:09 j_novak Exp $
* $Log: mtbl_cf_val_point.C,v $
* Revision 1.15 2014/10/13 08:53:09 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.14 2013/06/13 14:18:18 j_novak
* Inclusion of new bases R_LEG, R_LEGP and R_LEGI.
*
* Revision 1.13 2012/01/17 15:09:05 j_penner
* using MAX_BASE_2 for the phi coordinate
*
* Revision 1.12 2009/10/08 16:21:16 j_novak
* Addition of new bases T_COS and T_SIN.
*
* Revision 1.11 2007/12/20 09:11:08 jl_cornou
* Correction of an error in op_sxpun about Jacobi(0,2) polynomials
*
* Revision 1.10 2007/12/11 15:28:16 jl_cornou
* Jacobi(0,2) polynomials partially implemented
*
* Revision 1.9 2007/10/23 17:15:13 j_novak
* Added the bases T_COSSIN_C and T_COSSIN_S
*
* Revision 1.8 2006/06/06 14:57:01 j_novak
* Summation functions for angular coefficients at xi=+/-1.
*
* Revision 1.7 2006/05/30 08:30:15 n_vasset
* Implementation of sine-like bases (T_SIN_P, T_SIN_I, T_COSSIN_SI, etc...).
*
* Revision 1.6 2005/05/27 14:55:00 j_novak
* Added new bases T_COSSIN_CI and T_COS_I
*
* Revision 1.5 2005/02/16 15:10:39 m_forot
* Correct the case T_COSSIN_C
*
* Revision 1.4 2004/12/17 13:35:03 m_forot
* Add the case T_LEG
*
* Revision 1.3 2002/05/11 12:37:31 e_gourgoulhon
* Added basis T_COSSIN_SI.
*
* Revision 1.2 2002/05/05 16:22:33 e_gourgoulhon
* Added the case of the theta basis T_COSSIN_SP.
*
* Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
* LORENE
*
* Revision 1.7 2000/09/08 16:26:43 eric
* Ajout de la base T_SIN_I.
*
* Revision 1.6 2000/09/08 16:07:16 eric
* Ajout de la base P_COSSIN_I
*
* Revision 1.5 2000/09/06 14:00:19 eric
* Ajout de la base T_COS_I.
*
* Revision 1.4 1999/12/29 13:11:42 eric
* Ajout de la fonction val_point_jk.
*
* Revision 1.3 1999/12/07 15:10:45 eric
* *** empty log message ***
*
* Revision 1.2 1999/12/07 14:52:34 eric
* Changement ordre des arguments (phi,theta,xi) --> (xi,theta,phi)
*
* Revision 1.1 1999/12/06 16:47:39 eric
* Initial revision
*
*
* $Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_val_point.C,v 1.15 2014/10/13 08:53:09 j_novak Exp $
*
*/
// Headers Lorene
#include "mtbl_cf.h"
#include "proto.h"
//-------------------------------------------------------------//
namespace Lorene {
// version for an arbitrary point in (xi,theta',phi') //
//-------------------------------------------------------------//
double Mtbl_cf::val_point(int l, double x, double theta, double phi) const {
// Routines de sommation
static void (*som_r[MAX_BASE])
(double*, const int, const int, const int, const double, double*) ;
static void (*som_tet[MAX_BASE])
(double*, const int, const int, const double, double*) ;
static void (*som_phi[MAX_BASE_2])
(double*, const int, const double, double*) ;
static int premier_appel = 1 ;
// Initialisations au premier appel
// --------------------------------
if (premier_appel == 1) {
premier_appel = 0 ;
for (int i=0 ; i<MAX_BASE ; i++) {
if(i%2 == 0){
som_phi[i/2] = som_phi_pas_prevu ;
}
som_tet[i] = som_tet_pas_prevu ;
som_r[i] = som_r_pas_prevu ;
}
som_r[R_CHEB >> TRA_R] = som_r_cheb ;
som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
som_r[R_CHEBU >> TRA_R] = som_r_chebu ;
som_r[R_CHEBPIM_P >> TRA_R] = som_r_chebpim_p ;
som_r[R_CHEBPIM_I >> TRA_R] = som_r_chebpim_i ;
som_r[R_CHEBPI_P >> TRA_R] = som_r_chebpi_p ;
som_r[R_CHEBPI_I >> TRA_R] = som_r_chebpi_i ;
som_r[R_LEG >> TRA_R] = som_r_leg ;
som_r[R_LEGP >> TRA_R] = som_r_legp ;
som_r[R_LEGI >> TRA_R] = som_r_legi ;
som_r[R_JACO02 >> TRA_R] = som_r_jaco02 ;
som_tet[T_COS >> TRA_T] = som_tet_cos ;
som_tet[T_SIN >> TRA_T] = som_tet_sin ;
som_tet[T_COS_P >> TRA_T] = som_tet_cos_p ;
som_tet[T_COS_I >> TRA_T] = som_tet_cos_i ;
som_tet[T_SIN_P >> TRA_T] = som_tet_sin_p ;
som_tet[T_SIN_I >> TRA_T] = som_tet_sin_i ;
som_tet[T_COSSIN_CP >> TRA_T] = som_tet_cossin_cp ;
som_tet[T_COSSIN_CI >> TRA_T] = som_tet_cossin_ci ;
som_tet[T_COSSIN_SP >> TRA_T] = som_tet_cossin_sp ;
som_tet[T_COSSIN_SI >> TRA_T] = som_tet_cossin_si ;
som_tet[T_COSSIN_C >> TRA_T] = som_tet_cossin_c ;
som_tet[T_COSSIN_S >> TRA_T] = som_tet_cossin_s ;
som_phi[P_COSSIN >> TRA_P] = som_phi_cossin ;
som_phi[P_COSSIN_P >> TRA_P] = som_phi_cossin_p ;
som_phi[P_COSSIN_I >> TRA_P] = som_phi_cossin_i ;
} // fin des operations de premier appel
assert (etat != ETATNONDEF) ;
double resu ; // valeur de retour
// Cas ou tous les coefficients sont nuls :
if (etat == ETATZERO ) {
resu = 0 ;
return resu ;
}
// Nombre de points en phi, theta et r :
int np = mg->get_np(l) ;
int nt = mg->get_nt(l) ;
int nr = mg->get_nr(l) ;
// Bases de developpement :
int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
int base_t = (base.b[l] & MSQ_T) >> TRA_T ;
int base_p = (base.b[l] & MSQ_P) >> TRA_P ;
//--------------------------------------
// Calcul de la valeur au point demande
//--------------------------------------
// Pointeur sur le tableau contenant les coefficients:
assert(etat == ETATQCQ) ;
Tbl* tbcf = t[l] ;
if (tbcf->get_etat() == ETATZERO ) {
resu = 0 ;
return resu ;
}
assert(tbcf->get_etat() == ETATQCQ) ;
double* cf = tbcf->t ;
// Tableaux de travail
double* trp = new double [np+2] ;
double* trtp = new double [(np+2)*nt] ;
if (nr == 1) {
// Cas particulier nr = 1 (Fonction purement angulaire) :
// ----------------------
som_tet[base_t](cf, nt, np, theta, trp) ; // sommation sur theta
}
else {
// Cas general
// -----------
som_r[base_r](cf, nr, nt, np, x, trtp) ; // sommation sur r
som_tet[base_t](trtp, nt, np, theta, trp) ; // sommation sur theta
}
// Sommation sur phi
// -----------------
if (np == 1) {
resu = trp[0] ; // cas axisymetrique
}
else {
som_phi[base_p](trp, np, phi, &resu) ; // sommation sur phi
}
// Menage
delete [] trp ;
delete [] trtp ;
// Termine
return resu ;
}
//-------------------------------------------------------------//
// version for an arbitrary point in xi //
// but collocation point in (theta',phi') //
//-------------------------------------------------------------//
double Mtbl_cf::val_point_jk(int l, double x, int j0, int k0) const {
// Routines de sommation
static void (*som_r[MAX_BASE])
(double*, const int, const int, const int, const double, double*) ;
static int premier_appel = 1 ;
// Initialisations au premier appel
// --------------------------------
if (premier_appel == 1) {
premier_appel = 0 ;
for (int i=0 ; i<MAX_BASE ; i++) {
som_r[i] = som_r_pas_prevu ;
}
som_r[R_CHEB >> TRA_R] = som_r_cheb ;
som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
som_r[R_CHEBU >> TRA_R] = som_r_chebu ;
som_r[R_CHEBPIM_P >> TRA_R] = som_r_chebpim_p ;
som_r[R_CHEBPIM_I >> TRA_R] = som_r_chebpim_i ;
som_r[R_CHEBPI_P >> TRA_R] = som_r_chebpi_p ;
som_r[R_CHEBPI_I >> TRA_R] = som_r_chebpi_i ;
som_r[R_JACO02 >> TRA_R] = som_r_jaco02 ;
} // fin des operations de premier appel
assert (etat != ETATNONDEF) ;
double resu ; // valeur de retour
// Cas ou tous les coefficients sont nuls :
if (etat == ETATZERO ) {
resu = 0 ;
return resu ;
}
// Nombre de points en phi, theta et r :
int np = mg->get_np(l) ;
int nt = mg->get_nt(l) ;
int nr = mg->get_nr(l) ;
// Bases de developpement :
int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
//------------------------------------------------------------------------
// Valeurs des fonctions de base en phi aux points de collocation en phi
// et des fonctions de base en theta aux points de collocation en theta
//------------------------------------------------------------------------
Tbl tab_phi = base.phi_functions(l, np) ;
Tbl tab_theta = base.theta_functions(l, nt) ;
//--------------------------------------
// Calcul de la valeur au point demande
//--------------------------------------
// Pointeur sur le tableau contenant les coefficients:
assert(etat == ETATQCQ) ;
Tbl* tbcf = t[l] ;
if (tbcf->get_etat() == ETATZERO ) {
resu = 0 ;
return resu ;
}
assert(tbcf->get_etat() == ETATQCQ) ;
double* cf = tbcf->t ;
// Tableau de travail
double* coef_tp = new double [(np+2)*nt] ;
// 1/ Sommation sur r
// ------------------
som_r[base_r](cf, nr, nt, np, x, coef_tp) ;
// 2/ Sommation sur theta et phi
// -----------------------------
double* pi = coef_tp ; // pointeur courant sur les coef en theta et phi
// Sommation sur le premier phi, k=0
double somt = 0 ;
for (int j=0 ; j<nt ; j++) {
somt += (*pi) * tab_theta(0, j, j0) ;
pi++ ; // theta suivant
}
resu = somt * tab_phi(0, k0) ;
if (np > 1) { // sommation sur phi
// On saute le phi suivant (sin(0)), k=1
pi += nt ;
// Sommation sur le reste des phi (pour k=2,...,np)
int base_t = base.b[l] & MSQ_T ;
switch (base_t) {
case T_COS :
case T_SIN :
case T_SIN_P :
case T_SIN_I :
case T_COS_I :
case T_COS_P : {
for (int k=2 ; k<np+1 ; k++) {
somt = 0 ;
for (int j=0 ; j<nt ; j++) {
somt += (*pi) * tab_theta(0, j, j0) ;
pi++ ; // theta suivant
}
resu += somt * tab_phi(k, k0) ;
}
break ;
}
case T_COSSIN_C :
case T_COSSIN_S :
case T_COSSIN_SP :
case T_COSSIN_SI :
case T_COSSIN_CI :
case T_COSSIN_CP : {
for (int k=2 ; k<np+1 ; k++) {
int m_par = (k/2)%2 ; // 0 pour m pair, 1 pour m impair
somt = 0 ;
for (int j=0 ; j<nt ; j++) {
somt += (*pi) * tab_theta(m_par, j, j0) ;
pi++ ; // theta suivant
}
resu += somt * tab_phi(k, k0) ;
}
break ;
}
default: {
cout << "Mtbl_cf::val_point_jk: unknown theta basis ! " << endl ;
abort () ;
}
} // fin des cas sur base_t
} // fin du cas np > 1
// Menage
delete [] coef_tp ;
// Termine
return resu ;
}
//-------------------------------------------------------------//
// version for xi = 1 //
// and collocation point in (theta',phi') //
//-------------------------------------------------------------//
double Mtbl_cf::val_out_bound_jk(int l, int j0, int k0) const {
#ifndef NDEBUG
// Bases de developpement :
int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
assert((base_r == R_CHEB) || (base_r == R_CHEBU) || (base_r == R_CHEBP)
|| (base_r == R_CHEBI) || (base_r == R_CHEBPIM_P) || (base_r == R_CHEBPIM_I)
|| (base_r == R_CHEBPI_P) || (base_r == R_CHEBPI_I) || (base_r == R_JACO02)) ;
#endif
int nr = mg->get_nr(l) ;
double resu = 0 ;
for (int i=0; i<nr; i++)
resu += operator()(l, k0, j0, i) ;
return resu ;
}
//-------------------------------------------------------------//
// version for xi = -1 //
// and collocation point in (theta',phi') //
//-------------------------------------------------------------//
double Mtbl_cf::val_in_bound_jk(int l, int j0, int k0) const {
#ifndef NDEBUG
// Bases de developpement :
int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
assert((base_r == R_CHEB) || (base_r == R_CHEBU)) ;
#endif
int nr = mg->get_nr(l) ;
double resu = 0 ;
int pari = 1 ;
for (int i=0; i<nr; i++) {
resu += pari*operator()(l, k0, j0, i) ;
pari = - pari ;
}
return resu ;
}
}
Lorene::Mtbl_cf::val_point_jk
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...
Definition mtbl_cf_val_point.C:258
Lorene::Mtbl_cf::base
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
Lorene::Mtbl_cf::operator()
const Tbl & operator()(int l) const
Read-only of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:306
Lorene::Mtbl_cf::mg
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition mtbl_cf.h:192
Lorene::Mtbl_cf::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition mtbl_cf.h:196
Lorene::Mtbl_cf::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Lorene::Mtbl_cf::val_point
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 ...
Definition mtbl_cf_val_point.C:115
Lorene::Mtbl_cf::val_out_bound_jk
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...
Definition mtbl_cf_val_point.C:428
Lorene::Mtbl_cf::val_in_bound_jk
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...
Definition mtbl_cf_val_point.C:452
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
MAX_BASE_2
#define MAX_BASE_2
Smaller maximum bases used for phi (and higher dimensions for now).
Definition type_parite.h:146
TRA_T
#define TRA_T
Translation en Theta, used for a bitwise shift (in hex).
Definition type_parite.h:160
P_COSSIN_P
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition type_parite.h:247
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition type_parite.h:184
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition type_parite.h:180
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition type_parite.h:210
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition type_parite.h:202
MSQ_R
#define MSQ_R
Extraction de l'info sur R.
Definition type_parite.h:152
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition type_parite.h:194
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition type_parite.h:214
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
R_LEG
#define R_LEG
base de Legendre ordinaire (fin)
Definition type_parite.h:182
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition type_parite.h:200
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition type_parite.h:212
P_COSSIN
#define P_COSSIN
dev. standart
Definition type_parite.h:245
P_COSSIN_I
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition type_parite.h:249
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition type_parite.h:176
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex).
Definition type_parite.h:158
MSQ_T
#define MSQ_T
Extraction de l'info sur Theta.
Definition type_parite.h:154
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition type_parite.h:208
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition type_parite.h:206
T_COS
#define T_COS
dev. cos seulement
Definition type_parite.h:196
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
T_SIN
#define T_SIN
dev. sin seulement
Definition type_parite.h:198
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition type_parite.h:204
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition type_parite.h:192
TRA_P
#define TRA_P
Translation en Phi, used for a bitwise shift (in hex).
Definition type_parite.h:162
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:726