27char 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__ ;
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
void del_t() const
Logical destructor (deletes the Mtbl member *c ).
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
const double * get_beta() const
Returns the pointer on the array beta.
const double * get_alpha() const
Returns the pointer on the array alpha.
Radial mapping of rather general form.
const Valeur & get_ff() const
Returns a (constant) reference to the function .
const double * get_alpha() const
Returns a pointer on the array alpha (values of in each domain)
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Tbl ** ddaa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
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.
virtual void operator=(const Map_et &mp)
Assignment to another Map_et
virtual const Map_af & mp_angu(int) const
Returns the "angular" mapping for the outside of domain l_zone.
virtual ostream & operator>>(ostream &) const
Operator >>
Tbl ** ddbb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Tbl zaasx2
Values at the nr collocation points of in the outermost compactified domain.
Coord rsx2drdx
in the nucleus and the shells; \ in the outermost compactified domain.
virtual bool operator==(const Map &) const
Comparison operator (egality)
Tbl bbsx
Values at the nr collocation points of in the nucleus.
void set_gg(const Valeur &)
Assigns a given value to the function .
const Valeur & get_gg() const
Returns a (constant) reference to the function .
Map_et(const Mg3d &mgrille, const double *r_limits)
Standard Constructor.
Tbl ** dbb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
void fait_poly()
Construction of the polynomials and .
virtual void homothetie(double lambda)
Sets a new radial scale.
Tbl aasx
Values at the nr collocation points of in the nucleus.
const double * get_beta() const
Returns a pointer on the array beta (values of in each domain)
Tbl aasx2
Values at the nr collocation points of in the nucleus.
Tbl zaasx
Values at the nr collocation points of in the outermost compactified domain.
Tbl ** daa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
void set_coord()
Assignement of the building functions to the member Coords.
Coord rsxdxdr
in the nucleus; \ in the shells; \ in the outermost compactified domain.
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Tbl ** aa
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Tbl ** bb
Array (size: mg->nzone ) of Tbl which stores the values of in each domain.
virtual ~Map_et()
Destructor.
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Tbl bbsx2
Values at the nr collocation points of in the nucleus.
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
virtual void sauve(FILE *) const
Save in a file.
void set_ff(const Valeur &)
Assigns a given value to the function .
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
virtual void reset_coord()
Resets all the member Coords.
Base class for pure radial mappings.
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
virtual void reset_coord()
Resets all the member Coords.
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord drdt
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED).
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
virtual void sauve(FILE *) const
Save in a file.
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord sr2d2rdt2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Coord stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified external domain (CED).
Base class for coordinate mappings.
double get_ori_z() const
Returns the z coordinate of the origin.
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
double ori_x
Absolute coordinate x of the origin.
Coord ya
Absolute y coordinate.
void set_rot_phi(double phi0)
Sets a new rotation angle.
double ori_y
Absolute coordinate y of the origin.
Coord r
r coordinate centered on the grid
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
double get_ori_y() const
Returns the y coordinate of the origin.
Coord za
Absolute z coordinate.
Coord tet
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.
Base_vect_cart bvect_cart
Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Coord z
z coordinate centered on the grid
double get_ori_x() const
Returns the x coordinate of the origin.
Map_af * p_mp_angu
Pointer on the "angular" mapping.
Coord x
x coordinate centered on the grid
Coord phi
coordinate centered on the grid
double ori_z
Absolute coordinate z of the origin.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Coord xa
Absolute x coordinate.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
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.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
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_etat_zero()
Sets the logical state to ETATZERO (zero).
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 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.
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
void annule_hard()
Sets the Valeur to zero in a hard way.
void sauve(FILE *) const
Save in a file.
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
#define P_COSSIN
dev. standart
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
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.
Standard units of space, time and mass.