102 assert ((relax>0) && (relax<=1)) ;
104 cout <<
"Resolution LAPSE" << endl ;
113 for (
int i=0 ; i<3 ; i++) {
114 work.
set() = auxi(i, i) ;
129 Cmp soluce (source.poisson_dirichlet(limite, 0)) ;
130 soluce = soluce + 1 ;
133 n_auto.
set() = relax*soluce + (1-relax)*lapse_old ;
140 assert ((relax>0) && (relax<=1)) ;
142 cout <<
"Resolution PSI" << endl ;
150 for (
int i=0 ; i<3 ; i++) {
151 work.
set() = auxi(i, i) ;
164 for (
int k=0 ; k<np ; k++)
165 for (
int j=0 ; j<nt ; j++)
167 limite.std_base_scal() ;
170 soluce = soluce + 1 ;
173 psi_auto.
set() = relax*soluce + (1-relax)*psi_old ;
181 assert (precision > 0) ;
182 assert ((relax>0) && (relax<=1)) ;
184 cout <<
"resolution SHIFT" << endl ;
195 for (
int i=0 ; i<3 ; i++)
196 if (source(i).
get_etat() == ETATQCQ)
204 for (
int i=0 ; i<3 ; i++)
224 for (
int k=0 ; k<np ; k++)
225 for (
int j=0 ; j<nt ; j++)
227 lim_x.base = *bases[0] ;
231 for (
int k=0 ; k<np ; k++)
232 for (
int j=0 ; j<nt ; j++)
233 lim_y.
set(0, k, j, 0) = -
omega*x_mtbl(1, k, j, 0)-
boost[1] ;
234 lim_y.base = *bases[1] ;
238 for (
int k=0 ; k<np ; k++)
239 for (
int j=0 ; j<nt ; j++)
241 lim_z.base = *bases[2] ;
244 for (
int i=0 ; i<3 ; i++)
249 poisson_vect_frontiere(1./3., source,
shift_auto, lim_x, lim_y,
250 lim_z, 0, precision, 20) ;
261 for (
int i=0 ; i<3 ; i++) {
266 for (
int i=0 ; i<3 ; i++)
280 for (
int i=0 ; i<3 ; i++)
281 derive_r.set(i) = tbi(i).dsdr() ;
293 double r_0 =
mp.
val_r(1, -1, 0, 0) ;
294 double r_1 =
mp.
val_r(1, 1, 0, 0) ;
296 for (
int comp=0 ; comp<3 ; comp++) {
297 val_hor.annule_hard() ;
298 for (
int k=0 ; k<np ; k++)
299 for (
int j=0 ; j<nt ; j++)
300 for (
int i=0 ; i<nr ; i++)
301 val_hor.set(1, k, j, i) = derive_r(comp)(1, k, j, 0) ;
303 fonction_radiale =
pow(r_1-
mp.
r, 3.)* (
mp.
r-r_0)/
pow(r_1-r_0, 3.) ;
304 fonction_radiale.
annule(0) ;
305 fonction_radiale.annule(2, nz-1) ;
307 enleve = fonction_radiale*val_hor ;
318 if (norm(1) > 1e-5) {
333 Cmp lapse_non_sing (division_xpun(
n_auto(), 0)) ;
337 for (
int i=0 ; i<3 ; i++)
338 for (
int j=i ; j<3 ; j++) {
340 auxi = division_xpun (auxi, 0) ;
370 for (
int i=0 ; i<3 ; i++)
374 cout <<
"Calcul du moment non valable pour un boost != 0" << endl ;
383 for (
int i=0 ; i<3 ; i++)
386 Cmp integrant_un (vecteur_un(0)) ;
391 for (
int i=0 ; i<3 ; i++)
394 Cmp integrant_deux (vecteur_deux(0)) ;
417 for (
int i=0 ; i<3 ; i++)
421 cout <<
"Calcul du moment non valable pour un boost != 0" << endl ;
439 for (
int i=0 ; i<3 ; i++)
441 vecteur.set_std_base() ;
Coord xa
Absolute x coordinate.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
double moment_seul_hor() const
Calculates the angular momentum of the black hole using the formula on the horizon : where and H de...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void annule(int l)
Sets the Cmp to zero in a given domain.
const Valeur & mult_sp() const
Returns applied to *this.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
void coef_i() const
Computes the physical value of *this.
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
void solve_psi_seul(double relax)
Solves the equation for ~: with the condition that on the horizon, where f is the value of on the ...
void solve_lapse_seul(double relax)
Solves the equation for N ~: with the condition that N =0 on the horizon.
Tenseur psi_auto
Part of generated by the hole.
Values and coefficients of a (real-value) function.
void fait_taij_auto()
Calculates the part of generated by shift_auto .
Tenseur_sym taij_auto
Part of generated by the hole.
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Tenseur shift_auto
Part of generated by the hole.
void regularise_seul()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon...
void solve_shift_seul(double precis, double relax)
Solves the equation for ~: with on the horizon, being the boost and .
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
const Map * get_mp() const
Returns pointer on the mapping.
void fait_tkij()
Calculates the total .
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Map_af & mp
Affine mapping.
int get_nzone() const
Returns the number of domains.
const Valeur & mult_st() const
Returns applied to *this.
int get_etat() const
Returns the logical state.
double rayon
Radius of the horizon in LORENE's units.
Cmp pow(const Cmp &, int)
Power .
Tenseur_sym tkij_auto
Auto .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double * boost
Lapse on the horizon.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Coord ya
Absolute y coordinate.
Bases of the spectral expansions.
const Valeur & mult_cp() const
Returns applied to *this.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Coord y
y coordinate centered on the grid
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
double masse_komar_seul() const
Calculates the Komar mass of the black hole using : .
Coord x
x coordinate centered on the grid
Cmp poisson_neumann(const Valeur &, int) const
Idem as Cmp::poisson_dirichlet , the boundary condition being on the radial derivative of the solutio...
double omega
Angular velocity in LORENE's units.
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
double moment_seul_inf() const
Calculates the angular momentum of the black hole using the formula at infinity : where ...
Tenseur n_auto
Part of N generated by the hole.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double masse_adm_seul() const
Calculates the ADM mass of the black hole using : .
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Valeur va
The numerical value of the Cmp.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
double regul
Intensity of the correction on the shift vector.
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Coord r
r coordinate centered on the grid
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).