valeur_coef_i.C
/*
* Computations of the values at the collocation points from the spectral
* coefficients
*
*/
/*
* Copyright (c) 1999-2003 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 valeur_coef_i_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_coef_i.C,v 1.17 2014/10/13 08:53:49 j_novak Exp $" ;
/*
* $Id: valeur_coef_i.C,v 1.17 2014/10/13 08:53:49 j_novak Exp $
* $Log: valeur_coef_i.C,v $
* Revision 1.17 2014/10/13 08:53:49 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.16 2014/10/06 15:13:22 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.15 2013/06/06 15:31:33 j_novak
* Functions to compute Legendre coefficients (not fully tested yet).
*
* Revision 1.14 2012/01/17 15:08:02 j_penner
* using MAX_BASE_2 for the phi coordinate
*
* Revision 1.13 2009/10/23 12:56:29 j_novak
* New base T_LEG_MI
*
* Revision 1.12 2009/10/13 13:49:58 j_novak
* New base T_LEG_MP.
*
* Revision 1.11 2009/10/08 16:23:14 j_novak
* Addition of new bases T_COS and T_SIN.
*
* Revision 1.10 2008/10/07 15:01:58 j_novak
* The case nt=1 is now treated separately.
*
* Revision 1.9 2007/12/11 15:28:25 jl_cornou
* Jacobi(0,2) polynomials partially implemented
*
* Revision 1.8 2004/11/23 15:17:19 m_forot
* Added the bases for the cases without any equatorial symmetry
* (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
*
* Revision 1.7 2004/08/24 09:14:52 p_grandclement
* Addition of some new operators, like Poisson in 2d... It now requieres the
* GSL library to work.
*
* Also, the way a variable change is stored by a Param_elliptic is changed and
* no longer uses Change_var but rather 2 Scalars. The codes using that feature
* will requiere some modification. (It should concern only the ones about monopoles)
*
* Revision 1.6 2003/10/13 20:51:25 e_gourgoulhon
* Replaced malloc by new
*
* Revision 1.5 2003/09/17 12:30:22 j_novak
* New checks for changing to T_LEG* bases.
*
* Revision 1.4 2003/09/16 13:07:41 j_novak
* New files for coefficient trnasformation to/from the T_LEG_II base.
*
* Revision 1.3 2002/11/12 17:44:35 j_novak
* Added transformation function for T_COS basis.
*
* Revision 1.2 2002/10/16 14:37:15 j_novak
* Reorganization of #include instructions of standard C++, in order to
* use experimental version 3 of gcc.
*
* Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
* LORENE
*
* Revision 2.9 2000/10/04 14:41:40 eric
* Ajout des bases T_LEG_IP et T_LEG_PI
*
* Revision 2.8 2000/09/07 15:14:44 eric
* Ajout de la base P_COSSIN_I
*
* Revision 2.7 2000/08/16 10:33:04 eric
* Suppression de Mtbl::dzpuis.
*
* Revision 2.6 1999/12/16 16:41:48 phil
* *** empty log message ***
*
* Revision 2.5 1999/11/30 12:43:03 eric
* Valeur::base est desormais du type Base_val et non plus Base_val*.
*
* Revision 2.4 1999/10/28 07:44:20 eric
* Modif commentaires.
*
* Revision 2.3 1999/10/13 15:51:33 eric
* Anglisation.
*
* Revision 2.2 1999/06/22 15:10:02 phil
* ajout de dzpuis
*
* Revision 2.1 1999/02/22 15:40:25 hyc
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_coef_i.C,v 1.17 2014/10/13 08:53:49 j_novak Exp $
*
*/
#include <cmath>
// Header Lorene
#include "valeur.h"
#include "proto.h"
namespace Lorene {
void c_est_pas_fait(char * ) ;
void ipasprevu_r(const int*, const int*, double*, const int*, double*) ;
void ipasprevu_t(const int*, const int*, double*, const int*, double*) ;
void ipasprevu_p(const int* , const int* , const int* , double* , double* ) ;
void ibase_non_def_r(const int*, const int*, double*, const int*, double*) ;
void ibase_non_def_t(const int*, const int*, double*, const int*, double*) ;
void ibase_non_def_p(const int* , const int* , const int* , double* , double* ) ;
void Valeur::coef_i() const {
// Variables statiques
static void (*invcf_r[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
static void (*invcf_t[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
static void (*invcf_p[MAX_BASE_2])(const int* , const int* , const int*, double* , double* ) ;
static int premier_appel = 1 ;
// Premier appel
if (premier_appel) {
premier_appel = 0 ;
for (int i=0; i<MAX_BASE; i++) {
invcf_r[i] = ipasprevu_r ;
invcf_t[i] = ipasprevu_t ;
if(i%2==0){
invcf_p[i/2] = ipasprevu_p ; // saves a loop
}
}
invcf_r[NONDEF] = ibase_non_def_r ;
invcf_r[R_CHEB >> TRA_R] = circheb ;
invcf_r[R_CHEBU >> TRA_R] = circheb ;
invcf_r[R_CHEBP >> TRA_R] = circhebp ;
invcf_r[R_CHEBI >> TRA_R] = circhebi ;
invcf_r[R_CHEBPIM_P >> TRA_R] = circhebpimp ;
invcf_r[R_CHEBPIM_I >> TRA_R] = circhebpimi ;
invcf_r[R_CHEBPI_P >> TRA_R] = circhebpip ;
invcf_r[R_CHEBPI_I >> TRA_R] = circhebpii ;
invcf_r[R_LEG >> TRA_R] = cirleg ;
invcf_r[R_LEGP >> TRA_R] = cirlegp ;
invcf_r[R_LEGI >> TRA_R] = cirlegi ;
invcf_r[R_JACO02 >> TRA_R] = cirjaco02 ;
invcf_t[NONDEF] = ibase_non_def_t ;
invcf_t[T_COS >> TRA_T] = citcos ;
invcf_t[T_SIN >> TRA_T] = citsin ;
invcf_t[T_COS_P >> TRA_T] = citcosp ;
invcf_t[T_COS_I >> TRA_T] = citcosi ;
invcf_t[T_SIN_P >> TRA_T] = citsinp ;
invcf_t[T_SIN_I >> TRA_T] = citsini ;
invcf_t[T_COSSIN_CP >> TRA_T] = citcossincp ;
invcf_t[T_COSSIN_SI >> TRA_T] = citcossinsi ;
invcf_t[T_COSSIN_SP >> TRA_T] = citcossinsp ;
invcf_t[T_COSSIN_CI >> TRA_T] = citcossinci ;
invcf_t[T_COSSIN_S >> TRA_T] = citcossins ;
invcf_t[T_COSSIN_C >> TRA_T] = citcossinc ;
invcf_t[T_LEG_P >> TRA_T] = citlegp ;
invcf_t[T_LEG_PP >> TRA_T] = citlegpp ;
invcf_t[T_LEG_I >> TRA_T] = citlegi ;
invcf_t[T_LEG_IP >> TRA_T] = citlegip ;
invcf_t[T_LEG_PI >> TRA_T] = citlegpi ;
invcf_t[T_LEG_II >> TRA_T] = citlegii ;
invcf_t[T_LEG_MP >> TRA_T] = citlegmp ;
invcf_t[T_LEG_MI >> TRA_T] = citlegmi ;
invcf_t[T_LEG >> TRA_T] = citleg ;
invcf_p[NONDEF] = ibase_non_def_p ;
invcf_p[P_COSSIN >> TRA_P] = cipcossin ;
invcf_p[P_COSSIN_P >> TRA_P] = cipcossin ;
invcf_p[P_COSSIN_I >> TRA_P] = cipcossini ;
} // fin des operation de premier appel
//------------------//
// DEBUT DU CALCUL //
//------------------//
// Tout null ?
if (etat == ETATZERO) {
return ;
}
// Protection
assert(etat != ETATNONDEF) ;
// Peut-etre rien a faire
if (c != 0x0) {
return ;
}
// Il faut bosser
assert(c_cf != 0x0) ; // ..si on peut
assert(c_cf->base == base) ; // Consistence des bases
c = new Mtbl(mg) ;
c->set_etat_qcq() ;
// Boucles sur les zones
int nz = mg->get_nzone() ;
for (int l=0; l<nz; l++) {
// Initialisation des valeurs de this->c_cf avec celle de this->c :
Tbl* f = (c->t)[l] ;
const Tbl* cf = (c_cf->t)[l] ;
if (cf->get_etat() == ETATZERO) {
f->set_etat_zero() ;
continue ; // on ne fait rien si le tbl(cf) = 0
}
f->set_etat_qcq() ;
int np = f->get_dim(2) ;
int nt = f->get_dim(1) ;
int nr = f->get_dim(0) ;
int np_c = cf->get_dim(2) ;
int nt_c = cf->get_dim(1) ;
int nr_c = cf->get_dim(0) ;
// Attention a ce qui suit... (deg et dim)
int deg[3] ;
deg[0] = np ;
deg[1] = nt ;
deg[2] = nr ;
int dimc[3] ;
dimc[0] = np_c ;
dimc[1] = nt_c ;
dimc[2] = nr_c ;
// Allocation de l'espace memoire pour le tableau de travail trav
int ntot = cf->get_taille() ;
double* trav = new double[ntot] ;
// On recupere les bases en r, theta et phi :
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 ;
int vbase_t = base.b[l] & MSQ_T ;
int vbase_p = base.b[l] & MSQ_P ;
assert(base_r < MAX_BASE) ;
assert(base_t < MAX_BASE) ;
assert(base_p < MAX_BASE_2) ;
// Transformation inverse en r:
if ( nr == 1 ) {
for (int i=0; i<ntot; i++) {
trav[i] = cf->t[i] ; // simple recopie cf --> trav
}
}
else {
invcf_r[base_r]( deg, dimc, (cf->t), dimc, trav ) ;
}
// Partie angulaire
if ( np == 1) {
if (nt==1) {
for (int i=0 ; i<f->get_taille() ; i++)
f->t[i] = trav[i] ;
if ((vbase_t == T_LEG_PP) || (vbase_t == T_LEG_PI) ||
(vbase_t == T_LEG_IP) || (vbase_t == T_LEG_II) ||
(vbase_t == T_LEG_P) || (vbase_t == T_LEG_I) ||
(vbase_t == T_LEG_MP) || (vbase_t == T_LEG_MI) ||
(vbase_t == T_LEG) ) {
*f /=sqrt(2.) ;
}
}
else {
bool pair = ( (vbase_t == T_LEG_PP) || (vbase_t == T_LEG_IP)
|| (vbase_t == T_LEG_MP) ) ;
bool impair = ( (vbase_t == T_LEG_PI) || (vbase_t == T_LEG_II)
|| (vbase_t == T_LEG_MI)) ;
if ((pair && (vbase_p == P_COSSIN_I)) ||
(impair && (vbase_p == P_COSSIN_P)) )
ipasprevu_t(deg, dimc, trav, deg, (f->t) ) ;
else
invcf_t[base_t]( deg, dimc, trav, deg, (f->t) ) ;
}
}
else {
// Cas 3-D
// ...... Transformation inverse en theta:
bool pair = ( (vbase_t == T_LEG_PP) || (vbase_t == T_LEG_IP)
|| (vbase_t == T_LEG_MP) ) ;
bool impair = ( (vbase_t == T_LEG_PI) || (vbase_t == T_LEG_II)
|| (vbase_t == T_LEG_MI) ) ;
if ((pair && (vbase_p == P_COSSIN_I)) ||
(impair && (vbase_p == P_COSSIN_P)) )
ipasprevu_t(deg, dimc, trav, dimc, trav ) ;
else
invcf_t[base_t]( deg, dimc, trav, dimc, trav ) ;
// ...... Transformation inverse en phi:
invcf_p[base_p]( deg, dimc, deg, trav, (f->t) ) ;
}
// Menage
delete [] trav ;
} // fin de la boucle sur les differentes zones
}
//------------------------//
// Les machins pas prevus //
//------------------------//
void ipasprevu_r(const int*, const int*, double*, const int*, double*) {
cout << "Valeur::coef_i: the required expansion basis in r " << endl ;
cout << " is not implemented !" << endl ;
abort() ;
}
void ipasprevu_t(const int*, const int*, double*, const int*, double* ) {
cout << "Valeur::coef_i: the required expansion basis in theta " << endl ;
cout << " is not implemented !" << endl ;
abort() ;
}
void ipasprevu_p(const int*, const int*, const int*, double*, double* ) {
cout << "Valeur::coef_i: the required expansion basis in phi " << endl ;
cout << " is not implemented !" << endl ;
abort() ;
}
void ibase_non_def_r(const int*, const int*, double*, const int*, double*) {
cout << "Valeur::coef_i: the expansion basis in r is undefined !" << endl ;
abort() ;
}
void ibase_non_def_t(const int*, const int*, double*, const int*, double*) {
cout << "Valeur::coef_i: the expansion basis in theta is undefined !"
<< endl ;
abort() ;
}
void ibase_non_def_p(const int*, const int*, const int*, double*, double*) {
cout << "Valeur::coef_i: the expansion basis in phi is undefined !" << endl ;
abort() ;
}
}
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
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::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i]).
Definition tbl.h:403
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:299
Lorene::Valeur::mg
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition valeur_coef_i.C:140
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Lorene::Valeur::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
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
T_LEG_MP
#define T_LEG_MP
fct. de Legendre associees avec m pair
Definition type_parite.h:238
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
T_LEG_PI
#define T_LEG_PI
fct. de Legendre associees paires avec m impair
Definition type_parite.h:224
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_LEG
#define T_LEG
fct. de Legendre associees
Definition type_parite.h:236
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_LEG_P
#define T_LEG_P
fct. de Legendre associees paires
Definition type_parite.h:216
T_LEG_IP
#define T_LEG_IP
fct. de Legendre associees impaires avec m pair
Definition type_parite.h:222
NONDEF
#define NONDEF
base inconnue
Definition type_parite.h:150
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
T_LEG_MI
#define T_LEG_MI
fct. de Legendre associees avec m impair
Definition type_parite.h:240
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
T_LEG_II
#define T_LEG_II
fct. de Legendre associees impaires avec m impair
Definition type_parite.h:226
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
T_LEG_I
#define T_LEG_I
fct. de Legendre associees impaires
Definition type_parite.h:220
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_LEG_PP
#define T_LEG_PP
fct. de Legendre associees paires avec m pair
Definition type_parite.h:218
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::c_est_pas_fait
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Definition c_est_pas_fait.C:73
Lorene
Lorene prototypes.
Definition app_hor.h:64