map_af_primr.C
/*
* Method Map_af::primr
*
* (see file map.h for documentation).
*
*/
/*
* Copyright (c) 2004 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 version 2
* as published by the Free Software Foundation.
*
* 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 map_af_primr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_primr.C,v 1.7 2014/10/13 08:53:03 j_novak Exp $" ;
/*
* $Id: map_af_primr.C,v 1.7 2014/10/13 08:53:03 j_novak Exp $
* $Log: map_af_primr.C,v $
* Revision 1.7 2014/10/13 08:53:03 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2014/10/06 15:13:12 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.5 2013/04/25 15:46:05 j_novak
* Added special treatment in the case np = 1, for type_p = NONSYM.
*
* Revision 1.4 2007/12/20 09:11:05 jl_cornou
* Correction of an error in op_sxpun about Jacobi(0,2) polynomials
*
* Revision 1.3 2004/07/26 16:02:23 j_novak
* Added a flag to specify whether the primitive should be zero either at r=0
* or at r going to infinity.
*
* Revision 1.1 2004/06/14 15:25:34 e_gourgoulhon
* First version.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Map/map_af_primr.C,v 1.7 2014/10/13 08:53:03 j_novak Exp $
*
*/
// C headers
#include <cstdlib>
// Lorene headers
#include "map.h"
#include "tensor.h"
namespace Lorene {
void _primr_pas_prevu(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void _primr_r_cheb(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void _primr_r_chebp(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void _primr_r_chebi(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void _primr_r_chebpim_p(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void _primr_r_chebpim_i(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void _primr_r_jaco02(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
void Map_af::primr(const Scalar& uu, Scalar& resu, bool null_infty) const {
static void (*prim_domain[MAX_BASE])(const Tbl&, int bin, const Tbl&,
Tbl&, int&, Tbl& ) ;
static bool first_call = true ;
// Initialisation at first call of the array of primitive functions
// depending upon the basis in r
// ----------------------------------------------------------------
if (first_call) {
for (int i=0 ; i<MAX_BASE ; i++) prim_domain[i] = _primr_pas_prevu ;
prim_domain[R_CHEB >> TRA_R] = _primr_r_cheb ;
prim_domain[R_CHEBU >> TRA_R] = _primr_r_cheb ;
prim_domain[R_CHEBP >> TRA_R] = _primr_r_chebp ;
prim_domain[R_CHEBI >> TRA_R] = _primr_r_chebi ;
prim_domain[R_CHEBPIM_P >> TRA_R] = _primr_r_chebpim_p ;
prim_domain[R_CHEBPIM_I >> TRA_R] = _primr_r_chebpim_i ;
prim_domain[R_JACO02 >> TRA_R] = _primr_r_jaco02 ;
first_call = false ;
}
// End of first call operations
// ----------------------------
assert(uu.get_etat() != ETATNONDEF) ;
assert(uu.get_mp().get_mg() == mg) ;
assert(resu.get_mp().get_mg() == mg) ;
// Special case of vanishing input:
if (uu.get_etat() == ETATZERO) {
resu.set_etat_zero() ;
return ;
}
// General case
assert( (uu.get_etat() == ETATQCQ) || (uu.get_etat() == ETATUN) ) ;
assert(uu.check_dzpuis(2)) ;
int nz = mg->get_nzone() ;
int nzm1 = nz - 1 ;
int np = mg->get_np(0) ;
int nt = mg->get_nt(0) ;
#ifndef NDEBUG
for (int l=1; l<nz; l++) {
assert (mg->get_np(l) == np) ;
assert (mg->get_nt(l) == nt) ;
}
#endif
const Valeur& vuu = uu.get_spectral_va() ;
vuu.coef() ;
const Mtbl_cf& cuu = *(vuu.c_cf) ;
assert(cuu.t != 0x0) ;
const Base_val& buu = vuu.get_base() ; // spectral bases of the input
resu.set_etat_qcq() ; // result in ordinary state
Valeur& vprim = resu.set_spectral_va() ;
vprim.set_etat_cf_qcq() ; // allocates the Mtbl_cf for the coefficients
// of the result
Mtbl_cf& cprim = *(vprim.c_cf) ;
cprim.set_etat_qcq() ; // allocates the Tbl's to store the coefficients
// of the result in each domain
Base_val& bprim = cprim.base ; // spectral bases of the result
Tbl val_rmin(np+2,nt) ; // Values of primitive at the left boundary
// in the current domain
Tbl val_rmax(np+2,nt) ; // same but for the right boundary
val_rmin.set_etat_zero() ; // initialisation: primitive = 0 at r=0
int lmax = (mg->get_type_r(nzm1) == UNSURR) ? nz-2 : nzm1 ;
for (int l=0; l<=lmax; l++) {
assert(cuu.t[l] != 0x0) ;
assert(cprim.t[l] != 0x0) ;
const Tbl& cfuu = *(cuu.t[l]) ;
Tbl& cfprim = *(cprim.t[l]) ;
int buu_dom = buu.get_b(l) ;
int base_r = (buu_dom & MSQ_R) >> TRA_R ;
prim_domain[base_r](cfuu, buu_dom, val_rmin, cfprim, bprim.b[l],
val_rmax) ;
cfprim *= alpha[l] ;
val_rmin = alpha[l] * val_rmax / alpha[l+1] ; // for next domain
}
// Special case of compactified external domain (CED)
// --------------------------------------------------
if (mg->get_type_r(nzm1) == UNSURR) {
val_rmin = - val_rmin ;
const Tbl& cfuu = *(cuu.t[nzm1]) ;
Tbl& cfprim = *(cprim.t[nzm1]) ;
int buu_dom = buu.get_b(nzm1) ;
int base_r = (buu_dom & MSQ_R) >> TRA_R ;
assert(base_r == R_CHEBU) ;
prim_domain[base_r](cfuu, buu_dom, val_rmin, cfprim, bprim.b[nzm1],
val_rmax) ;
cfprim *= - alpha[nzm1] ;
}
if (null_infty)
for (int k=0; k<np; k++) //## not very elegant!
for(int j=0; j<nt; j++)
val_rmax.set(k,j) = cprim.val_out_bound_jk(nzm1, j, k) ;
// The output spectral bases (set on the Mtbl_cf) are copied to the Valeur:
vprim.set_base(bprim) ;
if (null_infty)
for (int l=0; l<nz; l++) //## not very elegant!
for (int k=0; k<np; k++)
for(int j=0; j<nt; j++)
for (int i=0; i<mg->get_nr(l); i++)
vprim.set(l, k, j, i) -= val_rmax(k,j) ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::get_b
int get_b(int l) const
Returns the code for the expansion basis in domain no. l.
Definition base_val.h:389
Lorene::Base_val::b
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition base_val.h:331
Lorene::Map_af::primr
virtual void primr(const Scalar &uu, Scalar &resu, bool null_infty) const
Computes the radial primitive which vanishes for .
Definition map_af_primr.C:74
Lorene::Map_af::alpha
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2033
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::base
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
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::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl_cf.C:300
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::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
Lorene::Scalar::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:873
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
Lorene::Valeur::get_base
const Base_val & get_base() const
Return the bases for spectral expansions (member base ).
Definition valeur.h:480
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base ).
Definition valeur.C:810
Lorene::Valeur::set
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
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
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
MSQ_R
#define MSQ_R
Extraction de l'info sur R.
Definition type_parite.h:152
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
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
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene
Lorene prototypes.
Definition app_hor.h:64