24char regularisation_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Isol_hor/regularisation.C,v 1.12 2014/10/13 08:53:01 j_novak Exp $" ;
89 const Vector& shift_comp_temp,
double om) {
91 Vector shift_auto(shift_auto_temp) ;
93 Vector shift_comp(shift_comp_temp) ;
95 Vector shift_old (shift_auto) ;
98 assert ((orientation==0) || (orientation == M_PI)) ;
100 assert ((orientation_autre==0) || (orientation_autre == M_PI)) ;
102 int alignes = (orientation == orientation_autre) ? 1 : -1 ;
110 shift_tot.
set(1).
import(alignes*shift_comp(1)) ;
111 shift_tot.
set(2).
import(alignes*shift_comp(2)) ;
114 shift_tot = shift_tot + shift_auto ;
116 double indic = (orientation == 0) ? 1 : -1 ;
120 for (
int i=1 ; i<=3 ; i++) {
125 tbi.
set(1) = *shift_tot(1).get_spectral_va().c - indic *om * shift_tot.
get_mp().
ya ;
126 tbi.
set(2) = *shift_tot(2).get_spectral_va().c + indic *om * shift_tot.
get_mp().
xa ;
133 for (
int i=1 ; i<=3 ; i++)
134 derive_r.
set(i) = tbi(i).dsdr() ;
144 for (
int comp=1 ; comp<=3 ; comp++) {
146 for (
int k=0 ; k<np ; k++)
147 for (
int j=0 ; j<nt ; j++)
148 for (
int i=0 ; i<nr ; i++)
149 val_hor.
set(1, k, j, i) = derive_r(comp).
150 val_grid_point(1, k, j, 0) ;
152 double r_0 = shift_auto.
get_mp().
val_r (1, -1, 0, 0) ;
153 double r_1 = shift_auto.
get_mp().
val_r (1, 1, 0, 0) ;
155 fonction_radiale =
pow(r_1-shift_auto.
get_mp().
r, 3.)*
156 (shift_auto.
get_mp().
r-r_0)/
pow(r_1-r_0, 3.) ;
157 fonction_radiale.
annule(0) ;
160 enleve = fonction_radiale * val_hor ;
162 get_spectral_va().get_base()) ;
164 if (
norme(enleve)(1) != 0)
165 shift_auto.
set(comp) = shift_auto(comp) - enleve ;
166 if (
norme(shift_auto(comp))(1) > 1e-5) {
168 if (erreur < diff(1))
191 for (
int i=1 ; i<=3 ; i++) {
196 for (
int i=1 ; i<=3 ; i++)
197 shift(i).get_spectral_va().coef_i() ;
200 tbi.
set(2) = *shift(2).get_spectral_va().c +
omega*
mp.
x ;
201 if (shift(3).get_etat() != ETATZERO)
202 tbi.
set(3) = *shift(3).get_spectral_va().c -
boost_z ;
213 for (
int i=1 ; i<=3 ; i++)
214 derive_r.
set(i) = tbi(i).dsdr() ;
225 double r_0 =
mp.
val_r(1, -1, 0, 0) ;
226 double r_1 =
mp.
val_r(1, 1, 0, 0) ;
228 for (
int comp=1 ; comp<=3 ; comp++) {
230 for (
int k=0 ; k<np ; k++)
231 for (
int j=0 ; j<nt ; j++)
232 for (
int i=0 ; i<nr ; i++)
233 val_hor.
set(1, k, j, i) = derive_r(comp).val_grid_point(1, k, j, 0) ;
235 fonction_radiale =
pow(r_1-
mp.
r, 3.)* (
mp.
r-r_0)/
pow(r_1-r_0, 3.) ;
236 fonction_radiale.
annule(0) ;
239 enleve = fonction_radiale*val_hor ;
242 Scalar copie (shift(comp)) ;
243 shift.
set(comp) = shift(comp)-enleve ;
246 assert (shift(comp).check_dzpuis(0)) ;
250 if (norm(1) > 1e-5) {
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
double regularise_one()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon.
double omega
Angular velocity in LORENE's units.
Evolution_std< Vector > beta_auto_evol
Values at successive time steps of the shift function .
double boost_z
Boost velocity in z-direction.
Map_af & mp
Affine mapping.
double regularisation(const Vector &shift_auto, const Vector &shift_comp, double ang_vel)
Corrects shift_auto in such a way that the total is equal to zero in the horizon,...
double boost_x
Boost velocity in x-direction.
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.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord y
y coordinate centered on the grid
Coord ya
Absolute y coordinate.
Coord r
r coordinate centered on the grid
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Coord x
x coordinate centered on the grid
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Coord xa
Absolute x coordinate.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tensor field of valence 0 (or component of a tensorial field).
Valeur & set_spectral_va()
Returns va (read/write version)
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
int jtime
Time step index of the latest slice.
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Evolution_std< double > the_time
Time label of each slice.
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Values and coefficients of a (real-value) function.
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
void annule(int l)
Sets the Valeur to zero in a given domain.
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
void coef_i() const
Computes the physical value of *this.
Base_val base
Bases on which the spectral expansion is performed.
void annule_hard()
Sets the Valeur to zero in a hard way.
Tensor field of valence 1.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
const Map & get_mp() const
Returns the mapping.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).