sym_tensor_tt_etamu.C
/*
* Methods of class Sym_tensor_tt related to eta and mu
*
* (see file sym_tensor.h for documentation)
*
*/
/*
* Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
*
* 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 sym_tensor_tt_etamu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $" ;
/*
* $Id: sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
* $Log: sym_tensor_tt_etamu.C,v $
* Revision 1.18 2014/10/13 08:53:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.17 2014/10/06 15:13:19 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.16 2006/10/24 13:03:19 j_novak
* New methods for the solution of the tensor wave equation. Perhaps, first
* operational version...
*
* Revision 1.15 2005/04/01 14:28:32 j_novak
* Members p_eta and p_mu are now defined in class Sym_tensor.
*
* Revision 1.14 2004/06/04 09:25:58 e_gourgoulhon
* Method eta(): eta is no longer computed from h^rr but from khi (in the
* case where khi is known).
*
* Revision 1.13 2004/05/25 15:08:44 f_limousin
* Add parameters in argument of the functions update, eta and mu for
* the case of a Map_et.
*
* Revision 1.12 2004/05/24 13:45:29 e_gourgoulhon
* Added parameter dzp to method Sym_tensor_tt::update.
*
* Revision 1.11 2004/05/05 14:24:54 e_gourgoulhon
* Corrected a bug in method set_khi_mu: the division of khi by r^2
* was ommitted in the case dzp=2 !!!
*
* Revision 1.10 2004/04/08 16:38:43 e_gourgoulhon
* Sym_tensor_tt::set_khi_mu: added argument dzp (dzpuis of resulting h^{ij}).
*
* Revision 1.9 2004/03/04 09:53:04 e_gourgoulhon
* Methods eta(), mu() and upate(): use of Scalar::mult_r_dzpuis and
* change of dzpuis behavior of eta and mu.
*
* Revision 1.8 2004/03/03 13:16:21 j_novak
* New potential khi (p_khi) and the functions manipulating it.
*
* Revision 1.7 2004/02/05 13:44:50 e_gourgoulhon
* Major modif. of methods eta(), mu() and update() to treat
* any value of dzpuis, thanks to the new definitions of
* Scalar::mult_r(), Scalar::dsdr(), etc...
*
* Revision 1.6 2004/01/28 13:25:41 j_novak
* The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
* In the div/mult _r_dzpuis, there is no more default value.
*
* Revision 1.5 2003/11/05 15:28:31 e_gourgoulhon
* Corrected error in update.
*
* Revision 1.4 2003/11/04 23:03:34 e_gourgoulhon
* First full version of method update().
* Add method set_rr_mu.
* Method set_eta_mu ---> set_rr_eta_mu.
*
* Revision 1.3 2003/11/04 09:35:27 e_gourgoulhon
* First operational version of update_tp().
*
* Revision 1.2 2003/11/03 22:33:36 e_gourgoulhon
* Added methods update_tp and set_eta_mu.
*
* Revision 1.1 2003/11/03 17:08:37 e_gourgoulhon
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
#include <cassert>
// Headers Lorene
#include "tensor.h"
//--------------//
// khi //
//--------------//
namespace Lorene {
const Scalar& Sym_tensor_tt::khi() const {
if (p_khi == 0x0) { // a new computation is necessary
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
// khi is computed from $h^{rr}$ component
p_khi = new Scalar(operator()(1,1)) ;
p_khi->mult_r() ;
p_khi->mult_r() ;
}
return *p_khi ;
}
//--------------//
// eta //
//--------------//
const Scalar& Sym_tensor_tt::eta(Param* par) const {
if (p_eta == 0x0) { // a new computation is necessary
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
// eta is computed from the divergence-free condition:
int dzp = operator()(1,1).get_dzpuis() ;
int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
Scalar source_eta(*mp) ;
if (p_khi == 0x0) { // eta is computed from h^rr
// -------------------------
source_eta = - operator()(1,1).dsdr() ;
// dhrr contains - dh^{rr}/dr in all domains but the CED,
// in the CED: - r^2 dh^{rr}/dr if dzp = 0 (1)
// - r^(dzp+1) dh^{rr}/dr if dzp > 0 (2)
// Multiplication by r of (-d h^{rr}/dr) (with the same dzpuis as h^{rr})
source_eta.mult_r_dzpuis( dzp ) ;
// Substraction of the h^rr part and multiplication by r :
source_eta -= 3. * operator()(1,1) ;
source_eta.mult_r_dzpuis(dzp_resu) ;
}
else { // eta is computed from khi
// ------------------------
source_eta = - p_khi->dsdr() ;
int diff_dzp = source_eta.get_dzpuis() - dzp_resu ;
assert( diff_dzp >= 0 ) ;
source_eta.dec_dzpuis(diff_dzp) ;
Scalar tmp(*p_khi) ;
tmp.div_r_dzpuis(dzp_resu) ;
source_eta -= tmp ;
}
// Resolution of the angular Poisson equation for eta
// --------------------------------------------------
if (dynamic_cast<const Map_af*>(mp) != 0x0) {
p_eta = new Scalar( source_eta.poisson_angu() ) ;
}
else {
Scalar resu (*mp) ;
resu = 0. ;
mp->poisson_angu(source_eta, *par, resu) ;
p_eta = new Scalar( resu ) ;
}
}
return *p_eta ;
}
//-------------------//
// set_rr_eta_mu //
//-------------------//
void Sym_tensor_tt::set_rr_eta_mu(const Scalar& hrr, const Scalar& eta_i,
const Scalar& mu_i) {
// All this has a meaning only for spherical components:
assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
set(1,1) = hrr ; // h^{rr}
// calls del_deriv() and therefore delete previous
// p_eta and p_mu
p_eta = new Scalar( eta_i ) ; // eta
p_mu = new Scalar( mu_i ) ; // mu
update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
}
//---------------//
// set_rr_mu //
//---------------//
void Sym_tensor_tt::set_rr_mu(const Scalar& hrr, const Scalar& mu_i) {
// All this has a meaning only for spherical components:
assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
set(1,1) = hrr ; // h^{rr}
// calls del_deriv() and therefore delete previous
// p_eta and p_mu
p_mu = new Scalar( mu_i ) ; // mu
eta() ; // computes eta form the divergence-free condition
update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
}
//-------------------//
// set_khi_eta_mu //
//-------------------//
void Sym_tensor_tt::set_khi_eta_mu(const Scalar& khi_i, const Scalar& eta_i,
const Scalar& mu_i) {
// All this has a meaning only for spherical components:
assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
set(1,1) = khi_i ;
set(1,1).div_r() ;
set(1,1).div_r() ; // h^{rr}
// calls del_deriv() and therefore delete previous
// p_khi, p_eta and p_mu
p_khi = new Scalar( khi_i ) ; // khi
p_eta = new Scalar( eta_i ) ; // eta
p_mu = new Scalar( mu_i ) ; // mu
update( khi_i.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
}
//---------------//
// set_khi_mu //
//---------------//
void Sym_tensor_tt::set_khi_mu(const Scalar& khi_i, const Scalar& mu_i,
int dzp, Param* par1, Param* par2, Param* par3) {
// All this has a meaning only for spherical components:
assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
set(1,1) = khi_i ;
// calls del_deriv() and therefore delete previous
// p_eta and p_mu
assert( khi_i.check_dzpuis(0) ) ;
if (dzp == 0) {
set(1,1).div_r() ;
set(1,1).div_r() ; // h^{rr}
}
else {
assert(dzp == 2) ; //## temporary: the other cases are not treated yet
set(1,1).div_r_dzpuis(1) ;
set(1,1).div_r_dzpuis(2) ;
}
p_khi = new Scalar ( khi_i ) ; // khi
p_mu = new Scalar( mu_i ) ; // mu
if (dynamic_cast<const Map_af*>(mp) != 0x0) {
eta() ; // computes eta form the divergence-free condition
// dzp = 0 ==> eta.dzpuis = 0
// dzp = 2 ==> eta.dzpuis = 1
update(dzp) ; // all h^{ij}, except for h^{rr}
}
else {
eta(par1) ; // computes eta form the divergence-free condition
// dzp = 0 ==> eta.dzpuis = 0
// dzp = 2 ==> eta.dzpuis = 1
update(dzp, par2, par3) ; // all h^{ij}, except for h^{rr}
}
}
//-------------//
// update //
//-------------//
void Sym_tensor_tt::update(int dzp, Param* par1, Param* par2) {
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
assert( (p_eta != 0x0) && (p_mu != 0x0) ) ;
Itbl idx(2) ;
idx.set(0) = 1 ; // r index
// h^{r theta} :
// ------------
idx.set(1) = 2 ; // theta index
*cmp[position(idx)] = p_eta->srdsdt() - p_mu->srstdsdp() ;
if (dzp == 0) {
assert( cmp[position(idx)]->check_dzpuis(2) ) ;
cmp[position(idx)]->dec_dzpuis(2) ;
}
assert( cmp[position(idx)]->check_dzpuis(dzp) ) ;
// h^{r phi} :
// ------------
idx.set(1) = 3 ; // phi index
*cmp[position(idx)] = p_eta->srstdsdp() + p_mu->srdsdt() ;
if (dzp == 0) {
assert( cmp[position(idx)]->check_dzpuis(2) ) ;
cmp[position(idx)]->dec_dzpuis(2) ;
}
assert( cmp[position(idx)]->check_dzpuis(dzp) ) ;
// h^{theta phi} and h^{phi phi}
// -----------------------------
//-------------- Computation of T^theta --> taut :
Scalar tautst = operator()(1,2).dsdr() ;
// dhrr contains dh^{rt}/dr in all domains but the CED,
// in the CED: r^2 dh^{rt}/dr if dzp = 0 (1)
// r^(dzp+1) dh^{rt}/dr if dzp > 0 (2)
// Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
tautst.mult_r_dzpuis( operator()(1,2).get_dzpuis() ) ;
// Addition of the remaining parts :
tautst += 3 * operator()(1,2) - operator()(1,1).dsdt() ;
tautst.mult_sint() ;
Scalar tmp = operator()(1,1) ;
tmp.mult_cost() ; // h^{rr} cos(th)
tautst -= tmp ; // T^th / sin(th)
Scalar taut = tautst ;
taut.mult_sint() ; // T^th
//----------- Computation of T^phi --> taup :
Scalar taupst = - operator()(1,3).dsdr() ;
// dhrr contains - dh^{rp}/dr in all domains but the CED,
// in the CED: - r^2 dh^{rp}/dr if dzp = 0 (3)
// - r^(dzp+1) dh^{rp}/dr if dzp > 0 (4)
// Multiplication by r of -dh^{rp}/dr (with the same dzpuis than h^{rp})
taupst.mult_r_dzpuis( operator()(1,3).get_dzpuis() ) ;
// Addition of the remaining part :
taupst -= 3 * operator()(1,3) ;
taupst.mult_sint() ; // T^ph / sin(th)
Scalar taup = taupst ;
taup.mult_sint() ; // T^ph
//------------------- Computation of F and h^[ph ph}
tmp = tautst ;
tmp.mult_cost() ; // T^th / tan(th)
// dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
tmp = taut.dsdt() + tmp + taup.stdsdp() ;
Scalar tmp2 (*mp) ;
if (dynamic_cast<const Map_af*>(mp) != 0x0) {
tmp2 = tmp.poisson_angu() ; // F
}
else {
tmp2 = 0. ;
mp->poisson_angu(tmp, *par1, tmp2) ; // F
}
tmp2.div_sint() ;
tmp2.div_sint() ; // h^{ph ph}
idx.set(0) = 3 ; // phi index
idx.set(1) = 3 ; // phi index
*cmp[position(idx)] = tmp2 ; // h^{ph ph} is updated
//------------------- Computation of G and h^[th ph}
tmp = taupst ;
tmp.mult_cost() ; // T^ph / tan(th)
// - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
tmp = - taut.stdsdp() + taup.dsdt() + tmp ;
if (dynamic_cast<const Map_af*>(mp) != 0x0) {
tmp2 = tmp.poisson_angu() ; // G
}
else {
tmp2 = 0. ;
mp->poisson_angu(tmp, *par2, tmp2) ; // G
}
tmp2.div_sint() ;
tmp2.div_sint() ; // h^{th ph}
idx.set(0) = 2 ; // theta index
idx.set(1) = 3 ; // phi index
*cmp[position(idx)] = tmp2 ; // h^{th ph} is updated
// h^{th th} (from the trace-free condition)
// ---------
idx.set(1) = 2 ; // theta index
*cmp[position(idx)] = - operator()(1,1) - operator()(3,3) ;
Sym_tensor_trans::del_deriv() ; //## in order not to delete p_eta and p_mu
}
//-----------------//
// set_A_tilde_B //
//-----------------//
void Sym_tensor_tt::set_A_tildeB(const Scalar& a_in, const Scalar& tb_in,
Param* par_bc, Param* par_mat) {
Scalar zero(*mp) ;
zero.set_etat_zero() ;
set_AtB_trace(a_in, tb_in, zero, par_bc, par_mat) ;
return ;
}
}
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::div_r_dzpuis
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:147
Lorene::Scalar::dsdt
const Scalar & dsdt() const
Returns of *this .
Definition scalar_deriv.C:208
Lorene::Scalar::div_sint
void div_sint()
Division by .
Definition scalar_th_manip.C:98
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Lorene::Scalar::mult_cost
void mult_cost()
Multiplication by .
Definition scalar_th_manip.C:58
Lorene::Scalar::poisson_angu
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition scalar_pde.C:200
Lorene::Scalar::stdsdp
const Scalar & stdsdp() const
Returns of *this .
Definition scalar_deriv.C:238
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::div_r
void div_r()
Division by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:133
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::mult_sint
void mult_sint()
Multiplication by .
Definition scalar_th_manip.C:84
Lorene::Scalar::mult_r_dzpuis
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:221
Lorene::Scalar::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:418
Lorene::Sym_tensor_trans::set_AtB_trace
void set_AtB_trace(const Scalar &a_in, const Scalar &tb_in, const Scalar &trace, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A , and the trace.
Definition sym_tensor_trans_aux.C:302
Lorene::Sym_tensor_trans::del_deriv
virtual void del_deriv() const
Deletes the derived quantities.
Definition sym_tensor_trans.C:164
Lorene::Sym_tensor_tt::set_khi_mu
void set_khi_mu(const Scalar &khi_i, const Scalar &mu_i, int dzp=0, Param *par1=0x0, Param *par2=0x0, Param *par3=0x0)
Sets the component , as well as the angular potential (see member p_khi and p_mu ).
Definition sym_tensor_tt_etamu.C:288
Lorene::Sym_tensor_tt::set_rr_mu
void set_rr_mu(const Scalar &hrr, const Scalar &mu_i)
Sets the component , as well as the angular potential (see member p_mu ).
Definition sym_tensor_tt_etamu.C:238
Lorene::Sym_tensor_tt::set_A_tildeB
void set_A_tildeB(const Scalar &a_in, const Scalar &tb_in, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A and .
Definition sym_tensor_tt_etamu.C:487
Lorene::Sym_tensor_tt::p_khi
Scalar * p_khi
Field such that the component .
Definition sym_tensor.h:946
Lorene::Sym_tensor_tt::set_khi_eta_mu
void set_khi_eta_mu(const Scalar &khi_i, const Scalar &eta_i, const Scalar &mu_i)
Sets the component , as well as the angular potentials and (see members p_khi , p_eta and p_mu ).
Definition sym_tensor_tt_etamu.C:260
Lorene::Sym_tensor_tt::update
void update(int dzp, Param *par1=0x0, Param *par2=0x0)
Computes the components , , , and , from and the potentials and .
Definition sym_tensor_tt_etamu.C:338
Lorene::Sym_tensor_tt::eta
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
Definition sym_tensor_tt_etamu.C:140
Lorene::Sym_tensor_tt::set_rr_eta_mu
void set_rr_eta_mu(const Scalar &hrr, const Scalar &eta_i, const Scalar &mu_i)
Sets the component , as well as the angular potentials and (see members p_eta and p_mu ).
Definition sym_tensor_tt_etamu.C:215
Lorene::Sym_tensor_tt::khi
const Scalar & khi() const
Gives the field such that the component .
Definition sym_tensor_tt_etamu.C:116
Lorene::Sym_tensor::p_mu
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:274
Lorene::Sym_tensor::p_eta
Scalar * p_eta
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:260
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::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
Lorene::Tensor_sym::position
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor_sym.C:245
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene
Lorene prototypes.
Definition app_hor.h:64