map_af_lap.C
/*
* Copyright (c) 1999-2001 Eric Gourgoulhon
* Copyright (c) 1999-2001 Philippe Grandclement
*
* 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 map_af_lap_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_lap.C,v 1.7 2014/10/13 08:53:02 j_novak Exp $" ;
/*
* $Id: map_af_lap.C,v 1.7 2014/10/13 08:53:02 j_novak Exp $
* $Log: map_af_lap.C,v $
* Revision 1.7 2014/10/13 08:53:02 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 2005/11/24 09:25:06 j_novak
* Added the Scalar version for the Laplacian
*
* Revision 1.4 2003/10/15 16:03:37 j_novak
* Added the angular Laplace operator for Scalar.
*
* Revision 1.3 2003/10/03 15:58:48 j_novak
* Cleaning of some headers
*
* Revision 1.2 2002/10/16 14:36:41 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.16 2000/08/16 10:35:41 eric
* Suppression de Mtbl::dzpuis.
*
* Revision 2.15 2000/02/25 09:01:47 eric
* Remplacement de (uu.get_dzpuis() == 0) par uu.check_dzpuis(0).
*
* Revision 2.14 2000/02/08 14:19:53 phil
* correction annulation derniere zone
*
* Revision 2.13 2000/01/27 17:52:16 phil
* corrections diverses et variees
*
* Revision 2.12 2000/01/26 13:10:34 eric
* Reprototypage complet :
* le resultat est desormais suppose alloue a l'exterieur de la routine
* et est passe en argument (Cmp& resu).
*
* Revision 2.11 1999/11/30 12:49:53 eric
* Valeur::base est desormais du type Base_val et non plus Base_val*.
*
* Revision 2.10 1999/11/29 12:57:42 eric
* REORGANISATION COMPLETE: nouveau prototype : Valeur --> Cmp
* Utilisation de la nouvelle arithmetique des Valeur's.
*
* Revision 2.9 1999/10/27 15:44:00 eric
* Suppression du membre Cmp::c.
*
* Revision 2.8 1999/10/14 14:27:35 eric
* Methode const.
*
* Revision 2.7 1999/09/06 16:26:03 phil
* ajout de la version Cmp
*
* Revision 2.6 1999/09/06 14:51:20 phil
* ajout du laplacien
*
* Revision 2.5 1999/05/03 15:22:00 phil
* Correction des bases
*
* Revision 2.4 1999/04/28 10:33:02 phil
* correction du cas zec_mult_r = 4
*
* Revision 2.3 1999/04/27 09:22:25 phil
* *** empty log message ***
*
* Revision 2.2 1999/04/27 09:17:27 phil
* corrections diverses et variees ....
*
* Revision 2.1 1999/04/26 17:24:21 phil
* correction de gestion de base
*
* Revision 2.0 1999/04/26 16:33:55 phil
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Map/map_af_lap.C,v 1.7 2014/10/13 08:53:02 j_novak Exp $
*
*/
// Fichiers include
// ----------------
#include <cstdlib>
#include <cmath>
#include "cmp.h"
#include "tensor.h"
//******************************************************************************
/*
* Calcul du laplacien d'un Scalar ou Cmp f dans le cas ou le mapping
* (coordonnees-grille) --> (coordonnees physiques) est affine.
* Le Laplacien est calcule suivant l'equation:
*
* Lap(f) = d^2 f/dr^2 + 1/r ( 2 df/dr + 1/r Lap_ang(f) ) (1)
*
* avec
*
* Lap_ang(f) := d^2 f/dtheta^2 + 1/tan(theta) df/dtheta
* + 1/sin^2(theta) d^2 f/dphi^2 (2)
*
* Le laplacien angulaire (2) est calcule par passage aux harmoniques
* spheriques, suivant la formule
*
* Lap_ang( f_{lm} Y_l^m ) = - l(l+1) f_{lm} Y_l^m (3)
*
* Dans la zone externe compactifiee (ZEC), la routine calcule soit Lap(f),
* soit r^2 Lap(f), soit r^4 Lap(f) suivant la valeur du drapeau zec_mult_r :
*
* -- pour zec_mult_r = 0, on calcule Lap(f) suivant la formule
*
* Lap(f) = u^2 [ u^2 d^2 f/du^2 + Lap_ang(f) ] , avec u = 1/r (4)
*
* -- pour zec_mult_r = 2, on calcule r^2 Lap(f) suivant la formule
*
* r^2 Lap(f) = u^2 d^2 f/du^2 + Lap_ang(f) (5)
*
* -- pour zec_mult_r = 4, on calcule 4^2 Lap(f) suivant la formule
*
* r^4 Lap(f) = d^2 f /du^2 + 1/u^2 Lap_ang(f) (6)
*
*
*
* Entree:
* ------
* const Scalar& / Cmp& uu : champ f dont on veut calculer le laplacien
*
*
* int zec_mult_r : drapeau indiquant la quantite calculee dans la ZEC:
* zec_mult_r = 0 : lapff = Lap(f)
* zec_mult_r = 2 : lapff = r^2 Lap(f)
* zec_mult_r = 4 : lapff = r^4 Lap(f)
* Sortie:
* ------
* Scalar& / Cmp& resu : Lap(f)
*
*/
namespace Lorene {
//**********************
// Scalar version
//**********************
void Map_af::laplacien(const Scalar& uu, int zec_mult_r, Scalar& resu) const {
assert (uu.get_etat() != ETATNONDEF) ;
assert (uu.get_mp().get_mg() == mg) ;
// Particular case of null value:
if ((uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN) ) {
resu.set_etat_zero() ;
return ;
}
assert( uu.get_etat() == ETATQCQ ) ;
assert( uu.check_dzpuis(0) ) ;
resu.set_etat_qcq() ;
int nz = get_mg()->get_nzone() ;
int nzm1 = nz - 1 ;
// On recupere les bases d'entree :
Base_val base_entree = uu.get_spectral_base() ;
// Separation zones internes / ZEC :
Scalar uuext(uu) ;
Valeur ff = uu.get_spectral_va() ;
if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
uuext.annule(0, nzm1-1) ;
ff.annule(nzm1) ;
}
else {
uuext.set_etat_zero() ; // pas de ZEC
}
//=========================================================================
// Calcul dans les zones internes (noyau + coquilles)
//=========================================================================
//----------------------------
// Derivations par rapport a x
//----------------------------
Valeur d2ff = ff.d2sdx2() ; // d^2/dx^2 partout
Valeur dff = ff.dsdx() ; // d/dx partout
//---------------------------
// Calcul de 1/x Lap_ang(f) ---> resultat dans ff
//---------------------------
//... Passage en harmoniques spheriques
ff.coef() ;
ff.ylm() ;
//... Multiplication par -l*(l+1)
ff = ff.lapang() ;
//... Division par x:
ff = ff.sx() ; // 1/x ds. le noyau
// Id ds. les coquilles
//----------------------------------------------
// On repasse dans l'espace des configurations
// pour effectuer les multiplications par
// les derivees du chgmt. de var.
//----------------------------------------------
d2ff.coef_i() ;
dff.coef_i() ;
ff.coef_i() ;
//-----------------------------------------------
// Calcul de 1/x ( 2 df/dr + 1/r Lap_ang(f) )
// Le resultat est mis dans ff
//-----------------------------------------------
Base_val sauve_base = dff.base ;
ff = double(2) * ( dxdr * dff) + xsr * ff ;
ff.base = sauve_base ;
ff = ff.sx() ; // 1/x ds. le noyau
// Id ds. les coquilles
ff.coef_i() ;
//---------------------------------------------
// Calcul de Lap(f) suivant l'Eq. (1)
// Le resultat est mis dans ff
//-----------------------------------------------
sauve_base = d2ff.base ;
ff = dxdr * dxdr * d2ff + xsr * ff ;
ff.base = sauve_base ;
if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
//=========================================================================
// Calcul dans la ZEC
//=========================================================================
Valeur& ffext = uuext.set_spectral_va() ;
//----------------------------
// Derivation par rapport a x
//----------------------------
d2ff = ffext.d2sdx2() ; // d^2/dx^2 partout
//---------------------------
// Calcul de Lap_ang(f) ---> resultat dans ffext
//---------------------------
//... Passage en harmoniques spheriques
ffext.coef() ;
ffext.ylm() ;
//... Multiplication par -l*(l+1)
ffext = ffext.lapang() ;
switch (zec_mult_r) {
case 0 : { // calcul de Lap(f) suivant l'Eq. (4)
d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
// par (x-1)^2
sauve_base = ffext.base ;
ffext = dxdr * dxdr / (xsr*xsr) * d2ff + ffext ;
ffext.base = sauve_base ;
ffext.mult2_xm1_zec() ; // multiplication par (x-1)^2
ffext.coef_i() ;
sauve_base = ffext.base ;
ffext = ffext / (xsr*xsr) ;
ffext.base = sauve_base ;
break ;
}
case 2 : { // calcul de r^2 Lap(f) suivant l'Eq. (5)
d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
// par (x-1)^2
sauve_base = ffext.base ;
ffext = dxdr*dxdr / (xsr*xsr) * d2ff + ffext ;
ffext.base = sauve_base ;
break ;
}
case 4 : { // calcul de r^4 Lap(f) suivant l'Eq. (6)
ffext = ffext.sx2() ; // division de Lap_ang(f) par (x-1)^2
sauve_base = ffext.base ;
ffext = dxdr*dxdr * d2ff + xsr*xsr * ffext ;
ffext.base = sauve_base ;
break ;
}
default : {
cout << "Map_af::laplacien : unknown value of zec_mult_r !"
<< endl << " zec_mult_r = " << zec_mult_r << endl ;
abort() ;
}
} // fin des differents cas pour zec_mult_r
// Resultat final
ff = ff + ffext ;
} // fin du cas ou il existe une ZEC
// Les bases de sorties sont egales aux bases d'entree:
ff.base = base_entree ;
resu = ff ;
resu.set_dzpuis(zec_mult_r) ;
}
//**********************
// Cmp version
//**********************
void Map_af::laplacien(const Cmp& uu, int zec_mult_r, Cmp& resu) const {
assert (uu.get_etat() != ETATNONDEF) ;
assert (uu.get_mp()->get_mg() == mg) ;
// Particular case of null value:
if (uu.get_etat() == ETATZERO) {
resu.set_etat_zero() ;
return ;
}
assert( uu.get_etat() == ETATQCQ ) ;
assert( uu.check_dzpuis(0) ) ;
resu.set_etat_qcq() ;
int nz = get_mg()->get_nzone() ;
int nzm1 = nz - 1 ;
// On recupere les bases d'entree :
Base_val base_entree = (uu.va).base ;
// Separation zones internes / ZEC :
Cmp uuext(uu) ;
Valeur ff = uu.va ;
if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
uuext.annule(0, nzm1-1) ;
ff.annule(nzm1) ;
}
else {
uuext.set_etat_zero() ; // pas de ZEC
}
//=========================================================================
// Calcul dans les zones internes (noyau + coquilles)
//=========================================================================
//----------------------------
// Derivations par rapport a x
//----------------------------
Valeur d2ff = ff.d2sdx2() ; // d^2/dx^2 partout
Valeur dff = ff.dsdx() ; // d/dx partout
//---------------------------
// Calcul de 1/x Lap_ang(f) ---> resultat dans ff
//---------------------------
//... Passage en harmoniques spheriques
ff.coef() ;
ff.ylm() ;
//... Multiplication par -l*(l+1)
ff = ff.lapang() ;
//... Division par x:
ff = ff.sx() ; // 1/x ds. le noyau
// Id ds. les coquilles
//----------------------------------------------
// On repasse dans l'espace des configurations
// pour effectuer les multiplications par
// les derivees du chgmt. de var.
//----------------------------------------------
d2ff.coef_i() ;
dff.coef_i() ;
ff.coef_i() ;
//-----------------------------------------------
// Calcul de 1/x ( 2 df/dr + 1/r Lap_ang(f) )
// Le resultat est mis dans ff
//-----------------------------------------------
Base_val sauve_base = dff.base ;
ff = double(2) * ( dxdr * dff) + xsr * ff ;
ff.base = sauve_base ;
ff = ff.sx() ; // 1/x ds. le noyau
// Id ds. les coquilles
ff.coef_i() ;
//---------------------------------------------
// Calcul de Lap(f) suivant l'Eq. (1)
// Le resultat est mis dans ff
//-----------------------------------------------
sauve_base = d2ff.base ;
ff = dxdr * dxdr * d2ff + xsr * ff ;
ff.base = sauve_base ;
if (get_mg()->get_type_r(nzm1) == UNSURR) { // il existe une ZEC
//=========================================================================
// Calcul dans la ZEC
//=========================================================================
Valeur& ffext = uuext.va ;
//----------------------------
// Derivation par rapport a x
//----------------------------
d2ff = ffext.d2sdx2() ; // d^2/dx^2 partout
//---------------------------
// Calcul de Lap_ang(f) ---> resultat dans ffext
//---------------------------
//... Passage en harmoniques spheriques
ffext.coef() ;
ffext.ylm() ;
//... Multiplication par -l*(l+1)
ffext = ffext.lapang() ;
switch (zec_mult_r) {
case 0 : { // calcul de Lap(f) suivant l'Eq. (4)
d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
// par (x-1)^2
sauve_base = ffext.base ;
ffext = dxdr * dxdr / (xsr*xsr) * d2ff + ffext ;
ffext.base = sauve_base ;
ffext.mult2_xm1_zec() ; // multiplication par (x-1)^2
ffext.coef_i() ;
sauve_base = ffext.base ;
ffext = ffext / (xsr*xsr) ;
ffext.base = sauve_base ;
break ;
}
case 2 : { // calcul de r^2 Lap(f) suivant l'Eq. (5)
d2ff.mult2_xm1_zec() ; // multiplication de d^2f/dx^2
// par (x-1)^2
sauve_base = ffext.base ;
ffext = dxdr*dxdr / (xsr*xsr) * d2ff + ffext ;
ffext.base = sauve_base ;
break ;
}
case 4 : { // calcul de r^4 Lap(f) suivant l'Eq. (6)
ffext = ffext.sx2() ; // division de Lap_ang(f) par (x-1)^2
sauve_base = ffext.base ;
ffext = dxdr*dxdr * d2ff + xsr*xsr * ffext ;
ffext.base = sauve_base ;
break ;
}
default : {
cout << "Map_af::laplacien : unknown value of zec_mult_r !"
<< endl << " zec_mult_r = " << zec_mult_r << endl ;
abort() ;
}
} // fin des differents cas pour zec_mult_r
// Resultat final
ff = ff + ffext ;
} // fin du cas ou il existe une ZEC
// Les bases de sorties sont egales aux bases d'entree:
ff.base = base_entree ;
resu = ff ;
resu.set_dzpuis(zec_mult_r) ;
}
void Map_af::lapang(const Scalar& uu, Scalar& resu) const {
assert (uu.get_etat() != ETATNONDEF) ;
assert (uu.get_mp().get_mg() == mg) ;
// Particular case of null result:
if ( (uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN) ) {
resu.set_etat_zero() ;
return ;
}
assert( uu.get_etat() == ETATQCQ ) ;
resu.set_etat_qcq() ;
Valeur ff = uu.get_spectral_va() ;
//... Passage en harmoniques spheriques
ff.ylm() ;
//... Multiplication par -l*(l+1)
resu = ff.lapang() ;
resu.set_dzpuis(uu.get_dzpuis()) ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::sx
void sx()
The basis is transformed as with a multiplication.
Definition base_val_manip.C:124
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::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:304
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
Lorene::Cmp::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 cmp.C:715
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Map_af::laplacien
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const
Computes the Laplacian of a scalar field.
Definition map_af_lap.C:179
Lorene::Map_af::lapang
virtual void lapang(const Scalar &uu, Scalar &lap) const
Computes the angular Laplacian of a scalar field.
Definition map_af_lap.C:549
Lorene::Map_radial::xsr
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
Lorene::Map_radial::dxdr
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1560
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::lapang
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition scalar_deriv.C:461
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
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::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
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::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_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::sx2
const Valeur & sx2() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ).
Definition valeur_sx2.C:114
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::dsdx
const Valeur & dsdx() const
Returns of *this.
Definition valeur_dsdx.C:111
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition valeur_coef_i.C:140
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Valeur::d2sdx2
const Valeur & d2sdx2() const
Returns of *this.
Definition valeur_d2sdx2.C:114
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Lorene::Valeur::lapang
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition valeur_lapang.C:72
Lorene::Valeur::mult2_xm1_zec
void mult2_xm1_zec()
Applies the following operator to *this : \ Id (r sampling = RARE, FIN ) \ (r -sampling = UNSURR ).
Definition valeur_mult2_xm1.C:76
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
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