27 char map_et_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Map/map_et.C,v 1.16 2014/10/13 08:53:03 j_novak Exp $" ;
138 #include "utilitaires.h"
151 aasx( mgrille.get_nr(0) ),
152 aasx2( mgrille.get_nr(0) ),
153 zaasx( mgrille.get_nr(mgrille.get_nzone()-1) ),
154 zaasx2( mgrille.get_nr(mgrille.get_nzone()-1) ),
155 bbsx( mgrille.get_nr(0) ),
156 bbsx2( mgrille.get_nr(0) ),
157 ff(mgrille.get_angu()) ,
158 gg(mgrille.get_angu())
172 alpha =
new double[nzone] ;
173 beta =
new double[nzone] ;
175 for (
int l=0 ; l<nzone ; l++) {
178 alpha[l] = bornes[l+1] - bornes[l] ;
179 beta[l] = bornes[l] ;
184 alpha[l] = (bornes[l+1] - bornes[l]) * .5 ;
185 beta[l] = (bornes[l+1] + bornes[l]) * .5 ;
190 double umax = double(1) / bornes[l] ;
191 double umin = double(1) /bornes[l+1] ;
192 alpha[l] = (umin - umax) *
double(0.5) ;
193 beta[l] = (umin + umax) *
double(0.5) ;
198 cout <<
"Map_et::Map_et: unkown type_r ! " << endl ;
225 aasx(grille.get_nr(0) ),
226 aasx2(grille.get_nr(0) ),
227 zaasx(grille.get_nr(grille.get_nzone()-1) ),
228 zaasx2(grille.get_nr(grille.get_nzone()-1) ),
229 bbsx(grille.get_nr(0) ),
230 bbsx2(grille.get_nr(0) ),
231 ff(grille.get_angu()) ,
232 gg(grille.get_angu()) {
238 Map_et mapping (grille, r_lim) ;
244 int np = grille.
get_np(0) ;
245 int nt = grille.
get_nt(0) ;
247 double * cf =
new double [nt*(np+2)] ;
248 for (
int k=0 ; k<np ; k++)
249 for (
int j=0 ; j<nt ; j++)
250 cf[k*nt+j] = S_0(k,j) - S_0(0,0) ;
252 int* deg =
new int [3] ;
257 int* dim =
new int [3] ;
262 Tbl ff_nucleus (np,nt) ;
265 Tbl gg_nucleus (np,nt) ;
274 double * coloc_even ;
278 cfpcossin (deg,dim,cf) ;
281 odd =
new double [nt*(np+2)] ;
282 even =
new double [nt*(np+2)] ;
285 for (
int k=0 ; k<np+2 ; k++)
286 if ((k%4 == 0) || (k%4==1))
287 for (
int j=0 ; j<nt ; j++) {
289 even[k*nt+j] = cf[k*nt+j] ;
292 if ((k%4 == 2) || (k%4 == 3))
293 for (
int j=0 ; j<nt ; j++) {
295 odd[k*nt+j] = cf[k*nt+j] ;
299 cout <<
"Erreur bizzare..." << endl ;
303 coloc_odd =
new double [nt*np] ;
304 coloc_even =
new double [nt*np] ;
306 cipcossin (deg,dim,deg,odd,coloc_odd) ;
307 cipcossin (deg,dim,deg,even,coloc_even) ;
308 for (
int k=0 ; k<np ; k++)
309 for (
int j=0 ; j<nt ; j++) {
310 gg_nucleus.
set(k,j) = coloc_even[k*nt+j] ;
311 ff_nucleus.
set(k,j) = coloc_odd[k*nt+j] ;
316 delete [] coloc_even ;
317 delete [] coloc_odd ;
324 cout <<
"Base_p != P_COSSIN not implemented in Map_et constructor" <<
329 double mu_nucleus = -
min(gg_nucleus) ;
330 double alpha_nucleus = S_0(0,0)-mu_nucleus ;
332 ff_nucleus /= alpha_nucleus ;
333 gg_nucleus += mu_nucleus ;
334 gg_nucleus /= alpha_nucleus ;
337 Tbl ff_shell (np,nt) ;
339 ff_shell = S_0 - S_0(0,0) ;
341 double lambda_shell = -
max(ff_shell) ;
343 double R_ext = r_lim[2] ;
345 double beta_shell = (R_ext+S_0(0,0)-lambda_shell)/2. ;
346 double alpha_shell = (R_ext-S_0(0,0)+lambda_shell)/2. ;
348 ff_shell += lambda_shell ;
349 ff_shell /= alpha_shell ;
357 for (
int k=0 ; k<np ; k++)
358 for (
int j=0 ; j<nt ; j++) {
359 ff.
set(0,k,j,0) = ff_nucleus(k,j) ;
360 gg.
set(0,k,j,0) = gg_nucleus(k,j) ;
361 ff.
set(1,k,j,0) = ff_shell(k,j) ;
370 alpha =
new double[nz] ;
371 alpha[0] = alpha_nucleus ;
372 alpha[1] = alpha_shell ;
374 beta =
new double[nz] ;
376 beta[1] = beta_shell ;
377 for (
int i=2 ; i<nz ; i++) {
392 zaasx2( mpi.zaasx2 ),
406 alpha =
new double[nzone] ;
407 beta =
new double[nzone] ;
409 for (
int l=0 ; l<nzone ; l++) {
427 assert(l<mg->get_nzone()) ;
438 assert(l<mg->get_nzone()) ;
451 aasx( mgi.get_nr(0) ),
452 aasx2( mgi.get_nr(0) ),
453 zaasx( mgi.get_nr(mgi.get_nzone()-1) ),
454 zaasx2( mgi.get_nr(mgi.get_nzone()-1) ),
455 bbsx( mgi.get_nr(0) ),
456 bbsx2( mgi.get_nr(0) ),
457 ff(*(mgi.get_angu()), fich) ,
458 gg(*(mgi.get_angu()), fich)
466 alpha =
new double[nz] ;
467 beta =
new double[nz] ;
587 r.
set(
this, map_et_fait_r) ;
588 tet.
set(
this, map_et_fait_tet) ;
589 phi.
set(
this, map_et_fait_phi) ;
590 sint.
set(
this, map_et_fait_sint) ;
591 cost.
set(
this, map_et_fait_cost) ;
592 sinp.
set(
this, map_et_fait_sinp) ;
593 cosp.
set(
this, map_et_fait_cosp) ;
595 x.
set(
this, map_et_fait_x) ;
596 y.
set(
this, map_et_fait_y) ;
597 z.
set(
this, map_et_fait_z) ;
599 xa.
set(
this, map_et_fait_xa) ;
600 ya.
set(
this, map_et_fait_ya) ;
601 za.
set(
this, map_et_fait_za) ;
604 xsr.
set(
this, map_et_fait_xsr) ;
605 dxdr.
set(
this, map_et_fait_dxdr) ;
606 drdt.
set(
this, map_et_fait_drdt) ;
649 aa =
new Tbl*[nzone] ;
652 bb =
new Tbl*[nzone] ;
656 for (
int l=0; l<nzone; l++) {
658 aa[l] =
new Tbl(nr) ;
661 bb[l] =
new Tbl(nr) ;
682 for (
int i=0; i<
mg->
get_nr(0); i++) {
685 double x2 = x1 * x1 ;
686 double x3 = x1 * x2 ;
697 aa[0]->
set(i) = x2 * x2 * (3. - 2.*x2) ;
698 daa[0]->
set(i) = 12. * x3 * (1. + x1) * (1. - x1) ;
699 ddaa[0]->
set(i) = 12. *x2 *(3. - 5.* x2) ;
700 aasx.
set(i) = x3 * (3. - 2.*x2) ;
705 bb[0]->
set(i) = 0.5 * x3 * (5. - 3.* x2) ;
706 dbb[0]->
set(i) = 7.5 * x2 * (1. + x1) * (1. - x1) ;
707 ddbb[0]->
set(i) = 15. * x1 * ( 1. - 2.*x2 ) ;
708 bbsx.
set(i) = 0.5 * x2 * (5. - 3.* x2) ;
709 bbsx2.
set(i) = 0.5 * x1 * (5. - 3.* x2) ;
715 for (
int l=1; l<nzone; l++) {
727 for (
int i=0; i<
mg->
get_nr(l); i++) {
730 double xm1 = x1 - 1. ;
731 double xp1 = x1 + 1. ;
735 aa[l]->
set(i) = 0.25* xm1 * xm1 * (x1 + 2.) ;
736 daa[l]->
set(i) = 0.75* xm1 * xp1 ;
741 bb[l]->
set(i) = 0.25* xp1 * xp1 * (2. - x1) ;
742 dbb[l]->
set(i) = - 0.75* xm1 * xp1 ;
751 int nzm1 = nzone - 1 ;
757 for (
int i=0; i<
mg->
get_nr(nzm1); i++) {
760 zaasx.
set(i) = 0.25 * (x1 - 1.) * (2. + x1) ;
803 "Radial mapping of form r = xi + A(xi)F(t,p) + B(xi)G(t,p) (class Map_et)"
806 for (
int l=0; l<nz; l++) {
807 ost <<
" Domain #" << l <<
" : alpha_l = " <<
alpha[l]
808 <<
" , beta_l = " <<
beta[l] << endl ;
810 ost << endl <<
"Function F(theta', phi') : " << endl ;
811 ost <<
"------------------------- " << endl ;
813 ost << endl <<
"Function G(theta', phi') : " << endl ;
814 ost <<
"------------------------- " << endl ;
821 <<
"Values of r at the outer boundary of each domain [km] :" << endl ;
822 ost <<
"------------------------------------------------------" << endl ;
823 ost <<
" 1/ for theta = Pi/2 and phi = 0 : " << endl ;
825 for (
int l=0; l<nz; l++) {
826 ost <<
" " <<
val_r(l, 1., M_PI/2, 0) / km ;
830 if ( type_t == SYM ) {
831 assert( (type_p == SYM) || (type_p == NONSYM) ) ;
832 ost <<
" Coord r : " ;
833 for (
int l=0; l<nz; l++) {
836 ost <<
" " << (+
r)(l, 0, ntm1, nrm1) / km ;
841 ost <<
" 2/ for theta = Pi/2 and phi = Pi/2 : " << endl ;
843 for (
int l=0; l<nz; l++) {
844 ost <<
" " <<
val_r(l, 1., M_PI/2, M_PI/2) / km ;
848 if ( type_t == SYM ) {
849 ost <<
" Coord r : " ;
850 for (
int l=0; l<nz; l++) {
854 if ( (type_p == NONSYM) && (np % 4 == 0) ) {
855 ost <<
" " << (+
r)(l, np/4, ntm1, nrm1) / km ;
857 if ( type_p == SYM ) {
858 ost <<
" " << (+
r)(l, np/2, ntm1, nrm1) / km ;
864 ost <<
" 3/ for theta = Pi/2 and phi = Pi : " << endl ;
866 for (
int l=0; l<nz; l++) {
867 ost <<
" " <<
val_r(l, 1., M_PI/2, M_PI) / km ;
871 if ( (type_t == SYM) && (type_p == NONSYM) ) {
872 ost <<
" Coord r : " ;
873 for (
int l=0; l<nz; l++) {
877 ost <<
" " << (+
r)(l, np/2, ntm1, nrm1) / km ;
882 ost <<
" 4/ for theta = 0 : " << endl ;
884 for (
int l=0; l<nz; l++) {
885 ost <<
" " <<
val_r(l, 1., 0., 0.) / km ;
889 ost <<
" Coord r : " ;
890 for (
int l=0; l<nz; l++) {
892 ost <<
" " << (+
r)(l, 0, 0, nrm1) / km ;
909 for (
int l=0; l<nz; l++) {
933 cout <<
"Map_et::resize can be applied only to a shell !" << endl ;
939 double n_alpha = 0.5 * ( (lambda + 1.) *
alpha[l] +
940 (lambda - 1.) *
beta[l] ) ;
942 double n_beta = 0.5 * ( (lambda - 1.) *
alpha[l] +
943 (lambda + 1.) *
beta[l] ) ;
953 assert(l<mg->get_nzone()-1) ;
958 assert(
ff(lp1).get_etat() == ETATZERO ) ;
959 assert(
gg(lp1).get_etat() == ETATZERO ) ;
974 alpha[lp1] = n_alpha ;
987 double precis = 1e-10 ;
1006 for (
int i=0 ; i<nz ; i++) {
1009 if ((i!=0) && (i!=nz-1))
1046 const char* f = __FILE__ ;
Radial mapping of rather general form.
Coord z
z coordinate centered on the grid
void set_coord()
Assignement of the building functions to the member Coords.
void set_rot_phi(double phi0)
Sets a new rotation angle.
#define P_COSSIN
dev. standart
virtual ~Map_et()
Destructor.
double & set(int i)
Read/write of a particular element (index i) (1D case)
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
double get_ori_x() const
Returns the x coordinate of the origin.
Values and coefficients of a (real-value) function.
Tbl ** dbb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Tbl aasx
Values at the nr collocation points of in the nucleus.
virtual bool operator==(const Map &) const
Comparison operator (egality)
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Map_af * p_mp_angu
Pointer on the "angular" mapping.
Coord sr2d2rdt2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
void sauve(FILE *) const
Save in a file.
void set_gg(const Valeur &)
Assigns a given value to the function .
virtual void reset_coord()
Resets all the member Coords.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double get_ori_y() const
Returns the y coordinate of the origin.
const double * get_alpha() const
Returns a pointer on the array alpha (values of in each domain)
virtual void sauve(FILE *) const
Save in a file.
const Valeur & get_ff() const
Returns a (constant) reference to the function .
const double * get_beta() const
Returns a pointer on the array beta (values of in each domain)
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Coord r
r coordinate centered on the grid
Tbl ** bb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Base class for pure radial mappings.
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
void annule(int l)
Sets the Valeur to zero in a given domain.
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Coord drdt
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED).
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.
Standard units of space, time and mass.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Tbl bbsx
Values at the nr collocation points of in the nucleus.
const Valeur & get_gg() const
Returns a (constant) reference to the function .
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Tbl ** daa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
void annule_hard()
Sets the Valeur to zero in a hard way.
void affiche_seuil(ostream &ostr, int type=0, int precision=4, double threshold=1.e-7) const
Prints only the values greater than a given threshold.
Tbl ** ddbb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Coord stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED).
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
double ori_y
Absolute coordinate y of the origin.
Tbl aasx2
Values at the nr collocation points of in the nucleus.
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
double ori_x
Absolute coordinate x of the origin.
int get_nzone() const
Returns the number of domains.
Tbl bbsx2
Values at the nr collocation points of in the nucleus.
Coord tet
coordinate centered on the grid
virtual ostream & operator>>(ostream &) const
Operator >>
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord phi
coordinate centered on the grid
Base_vect_cart bvect_cart
Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Map_et(const Mg3d &mgrille, const double *r_limits)
Standard Constructor.
void set_ff(const Valeur &)
Assigns a given value to the function .
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Coord x
x coordinate centered on the grid
Tbl ** aa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Coord xa
Absolute x coordinate.
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
double get_ori_z() const
Returns the z coordinate of the origin.
const double * get_beta() const
Returns the pointer on the array beta.
Coord rsx2drdx
in the nucleus and the shells; \ in the outermost compactified domain.
virtual void sauve(FILE *) const
Save in a file.
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
virtual void homothetie(double lambda)
Sets a new radial scale.
double ori_z
Absolute coordinate z of the origin.
Tbl zaasx
Values at the nr collocation points of in the outermost compactified domain.
Tbl ** ddaa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
void del_t() const
Logical destructor (deletes the Mtbl member *c ).
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord y
y coordinate centered on the grid
Tbl zaasx2
Values at the nr collocation points of in the outermost compactified domain.
Coord za
Absolute z coordinate.
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
const double * get_alpha() const
Returns the pointer on the array alpha.
void fait_poly()
Construction of the polynomials and .
Base class for coordinate mappings.
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
Coord rsxdxdr
in the nucleus; \ in the shells; \ in the outermost compactified domain.
virtual const Map_af & mp_angu(int) const
Returns the "angular" mapping for the outside of domain l_zone.
virtual void reset_coord()
Resets all the member Coords.
double * beta
Array (size: mg->nzone ) of the values of in each domain.
virtual void operator=(const Map_et &mp)
Assignment to another Map_et
Coord ya
Absolute y coordinate.
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined