scalar_manip.C
/*
* Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 2000-2001 Philippe Grandclement (Cmp version)
*
* 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 scalar_manip_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_manip.C,v 1.19 2014/10/13 08:53:46 j_novak Exp $" ;
/*
* $Id: scalar_manip.C,v 1.19 2014/10/13 08:53:46 j_novak Exp $
* $Log: scalar_manip.C,v $
* Revision 1.19 2014/10/13 08:53:46 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.18 2014/10/06 15:16:15 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.17 2008/10/03 09:03:52 j_novak
* Correction of yet another mistake (the array values in physical space was not
* destroyed).
*
* Revision 1.16 2008/09/30 08:35:18 j_novak
* Correction of forgotten call to coef()
*
* Revision 1.15 2008/09/29 13:23:51 j_novak
* Implementation of the angular mapping associated with an affine
* mapping. Things must be improved to take into account the domain index.
*
* Revision 1.14 2008/09/22 19:08:01 j_novak
* New methods to deal with boundary conditions
*
* Revision 1.13 2007/06/21 20:00:00 k_taniguchi
* Addition of the method filtre_r (int n, int nz)
*
* Revision 1.12 2006/06/28 07:46:41 j_novak
* Better treatment in the case of a domain set to zero.
*
* Revision 1.11 2005/09/07 13:39:10 j_novak
* *** empty log message ***
*
* Revision 1.10 2005/09/07 13:10:48 j_novak
* Added a filter setting to zero all mulitpoles in a given range.
*
* Revision 1.9 2004/11/23 12:47:44 f_limousin
* Add function filtre_tp(int nn, int nz1, int nz2).
*
* Revision 1.8 2004/05/07 11:24:54 f_limousin
* Implement new method filtre_r (int* nn)
*
* Revision 1.7 2004/02/27 09:47:26 f_limousin
* New methods filtre_phi(int) and filtre_theta(int).
*
* Revision 1.6 2004/01/23 13:26:28 e_gourgoulhon
* Added methods set_inner_boundary and set_outer_boundary.
* Methods set_val_inf and set_val_hor, which are particular cases of
* the above, have been suppressed.
*
* Revision 1.5 2003/11/13 13:43:55 p_grandclement
* Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
*
* Revision 1.4 2003/10/11 14:47:17 e_gourgoulhon
* Lines 112 and 145 : replaced "0" by "double(0)" in the comparison
* statement.
*
* Revision 1.3 2003/10/10 15:57:29 j_novak
* Added the state one (ETATUN) to the class Scalar
*
* Revision 1.2 2003/10/08 14:24:09 j_novak
* replaced mult_r_zec with mult_r_ced
*
* Revision 1.1 2003/09/25 09:33:36 j_novak
* Added methods for integral calculation and various manipulations
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_manip.C,v 1.19 2014/10/13 08:53:46 j_novak Exp $
*
*/
//standard
#include <cstdlib>
#include <cmath>
// Lorene
#include "tensor.h"
#include "proto.h"
#include "utilitaires.h"
/*
* Annule tous les l entre l_min et l_max (compris)
*/
namespace Lorene {
void Scalar::annule_l (int l_min, int l_max, bool ylm_output) {
assert (etat != ETATNONDEF) ;
assert (l_min <= l_max) ;
assert (l_min >= 0) ;
if (etat == ETATZERO )
return ;
if (etat == ETATUN) {
if (l_min == 0) set_etat_zero() ;
else return ;
}
va.ylm() ;
Mtbl_cf& m_coef = *va.c_cf ;
const Base_val& base = va.base ;
int l_q, m_q, base_r ;
for (int lz=0; lz<mp->get_mg()->get_nzone(); lz++)
if (m_coef(lz).get_etat() != ETATZERO)
for (int k=0; k<mp->get_mg()->get_np(lz)+1; k++)
for (int j=0; j<mp->get_mg()->get_nt(lz); j++)
for (int i=0; i<mp->get_mg()->get_nr(lz); i++) {
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if ((l_min <= l_q) && (l_q<= l_max))
m_coef.set(lz, k, j, i) = 0 ;
}
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
if (!ylm_output)
va.ylm_i() ;
return ;
}
/*
* Annule les n derniers coefficients en r dans la derniere zone
*/
void Scalar::filtre (int n) {
assert (etat != ETATNONDEF) ;
if ( (etat == ETATZERO) || (etat == ETATUN) )
return ;
int nz = mp->get_mg()->get_nzone() ;
int np = mp->get_mg()->get_np(nz-1) ;
int nt = mp->get_mg()->get_nt(nz-1) ;
int nr = mp->get_mg()->get_nr(nz-1) ;
del_deriv() ;
va.coef() ;
va.set_etat_cf_qcq() ;
for (int k=0 ; k<np+1 ; k++)
if (k!=1)
for (int j=0 ; j<nt ; j++)
for (int i=nr-1 ; i>nr-1-n ; i--)
va.c_cf->set(nz-1, k, j, i) = 0 ;
}
/*
* Annule les n derniers coefficients en r dans toutes les zones
*/
void Scalar::filtre_r (int* nn) {
assert (etat != ETATNONDEF) ;
if ( (etat == ETATZERO) || (etat == ETATUN) )
return ;
del_deriv() ;
va.coef() ;
va.set_etat_cf_qcq() ;
int nz = mp->get_mg()->get_nzone() ;
int* nr = new int[nz];
int* nt = new int[nz];
int* np = new int[nz];
for (int l=0; l<=nz-1; l++) {
nr[l] = mp->get_mg()->get_nr(l) ;
nt[l] = mp->get_mg()->get_nt(l) ;
np[l] = mp->get_mg()->get_np(l) ;
}
for (int l=0; l<=nz-1; l++)
for (int k=0 ; k<np[l]+1 ; k++)
if (k!=1)
for (int j=0 ; j<nt[l] ; j++)
for (int i=nr[l]-1; i>nr[l]-1-nn[l] ; i--)
va.c_cf->set(l, k, j, i) = 0 ;
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
}
/*
* Annule les n derniers coefficients en r dans zone nz
*/
void Scalar::filtre_r (int n, int nz) {
assert (etat != ETATNONDEF) ;
if ( (etat == ETATZERO) || (etat == ETATUN) )
return ;
del_deriv() ;
va.coef() ;
va.set_etat_cf_qcq() ;
int nr = mp->get_mg()->get_nr(nz) ;
int nt = mp->get_mg()->get_nt(nz) ;
int np = mp->get_mg()->get_np(nz) ;
for (int k=0 ; k<np+1 ; k++)
if (k!=1)
for (int j=0 ; j<nt ; j++)
for (int i=nr-1; i>nr-1-n ; i--)
va.c_cf->set(nz, k, j, i) = 0 ;
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
}
/*
* Annule les n derniers coefficients en phi dans zone nz
*/
void Scalar::filtre_phi (int n, int nz) {
assert (etat != ETATNONDEF) ;
if ( (etat == ETATZERO) || (etat == ETATUN) )
return ;
del_deriv() ;
va.coef() ;
va.set_etat_cf_qcq() ;
int np = mp->get_mg()->get_np(nz) ;
int nt = mp->get_mg()->get_nt(nz) ;
int nr = mp->get_mg()->get_nr(nz) ;
for (int k=np+1-n ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++)
for (int i=0 ; i<nr ; i++)
va.c_cf->set(nz, k, j, i) = 0 ;
}
void Scalar::filtre_tp(int nn, int nz1, int nz2) {
va.filtre_tp(nn, nz1, nz2) ;
}
/* Sets the value of the {\tt Scalar} at the inner boundary of a given
* domain.
* @param l [input] domain index
* @param x [input] (constant) value at the inner boundary of domain no. {\tt l}
*/
void Scalar::set_inner_boundary(int l0, double x0) {
assert (etat != ETATNONDEF) ;
if (etat == ETATZERO) {
if (x0 == double(0)) return ;
else annule_hard() ;
}
if (etat == ETATUN) {
if (x0 == double(1)) return ;
else etat = ETATQCQ ;
}
del_deriv() ;
int nt = mp->get_mg()->get_nt(l0) ;
int np = mp->get_mg()->get_np(l0) ;
va.coef_i() ;
va.set_etat_c_qcq() ;
for (int k=0 ; k<np ; k++)
for (int j=0 ; j<nt ; j++)
va.set(l0, k, j, 0) = x0 ;
}
/* Sets the value of the {\tt Scalar} at the outer boundary of a given
* domain.
* @param l [input] domain index
* @param x [input] (constant) value at the outer boundary of domain no. {\tt l}
*/
void Scalar::set_outer_boundary(int l0, double x0) {
assert (etat != ETATNONDEF) ;
if (etat == ETATZERO) {
if (x0 == double(0)) return ;
else annule_hard() ;
}
if (etat == ETATUN) {
if (x0 == double(1)) return ;
else etat = ETATQCQ ;
}
del_deriv() ;
int nrm1 = mp->get_mg()->get_nr(l0) - 1 ;
int nt = mp->get_mg()->get_nt(l0) ;
int np = mp->get_mg()->get_np(l0) ;
va.coef_i() ;
va.set_etat_c_qcq() ;
for (int k=0 ; k<np ; k++)
for (int j=0 ; j<nt ; j++)
va.set(l0, k, j, nrm1) = x0 ;
}
/*
* Permet de fixer la decroissance du cmp a l infini en viurant les
* termes en 1/r^n
*/
void Scalar::fixe_decroissance (int puis) {
if (puis<dzpuis)
return ;
else {
int nbre = puis-dzpuis ;
// le confort avant tout ! (c'est bien le confort ...)
int nz = mp->get_mg()->get_nzone() ;
int np = mp->get_mg()->get_np(nz-1) ;
int nt = mp->get_mg()->get_nt(nz-1) ;
int nr = mp->get_mg()->get_nr(nz-1) ;
const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
if (map == 0x0) {
cout << "Le mapping doit etre affine" << endl ;
abort() ;
}
double alpha = map->get_alpha()[nz-1] ;
Scalar courant (*this) ;
va.coef() ;
va.set_etat_cf_qcq() ;
for (int conte=0 ; conte<nbre ; conte++) {
int base_r = courant.va.base.get_base_r(nz-1) ;
courant.va.coef() ;
// On calcul les coefficients de 1/r^conte
double* coloc = new double [nr] ;
int * deg = new int[3] ;
deg[0] = 1 ;
deg[1] = 1 ;
deg[2] = nr ;
for (int i=0 ; i<nr ; i++)
coloc[i] =pow(alpha, double(conte))*
pow(-1-cos(M_PI*i/(nr-1)), double(conte)) ;
cfrcheb(deg, deg, coloc, deg, coloc) ;
for (int k=0 ; k<np+1 ; k++)
if (k != 1)
for (int j=0 ; j<nt ; j++) {
// On doit determiner le coefficient du truc courant :
double* coef = new double [nr] ;
double* auxi = new double[1] ;
for (int i=0 ; i<nr ; i++)
coef[i] = (*courant.va.c_cf)(nz-1, k, j, i) ;
switch (base_r) {
case R_CHEBU :
som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
break ;
default :
som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
break ;
}
// On modifie le cmp courant :
courant.va.coef() ;
courant.va.set_etat_cf_qcq() ;
courant.va.c_cf->set(nz-1, k, j, 0) -= *auxi ;
for (int i=0 ; i<nr ; i++)
this->va.c_cf->set(nz-1, k, j, i) -= *auxi * coloc[i] ;
delete [] coef ;
delete [] auxi ;
}
delete [] coloc ;
delete [] deg ;
courant.mult_r_ced() ;
}
}
}
Tbl Scalar::tbl_out_bound(int l_zone, bool output_ylm) {
va.coef() ;
if (output_ylm) va.ylm() ;
int np = mp->get_mg()->get_np(l_zone) ;
int nt = mp->get_mg()->get_nt(l_zone) ;
Tbl resu(np+2, nt) ;
if (etat == ETATZERO) resu.set_etat_zero() ;
else {
assert(etat == ETATQCQ) ;
resu.set_etat_qcq() ;
for (int k=0; k<np+2; k++)
for (int j=0; j<nt; j++)
resu.set(k, j) = va.c_cf->val_out_bound_jk(l_zone, j, k) ;
}
return resu ;
}
Tbl Scalar::tbl_in_bound(int l_zone, bool output_ylm) {
assert(mp->get_mg()->get_type_r(l_zone) != RARE) ;
va.coef() ;
if (output_ylm) va.ylm() ;
int np = mp->get_mg()->get_np(l_zone) ;
int nt = mp->get_mg()->get_nt(l_zone) ;
Tbl resu(np+2, nt) ;
if (etat == ETATZERO) resu.set_etat_zero() ;
else {
assert(etat == ETATQCQ) ;
resu.set_etat_qcq() ;
for (int k=0; k<np+2; k++)
for (int j=0; j<nt; j++)
resu.set(k, j) = va.c_cf->val_in_bound_jk(l_zone, j, k) ;
}
return resu ;
}
Scalar Scalar::scalar_out_bound(int l_zone, bool output_ylm) {
va.coef() ;
if (output_ylm) va.ylm() ;
Scalar resu(mp->mp_angu(l_zone)) ;
resu.std_spectral_base() ;
Base_val base = resu.get_spectral_base() ;
base.set_base_t(va.base.get_base_t(l_zone)) ;
resu.set_spectral_base(base) ;
if (etat == ETATZERO) resu.set_etat_zero() ;
else {
assert(etat == ETATQCQ) ;
resu.annule_hard() ;
int np = mp->get_mg()->get_np(l_zone) ;
int nt = mp->get_mg()->get_nt(l_zone) ;
for (int k=0; k<np+2; k++)
for (int j=0; j<nt; j++)
resu.set_spectral_va().c_cf->set(0, k, j, 0)
= va.c_cf->val_out_bound_jk(l_zone, j, k) ;
delete resu.set_spectral_va().c ;
resu.set_spectral_va().c = 0x0 ;
}
return resu ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::set_base_t
void set_base_t(int base_t)
Sets the expansion basis for functions in all domains.
Definition base_val.C:195
Lorene::Base_val::get_base_r
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition base_val.h:400
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition base_val_quantum.C:81
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Lorene::Scalar::annule_l
void annule_l(int l_min, int l_max, bool ylm_output=false)
Sets all the multipolar components between l_min and l_max to zero.
Definition scalar_manip.C:111
Lorene::Scalar::del_deriv
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition scalar.C:287
Lorene::Scalar::scalar_out_bound
Scalar scalar_out_bound(int n, bool leave_ylm=false)
Returns the Scalar containing the values of angular coefficients at the outer boundary.
Definition scalar_manip.C:467
Lorene::Scalar::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition scalar_manip.C:151
Lorene::Scalar::fixe_decroissance
void fixe_decroissance(int puis)
Substracts all the components behaving like in the external domain, with n strictly lower than puis ...
Definition scalar_manip.C:346
Lorene::Scalar::tbl_in_bound
Tbl tbl_in_bound(int n, bool leave_ylm=false)
Returns the Tbl containing the values of angular coefficients at the inner boundary.
Definition scalar_manip.C:448
Lorene::Scalar::cos
friend Scalar cos(const Scalar &)
Cosine.
Definition scalar_math.C:104
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
Lorene::Scalar::tbl_out_bound
Tbl tbl_out_bound(int l_dom, bool leave_ylm=false)
Returns the Tbl containing the values of angular coefficients at the outer boundary.
Definition scalar_manip.C:430
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::Scalar
Scalar(const Map &mpi)
Constructor from mapping.
Definition scalar.C:204
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::mult_r_ced
void mult_r_ced()
Multiplication by r in the compactified external domain (CED), the dzpuis flag is not changed.
Definition scalar_r_manip.C:269
Lorene::Scalar::filtre_tp
void filtre_tp(int nn, int nz1, int nz2)
Sets the n lasts coefficients in to 0 in the domain nz1 to nz2 when expressed in spherical harmonics...
Definition scalar_manip.C:270
Lorene::Scalar::etat
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition scalar.h:1294
Lorene::Scalar::set_inner_boundary
void set_inner_boundary(int l, double x)
Sets the value of the Scalar at the inner boundary of a given domain.
Definition scalar_manip.C:283
Lorene::Scalar::filtre_r
void filtre_r(int *nn)
Sets the n lasts coefficients in r to 0 in all domains.
Definition scalar_manip.C:180
Lorene::Scalar::filtre_phi
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
Definition scalar_manip.C:249
Lorene::Scalar::set_outer_boundary
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition scalar_manip.C:315
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:405
Lorene::Scalar::pow
friend Scalar pow(const Scalar &, int)
Power .
Definition scalar_math.C:451
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition scalar.C:797
Lorene::Scalar::dzpuis
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external d...
Definition scalar.h:403
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_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
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::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:299
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
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition type_parite.h:180
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene
Lorene prototypes.
Definition app_hor.h:64