82 #include "utilitaires.h" 83 #include "graphique.h" 89 assert ((relax >0) && (relax<=1)) ;
91 cout <<
"-----------------------------------------------" << endl ;
92 cout <<
"Resolution LAPSE" << endl ;
98 Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
101 Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
119 dirichlet_binaire (source_un, source_deux, -1., -1.,
hole1.
n_auto.
set(),
139 assert ((relax>0) && (relax<=1)) ;
141 cout <<
"-----------------------------------------------" << endl ;
142 cout <<
"Resolution PSI" << endl ;
148 Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
151 Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
165 for (
int k=0 ; k<np_un ; k++)
166 for (
int j=0 ; j<nt_un ; j++)
168 lim_un.std_base_scal() ;
174 for (
int k=0 ; k<np_deux ; k++)
175 for (
int j=0 ; j<nt_deux ; j++)
177 lim_deux.std_base_scal() ;
180 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
201 cout <<
"------------------------------------------------" << endl ;
202 cout <<
"Resolution shift : Omega = " <<
omega << endl ;
208 if (source_un.
get_etat() == ETATZERO) {
210 for (
int i=0 ; i<3 ; i++)
218 if (source_deux.
get_etat() == ETATZERO) {
220 for (
int i=0 ; i<3 ; i++)
226 for (
int i=0 ; i<3 ; i++) {
227 if (source_un(i).get_etat() != ETATZERO)
229 if (source_deux(i).get_etat() != ETATZERO)
235 assert ((orientation_un==0) || (orientation_un == M_PI)) ;
238 assert ((orientation_deux==0) || (orientation_deux == M_PI)) ;
240 int aligne_un = (orientation_un == 0) ? 1 : -1 ;
241 int aligne_deux = (orientation_deux == 0) ? 1 : -1 ;
246 for (
int i=0 ; i<3 ; i++) {
247 if (source_un(i).
get_etat() == ETATQCQ)
249 if (source_deux(i).
get_etat() == ETATQCQ)
267 xa_mtbl_un = source_un.
get_mp()->
xa ;
268 ya_mtbl_un = source_un.
get_mp()->
ya ;
275 x_mtbl_un = source_un.
get_mp()->
x ;
276 y_mtbl_un = source_un.
get_mp()->
y ;
285 for (
int k=0 ; k<np_un ; k++)
286 for (
int j=0 ; j<nt_un ; j++)
287 lim_x_un.set(0, k, j, 0) = aligne_un*
omega*ya_mtbl_un(0, 0, 0, 0) + aligne_un*
hole1.
omega_local*y_mtbl_un(1,k,j,0) ;
288 lim_x_un.base = *bases_un[0] ;
293 for (
int k=0 ; k<np_un ; k++)
294 for (
int j=0 ; j<nt_un ; j++)
295 lim_y_un.set(0, k, j, 0) = -aligne_un*
omega*xa_mtbl_un(0, 0, 0, 0) - aligne_un*
hole1.
omega_local*x_mtbl_un(1,k,j,0) ;;
296 lim_y_un.base = *bases_un[1] ;
300 for (
int k=0 ; k<np_un ; k++)
301 for (
int j=0 ; j<nt_un ; j++)
302 lim_z_un.
set(0, k, j, 0) = 0 ;
303 lim_z_un.base = *bases_un[2] ;
314 xa_mtbl_deux = source_deux.
get_mp()->
xa ;
315 ya_mtbl_deux = source_deux.
get_mp()->
ya ;
322 x_mtbl_deux = source_deux.
get_mp()->
x ;
323 y_mtbl_deux = source_deux.
get_mp()->
y ;
328 for (
int k=0 ; k<np_deux ; k++)
329 for (
int j=0 ; j<nt_deux ; j++)
330 lim_x_deux.set(0, k, j, 0) = aligne_deux*
omega*ya_mtbl_deux(0, 0, 0, 0) + aligne_deux*
hole2.
omega_local*y_mtbl_deux(1,k,j,0) ;
331 lim_x_deux.base = *bases_deux[0] ;
336 for (
int k=0 ; k<np_deux ; k++)
337 for (
int j=0 ; j<nt_deux ; j++)
338 lim_y_deux.set(0, k, j, 0) = -aligne_deux*
omega*xa_mtbl_deux(0, 0, 0, 0) - aligne_deux*
hole2.
omega_local*x_mtbl_deux(1,k,j,0) ;
339 lim_y_deux.base = *bases_deux[1] ;
343 for (
int k=0 ; k<np_deux ; k++)
344 for (
int j=0 ; j<nt_deux ; j++)
345 lim_z_deux.
set(0, k, j, 0) = 0 ;
346 lim_z_deux.base = *bases_deux[2] ;
348 for (
int i=0 ; i<3 ; i++) {
350 delete bases_deux[i] ;
353 delete [] bases_deux ;
359 poisson_vect_binaire (1./3., source_un, source_deux,
360 lim_x_un, lim_y_un, lim_z_un,
361 lim_x_deux, lim_y_deux, lim_z_deux,
364 for (
int i=0 ; i<3 ; i++) {
371 (1-relax)*shift_un_old ;
373 (1-relax)*shift_deux_old ;
380 cout <<
"Difference relatives : " << diff_un <<
" " << diff_deux << endl ;
Coord xa
Absolute x coordinate.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Tenseur grad_n_tot
Total gradient of N .
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void solve_shift(double precis, double relax)
Solves the equation for the shift, using the Oohara-Nakarmure scheme~: using the current for the ...
void set_std_base()
Set the standard spectal basis of decomposition for each component.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
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...
double omega_local
local angular velocity
Tenseur psi_auto
Part of generated by the hole.
Values and coefficients of a (real-value) function.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Tenseur_sym taij_tot
Total , which must be zero on the horizon of the regularisation on the shift has been done...
void fait_n_comp(const Bhole &comp)
Imports the part of N due to the companion hole comp .
Tenseur shift_auto
Part of generated by the hole.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
double omega
Position of the axis of rotation.
const Map * get_mp() const
Returns pointer on the mapping.
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.
Bhole hole1
Black hole one.
void solve_psi(double precis, double relax)
Solves the equation for the conformal factor~: The fields are the total values excpet those with s...
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).
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Tenseur_sym tkij_tot
Total .
Coord ya
Absolute y coordinate.
Bases of the spectral expansions.
void solve_lapse(double precis, double relax)
Solves the equation for the lapse~: The fields are the total values excpet those with subscript ...
Coord y
y coordinate centered on the grid
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Coord x
x coordinate centered on the grid
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
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).
void fait_psi_comp(const Bhole &comp)
Imports the part of due to the companion hole comp .
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Tensor handling *** DEPRECATED : use class Tensor instead ***.
double regul
Intensity of the correction on the shift vector.
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Bhole hole2
Black hole two.