26 char gval_from_spectral_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Valencia/gval_from_spectral.C,v 1.14 2014/10/13 08:53:48 j_novak Exp $" ;
81 #include "grille_val.h"
82 #include "proto_f77.h"
93 if (taille != taille_in) {
94 cout <<
"Gval_spher::somme_spectral2():\n" ;
95 cout <<
"grid size incompatible with array size... exiting!" << endl ;
102 for (
int i=0; i<
nfantome; i++) resu[i] = 0 ;
107 for (
int i=nrv; i<taille; i++) resu[i] = 0 ;
115 int taille = nxv2*nzv2 ;
116 if (taille != taille_in) {
117 cout <<
"Gval_spher::somme_spectral2():\n" ;
118 cout <<
"grid size incompatible with array size... exiting!" << endl ;
123 double xi0, rr, theta ;
127 for (
int iz=0; iz<nzv2; iz++) {
132 for (
int ix=
nfantome; ix<nxv; ix++) {
137 double xx2 = (
x->
t[ix])*(
x->
t[ix]) ;
138 for (
int iz=
nfantome; iz<nzv; iz++) {
140 theta = (rr != 0. ?
acos((
zr->
t[iz])/rr) : 0) ;
141 mp.
val_lx(rr, theta, phi, l, xi0) ;
145 for (
int iz=nzv; iz<nzv2; iz++) {
150 for (
int ix=nxv; ix<nxv2; ix++) {
151 for (
int iz=0; iz<nzv2; iz++) {
165 int taille = nyv2*nxv2*nzv2 ;
166 if (taille != taille_in) {
167 cout <<
"Gval_spher::somme_spectral2():\n" ;
168 cout <<
"grid size incompatible with array size... exiting!" << endl ;
173 double xi0, rr, theta, phi ;
176 for (
int ix=0; ix<nxv2; ix++) {
177 for (
int iz=0; iz<nzv2; iz++){
183 for (
int iy=
nfantome; iy<nyv; iy++) {
184 double yy =
x->
t[iy] ;
187 for (
int iz=0; iz<nzv2; iz++) {
192 for (
int ix=
nfantome; ix<nxv; ix++) {
197 double xx =
x->
t[ix] ;
199 for (
int iz=
nfantome; iz<nzv; iz++) {
200 rr =
sqrt((
zr->
t[iz])*(
zr->
t[iz]) + xx2 + yy2) ;
201 theta = (rr != 0. ?
acos((
zr->
t[iz])/rr) : 0. );
202 phi = (rr != 0. ? atan2(yy, xx) : 0. ) ;
203 mp.
val_lx(rr, theta, phi, l, xi0) ;
207 for (
int iz=nzv; iz<nzv2; iz++) {
212 for (
int ix=nxv; ix<nxv2; ix++) {
213 for (
int iz=0; iz<nzv2; iz++) {
219 for (
int iy=nyv; iy<nyv2; iy++) {
220 for (
int ix=0; ix<nxv2; ix++) {
221 for (
int iz=0; iz<nzv2; iz++){
236 int taille = ntv2*nrv2 ;
237 if (taille != taille_in) {
238 cout <<
"Gval_spher::somme_spectral2():\n" ;
239 cout <<
"grid size incompatible with array size... exiting!" << endl ;
244 double xi, rr, theta ;
248 for (
int ir=0; ir<nrv2; ir++) {
253 for (
int it=
nfantome; it<ntv; it++) {
259 for (
int ir=
nfantome; ir<nrv; ir++) {
261 mp.
val_lx(rr, theta, phi0, l, xi) ;
265 for (
int ir=nrv; ir<nrv2; ir++) {
270 for (
int it=ntv; it<ntv2; it++) {
271 for (
int ir=0; ir<nrv2; ir++) {
283 int taille = ntv2*nrv2 ;
285 double* resu =
new double[taille] ;
287 double xi, rr, theta ;
291 for (
int ir=0; ir<nrv2; ir++) {
296 for (
int it=
nfantome; it<ntv; it++) {
302 for (
int ir=
nfantome; ir<nrv; ir++) {
304 mp.
val_lx(rr, theta, phi0, l, xi) ;
308 for (
int ir=nrv; ir<nrv2; ir++) {
313 for (
int it=ntv; it<ntv2; it++) {
314 for (
int ir=0; ir<nrv2; ir++) {
327 int taille = ntv2*nrv2 ;
329 double* resu =
new double[taille] ;
331 double xi, rr, theta ;
335 for (
int ir=0; ir<nrv2; ir++) {
340 for (
int it=
nfantome; it<ntv; it++) {
345 theta =
teti->
t[it] ;
346 for (
int ir=
nfantome; ir<nrv; ir++) {
348 mp.
val_lx(rr, theta, phi0, l, xi) ;
352 for (
int ir=nrv; ir<nrv2; ir++) {
357 for (
int it=ntv; it<ntv2; it++) {
358 for (
int ir=0; ir<nrv2; ir++) {
368 assert(meudon.
get_etat() == ETATQCQ) ;
381 int taille = npv2*ntv2*nrv2 ;
382 if (taille != taille_in) {
383 cout <<
"Gval_spher::somme_spectral3():\n" ;
384 cout <<
"grid size incompatible with array size... exiting!" << endl ;
389 const Map_af* mpaff =
dynamic_cast<const Map_af*
>(&mp) ;
390 assert(mpaff != 0x0) ;
397 for (
int lz=1; lz<nz; lz++) {
398 assert (ntm == mg->
get_nt(lz)) ;
399 assert (npm == mg->
get_np(lz)) ;
405 double* alpha =
new double[nrv0*(npm+2)*ntm] ;
406 double* p_coef = alpha ;
407 double* chebnri = 0x0 ;
410 double* p_func = chebnri ;
412 double** coefm =
new double*[nz] ;
413 for (
int lz=0; lz<nz; lz++) {
414 assert((mtbcf.
t[lz])->get_etat() != ETATNONDEF) ;
415 coefm[lz] = (mtbcf.
t[lz])->t ;
416 if (coefm[lz] == 0x0) {
417 int sizem = mg->
get_nr(lz)*ntm*(npm+2) ;
418 coefm[lz] =
new double[sizem] ;
419 double* pcf = coefm[lz] ;
420 for (
int i=0; i<sizem; i++)
427 for (
int irv=0; irv<nrv0; irv++) {
429 double* tbcf = coefm[lz] ;
430 int nrm = mg->
get_nr(lz) ;
431 for (
int mpm=0; mpm<npm+2; mpm++) {
432 for (
int ltm=0; ltm<ntm; ltm++) {
434 for (
int irm=0; irm<nrm; irm++) {
435 *p_coef += (*tbcf)*(*p_func) ;
445 for (
int lz=0; lz<nz; lz++) {
446 if ((mtbcf.
t[lz])->t == 0x0)
delete [] coefm[lz] ;
452 double* beta =
new double[ntv0*nrv0*(npm+2)] ;
454 double* tetlj = 0x0 ;
457 double* p_interm = alpha ;
461 for (
int jtv=0; jtv<ntv0; jtv++) {
462 for (
int irv=0; irv<nrv0; irv++) {
463 for (
int mpm=0; mpm<npm+2; mpm++) {
465 for (
int ltm=0; ltm<ntm; ltm++) {
466 *p_coef += (*p_interm) * (*p_func) ;
472 p_func -= (npm+2)*ntm ;
475 p_func += (npm+2)*ntm ;
486 double* expmk = 0x0 ;
491 for (
int it=0; it<ntv2; it++) {
492 for (
int ir=0; ir<nrv2; ir++){
498 for (
int ip=
nfantome; ip<npv; ip++) {
500 for (
int ir=0; ir<nrv2; ir++) {
505 for (
int it=
nfantome; it<ntv; it++) {
510 for (
int ir=
nfantome; ir<nrv; ir++) {
512 for (
int mpm=0; mpm<npm+2; mpm++) {
513 *p_coef += (*p_interm) * (*p_func) ;
520 for (
int ir=nrv; ir<nrv2; ir++) {
525 for (
int it=ntv; it<ntv2; it++) {
526 for (
int ir=0; ir<nrv2; ir++) {
534 for (
int ip=npv; ip<npv2; ip++) {
535 for (
int it=0; it<ntv2; it++) {
536 for (
int ir=0; ir<nrv2; ir++){
547 void Gval_spher::initialize_spectral_r(
const Map& mp,
const Base_val& base,
548 int*& idom,
double*& chebnri)
const {
555 assert (idom == 0x0) ;
556 idom =
new int[nrv0] ;
557 double* xi =
new double[nrv0] ;
560 for (
int i=0; i<nrv0; i++) {
562 nrmax += mg->
get_nr(idom[i]) ;
565 assert (chebnri == 0x0) ;
566 chebnri =
new double[(npm+2)*ntm*nrmax] ;
567 double* p_out = chebnri ;
568 for (
int irv=0; irv<nrv0; irv++) {
569 bool nucleus = (mg->
get_type_r(idom[irv]) == RARE) ;
570 int nmax = (nucleus ? 2*mg->
get_nr(idom[irv]) + 1
571 : mg->
get_nr(idom[irv])) ;
572 double* cheb =
new double[nmax] ;
575 for (
int ir=2; ir<nmax; ir++) {
576 cheb[ir] = 2*xi[irv]*cheb[ir-1] - cheb[ir-2] ;
581 for (
int ip=0; ip<npm+2; ip++) {
582 for (
int it=0; it<ntm; it++) {
613 par = 1 - ((ip/2) % 2) ;
618 cout <<
"Gval_spher::initialize_spectral_r : " <<
'\n'
619 <<
"Unexpected radial base !" <<
'\n'
620 <<
"Base : " << base_r << endl ;
626 for (
int ir=0; ir<mg->
get_nr(idom[irv]); ir++) {
627 *p_out = cheb[fact*ir+par] ;
641 void Gval_spher::initialize_spectral_theta(
const Map& mp,
const Base_val& base,
642 double*& tetlj)
const {
645 const Mg3d* mg = mp.get_mg() ;
647 int ntm = mg->get_nt(0) ;
648 int base_t = base.get_base_t(0) ;
650 assert (tetlj == 0x0) ;
651 tetlj =
new double[(npm+2)*ntv0*ntm] ;
652 double* p_out = tetlj ;
654 for (
int jtv=0; jtv<ntv0; jtv++) {
656 for (
int mpm=0; mpm < npm+2; mpm++) {
657 for (
int ltm=0; ltm<ntm; ltm++) {
660 *p_out =
cos(ltm*teta) ;
664 *p_out =
sin(ltm*teta) ;
668 *p_out =
cos(2*ltm*teta) ;
672 *p_out =
cos((2*ltm+1)*teta) ;
676 *p_out =
sin(2*ltm*teta) ;
680 *p_out =
sin((2*ltm+1)*teta) ;
684 *p_out = ( ((mpm/2) % 2 == 0) ?
cos(2*ltm*teta)
685 :
sin((2*ltm+1)*teta)) ;
689 *p_out = ( ((mpm/2) % 2 == 0) ?
cos((2*ltm+1)*teta)
694 *p_out = ( ((mpm/2) % 2 == 0) ?
sin(2*ltm*teta)
695 :
cos((2*ltm+1)*teta)) ;
699 *p_out = ( ((mpm/2) % 2 == 0) ?
sin((2*ltm+1)*teta)
704 *p_out = ( ((mpm/2) % 2 == 0) ?
cos(ltm*teta)
709 *p_out = ( ((mpm/2) % 2 == 0) ?
sin(ltm*teta)
714 cout <<
"Gval_spher::initialize_spectral_theta : " <<
'\n'
715 <<
"Unexpected theta base !" <<
'\n'
716 <<
"Base : " << base_t << endl ;
735 void Gval_spher::initialize_spectral_phi(
const Map& mp,
const Base_val& base,
736 double*& expmk)
const {
739 const Mg3d* mg = mp.get_mg() ;
740 int npm = mg->get_np(0) ;
741 int base_p = base.get_base_p(0) ;
743 assert (expmk == 0x0) ;
744 expmk =
new double[(npm+2)*npv0] ;
745 double* p_out = expmk ;
747 for (
int kpv=0; kpv<npv0; kpv++) {
749 for (
int mpm=0; mpm < npm+2; mpm++) {
753 *p_out = ( (mpm%2 == 0) ?
cos(m*fi) :
sin(m*fi) ) ;
758 *p_out = ( (mpm%2 == 0) ?
cos(2*m*fi) :
sin(2*m*fi) ) ;
763 *p_out = ( (mpm%2 == 0) ?
cos((2*m+1)*fi) :
sin((2*m+1)*fi) ) ;
767 cout <<
"Gval_spher::initialize_spectral_phi : " <<
'\n'
768 <<
"Unexpected phi base !" <<
'\n'
769 <<
"Base : " << base_p << endl ;
Bases of the spectral expansions.
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
int * dim
Array of dimensions (size: ndim).
int ndim
Number of dimensions of the Tbl: can be 1, 2 or 3.
int nfantome
The number of hidden cells (same on each side)
void somme_spectrale1(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...
Dim_tbl dim
The dimensions of the grid.
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes
Tbl * zri
Arrays containing the values of coordinate z (or r) on the interfaces.
virtual void somme_spectrale2(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
Tbl * x
Arrays containing the values of coordinate x on the nodes.
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
virtual void somme_spectrale2(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...
Tbl * phi
Arrays containing the values of coordinate on the nodes.
Tbl * teti
Arrays containing the values of coordinate on the interfaces.
double * somme_spectrale2ri(const Scalar &meudon) const
Same as before but at radial grid interfaces.
double * somme_spectrale2ti(const Scalar &meudon) const
Same as before but at angular grid interfaces.
Tbl * tet
Arrays containing the values of coordinate on the nodes.
Base class for coordinate mappings.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point given by its physical coordina...
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.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Coefficients storage for the multi-domain spectral method.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Tensor field of valence 0 (or component of a tensorial field).
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
const Valeur & get_spectral_va() const
Returns va (read only version)
double * t
The array of double.
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Cmp sqrt(const Cmp &)
Square root.
Cmp sin(const Cmp &)
Sine.
Cmp acos(const Cmp &)
Arccosine.
Cmp cos(const Cmp &)
Cosine.
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define T_SIN_P
dev. sin seulement, harmoniques paires
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
#define T_COS_P
dev. cos seulement, harmoniques paires
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
#define P_COSSIN
dev. standart
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
#define T_SIN_I
dev. sin seulement, harmoniques impaires
#define T_COS
dev. cos seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define T_SIN
dev. sin seulement
#define T_COS_I
dev. cos seulement, harmoniques impaires
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
const Map & get_mp() const
Returns the mapping.