spheroid.C
/*
* Definition of methods for the class Spheroid and its subclass App_hor
*
*/
/*
* Copyright (c) 2006 Nicolas Vasset, Jerome Novak & Jose-Luis Jaramillo
*
* 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 spheroid_C[] = "$Header: /cvsroot/Lorene/C++/Source/App_hor/spheroid.C,v 1.21 2014/10/13 08:52:38 j_novak Exp $" ;
/*
* $Id: spheroid.C,v 1.21 2014/10/13 08:52:38 j_novak Exp $
* $Log: spheroid.C,v $
* Revision 1.21 2014/10/13 08:52:38 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.20 2014/10/06 15:12:56 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.19 2012/05/18 16:27:05 j_novak
* ... and another memory leak
*
* Revision 1.18 2012/05/18 15:19:14 j_novak
* Corrected a memory leak
*
* Revision 1.17 2012/05/11 14:11:24 j_novak
* Modifications to deal with the accretion of a scalar field
*
* Revision 1.16 2009/12/01 08:44:24 j_novak
* Added the missing operator=().
*
* Revision 1.15 2008/12/10 13:55:55 jl_jaramillo
* versions developed at Meudon in November 2008
*
* Revision 1.14 2008/11/17 08:30:01 jl_jaramillo
* Nicolas's changes on multipoles. Sign correction in outgoing shear expression
*
* Revision 1.13 2008/11/12 15:17:47 n_vasset
* Bug-hunting, and new definition for computation of the Ricci scalar
* (instead of Ricci tensor previously)
*
* Revision 1.12 2008/07/09 08:47:33 n_vasset
* new version for multipole calculation. Function zeta implemented.
* Revision 1.11 2008/06/04 12:31:23 n_vasset
* New functions multipole_mass and multipole_angu. first version.
*
* Revision 1.9 2006/09/07 08:39:45 j_novak
* Minor changes.
*
* Revision 1.8 2006/07/03 10:13:48 n_vasset
* More efficient method for calculation of ricci tensor. Adding of flag issphere
*
* Revision 1.4 2006/06/01 11:47:50 j_novak
* Memory error hunt.
*
* Revision 1.3 2006/06/01 08:18:16 n_vasset
* Further implementation of Spheroid class definitions
*
* $Header: /cvsroot/Lorene/C++/Source/App_hor/spheroid.C,v 1.21 2014/10/13 08:52:38 j_novak Exp $
*
*/
// C headers
#include <cmath>
// Lorene headers
#include "metric.h"
#include "spheroid.h"
#include "proto.h"
#include "param.h"
//---------------//
// Constructors //
//--------------//
namespace Lorene {
Spheroid::Spheroid(const Map_af& map, double radius):
h_surf(map),
jac2d(map, 2, COV, map.get_bvect_spher()),
proj(map, 2, COV, map.get_bvect_spher()),
qq(map, COV, map.get_bvect_spher()),
ss (map, COV, map.get_bvect_spher()),
ephi(map, CON, map.get_bvect_spher()),
qab(map.flat_met_spher()),
ricci(map),
hh(map, COV, map.get_bvect_spher()),
trk(map),
ll(map, COV, map.get_bvect_spher()),
jj(map, COV, map.get_bvect_spher()),
fff(map),
ggg(map),
zeta(map),
issphere(true)
{
// Itbl type (2) ;
// type.set(1) = CON ;
// type.set(2) = COV ;
// Tensor proj(map, 2, type, map.get_bvect_spher());
assert(radius > 1.e-15) ;
assert(map.get_mg()->get_nzone() == 1) ; // one domain
assert(map.get_mg()->get_nr(0) == 1) ; // angular grid
assert(map.get_mg()->get_type_r(0) == FIN) ; //considered as a shell
// Setting of real index types for jacobian and projector (first contravariant, second covariant)
jac2d.set_index_type(0) = CON ;
proj.set_index_type(0) = CON ;
jac2d.set_etat_zero() ;
h_surf = radius ;
ss.set_etat_zero();
ephi.set_etat_zero();
proj.set_etat_zero();
hh.set_etat_zero() ;
for (int i=1; i<=3; i++)
hh.set(i,i) = 2./radius ;
trk.set_etat_zero() ;
ll.set_etat_zero() ;
jj.set_etat_zero() ;
fff.set_etat_zero();
ggg.set_etat_zero();
zeta.set_etat_zero();
set_der_0x0() ;
}
Spheroid::Spheroid(const Scalar& h_in, const Metric& gamij, const Sym_tensor& Kij):
h_surf(h_in),
jac2d(h_in.get_mp(),2, COV, h_in.get_mp().get_bvect_spher()),
proj(h_in.get_mp(),2, COV, h_in.get_mp().get_bvect_spher()),
qq(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
ss (h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
ephi(h_in.get_mp(), CON, h_in.get_mp().get_bvect_spher()),
qab(h_in.get_mp().flat_met_spher()),
ricci(h_in.get_mp()),
hh(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
trk(h_in.get_mp()),
ll(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
jj(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
fff(h_in.get_mp()),
ggg(h_in.get_mp()),
zeta(h_in.get_mp()),
issphere(true) {
set_der_0x0() ;
const Map& map = h_in.get_mp() ; // The 2-d 1-domain map
const Map& map3 = Kij.get_mp(); // The 3-d map
const Mg3d& gri2d = *map.get_mg() ;
const Map_radial* mp_rad = dynamic_cast<const Map_radial*>(&Kij.get_mp()) ;
assert(mp_rad != 0x0) ;
const Mg3d& gri3d = *map3.get_mg();
assert(&gri2d == Kij.get_mp().get_mg()->get_angu_1dom()) ;
int np = gri2d.get_np(0) ;
int nt = gri2d.get_nt(0) ;
int nr3 = gri3d.get_nr(0) ;
int nt3 = gri3d.get_nt(0) ;
int np3 = gri3d.get_np(0) ;
int nz = gri3d.get_nzone() ;
assert( nt == nt3 ) ;
assert ( np == np3 );
assert ( nr3 != 1 );
Param pipo ;
double xi = 0. ;
int lz = 0 ;
if(nz >2){
lz =1;
}
// Setting of real index types forjacobian and projector(first contravariant, other covariant)
proj.set_index_type(0) = CON;
jac2d.set_index_type(0) = CON;
// Copy of the 2-dimensional h_surf to a 3_d h_surf (calculus commodity, in order to match grids)
Scalar h_surf3 (map3);
h_surf3.allocate_all();
h_surf3.std_spectral_base();
for (int f= 0; f<nz; f++)
for (int k=0; k<np3; k++)
for (int j=0; j<nt3; j++) {
for (int l=0; l<nr3; l++) {
h_surf3.set_grid_point(f,k,j,l) = h_surf.val_grid_point(0,k,j,0);
}
}
if (nz >2){
h_surf3.annule_domain(0);
h_surf3.annule_domain(nz - 1);
}
// h_surf3.std_spectral_base();
/* Computation of the jacobian projector linked to the mapping from the
spheroid to a coordinate sphere. All quantities will then be calculated
as from a real coordinate sphere
*/
Itbl type (2);
type.set(0) = CON ;
type.set(1) = COV ;
Tensor jac (Kij.get_mp(),2,type, Kij.get_mp().get_bvect_spher());
jac.set_etat_zero();
jac.std_spectral_base();
jac.set(1,1) = 1. ;
jac.set(2,2)= 1. ;
jac.set(3,3) = 1. ;
jac.set(1,2) = -h_surf3.srdsdt() ;
jac.set(1,3) = -h_surf3.srstdsdp() ;
jac.std_spectral_base() ;
// Copy on the 2-d grid
jac2d.allocate_all() ;
jac2d.std_spectral_base();
for (int l=1; l<4; l++)
for (int m=1; m<4; m++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk((h_surf.val_grid_point(0, k, j, 0))*1.0000000000001,
j, k, pipo, lz, xi) ;
jac2d.set(l,m).set_grid_point(0, k, j, 0) =
jac(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Inverse jacobian (on 3-d grid)
Tensor jac_inv = jac ;
jac_inv.set(1,2) = - jac_inv(1,2);
jac_inv.set(1,3) = - jac_inv(1,3) ;
// Scalar field which annulation characterizes the 2-surface
Scalar carac = Kij.get_mp().r - h_surf3;
carac.std_spectral_base();
// Computation of the normal vector (covariant form) on both grids
Vector ss3d= carac.derive_cov(gamij) ;
Vector ss3dcon= carac.derive_con(gamij) ;
Scalar ssnorm = contract (ss3d.up(0, gamij), 0, ss3d, 0);
ssnorm.std_spectral_base() ;
ss3d = ss3d / sqrt (ssnorm) ;
ss3dcon = ss3dcon / sqrt (ssnorm) ;
if (nz >2){
ss3d.annule_domain(0);
ss3dcon.annule_domain(0);
ss3d.annule_domain(nz-1);
ss3dcon.annule_domain(nz -1);
}
ss3d.std_spectral_base();
ss3dcon.std_spectral_base();
// Provisory handling of dzpuis problems
if (nz >2){
h_surf3.annule_domain(nz-1);
}
for (int ii=1; ii <=3; ii++){
ss3d.set(ii).dec_dzpuis(ss3d(ii).get_dzpuis());
ss3dcon.set(ii).dec_dzpuis(ss3dcon(ii).get_dzpuis());
}
for (int ii=1; ii <=3; ii++)
for (int jjy = 1; jjy <=3; jjy ++){
jac_inv.set(ii, jjy).dec_dzpuis(jac_inv(ii, jjy).get_dzpuis());
jac.set(ii, jjy).dec_dzpuis(jac(ii, jjy).get_dzpuis());
}
// End of dzpuis handling.
Sym_tensor sxss3d = ss3d * ss3d ;
Sym_tensor sxss3dcon = ss3dcon * ss3dcon ;
Vector ss3 (Kij.get_mp(), COV, Kij.get_mp().get_bvect_spher()) ;
Vector ss3con(Kij.get_mp(), CON, Kij.get_mp().get_bvect_spher()) ;
ss3 = contract(jac_inv, 0, ss3d, 0);
ss.allocate_all() ;
ss.std_spectral_base();
for (int l=1; l<4; l++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.0000000000001, j, k, pipo,
lz, xi) ;
ss.set(l).set_grid_point(0, k, j, 0) =
ss3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// The vector field associated with Kiling conformal symmetry
Vector ephi3(gamij.get_mp(), CON, gamij.get_mp().get_bvect_spher());
ephi3.set(1).annule_hard();
ephi3.set(2).annule_hard();
Scalar ephi33(gamij.get_mp()); ephi33 = 1.; ephi33.std_spectral_base();
ephi33.mult_r(); ephi33.mult_sint();
ephi3.set(3) = ephi33;
ephi3.std_spectral_base();
ephi3 = contract (jac, 1, ephi3,0);
ephi.allocate_all() ;
ephi.std_spectral_base();
for (int l=1; l<4; l++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
j, k, pipo, lz, xi) ;
ephi.set(l).set_grid_point(0, k, j, 0) =
ephi3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Computation of the 3-d projector on the 2-sphere
Tensor proj3d(Kij.get_mp(),2, type, Kij.get_mp().get_bvect_spher());
Tensor proj_prov = gamij.con().down(1, gamij) - ss3dcon*ss3d;
proj.allocate_all();
proj.std_spectral_base();
proj3d.allocate_all();
proj3d = contract(jac, 1, contract( jac_inv, 0, proj_prov , 1), 1 );
proj3d.std_spectral_base();
for (int l=1; l<4; l++)
for (int m=1; m<4; m++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.00000000000001,
j, k, pipo, lz, xi) ;
proj.set(l,m).set_grid_point(0, k, j, 0) =
proj3d(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
/* Computation of the metric linked to the 2-surface (linked to covariant form
of the degenerated 2-metric).
It will be designed as a block-diagonal 3-metric, with 1 for
the first coordinate and the concerned 2-d metric as a
second block */
Sym_tensor qq3d (Kij.get_mp(), COV, Kij.get_mp().get_bvect_spher());
Sym_tensor qab2 (h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher());
qq3d.set_etat_zero();
qq3d.set(1,1) = 1;
qq3d.set(2,2)= gamij.cov()(1,1) * (h_surf3.srdsdt())* (h_surf3.srdsdt())
+ 2* gamij.cov()(1,2)* (h_surf3.srdsdt()) + gamij.cov()(2,2);
qq3d.set(3,3)= gamij.cov()(1,1)* (h_surf3.srstdsdp())*(h_surf3.srstdsdp())
+2*gamij.cov()(1,3)* (h_surf3.srstdsdp()) +gamij.cov()(3,3);
qq3d.set(2,3)= gamij.cov()(1,1)* (h_surf3.srdsdt())* (h_surf3.srstdsdp())+
gamij.cov()(1,2)* (h_surf3.srstdsdp())+gamij.cov()(1,3)* (h_surf3.srdsdt())
+ gamij.cov()(2,3) ;
qq3d.set(3,2)= gamij.cov()(1,1)* (h_surf3.srdsdt())* (h_surf3.srstdsdp())+
gamij.cov()(1,2)* (h_surf3.srstdsdp())+gamij.cov()(1,3)* (h_surf3.srdsdt())
+ gamij.cov()(2,3) ;
qq3d.std_spectral_base();
qab2.allocate_all() ;
qab2.std_spectral_base();
for (int l=1; l<4; l++)
for (int m=1; m<4; m++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
j, k, pipo, lz, xi) ;
qab2.set(l,m).set_grid_point(0, k, j, 0) =
qq3d(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
qab= qab2;
// Computation of the degenerated 3d degenerated covariant metric on the 2-surface
Sym_tensor qqq = contract(jac_inv, 0, contract( jac_inv, 0, (gamij.cov()
- ss3d * ss3d) , 0), 1) ;
qq.allocate_all() ;
qq.std_spectral_base();
for (int l=1; l<4; l++)
for (int m=1; m<4; m++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
j, k, pipo, lz, xi) ;
qq.set(l,m).set_grid_point(0, k, j, 0) =
qqq(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Computation of the trace of the extrinsic curvature of 3-slice
Scalar trk3d = Kij.trace(gamij) ;
trk.allocate_all() ;
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.00000000000001,
j, k, pipo, lz, xi) ;
trk.set_grid_point(0, k, j, 0) =
trk3d.get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Computation of the normalization factor of the outgoing null vector.
Scalar fff3d (map3);
fff3d = 1. ;
fff.allocate_all() ;
fff3d.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001, j, k, pipo,
lz, xi) ;
fff.set_grid_point(0, k, j, 0) =
fff3d.get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Computation of the normalization factor of the ingoing null vector.
Scalar ggg3d (map3);
ggg3d = 1. ;
ggg.allocate_all() ;
ggg3d.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.00000000000001, j, k, pipo,
lz, xi) ;
ggg.set_grid_point(0, k, j, 0) =
ggg3d.get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
//Computation of function zeta, changing spheroid coordinates to get a round measure.
Scalar ftilde = sqrt_q();
double rayon = sqrt(area()/(4.*M_PI));
ftilde = -ftilde/(rayon*rayon);
ftilde = ftilde*h_surf*h_surf;
ftilde.set_spectral_va().ylm();
Base_val base = ftilde.get_spectral_base() ;
Mtbl_cf *coefftilde = ftilde.get_spectral_va().c_cf;
int nombre = 2*nt; // ### Doubled in SYM base!!
double *a_tilde = new double[nombre];
lz = 0; // Now we work with 2d map associated with sqrt(q)
for (int k=0; k<np+1; k++)
for (int j=0; j<nt; j++) {
int l_q, m_q, base_r ;
base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
if (nullite_plm(j, nt, k, np, base) == 1) {
donne_lm(1, lz, j, k, base, m_q, l_q,base_r) ;
if (m_q ==0) {
a_tilde[l_q] = (*coefftilde)(0, k, j, 0);
}
}
}
Scalar zeta2(map); zeta2= 3.; zeta2.std_spectral_base(); zeta2.mult_cost();
zeta2.set_spectral_va().ylm();
Base_val base2 = zeta2.get_spectral_base() ;
Mtbl_cf *dzeta = zeta2.set_spectral_va().c_cf;
for (int k=0; k<np+1; k++)
for (int j=0; j<nt; j++) {
int l_q2, m_q2, base_r2 ;
base.give_quant_numbers(lz, k, j, m_q2, l_q2, base_r2) ;
if (nullite_plm(j, nt, k, np, base2) == 1) {
donne_lm(1, lz, j, k, base2, m_q2, l_q2,base_r2) ;
if (m_q2 ==0){
if(l_q2 ==0){
(*dzeta).set(0,k,j,0) = a_tilde[0] + (a_tilde[1]/3.)*sqrt(2.*1. + 1.);
}
if (l_q2 >0){
(*dzeta).set(0,k,j,0) =
(a_tilde[l_q2 +1]/(2.*l_q2 + 3.))*sqrt((2.*(l_q2 +1.)+1.)/(2.*l_q2 + 1.))
- (a_tilde[l_q2 -1]/(2.*l_q2 - 1.))
*sqrt((2.*(l_q2 - 1.)+1.)/(2.*l_q2 + 1.));
}
}
}
}
zeta2.set_spectral_va().coef();
zeta = zeta2;
/* Computation of the tangent part of the extrinsic curvature of
* the 2 surface embedded in the 3 slice. The reference vector
used is the vector field s */
Vector ll3d = contract( proj_prov, 0, contract(Kij, 1, ss3dcon, 0), 0) ;
Vector ll3 = contract( jac_inv, 0, ll3d, 0) ;
ll.allocate_all() ;
ll.std_spectral_base();
for (int l=1; l<4; l++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000001, j, k, pipo,
lz, xi) ;
ll.set(l).set_grid_point(0, k, j, 0) =
ll3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
/* computation of the tangent components of the extrinsic curvature
*of the 3-slice
*(extracted from the curvature of the timeslice)
* Note: this is not the actual 2d_ extrinsic curvature, but the
*tangent part of the time-slice extrinsic curvature */
Tensor jj3d = Kij - ss3d*ll3d - ll3d*ss3d
- contract(Kij, 0 , 1, sxss3dcon , 0, 1)* sxss3d ;
Tensor jj3 =contract(jac_inv, 0 , contract(jac_inv,0 , jj3d,1),1);
jj.allocate_all() ;
jj.std_spectral_base();
for (int l=1; l<4; l++)
for (int m=1; m<4; m++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
j, k, pipo, lz, xi) ;
jj.set(l,m).set_grid_point(0, k, j, 0) =
jj3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Computation of 2d extrinsic curvature in the 3-slice
Tensor hh3d = contract(proj_prov, 0,
contract(proj_prov, 0,ss3d.derive_cov(gamij),1), 1 ) ;
Tensor hh3 =contract(jac_inv, 0 , contract(jac_inv,0 , hh3d,1),1);
hh.allocate_all() ;
hh.std_spectral_base();
for (int l=1; l<4; l++)
for (int m=1; m<4; m++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
j, k, pipo, lz, xi) ;
hh.set(l,m).set_grid_point(0, k, j, 0) =
hh3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
// Computation of 2d ricci scalar
Tensor hh3dupdown = hh3d.up(0, gamij);
Scalar ricciscal3 = gamij.ricci().trace(gamij);
if (nz>2){
// Ricci scalar on the 3-surface.
ricciscal3.annule_domain(nz-1); ricciscal3.std_spectral_base();
}
Tensor hh3dupup = hh3dupdown.up(1,gamij);
if (nz>2){
hh3dupup.annule_domain(nz-1); hh3dupup.std_spectral_base();
}
Scalar ricci22 = ricciscal3
- 2.*contract(contract(gamij.ricci(), 1, ss3dcon, 0),0, ss3dcon, 0);
if (nz >2){
ricci22.annule_domain(nz-1);
ricci22.std_spectral_base();
}
ricci22 += (hh3d.trace(gamij)*hh3d.trace(gamij))
- contract(contract(hh3dupup,0, hh3d,0),0,1);
if (nz >2){
ricci22.annule_domain(nz-1);
}
ricci22.std_spectral_base();
ricci.allocate_all();
ricci.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001, j, k, pipo,
lz, xi) ;
ricci.set_grid_point(0, k, j, 0) =
ricci22.get_spectral_va().val_point_jk(lz, xi, j, k) ;
}
del_deriv() ; //### to be checked...
set_der_0x0() ;
delete [] a_tilde ;
}
//Copy constructor//
Spheroid::Spheroid (const Spheroid &sph_in) :h_surf(sph_in.h_surf),
jac2d(sph_in.jac2d),
proj(sph_in.proj),
qq(sph_in.qq),
ss (sph_in.ss),
ephi (sph_in.ephi),
qab(sph_in.qab),
ricci(sph_in.ricci),
hh(sph_in.hh),
trk(sph_in.trk),
ll(sph_in.ll),
jj(sph_in.jj),
fff(sph_in.fff),
ggg(sph_in.ggg),
zeta(sph_in.zeta),
issphere(sph_in.issphere)
{
set_der_0x0() ;
}
// Assignment to another Spheroid
void Spheroid::operator=(const Spheroid& sph_in)
{
h_surf = sph_in.h_surf ;
jac2d = sph_in.jac2d ;
proj = sph_in.proj ;
qq = sph_in.qq ;
ss = sph_in.ss ;
ephi = sph_in.ephi ;
qab = sph_in.qab ;
ricci = sph_in.ricci ;
hh = sph_in.hh ;
trk = sph_in.trk ;
ll = sph_in.ll ;
jj = sph_in.jj ;
fff = sph_in.fff ;
ggg = sph_in.ggg ;
zeta = sph_in.zeta ;
issphere = sph_in.issphere ;
del_deriv() ; // Deletes all derived quantities
}
//------------//
//Destructor //
//-----------//
Spheroid::~Spheroid()
{
del_deriv() ;
}
// -----------------//
// Memory management//
//------------------//
void Spheroid::del_deriv() const {
if (p_sqrt_q != 0x0) delete p_sqrt_q ;
if (p_area != 0x0) delete p_area ;
if (p_angu_mom != 0x0) delete p_angu_mom ;
if (p_mass != 0x0) delete p_mass ;
if (p_multipole_mass != 0x0) delete p_multipole_mass ;
if (p_multipole_angu != 0x0) delete p_multipole_angu ;
if (p_epsilon_A_minus_one != 0x0) delete p_epsilon_A_minus_one ;
if (p_epsilon_P_minus_one != 0x0) delete p_epsilon_P_minus_one ;
if (p_theta_plus != 0x0) delete p_theta_plus ;
if (p_theta_minus != 0x0) delete p_theta_minus ;
if (p_shear != 0x0) delete p_shear ;
if (p_delta != 0x0) delete p_delta ;
set_der_0x0() ;
}
void Spheroid::set_der_0x0() const {
p_sqrt_q = 0x0 ;
p_area = 0x0 ;
p_angu_mom = 0x0 ;
p_mass = 0x0 ;
p_multipole_mass = 0x0;
p_multipole_angu = 0x0;
p_epsilon_A_minus_one = 0x0;
p_epsilon_P_minus_one = 0x0;
p_theta_plus = 0x0 ;
p_theta_minus = 0x0 ;
p_shear = 0x0 ;
p_delta = 0x0;
}
//---------//
//Accessors//
//---------//
// Computation of the 2-dimensional Jacobian amplitude for the surface
const Scalar& Spheroid::sqrt_q() const {
if (p_sqrt_q == 0x0) {
p_sqrt_q = new Scalar(sqrt((get_qq()(2,2)*get_qq()(3,3))- (get_qq()(2,3)*get_qq()(2,3)))) ;
p_sqrt_q->std_spectral_base() ;
}
return *p_sqrt_q ;
}
// Computation of the 2-dimensional area of the surface
double Spheroid::area() const {
if (p_area == 0x0) {
const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
p_area = new double (mp_ang.integrale_surface((sqrt_q()) * h_surf *h_surf, 1)) ;
}
return *p_area;
}
// Computation of the angular momentum of the surface (G is set to be identically one)
double Spheroid::angu_mom() const {
if (p_angu_mom == 0x0) {
const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
Vector phi(mp_ang, CON, mp_ang.get_bvect_spher());
phi = get_ephi();
Scalar tmp = contract(ll,0, contract (jac2d, 1,phi,0), 0 );
p_angu_mom = new double (mp_ang.integrale_surface((sqrt_q()*h_surf*h_surf*tmp),1)) ;
*p_angu_mom = *p_angu_mom /(8. * M_PI) ;
}
return *p_angu_mom;
}
double Spheroid::mass() const {
if (p_mass == 0x0) {
double rayon = sqrt(area()/(4.*M_PI));
p_mass = new double ((1/(2.*rayon))*sqrt(rayon*rayon*rayon*rayon + 4.*angu_mom()*angu_mom()));
}
return *p_mass;
}
double Spheroid::multipole_mass(const int order) const{
const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
double rayon = sqrt(area()/(4.*M_PI));
// Multiplicative factor before integral.
double factor = mass()/(8. * M_PI); // To check later
if (order >0)
{ for (int compte=0; compte <=order -1; compte++)
factor = factor*rayon;
}
// Calculus of legendre polynomial of order n, as function of cos(theta)
Scalar Pn(mp_ang); Pn=1; Pn.std_spectral_base(); Pn.set_spectral_va().ylm();
if (order >0)
{ Pn = Pn*zeta;
}
if (order >1)
{ Scalar Pnold(mp_ang); Pnold = 1; Pnold.std_spectral_base(); Pnold.set_spectral_va().ylm();
for (int nn=1; nn<order; nn++){
Scalar Pnnew = (2.*nn +1.)*Pn;
Pnnew = Pnnew*zeta;
Pnnew = Pnnew - nn*Pnold;
Pnnew = Pnnew/(double(nn) + 1.);
Pnold = Pn;
Pn = Pnnew;
}
}
// Calculus of Ricci Scalar over the surface
Scalar ricciscal(mp_ang);
ricciscal = get_ricci();
ricciscal.set_spectral_va().ylm();
Scalar rayyon = h_surf;
rayyon.std_spectral_base();
rayyon.set_spectral_va().ylm();
Scalar sqq = sqrt_q();
Scalar integrande = sqq * rayyon *rayyon*ricciscal*Pn; integrande.std_spectral_base();
p_multipole_mass = new double (factor*mp_ang.integrale_surface(integrande, 1)) ;
return *p_multipole_mass;
}
double Spheroid::multipole_angu(const int order) const{
assert (order >=1) ;
const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
Vector phi(mp_ang, CON, mp_ang.get_bvect_spher());
phi = get_ephi();
double rayon = sqrt(area()/(4.*M_PI));
double factor = 1./(8. * M_PI);
if (order >1)
{ for (int compte=0; compte <=order -2; compte++)
factor = factor*rayon;
}
// Calculus of legendre polynomial of order n, as function of cos(theta)
Scalar Pn(mp_ang); Pn=1; Pn.std_spectral_base(); Pn.set_spectral_va().ylm();
Scalar dPn = Pn;
Pn = Pn*zeta;
if (order >1)
{ Scalar Pnold(mp_ang); Pnold = 1; Pnold.std_spectral_base(); Pnold.set_spectral_va().ylm();
for (int nn=1; nn<order; nn++){
Scalar Pnnew = (2.*nn +1.)*Pn;
Pnnew = Pnnew*zeta;
Pnnew = Pnnew - nn*Pnold;
Pnnew = Pnnew/(double(nn) + 1.);
Pnold = Pn;
Pn = Pnnew; // Pn is now P(n+1)
}
// Calculus of functional derivative of order N legendre polynomial.
dPn = Pn* zeta; dPn = dPn - Pnold; dPn = double(order)*dPn;
Scalar quotient(mp_ang); quotient = 1.; quotient.std_spectral_base();
quotient = quotient*zeta*zeta; quotient = quotient -1.;
dPn = dPn/quotient;
}
// Computation of the multipole;
Scalar tmp = contract(ll,0, contract (jac2d, 1,phi,0), 0 ); tmp.std_spectral_base();
Scalar tmp2 = (sqrt_q()) * h_surf *h_surf*tmp*dPn; tmp2.std_spectral_base();
p_multipole_angu = new double (factor*mp_ang.integrale_surface(tmp2, 1)) ;
return *p_multipole_angu;
}
// Computation of the refined Penrose parameter for axisymmetric spacetimes, and its difference wrt one.
double Spheroid::epsilon_A_minus_one() const {
if (p_epsilon_A_minus_one == 0x0) {
assert (pow(mass(),4) - pow (angu_mom(),2) > 0.);
p_epsilon_A_minus_one = new double(area()/(8.*M_PI*(mass()*mass() + sqrt(pow(mass(),4) - pow (angu_mom(),2)))) - 1.);
}
return *p_epsilon_A_minus_one;
}
// Computation of the classical Penrose parameter, and its difference wrt one.
// To use in replacement of epsilon_A_minus_one when the computed spacetime is not axisymmetric.
double Spheroid::epsilon_P_minus_one() const {
if (p_epsilon_P_minus_one == 0x0) {
assert (pow(mass(),4) - pow (angu_mom(),2) > 0.);
p_epsilon_P_minus_one = new double(area()/(16.*M_PI*mass()*mass()) - 1.);
}
return *p_epsilon_P_minus_one;
}
// Outgoing null expansion of 2-surface
const Scalar &Spheroid::theta_plus() const {
if (p_theta_plus == 0x0) {
p_theta_plus = new Scalar(fff*(hh.trace(qab) - jj.trace(qab))) ;
p_theta_plus->std_spectral_base() ;
p_theta_plus->set_spectral_va().ylm();
}
return *p_theta_plus;
}
// ingoing null expansion of 2-surface
const Scalar& Spheroid::theta_minus() const {
if (p_theta_minus == 0x0) {
p_theta_minus = new Scalar(ggg*(-hh.trace(qab) - jj.trace(qab))) ;
p_theta_minus->std_spectral_base() ;
}
return *p_theta_minus;
}
//outer null-oriented shear of 2-surface
const Sym_tensor& Spheroid::shear() const {
if (p_shear == 0x0) {
p_shear = new Sym_tensor( fff*(hh - jj) - (qab.cov()/2) *(hh.trace(qab) - jj.trace(qab))) ;
// This is associated with the null vector: "l = n + s";
// For a null vector, "l = f (n + s)", multiply by the global factor "f" (e.g. the lapse "N" for "l=N(n+s)".
p_shear->std_spectral_base() ;
}
return *p_shear;
}
//-------------------------------------------//
// Covariant flat derivative, returning a pointer.//
//-------------------------------------------//
Tensor Spheroid::derive_cov2dflat(const Tensor& uu) const{
// Notations: suffix 0 in name <=> input tensor
// suffix 1 in name <=> output tensor
int valence0 = uu.get_valence() ;
int valence1 = valence0 + 1 ;
int valence1m1 = valence1 - 1 ; // same as valence0, but introduced for
// the sake of clarity
// Protections
// -----------
if (valence0 >= 1) {
}
// Creation of the result (pointer)
// --------------------------------
Tensor *resu ;
// If uu is a Scalar, the result is a vector
if (valence0 == 0) {
resu = new Vector(uu.get_mp(), COV, uu.get_mp().get_bvect_spher()) ;
}
else {
// Type of indices of the result :
Itbl tipe(valence1) ;
const Itbl& tipeuu = uu.get_index_type() ;
for (int id = 0; id<valence0; id++) {
tipe.set(id) = tipeuu(id) ; // First indices = same as uu
}
tipe.set(valence1m1) = COV ; // last index is the derivation index
// if uu is a Tensor_sym, the result is also a Tensor_sym:
const Tensor* puu = &uu ;
const Tensor_sym* puus = dynamic_cast<const Tensor_sym*>(puu) ;
if ( puus != 0x0 ) { // the input tensor is symmetric
resu = new Tensor_sym(uu.get_mp(), valence1, tipe, *uu.get_triad(),
puus->sym_index1(), puus->sym_index2()) ;
}
else {
resu = new Tensor(uu.get_mp(), valence1, tipe, *uu.get_triad()) ; // no symmetry
}
}
int ncomp1 = resu->get_n_comp() ;
Itbl ind1(valence1) ; // working Itbl to store the indices of resu
Itbl ind0(valence0) ; // working Itbl to store the indices of uu
Itbl ind(valence0) ; // working Itbl to store the indices of uu
Scalar tmp(uu.get_mp()) ; // working scalar
// Determination of the dzpuis parameter of the result --> dz_resu
// --------------------------------------------------
int dz_resu = 0;
// (We only work here on a non-compactified shell) //
// Loop on all the components of the output tensor
// -----------------------------------------------
/* Note: we have here preserved all the non-useful terms in this case(typically christoffel symbols) for the sake of understandng what's going on...
*/
for (int ic=0; ic<ncomp1; ic++) {
// indices corresponding to the component no. ic in the output tensor
ind1 = resu->indices(ic) ;
// Component no. ic:
Scalar& cresu = resu->set(ind1) ;
// Indices of the input tensor
for (int id = 0; id < valence0; id++) {
ind0.set(id) = ind1(id) ;
}
// Value of last index (derivation index)
int k = ind1(valence1m1) ;
switch (k) {
case 1 : { // Derivation index = r
//---------------------
cresu = 0; //(uu(ind0)).dsdr() ; // d/dr
// all the connection coefficients Gamma^i_{jk} are zero for k=1
break ;
}
case 2 : { // Derivation index = theta
//-------------------------
cresu = (uu(ind0)).srdsdt() ; // 1/r d/dtheta
// Loop on all the indices of uu
for (int id=0; id<valence0; id++) {
switch ( ind0(id) ) {
case 1 : { // Gamma^r_{l theta} V^l
// or -Gamma^l_{r theta} V_l
/* ind = ind0 ;
ind.set(id) = 2 ; // l = theta
// Division by r :
tmp = uu(ind) ;
tmp.div_r_dzpuis(dz_resu) ;
cresu -= tmp ; */
break ;
}
case 2 : { // Gamma^theta_{l theta} V^l
// or -Gamma^l_{theta theta} V_l
/* ind = ind0 ;
ind.set(id) = 1 ; // l = r
tmp = uu(ind) ;
tmp.div_r_dzpuis(dz_resu) ;
cresu += tmp ; */
break ;
}
case 3 : { // Gamma^phi_{l theta} V^l
// or -Gamma^l_{phi theta} V_l
break ;
}
default : {
cerr << "Connection_fspher::p_derive_cov : index problem ! "
<< endl ;
abort() ;
}
}
}
break ;
}
case 3 : { // Derivation index = phi
//-----------------------
cresu = (uu(ind0)).srstdsdp() ; // 1/(r sin(theta)) d/dphi
// Loop on all the indices of uu
for (int id=0; id<valence0; id++) {
switch ( ind0(id) ) {
case 1 : { // Gamma^r_{l phi} V^l
// or -Gamma^l_{r phi} V_l
/* ind = ind0 ;
ind.set(id) = 3 ; // l = phi
tmp = uu(ind) ;
tmp.div_r_dzpuis(dz_resu) ;
cresu -= tmp ; */
break ;
}
case 2 : { // Gamma^theta_{l phi} V^l
// or -Gamma^l_{theta phi} V_l
ind = ind0 ;
ind.set(id) = 3 ; // l = phi
tmp = uu(ind) ;
tmp.div_r_dzpuis(dz_resu) ;
tmp.div_tant() ; // division by tan(theta)
cresu -= tmp ;
break ;
}
case 3 : { // Gamma^phi_{l phi} V^l
// or -Gamma^l_{phi phi} V_l
ind = ind0 ;
// ind.set(id) = 1 ; // l = r
// tmp = uu(ind) ;
// tmp.div_r_dzpuis(dz_resu) ;
// cresu += tmp ;
ind.set(id) = 2 ; // l = theta
tmp = uu(ind) ;
tmp.div_r_dzpuis(dz_resu) ;
tmp.div_tant() ; // division by tan(theta)
cresu += tmp ;
break ;
}
default : {
cerr << "Connection_fspher::p_derive_cov : index problem ! \n"
<< endl ;
abort() ;
}
}
}
break ;
}
default : {
cerr << "Connection_fspher::p_derive_cov : index problem ! \n" ;
abort() ;
}
} // End of switch on the derivation index
} // End of loop on all the components of the output tensor
// C'est fini !
// -----------
return *resu ;
}
void Spheroid::sauve(FILE* ) const {
cout << "c'est pas fait!" << endl ;
return ;
}
// Computation of the delta coefficients
const Tensor& Spheroid::delta() const {
if (p_delta == 0x0) {
Tensor christoflat(qab.get_mp(), 3, COV, qab.get_mp().get_bvect_spher());
christoflat.set_index_type(0) = CON;
christoflat.set_etat_zero() ;
// assert(flat_met != 0x0) ;
Tensor dgam = derive_cov2dflat(qab.cov()) ;
for (int k=1; k<=3; k++) {
for (int i=1; i<=3; i++) {
for (int j=1; j<=i; j++) {
Scalar& cc= christoflat.set(k,i,j);
for (int l=1; l<=3; l++) {
cc += qab.con()(k,l) * (
dgam(l,j,i) + dgam(i,l,j) - dgam(i,j,l) ) ;
}
cc = 0.5 * cc ;
}
}
}
p_delta = new Tensor (christoflat) ;
}
return *p_delta;
}
// Computation of global derivative on 2-sphere
Tensor Spheroid::derive_cov2d(const Tensor& uu) const {
if(uu.get_valence()>=1){
int nbboucle = uu.get_valence();
Tensor resu = derive_cov2dflat(uu);
for (int y=1; y<=nbboucle; y++){
int df = uu.get_index_type(y-1);
if (df == COV) {
resu -= contract(delta(),0, uu, y-1);
}
else {resu += contract(delta(),1, uu, y-1);}
return resu;
}
}
else return derive_cov2dflat(uu);
return derive_cov2dflat(uu); // to avoid warnings...
}
// // COmputation of the ricci tensor
// const Sym_tensor& Spheroid::ricci() const {
// if (p_ricci == 0x0) { // a new computation is necessary
// assert( issphere == true ) ;
// Sym_tensor riccia(h_surf.get_mp(), CON, h_surf.get_mp().get_bvect_spher()) ;
// riccia.set_etat_zero();
// const Tensor& d_delta = derive_cov2dflat(delta()) ;
// for (int i=1; i<=3; i++) {
// int jmax = 3 ;
// for (int j=1; j<=jmax; j++) {
// Scalar tmp1(h_surf.get_mp()) ;
// tmp1.set_etat_zero() ;
// for (int k=1; k<=3; k++) {
// tmp1 += d_delta(k,i,j,k) ;
// }
// Scalar tmp2(h_surf.get_mp()) ;
// tmp2.set_etat_zero() ;
// for (int k=1; k<=3; k++) {
// tmp2 += d_delta(k,i,k,j) ;
// }
// Scalar tmp3(h_surf.get_mp()) ;
// tmp3.set_etat_zero() ;
// for (int k=1; k<=3; k++) {
// for (int m=1; m<=3; m++) {
// tmp3 += delta()(k,k,m) * delta()(m,i,j) ;
// }
// }
// tmp3.dec_dzpuis() ; // dzpuis 4 -> 3
// Scalar tmp4(h_surf.get_mp()) ;
// tmp4.set_etat_zero() ;
// for (int k=1; k<=3; k++) {
// for (int m=1; m<=3; m++) {
// tmp4 += delta()(k,j,m) * delta()(m,i,k) ;
// }
// }
// tmp4.dec_dzpuis() ; // dzpuis 4 -> 3
// riccia.set(i,j) = tmp1 - tmp2 + tmp3 - tmp4 ;
// }
// }
// /* Note: Here we must take into account the fact that a round metric on a spheroid doesn't give zero as "flat" ricci part. Then a diagonal scalar term must be added.
// WARNING: this only works with "round" horizons!! */
// double rayon = sqrt(area()/(4.*M_PI));
// Scalar rayon2 = h_surf;
// rayon2 = rayon;
// rayon2.std_spectral_base();
// for (int hi=1; hi<=3; hi++){
// riccia.set(hi,hi) += 2/(rayon2 * rayon2) ; // Plutot 1/hsurf^2, non?
// }
// p_ricci = new Sym_tensor(riccia);
// }
// return *p_ricci ;
// }
// COmputation of the ricci tensor
// const Sym_tensor& Spheroid::ricci() const {
// if (p_ricci == 0x0) { // a new computation is necessary
// Sym_tensor riccia(h_surf.get_mp(), CON, h_surf.get_mp().get_bvect_spher()) ;
// Sym_tensor ricci3 = gamij.ricci();
// Sym_tensor ricci3-2d(h_surf.get_mp(), COV, h_surf.get_mp().get_bvect_spher());
// ricci3-2d.allocate_all() ;
// ricci3-2d.std_spectral_base();
// for (int l=1; l<4; l++)
// for (int m=1; m<4; m++)
// for (int k=0; k<np; k++)
// for (int j=0; j<nt; j++)
// {
// mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000001, j, k, pipo,
// lz, xi) ;
// ricci3-2d.set(l,m).set_grid_point(0, k, j, 0) =
// ricci3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
// }
// }
// return *p_ricci ;
// }
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
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::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::Map_af::integrale_surface
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Definition map_af_integ_surf.C:90
Lorene::Map_radial
Base class for pure radial mappings.
Definition map.h:1536
Lorene::Map_radial::val_lx_jk
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
Lorene::Metric
Metric for tensor calculation.
Definition metric.h:90
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
Lorene::Metric::get_mp
const Map & get_mp() const
Returns the mapping.
Definition metric.h:202
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Metric::ricci
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect ).
Definition metric.C:338
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::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
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::srdsdt
const Scalar & srdsdt() const
Returns of *this .
Definition scalar_deriv.C:145
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition scalar_deriv.C:390
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::mult_cost
void mult_cost()
Multiplication by .
Definition scalar_th_manip.C:58
Lorene::Scalar::srstdsdp
const Scalar & srstdsdp() const
Returns of *this .
Definition scalar_deriv.C:177
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::div_tant
void div_tant()
Division by .
Definition scalar_th_manip.C:111
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::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::mult_sint
void mult_sint()
Multiplication by .
Definition scalar_th_manip.C:84
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::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_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
Lorene::Scalar::derive_con
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition scalar_deriv.C:402
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::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:208
Lorene::Spheroid::derive_cov2dflat
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
Definition spheroid.C:954
Lorene::Spheroid::~Spheroid
virtual ~Spheroid()
Destructor.
Definition spheroid.C:642
Lorene::Spheroid::sqrt_q
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Definition spheroid.C:693
Lorene::Spheroid::proj
Tensor proj
The 3-d projector on the 2-surface (contravariant-covariant form).
Definition spheroid.h:100
Lorene::Spheroid::qq
Sym_tensor qq
The 3-d covariant degenerated 2-metric on the surface.
Definition spheroid.h:104
Lorene::Spheroid::theta_minus
const Scalar & theta_minus() const
Computes the ingoing null expansion .
Definition spheroid.C:909
Lorene::Spheroid::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition spheroid.C:1186
Lorene::Spheroid::ll
Vector ll
Normal-tangent component of the extrinsic curvature of the 3-slice.
Definition spheroid.h:129
Lorene::Spheroid::p_multipole_angu
double * p_multipole_angu
Angular momentum multipole for the spheroid.
Definition spheroid.h:162
Lorene::Spheroid::p_mass
double * p_mass
Mass defined from angular momentum.
Definition spheroid.h:160
Lorene::Spheroid::trk
Scalar trk
Trace K of the extrinsic curvature of the 3-slice.
Definition spheroid.h:124
Lorene::Spheroid::p_theta_plus
Scalar * p_theta_plus
Classical Penrose parameter, difference wrt one.
Definition spheroid.h:165
Lorene::Spheroid::ggg
Scalar ggg
Normalization function for the ingoing null vector k.
Definition spheroid.h:143
Lorene::Spheroid::operator=
void operator=(const Spheroid &)
Assignment to another Spheroid.
Definition spheroid.C:614
Lorene::Spheroid::epsilon_A_minus_one
double epsilon_A_minus_one() const
Computation of the refined Penrose parameter for axisymmetric spacetimes, and its difference wrt one.
Definition spheroid.C:864
Lorene::Spheroid::set_der_0x0
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition spheroid.C:666
Lorene::Spheroid::Spheroid
Spheroid(const Map_af &map, double radius)
The delta tensorial fields linked to Christoffel symbols.
Definition spheroid.C:94
Lorene::Spheroid::delta
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
Definition spheroid.C:1203
Lorene::Spheroid::multipole_angu
double multipole_angu(const int order) const
Computes the angular multipole of a given order for the spheroid, assumed to be spherical.
Definition spheroid.C:801
Lorene::Spheroid::issphere
bool issphere
Flag to know whether the horizon is geometrically round or distorted.
Definition spheroid.h:151
Lorene::Spheroid::jac2d
Tensor jac2d
The jacobian of the adaptation of coordinates (contravariant/covariant representation).
Definition spheroid.h:96
Lorene::Spheroid::mass
double mass() const
Computes the mass as defined from the calculus of angular momentum, done with respect to a divergence...
Definition spheroid.C:736
Lorene::Spheroid::jj
Sym_tensor jj
Tangent components of the extrinsic curvature of the 3-slice.
Definition spheroid.h:134
Lorene::Spheroid::p_area
double * p_area
The area of the 2-surface.
Definition spheroid.h:158
Lorene::Spheroid::fff
Scalar fff
Normalization function for the outgoing null vector l.
Definition spheroid.h:138
Lorene::Spheroid::area
double area() const
Computes the area of the 2-surface.
Definition spheroid.C:708
Lorene::Spheroid::multipole_mass
double multipole_mass(const int order) const
Computes the mass multipole of a given order for the spheroid, assumed to be spherical.
Definition spheroid.C:748
Lorene::Spheroid::h_surf
Scalar h_surf
The location of the 2-surface as r = h_surf .
Definition spheroid.h:91
Lorene::Spheroid::hh
Sym_tensor hh
The ricci scalar on the 2-surface.
Definition spheroid.h:122
Lorene::Spheroid::p_epsilon_P_minus_one
double * p_epsilon_P_minus_one
Refined Penrose parameter, difference wrt one.
Definition spheroid.h:164
Lorene::Spheroid::p_sqrt_q
Scalar * p_sqrt_q
Surface element .
Definition spheroid.h:157
Lorene::Spheroid::ss
Vector ss
The adapted normal vector field to spheroid in the 3-slice.
Definition spheroid.h:108
Lorene::Spheroid::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition spheroid.C:650
Lorene::Spheroid::get_ricci
const Scalar & get_ricci() const
Returns the 2-ricci scalar .
Definition spheroid.h:229
Lorene::Spheroid::shear
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
Definition spheroid.C:925
Lorene::Spheroid::p_theta_minus
Scalar * p_theta_minus
Null ingoing expansion.
Definition spheroid.h:166
Lorene::Spheroid::angu_mom
double angu_mom() const
Computes the angular momentum with respect to a divergence-free vector field tangent to the 2-surface...
Definition spheroid.C:721
Lorene::Spheroid::get_qq
const Sym_tensor & get_qq() const
returns the 3-d degenerate 2-metric
Definition spheroid.h:235
Lorene::Spheroid::ephi
Vector ephi
The conformal Killing vector field on the 2-surface (set to by default to the axial vector associated...
Definition spheroid.h:113
Lorene::Spheroid::get_ephi
const Vector & get_ephi() const
Returns the conformal Killing symmetry vector on the 2-surface.
Definition spheroid.h:253
Lorene::Spheroid::p_multipole_mass
double * p_multipole_mass
Mass multipole for the spheroid.
Definition spheroid.h:161
Lorene::Spheroid::p_angu_mom
double * p_angu_mom
The angular momentum.
Definition spheroid.h:159
Lorene::Spheroid::derive_cov2d
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
Definition spheroid.C:1240
Lorene::Spheroid::ricci
Scalar ricci
Induced metric on the 2-surface .
Definition spheroid.h:117
Lorene::Spheroid::p_shear
Sym_tensor * p_shear
The shear tensor.
Definition spheroid.h:167
Lorene::Spheroid::epsilon_P_minus_one
double epsilon_P_minus_one() const
Computation of the classical Penrose parameter, and its difference wrt one.
Definition spheroid.C:876
Lorene::Spheroid::theta_plus
const Scalar & theta_plus() const
Computes the outgoing null expansion .
Definition spheroid.C:889
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Tensor_sym
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Valeur::val_point_jk
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Definition valeur.C:900
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::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene::Tensor_sym::sym_index1
int sym_index1() const
Number of the first symmetric index (0<= id_sym1 < valence ).
Definition tensor.h:1149
Lorene::Tensor::get_index_type
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition tensor.h:886
Lorene::Tensor::allocate_all
virtual void allocate_all()
Performs the memory allocation of all the elements, down to the double arrays of the Tbl s.
Definition tensor.C:508
Lorene::Tensor_sym::sym_index2
int sym_index2() const
Number of the second symmetric index (id_sym1 < id_sym2 < valence ).
Definition tensor.h:1154
Lorene::Tensor::get_valence
int get_valence() const
Returns the valence.
Definition tensor.h:869
Lorene::Tensor::indices
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition tensor.C:539
Lorene::Tensor::set_index_type
int & set_index_type(int i)
Sets the type of the index number i .
Definition tensor.h:909
Lorene::Tensor::get_n_comp
int get_n_comp() const
Returns the number of stored components.
Definition tensor.h:872
Lorene::Tensor::down
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
Definition tensor_calculus.C:265
Lorene::Tensor::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Lorene::Tensor::up
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Definition tensor_calculus.C:225
Lorene::Tensor::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition tensor_calculus.C:94
Lorene::Tensor::set_etat_zero
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::srstdsdp
virtual void srstdsdp(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .
Lorene::get_bvect_spher
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Lorene::y
Coord y
y coordinate centered on the grid
Definition map.h:727
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Lorene::flat_met_spher
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Lorene::srdsdt
virtual void srdsdt(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .