scalar.C
/*
* Methods of class Scalar
*
* (see file scalar.h for documentation)
*
*/
/*
* Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Cmp)
*
* 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_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar.C,v 1.40 2015/12/18 15:52:52 j_novak Exp $" ;
/*
* $Id: scalar.C,v 1.40 2015/12/18 15:52:52 j_novak Exp $
* $Log: scalar.C,v $
* Revision 1.40 2015/12/18 15:52:52 j_novak
* New method is_nan() for class Scalar.
*
* Revision 1.39 2014/10/13 08:53:45 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.38 2014/10/06 15:16:15 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.37 2013/06/05 15:06:11 j_novak
* Legendre bases are treated as standard bases, when the multi-grid
* (Mg3d) is built with BASE_LEG.
*
* Revision 1.36 2013/03/27 09:51:42 e_gourgoulhon
* Restored log
*
*
* revision 1.35
* date: 2013/01/11 15:44:54; author: j_novak; state: Exp; lines: +15 -4
* Addition of Legendre bases (part 2).
*
* revision 1.34
* date: 2012/01/17 15:10:12; author: j_penner; state: Exp; lines: +2 -7
*
* revision 1.33
* date: 2012/01/17 10:26:48; author: j_penner; state: Exp; lines: +11 -3
* added a derivative with respect to the computational coordinate xi
*
* Revision 1.32 2011/04/08 12:55:50 e_gourgoulhon
* Scalar::val_point : division by r^dzpuis to return the actual
* (i.e. physical) value of the field in the external compactified
* domain.
*
* Revision 1.31 2005/10/25 08:56:39 p_grandclement
* addition of std_spectral_base in the case of odd functions near the origin
*
* Revision 1.30 2005/03/11 13:16:37 j_novak
* Corrected an error in multipole_spectrum().
*
* Revision 1.29 2004/10/11 15:09:04 j_novak
* The radial manipulation functions take Scalar as arguments, instead of Cmp.
* Added a conversion operator from Scalar to Cmp.
* The Cmp radial manipulation function make conversion to Scalar, call to the
* Map_radial version with a Scalar argument and back.
*
* Revision 1.28 2004/08/24 09:14:51 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.27 2004/06/22 08:50:00 p_grandclement
* Addition of everything needed for using the logarithmic mapping
*
* Revision 1.26 2004/06/11 14:29:58 j_novak
* Scalar::multipole_spectrum() is now a const method.
*
* Revision 1.25 2004/03/04 09:51:36 e_gourgoulhon
* Method dz_nonzero() : treated the case ETATUN.
*
* Revision 1.24 2004/02/19 22:11:50 e_gourgoulhon
* Added argument "comment" in method spectral_display.
*
* Revision 1.23 2004/01/12 15:32:25 j_novak
* Yet another problem with ETATUN
*
* Revision 1.22 2003/11/05 15:31:13 e_gourgoulhon
* Method set_etat_qcq(): del_t() is not invoqued for etat == ETATUN.
*
* Revision 1.21 2003/11/03 10:25:27 j_novak
* Suppressed the call to del_deriv() in set_etat_qcq() method.
*
* Revision 1.20 2003/10/28 21:33:13 e_gourgoulhon
* In operator=(const Scalar& ci) changed the place of va.del_t()
* in order to correct an error in the case where both this and ci
* are in the state ETATUN.
*
* Revision 1.19 2003/10/28 12:36:10 e_gourgoulhon
* Corrected operator=(const Scalar& ci) for the case ETATUN: the
* spectral bases are now copied.
*
* Revision 1.18 2003/10/27 14:51:25 e_gourgoulhon
* Correction of errors in constructor from a Tensor.
*
* Revision 1.17 2003/10/22 08:35:30 j_novak
* Error corrected
*
* Revision 1.16 2003/10/19 19:54:37 e_gourgoulhon
* -- Modified method spectral_display: now calling Valeur::display_coef.
*
* Revision 1.15 2003/10/15 16:03:38 j_novak
* Added the angular Laplace operator for Scalar.
*
* Revision 1.14 2003/10/15 10:43:06 e_gourgoulhon
* Added new members p_dsdt and p_stdsdp.
*
* Revision 1.13 2003/10/13 13:52:40 j_novak
* Better managment of derived quantities.
*
* Revision 1.12 2003/10/11 14:42:16 e_gourgoulhon
* Suppressed unusued argument new_triad in method change_triad.
*
* Revision 1.11 2003/10/10 15:57:29 j_novak
* Added the state one (ETATUN) to the class Scalar
*
* Revision 1.10 2003/10/07 08:05:03 j_novak
* Added an assert for the constructor from a Tensor.
*
* Revision 1.9 2003/10/06 16:16:03 j_novak
* New constructor from a Tensor.
*
* Revision 1.8 2003/10/06 13:58:48 j_novak
* The memory management has been improved.
* Implementation of the covariant derivative with respect to the exact Tensor
* type.
*
* Revision 1.7 2003/10/05 21:15:42 e_gourgoulhon
* - Suppressed method std_spectral_base_scal().
* - Added method std_spectral_base().
*
* Revision 1.6 2003/10/01 13:04:44 e_gourgoulhon
* The method Tensor::get_mp() returns now a reference (and not
* a pointer) onto a mapping.
*
* Revision 1.5 2003/09/29 12:52:58 j_novak
* Methods for changing the triad are implemented.
*
* Revision 1.4 2003/09/24 20:55:27 e_gourgoulhon
* Added -- constructor by conversion of a Cmp
* -- operator=(const Cmp&)
*
* Revision 1.3 2003/09/24 15:10:54 j_novak
* Suppression of the etat flag in class Tensor (still present in Scalar)
*
* Revision 1.2 2003/09/24 10:21:07 e_gourgoulhon
* added more methods
*
* Revision 1.1 2003/09/23 08:52:17 e_gourgoulhon
* new version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar.C,v 1.40 2015/12/18 15:52:52 j_novak Exp $
*
*/
// headers C
#include <cassert>
#include <cstdlib>
#include <cmath>
// headers Lorene
#include "tensor.h"
#include "type_parite.h"
#include "utilitaires.h"
#include "proto.h"
#include "cmp.h"
//---------------//
// Constructors //
//---------------//
namespace Lorene {
Scalar::Scalar(const Map& mpi) : Tensor(mpi), etat(ETATNONDEF), dzpuis(0),
va(mpi.get_mg()) {
cmp[0] = this ;
set_der_0x0() ;
}
// Constructor from a Tensor
// -------------------------
Scalar::Scalar(const Tensor& ti) : Tensor(*(ti.mp)), etat(ti.cmp[0]->etat),
dzpuis(ti.cmp[0]->dzpuis), va(ti.cmp[0]->va) {
assert(ti.valence == 0) ;
cmp[0] = this ;
set_der_0x0() ;
}
// Copy constructor
// ----------------
Scalar::Scalar(const Scalar& sci) : Tensor(*(sci.mp)), etat(sci.etat),
dzpuis(sci.dzpuis), va(sci.va) {
cmp[0] = this ;
set_der_0x0() ; // On ne recopie pas les derivees
}
// Conversion from a Cmp
//----------------------
Scalar::Scalar(const Cmp& ci) : Tensor(*(ci.get_mp())),
etat(ci.get_etat()),
dzpuis(ci.get_dzpuis()),
va(ci.va) {
cmp[0] = this ;
set_der_0x0() ;
}
// From file
// ---------
Scalar::Scalar(const Map& mpi, const Mg3d& mgi, FILE* fd) : Tensor(mpi),
va(mgi, fd) {
assert( mpi.get_mg() == &mgi ) ;
fread_be(&etat, sizeof(int), 1, fd) ; // L'etat
fread_be(&dzpuis, sizeof(int), 1, fd) ; // dzpuis
cmp[0] = this ;
set_der_0x0() ; // Les derivees sont initialisees a zero
}
//--------------//
// Destructor //
//--------------//
// Destructeur
Scalar::~Scalar() {
del_t() ;
cmp[0] = 0x0 ; //cmp[0] was set to 'this' before (now 0x0 to break a
// in loop in ~Tensor!)
}
//-----------------------//
// Memory management //
//-----------------------//
// Destructeur logique
void Scalar::del_t() {
va.del_t() ;
va.set_etat_nondef() ;
Scalar::del_deriv() ;
}
void Scalar::del_deriv() const{
if (p_dsdr != 0x0) delete p_dsdr ;
if (p_srdsdt != 0x0) delete p_srdsdt ;
if (p_srstdsdp != 0x0) delete p_srstdsdp ;
if (p_dsdt != 0x0) delete p_dsdt ;
if (p_stdsdp != 0x0) delete p_stdsdp ;
if (p_dsdx != 0x0) delete p_dsdx ;
if (p_dsdy != 0x0) delete p_dsdy ;
if (p_dsdz != 0x0) delete p_dsdz ;
if (p_lap != 0x0) delete p_lap ;
if (p_lapang != 0x0) delete p_lapang ;
if (p_integ != 0x0) delete p_integ ;
if (p_dsdradial != 0x0) delete p_dsdradial ;
if (p_dsdrho != 0x0) delete p_dsdrho ;
set_der_0x0() ;
Tensor::del_deriv() ;
}
void Scalar::set_der_0x0() const {
p_dsdr = 0x0 ;
p_srdsdt = 0x0 ;
p_srstdsdp = 0x0 ;
p_dsdt = 0x0 ;
p_stdsdp = 0x0 ;
p_dsdx = 0x0 ;
p_dsdy = 0x0 ;
p_dsdz = 0x0 ;
p_lap = 0x0 ;
p_lapang = 0x0 ;
ind_lap = - 1 ;
p_integ = 0x0 ;
p_dsdradial = 0x0 ;
p_dsdrho = 0x0 ;
}
// ETATZERO
void Scalar::set_etat_zero() {
if (etat == ETATZERO) return ;
else {
del_deriv() ;
va.set_etat_zero() ;
etat = ETATZERO ;
}
}
// ETATUN
void Scalar::set_etat_one() {
if (etat == ETATUN) return ;
else {
del_deriv() ;
va = 1 ;
etat = ETATUN ;
}
}
// ETATNONDEF
void Scalar::set_etat_nondef() {
if (etat == ETATNONDEF) return ;
else {
del_t() ;
etat = ETATNONDEF ;
}
}
// ETATQCQ
void Scalar::set_etat_qcq() {
if (etat == ETATQCQ) return ;
else {
if (etat != ETATUN) del_t() ;
etat = ETATQCQ ;
}
}
// Allocates everything
// --------------------
void Scalar::allocate_all() {
set_etat_qcq() ;
va.set_etat_c_qcq() ; // allocation in configuration space
Mtbl* mt = va.c ;
mt->set_etat_qcq() ;
for (int l=0; l<mt->get_nzone(); l++) {
mt->t[l]->set_etat_qcq() ;
}
}
// ZERO hard
void Scalar::annule_hard() {
va.annule_hard() ;
del_deriv() ;
etat = ETATQCQ ;
}
// Sets the Scalar to zero in several domains
// ---------------------------------------
void Scalar::annule(int l_min, int l_max) {
// Cas particulier: annulation globale :
if ( (l_min == 0) && (l_max == va.mg->get_nzone()-1) ) {
set_etat_zero() ;
return ;
}
assert( etat != ETATNONDEF ) ;
if ( etat == ETATZERO ) {
return ; // rien n'a faire si c'est deja zero
}
else {
assert( (etat == ETATQCQ) || (etat == ETATUN) ) ; // sinon...
etat = ETATQCQ ;
va.annule(l_min, l_max) ; // Annule la Valeur
// Annulation des membres derives
if (p_dsdr != 0x0) p_dsdr->annule(l_min, l_max) ;
if (p_srdsdt != 0x0) p_srdsdt->annule(l_min, l_max) ;
if (p_srstdsdp != 0x0) p_srstdsdp->annule(l_min, l_max) ;
if (p_dsdt != 0x0) p_dsdt->annule(l_min, l_max) ;
if (p_stdsdp != 0x0) p_stdsdp->annule(l_min, l_max) ;
if (p_dsdx != 0x0) p_dsdx->annule(l_min, l_max) ;
if (p_dsdy != 0x0) p_dsdy->annule(l_min, l_max) ;
if (p_dsdz != 0x0) p_dsdz->annule(l_min, l_max) ;
if (p_lap != 0x0) p_lap->annule(l_min, l_max) ;
if (p_lapang != 0x0) p_lapang->annule(l_min, l_max) ;
if (p_integ != 0x0) delete p_integ ;
if (p_dsdradial != 0x0) p_dsdradial->annule(l_min, l_max) ;
}
}
//------------//
// Assignment //
//------------//
// From tensor
// -----------
void Scalar::operator=(const Tensor& uu) {
assert(uu.valence == 0) ;
operator=(*(uu.cmp[0])) ;
}
// From Scalar
// ----------
void Scalar::operator=(const Scalar& ci) {
assert(&ci != this) ; // pour eviter l'auto-affectation
// Les elements fixes
assert( mp == ci.mp ) ;
dzpuis = ci.dzpuis ;
// La valeur eventuelle
switch(ci.etat) {
case ETATNONDEF: {
set_etat_nondef() ;
break ; // valeur par defaut
}
case ETATZERO: {
set_etat_zero() ;
break ;
}
case ETATUN: {
set_etat_one() ;
va.set_base( ci.va.get_base() ) ;
break ;
}
case ETATQCQ: {
// Menage general de la Valeur, mais pas des quantites derivees !
va.del_t() ;
set_etat_qcq() ; // set_etat_qcq n'appelle pas del_deriv()
va = ci.va ;
// On detruit les quantites derivees (seulement lorsque tout est fini !)
del_deriv() ;
break ;
}
default: {
cout << "Unkwown state in Scalar::operator=(const Scalar&) !"
<< endl ;
abort() ;
break ;
}
}
}
// From Cmp
// --------
void Scalar::operator=(const Cmp& ci) {
// Menage general de la Valeur, mais pas des quantites derivees !
va.del_t() ;
// Les elements fixes
assert( mp == ci.get_mp() ) ;
dzpuis = ci.get_dzpuis() ;
// La valeur eventuelle
switch(ci.get_etat()) {
case ETATNONDEF: {
set_etat_nondef() ;
break ; // valeur par defaut
}
case ETATZERO: {
set_etat_zero() ;
break ;
}
case ETATQCQ: {
set_etat_qcq() ;
va = ci.va ;
// On detruit les quantites derivees (seulement lorsque tout est fini !)
del_deriv() ;
break ;
}
default: {
cout << "Unkwown state in Scalar::operator=(const Cmp&) !"
<< endl ;
abort() ;
break ;
}
}
}
// From Valeur
// -----------
void Scalar::operator=(const Valeur& vi) {
// Traitement de l'auto-affectation :
if (&vi == &va) {
return ;
}
// Protection
assert(vi.get_etat() != ETATNONDEF) ;
// Menage general de la Valeur, mais pas des quantites derivees !
va.del_t() ;
// La valeure eventuelle
switch(vi.get_etat()) {
case ETATZERO: {
set_etat_zero() ;
break ;
}
case ETATQCQ: {
set_etat_qcq() ;
va = vi ;
// On detruit les quantites derivees (seulement lorsque tout est fini !)
del_deriv() ;
break ;
}
default: {
cout << "Unkwown state in Scalar::operator=(const Valeur&) !" << endl ;
abort() ;
break ;
}
}
}
// From Mtbl
// ---------
void Scalar::operator=(const Mtbl& mi) {
// Protection
assert(mi.get_etat() != ETATNONDEF) ;
assert(&mi != va.c) ; // pour eviter l'auto-affectation
// Menage general de la Valeur, mais pas des quantites derivees !
va.del_t() ;
// La valeure eventuelle
switch(mi.get_etat()) {
case ETATZERO: {
set_etat_zero() ;
break ;
}
case ETATQCQ: {
set_etat_qcq() ;
va = mi ;
// On detruit les quantites derivees (seulement lorsque tout est fini !)
del_deriv() ;
break ;
}
default: {
cout << "Unkwown state in Scalar::operator=(const Mtbl&) !" << endl ;
abort() ;
break ;
}
}
}
// From double
// -----------
void Scalar::operator=(double x) {
if (x == double(0)) {
set_etat_zero() ;
}
else {
if (x == double(1)) {
set_etat_one() ;
}
else {
set_etat_qcq() ;
del_deriv() ;
va = x ;
}
}
dzpuis = 0 ;
}
// From int
// --------
void Scalar::operator=(int n) {
if (n == 0) {
set_etat_zero() ;
}
else {
if (n == 1) {
set_etat_one() ;
}
else {
set_etat_qcq() ;
del_deriv() ;
va = n ;
}
}
dzpuis = 0 ;
}
// Conversion to a Cmp
//----------------------
Scalar::operator Cmp() const {
Cmp resu(mp) ;
resu = va ;
resu.set_dzpuis(dzpuis) ;
return resu ;
}
//------------//
// Sauvegarde //
//------------//
void Scalar::sauve(FILE* fd) const {
va.sauve(fd) ; // la valeur (en premier pour la construction
// lors de la lecture du fichier)
fwrite_be(&etat, sizeof(int), 1, fd) ; // l'etat
fwrite_be(&dzpuis, sizeof(int), 1, fd) ; // dzpuis
}
//------------//
// Impression //
//------------//
// Operator <<
// -----------
ostream& operator<<(ostream& ost, const Scalar& ci) {
switch(ci.etat) {
case ETATNONDEF: {
ost << "*** UNDEFINED STATE *** " << endl ;
break ;
}
case ETATZERO: {
ost << "*** identically ZERO ***" << endl ;
break ;
}
case ETATUN: {
ost << "*** identically ONE ***" << endl ;
break ;
}
case ETATQCQ: {
ost << "*** dzpuis = " << ci.get_dzpuis() << endl ;
ost << ci.va << endl ;
break ;
}
default: {
cout << "operator<<(ostream&, const Scalar&) : unknown state !"
<< endl ;
abort() ;
break ;
}
}
// Termine
return ost ;
}
// spectral_display
//-----------------
void Scalar::spectral_display(const char* comment,
double thres, int precis, ostream& ost) const {
if (comment != 0x0) {
ost << comment << " : " << endl ;
}
// Cas particuliers
//-----------------
if (etat == ETATNONDEF) {
ost << "*** UNDEFINED ***" << endl << endl ;
return ;
}
if (etat == ETATZERO) {
ost << "*** identically ZERO ***" << endl ;
return ;
}
if (etat == ETATUN) {
ost << "*** identically ONE ***" << endl ;
return ;
}
// Cas general : on affiche la Valeur
//------------
if (dzpuis != 0) {
ost << "*** dzpuis = " << dzpuis << endl ;
}
va.display_coef(thres, precis, ost) ;
}
//------------------------------------//
// Spectral bases of the Valeur va //
//------------------------------------//
void Scalar::std_spectral_base() {
va.std_base_scal() ;
}
void Scalar::std_spectral_base_odd() {
va.std_base_scal_odd() ;
}
void Scalar::set_spectral_base(const Base_val& bi) {
va.set_base(bi) ;
}
//--------------------------//
// dzpuis manipulations //
//--------------------------//
void Scalar::set_dzpuis(int dzi) {
dzpuis = dzi ;
}
bool Scalar::dz_nonzero() const {
assert(etat != ETATNONDEF) ;
const Mg3d* mg = mp->get_mg() ;
int nzm1 = mg->get_nzone() - 1;
if (mg->get_type_r(nzm1) != UNSURR) {
return false ;
}
if (etat == ETATZERO) {
return false ;
}
if (etat == ETATUN) {
return true ;
}
assert( etat == ETATQCQ ) ;
if (va.etat == ETATZERO) {
return false ;
}
assert(va.etat == ETATQCQ) ;
if (va.c != 0x0) {
if ( (va.c)->get_etat() == ETATZERO ) {
return false ;
}
assert( (va.c)->get_etat() == ETATQCQ ) ;
if ( (va.c)->t[nzm1]->get_etat() == ETATZERO ) {
return false ;
}
else {
assert( (va.c)->t[nzm1]->get_etat() == ETATQCQ ) ;
return true ;
}
}
else{
assert(va.c_cf != 0x0) ;
if ( (va.c_cf)->get_etat() == ETATZERO ) {
return false ;
}
assert( (va.c_cf)->get_etat() == ETATQCQ ) ;
if ( (va.c_cf)->t[nzm1]->get_etat() == ETATZERO ) {
return false ;
}
else {
assert( (va.c_cf)->t[nzm1]->get_etat() == ETATQCQ ) ;
return true ;
}
}
}
bool Scalar::check_dzpuis(int dzi) const {
if (dz_nonzero()) { // the check must be done
return (dzpuis == dzi) ;
}
else{
return true ;
}
}
//---------------------------------------//
// Value at an arbitrary point //
//---------------------------------------//
double Scalar::val_point(double r, double theta, double phi) const {
assert(etat != ETATNONDEF) ;
if (etat == ETATZERO) {
return double(0) ;
}
if (etat == ETATUN) {
return double(1) ;
}
assert(etat == ETATQCQ) ;
// 1/ Search for the domain and the grid coordinates (xi,theta',phi')
// which corresponds to the point (r,theta,phi)
int l ;
double xi ;
mp->val_lx(r, theta, phi, l, xi) ; // call of val_lx with default
// accuracy parameters
// 2/ Call to the Valeur version
if (l == mp->get_mg()->get_nzone() - 1){ // in the last domain, one must take into account dzpuis
switch (dzpuis) {
case 0:
{
return va.val_point(l, xi, theta, phi);
break;
}
case 1:
{
return va.val_point(l, xi, theta, phi)/r ;
break;
}
case 2:
{
return va.val_point(l, xi, theta, phi)/(r*r) ;
break;
}
case 3:
{
return va.val_point(l, xi, theta, phi)/(r*r*r) ;
break;
}
case 4:
{
return va.val_point(l, xi, theta, phi)/(r*r*r*r) ;
break;
}
default:
{
cout << "scalar::val_point unexpected value of dzpuis !" << endl;
abort();
}
}
}
else{
return va.val_point(l, xi, theta, phi);
}
}
//---------------------------------//
// Multipolar spectrum //
//---------------------------------//
Tbl Scalar::multipole_spectrum() const {
assert (etat != ETATNONDEF) ;
const Mg3d* mg = mp->get_mg() ;
int nzone = mg->get_nzone() ;
int lmax = 0 ;
for (int lz=0; lz<nzone; lz++)
lmax = (lmax < 2*mg->get_nt(lz) - 1 ? 2*mg->get_nt(lz) - 1 : lmax) ;
Tbl resu(nzone, lmax) ;
if (etat == ETATZERO) {
resu.set_etat_zero() ;
return resu ;
}
assert((etat == ETATQCQ) || (etat == ETATUN));
Valeur va_copy = va ;
va_copy.coef() ;
va_copy.ylm() ;
resu.annule_hard() ;
const Base_val& base = va_copy.c_cf->base ;
int m_quant, l_quant, base_r ;
for (int lz=0; lz<nzone; lz++)
for (int k=0 ; k<mg->get_np(lz) ; k++)
for (int j=0 ; j<mg->get_nt(lz) ; j++) {
if (nullite_plm(j, mg->get_nt(lz), k, mg->get_np(lz), base) == 1)
{
// quantic numbers and spectral bases
donne_lm(nzone, lz, j, k, base, m_quant, l_quant, base_r) ;
for (int i=0; i<mg->get_nr(lz); i++) resu.set(lz, l_quant)
+= fabs((*va_copy.c_cf)(lz, k, j, i)) ;
}
}
return resu ;
}
void Scalar::change_triad(const Base_vect& ) {
cout << "WARNING: Scalar::change_triad : "<< endl ;
cout << "This method does nothing ... " << endl ;
}
//---------------------------------//
// Check for NaNs //
//---------------------------------//
bool Scalar::is_nan(bool verb) const {
bool resu = false ;
if (etat == ETATQCQ) {
bool phys = false ;
bool coef = false ;
if (va.c != 0x0)
phys = true ;
if (va.c_cf != 0x0)
coef = true ;
int nz = mp->get_mg()->get_nzone() ;
for (int lz=0; lz<nz; lz++) {
bool domain_p = false ;
bool domain_c = false ;
if (phys) domain_p = (va.c->get_etat() == ETATQCQ) ;
if (coef) domain_c = (va.c_cf->get_etat() == ETATQCQ) ;
if (domain_p) {
int np = mp->get_mg()->get_np(lz) ;
int nt = mp->get_mg()->get_nt(lz) ;
int nr = mp->get_mg()->get_nr(lz) ;
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
if (isnan( (*va.c)(lz, k, j, i) ) ) {
resu = true ;
if (verb) {
cout << "NaN found at physical grid point domain = " << lz
<< ", i= " << i << ", j= " << j << ", k= " << k << endl ;
}
}
} // i loop
} // j loop
} // k loop
}
if (domain_c) {
int np = mp->get_mg()->get_np(lz) + 2 ;
int nt = mp->get_mg()->get_nt(lz) ;
int nr = mp->get_mg()->get_nr(lz) ;
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
if (isnan( (*va.c_cf)(lz, k, j, i) ) ) {
resu = true ;
if (verb) {
cout << "NaN found at coefficient, domain = " << lz
<< ", i= " << i << ", j= " << j << ", k= " << k << endl ;
}
}
} // i loop
} // j loop
} // k loop
}
} // lz loop
} // etat condition
return resu ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition cmp.h:903
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Lorene::Mtbl_cf::base
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Mtbl::get_etat
int get_etat() const
Gives the logical state.
Definition mtbl.h:277
Lorene::Mtbl::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
Lorene::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Lorene::Mtbl::get_nzone
int get_nzone() const
Gives the number of zones (domains).
Definition mtbl.h:280
Lorene::Scalar::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
Lorene::Scalar::set_etat_one
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition scalar.C:334
Lorene::Scalar::del_deriv
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition scalar.C:287
Lorene::Scalar::p_lap
Scalar * p_lap
Pointer on the Laplacian of *this (0x0 if not up to date).
Definition scalar.h:448
Lorene::Scalar::p_dsdy
Scalar * p_dsdy
Pointer on of *this , where (0x0 if not up to date).
Definition scalar.h:439
Lorene::Scalar::~Scalar
virtual ~Scalar()
Destructor.
Definition scalar.C:267
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Lorene::Scalar::p_dsdrho
Scalar * p_dsdrho
Pointer on of *this.
Definition scalar.h:458
Lorene::Scalar::operator=
void operator=(const Scalar &a)
Assignment to another Scalar defined on the same mapping.
Definition scalar.C:446
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
Lorene::Scalar::dz_nonzero
bool dz_nonzero() const
Returns true if the last domain is compactified and *this is not zero in this domain.
Definition scalar.C:814
Lorene::Scalar::operator<<
friend ostream & operator<<(ostream &, const Scalar &)
Display.
Definition scalar.C:702
Lorene::Scalar::p_dsdr
Scalar * p_dsdr
Pointer on of *this (0x0 if not up to date).
Definition scalar.h:411
Lorene::Scalar::p_stdsdp
Scalar * p_stdsdp
Pointer on of *this (0x0 if not up to date).
Definition scalar.h:429
Lorene::Scalar::p_dsdx
Scalar * p_dsdx
Pointer on of *this , where (0x0 if not up to date).
Definition scalar.h:434
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::is_nan
bool is_nan(bool verbose=false) const
Looks for NaNs (not a number) in the scalar field.
Definition scalar.C:1008
Lorene::Scalar::p_dsdz
Scalar * p_dsdz
Pointer on of *this , where (0x0 if not up to date).
Definition scalar.h:444
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::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition scalar.C:997
Lorene::Scalar::ind_lap
int ind_lap
Power of r by which the last computed Laplacian has been multiplied in the compactified external doma...
Definition scalar.h:463
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Lorene::Scalar::Scalar
Scalar(const Map &mpi)
Constructor from mapping.
Definition scalar.C:204
Lorene::Scalar::p_srstdsdp
Scalar * p_srstdsdp
Pointer on of *this (0x0 if not up to date).
Definition scalar.h:421
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::del_t
void del_t()
Logical destructor.
Definition scalar.C:279
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::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
Lorene::Scalar::set_etat_nondef
virtual void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition scalar.C:344
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Scalar::p_lapang
Scalar * p_lapang
Pointer on the Laplacian of *this (0x0 if not up to date).
Definition scalar.h:452
Lorene::Scalar::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition scalar.C:890
Lorene::Scalar::multipole_spectrum
Tbl multipole_spectrum() const
Gives the spectrum in terms of multipolar modes l .
Definition scalar.C:956
Lorene::Scalar::std_spectral_base_odd
virtual void std_spectral_base_odd()
Sets the spectral bases of the Valeur va to the standard odd ones for a scalar field.
Definition scalar.C:791
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:405
Lorene::Scalar::spectral_display
virtual void spectral_display(const char *comment=0x0, double threshold=1.e-7, int precision=4, ostream &ostr=cout) const
Displays the spectral coefficients and the associated basis functions.
Definition scalar.C:741
Lorene::Scalar::p_dsdradial
Scalar * p_dsdradial
Pointer on of *this.
Definition scalar.h:455
Lorene::Scalar::p_integ
Tbl * p_integ
Pointer on the space integral of *this (values in each domain) (0x0 if not up to date).
Definition scalar.h:468
Lorene::Scalar::p_dsdt
Scalar * p_dsdt
Pointer on of *this (0x0 if not up to date).
Definition scalar.h:424
Lorene::Scalar::set_der_0x0
void set_der_0x0() const
Sets the pointers for derivatives to 0x0.
Definition scalar.C:306
Lorene::Scalar::p_srdsdt
Scalar * p_srdsdt
Pointer on of *this (0x0 if not up to date).
Definition scalar.h:416
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::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
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
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::get_base
const Base_val & get_base() const
Return the bases for spectral expansions (member base ).
Definition valeur.h:480
Lorene::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition valeur.h:726
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
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::fread_be
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
Lorene::fwrite_be
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene::Tensor::Tensor
Tensor(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i)
Standard constructor.
Definition tensor.C:202
Lorene::Tensor::valence
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...).
Definition tensor.h:298
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
Lorene::Tensor::del_deriv
virtual void del_deriv() const
Deletes the derived quantities.
Definition tensor.C:398
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
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
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:718