gval_from_spectral.C
/*
* Functions for spectral summation to a Valencia-type grid (see grille_val.h)
*
*/
/*
* Copyright (c) 2001 and 2004 Jerome Novak
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2
* as published by the Free Software Foundation.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: gval_from_spectral.C,v 1.14 2014/10/13 08:53:48 j_novak Exp $
* $Log: gval_from_spectral.C,v $
* Revision 1.14 2014/10/13 08:53:48 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.13 2014/10/06 15:13:22 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.12 2009/10/28 13:40:23 j_novak
* General case for the theta symmetry (now should work).
*
* Revision 1.11 2009/10/21 13:19:04 j_novak
* Going back (temporary) to previous version.
*
* Revision 1.9 2007/12/21 10:46:29 j_novak
* In "from_spectral..." functions: better treatment of ETATZERO case.
*
* Revision 1.8 2007/11/02 16:49:12 j_novak
* Suppression of intermediate array for spectral summation.
*
* Revision 1.7 2006/10/02 07:41:03 j_novak
* Corrected an error in the case r=0, when exporting to a cartesian grid.
*
* Revision 1.6 2005/06/23 13:44:18 j_novak
* Removed some old comments.
*
* Revision 1.5 2005/06/23 13:40:08 j_novak
* The tests on the number of dimensions have been changed to handle better the
* axisymmetric case.
*
* Revision 1.4 2005/06/22 09:11:17 lm_lin
*
* Grid wedding: convert from the old C++ object "Cmp" to "Scalar".
*
* Revision 1.3 2004/12/17 13:35:04 m_forot
* Add the case T_LEG
*
* Revision 1.2 2004/05/07 13:19:24 j_novak
* Prevention of warnings
*
* Revision 1.1 2004/05/07 12:32:13 j_novak
* New summation from spectral to FD grid. Much faster!
*
*
* $Header: /cvsroot/Lorene/C++/Source/Valencia/gval_from_spectral.C,v 1.14 2014/10/13 08:53:48 j_novak Exp $
*
*/
#include <cmath>
// Lorene headers
#include "grille_val.h"
#include "proto_f77.h"
//--------------------------------------
// Sommation depuis une grille spectrale
//--------------------------------------
namespace Lorene {
void Grille_val::somme_spectrale1(const Scalar& meudon, double* resu, int taille_in) const {
int taille = dim.dim[0]+2*nfantome ;
if (taille != taille_in) {
cout << "Gval_spher::somme_spectral2():\n" ;
cout << "grid size incompatible with array size... exiting!" << endl ;
abort() ;
}
int nrv = dim.dim[0]+nfantome ;
const Map& mp = meudon.get_mp() ;
int l ;
double xi ;
for (int i=0; i<nfantome; i++) resu[i] = 0 ;
for (int i=nfantome; i<nrv; i++) {
mp.val_lx(zr->t[i],0.,0.,l,xi) ;
resu[i] = meudon.get_spectral_va().val_point_jk(l, xi, 0, 0) ;
}
for (int i=nrv; i<taille; i++) resu[i] = 0 ;
}
void Gval_cart::somme_spectrale2(const Scalar& meudon, double* resu, int taille_in) const {
int nzv = dim.dim[0] + nfantome ;
int nxv = dim.dim[1] + nfantome ;
int nzv2 = dim.dim[0] + 2*nfantome ;
int nxv2 = dim.dim[1] + 2*nfantome ;
int taille = nxv2*nzv2 ;
if (taille != taille_in) {
cout << "Gval_spher::somme_spectral2():\n" ;
cout << "grid size incompatible with array size... exiting!" << endl ;
abort() ;
}
const Map& mp = meudon.get_mp() ;
int l ;
double xi0, rr, theta ;
double phi = 0 ;
int inum = 0 ;
for (int ix=0; ix<nfantome; ix++) {
for (int iz=0; iz<nzv2; iz++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int ix=nfantome; ix<nxv; ix++) {
for (int iz=0; iz<nfantome; iz++) {
resu[inum] = 0. ;
inum++ ;
}
double xx2 = (x->t[ix])*(x->t[ix]) ;
for (int iz=nfantome; iz<nzv; iz++) {
rr = sqrt((zr->t[iz])*(zr->t[iz]) + xx2) ;
theta = (rr != 0. ? acos((zr->t[iz])/rr) : 0) ;
mp.val_lx(rr, theta, phi, l, xi0) ;
resu[inum] = meudon.get_spectral_va().val_point(l, xi0, theta, phi) ;
inum++ ;
}
for (int iz=nzv; iz<nzv2; iz++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int ix=nxv; ix<nxv2; ix++) {
for (int iz=0; iz<nzv2; iz++) {
resu[inum] = 0. ;
inum++ ;
}
}
}
void Gval_cart::somme_spectrale3(const Scalar& meudon, double* resu, int taille_in) const{
int nzv = dim.dim[0] + nfantome ;
int nxv = dim.dim[1] + nfantome ;
int nyv = dim.dim[2] + nfantome ;
int nzv2 = dim.dim[0] + 2*nfantome ;
int nxv2 = dim.dim[1] + 2*nfantome ;
int nyv2 = dim.dim[2] + 2*nfantome ;
int taille = nyv2*nxv2*nzv2 ;
if (taille != taille_in) {
cout << "Gval_spher::somme_spectral2():\n" ;
cout << "grid size incompatible with array size... exiting!" << endl ;
abort() ;
}
const Map& mp = meudon.get_mp() ;
int l ;
double xi0, rr, theta, phi ;
int inum = 0 ;
for (int iy=0; iy<nfantome; iy++) {
for (int ix=0; ix<nxv2; ix++) {
for (int iz=0; iz<nzv2; iz++){
resu[inum] = 0. ;
inum++ ;
}
}
}
for (int iy=nfantome; iy<nyv; iy++) {
double yy = x->t[iy] ;
double yy2 = yy*yy ;
for (int ix=0; ix<nfantome; ix++) {
for (int iz=0; iz<nzv2; iz++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int ix=nfantome; ix<nxv; ix++) {
for (int iz=0; iz<nfantome; iz++) {
resu[inum] = 0. ;
inum++ ;
}
double xx = x->t[ix] ;
double xx2 = xx*xx ;
for (int iz=nfantome; iz<nzv; iz++) {
rr = sqrt((zr->t[iz])*(zr->t[iz]) + xx2 + yy2) ;
theta = (rr != 0. ? acos((zr->t[iz])/rr) : 0. );
phi = (rr != 0. ? atan2(yy, xx) : 0. ) ; // return value in [-M_PI,M_PI], should work
mp.val_lx(rr, theta, phi, l, xi0) ;
resu[inum] = meudon.get_spectral_va().val_point(l, xi0, theta, phi) ;
inum++ ;
}
for (int iz=nzv; iz<nzv2; iz++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int ix=nxv; ix<nxv2; ix++) {
for (int iz=0; iz<nzv2; iz++) {
resu[inum] = 0. ;
inum++ ;
}
}
}
for (int iy=nyv; iy<nyv2; iy++) {
for (int ix=0; ix<nxv2; ix++) {
for (int iz=0; iz<nzv2; iz++){
resu[inum] = 0. ;
inum++ ;
}
}
}
}
void Gval_spher::somme_spectrale2(const Scalar& meudon, double* resu, int taille_in) const {
assert (dim.ndim >=2) ;
int nrv = dim.dim[0] + nfantome ;
int ntv = dim.dim[1] + nfantome ;
int nrv2 = dim.dim[0] + 2*nfantome ;
int ntv2 = dim.dim[1] + 2*nfantome ;
int taille = ntv2*nrv2 ;
if (taille != taille_in) {
cout << "Gval_spher::somme_spectral2():\n" ;
cout << "grid size incompatible with array size... exiting!" << endl ;
abort() ;
}
const Map& mp = meudon.get_mp() ;
int l ;
double xi, rr, theta ;
double phi0 = 0 ;
int inum = 0 ;
for (int it=0; it<nfantome; it++) {
for (int ir=0; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int it=nfantome; it<ntv; it++) {
for (int ir=0; ir<nfantome; ir++) {
resu[inum] = 0. ;
inum++ ;
}
theta = tet->t[it] ;
for (int ir=nfantome; ir<nrv; ir++) {
rr = zr->t[ir] ;
mp.val_lx(rr, theta, phi0, l, xi) ;
resu[inum] = meudon.get_spectral_va().val_point(l, xi, theta, phi0) ;
inum++ ;
}
for (int ir=nrv; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int it=ntv; it<ntv2; it++) {
for (int ir=0; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
}
double* Gval_spher::somme_spectrale2ri(const Scalar& meudon) const {
int nrv = dim.dim[0] + 1 + nfantome ;
int ntv = dim.dim[1] + nfantome ;
int nrv2 = dim.dim[0] + 1 + 2*nfantome ;
int ntv2 = dim.dim[1] + 2*nfantome ;
int taille = ntv2*nrv2 ;
const Map& mp = meudon.get_mp() ;
double* resu = new double[taille] ;
int l ;
double xi, rr, theta ;
double phi0 = 0 ;
int inum = 0 ;
for (int it=0; it<nfantome; it++) {
for (int ir=0; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int it=nfantome; it<ntv; it++) {
for (int ir=0; ir<nfantome; ir++) {
resu[inum] = 0. ;
inum++ ;
}
theta = tet->t[it] ;
for (int ir=nfantome; ir<nrv; ir++) {
rr = zri->t[ir] ;
mp.val_lx(rr, theta, phi0, l, xi) ;
resu[inum] = meudon.get_spectral_va().val_point(l, xi, theta, phi0) ;
inum++ ;
}
for (int ir=nrv; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int it=ntv; it<ntv2; it++) {
for (int ir=0; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
return resu ;
}
double* Gval_spher::somme_spectrale2ti(const Scalar& meudon) const {
int nrv = dim.dim[0] + nfantome ;
int ntv = dim.dim[1] + 1 + nfantome ;
int nrv2 = dim.dim[0] + 2*nfantome ;
int ntv2 = dim.dim[1] + 1 + 2*nfantome ;
int taille = ntv2*nrv2 ;
const Map& mp = meudon.get_mp() ;
double* resu = new double[taille] ;
int l ;
double xi, rr, theta ;
double phi0 = 0 ;
int inum = 0 ;
for (int it=0; it<nfantome; it++) {
for (int ir=0; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int it=nfantome; it<ntv; it++) {
for (int ir=0; ir<nfantome; ir++) {
resu[inum] = 0. ;
inum++ ;
}
theta = teti->t[it] ;
for (int ir=nfantome; ir<nrv; ir++) {
rr = zr->t[ir] ;
mp.val_lx(rr, theta, phi0, l, xi) ;
resu[inum] = meudon.get_spectral_va().val_point(l, xi, theta, phi0) ;
inum++ ;
}
for (int ir=nrv; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
for (int it=ntv; it<ntv2; it++) {
for (int ir=0; ir<nrv2; ir++) {
resu[inum] = 0. ;
inum++ ;
}
}
return resu ;
}
void Gval_spher::somme_spectrale3(const Scalar& meudon, double* resu, int taille_in) const{
assert(meudon.get_etat() == ETATQCQ) ;
meudon.get_spectral_va().coef() ;
//Sizes of both grids
//-------------------
int nrv0 = dim.dim[0] ;
int ntv0 = dim.dim[1] ;
int nrv = dim.dim[0] + nfantome ;
int ntv = dim.dim[1] + nfantome ;
int npv = dim.dim[2] + nfantome ;
int nrv2 = dim.dim[0] + 2*nfantome ;
int ntv2 = dim.dim[1] + 2*nfantome ;
int npv2 = dim.dim[2] + 2*nfantome ;
int taille = npv2*ntv2*nrv2 ;
if (taille != taille_in) {
cout << "Gval_spher::somme_spectral3():\n" ;
cout << "grid size incompatible with array size... exiting!" << endl ;
abort() ;
}
const Map& mp = meudon.get_mp() ;
#ifndef NDEBUG
const Map_af* mpaff = dynamic_cast<const Map_af*>(&mp) ;
assert(mpaff != 0x0) ;
#endif
const Mg3d* mg = mp.get_mg() ;
int ntm = mg->get_nt(0) ;
int npm = mg->get_np(0) ;
int nz = mg->get_nzone() ;
#ifndef NDEBUG
for (int lz=1; lz<nz; lz++) {
assert (ntm == mg->get_nt(lz)) ; //Same angular grids in all domains...
assert (npm == mg->get_np(lz)) ;
}
#endif
//Intermediate quantities
//-----------------------
double* alpha = new double[nrv0*(npm+2)*ntm] ;
double* p_coef = alpha ;
double* chebnri = 0x0 ; //size ~ nrv0 * (npm+2) * nr ...
int* idom = 0x0 ;
initialize_spectral_r(mp, meudon.get_spectral_va().get_base(), idom, chebnri) ;
double* p_func = chebnri ;
Mtbl_cf& mtbcf = *meudon.get_spectral_va().c_cf ;
double** coefm = new double*[nz] ;
for (int lz=0; lz<nz; lz++) {
assert((mtbcf.t[lz])->get_etat() != ETATNONDEF) ;
coefm[lz] = (mtbcf.t[lz])->t ;
if (coefm[lz] == 0x0) {
int sizem = mg->get_nr(lz)*ntm*(npm+2) ;
coefm[lz] = new double[sizem] ;
double* pcf = coefm[lz] ;
for (int i=0; i<sizem; i++)
pcf[i] = 0. ;
}
}
//First partial summation
//-----------------------
for (int irv=0; irv<nrv0; irv++) {
int lz = idom[irv] ;
double* tbcf = coefm[lz] ;
int nrm = mg->get_nr(lz) ;
for (int mpm=0; mpm<npm+2; mpm++) {
for (int ltm=0; ltm<ntm; ltm++) {
*p_coef = 0 ;
for (int irm=0; irm<nrm; irm++) {
*p_coef += (*tbcf)*(*p_func) ;
tbcf++ ;
p_func++ ;
// cout << *p_func << ", " << *tbcf << ", " << *p_coef << endl ;
}
p_coef++ ;
}
}
}
for (int lz=0; lz<nz; lz++) {
if ((mtbcf.t[lz])->t == 0x0) delete [] coefm[lz] ;
}
delete [] coefm ;
delete [] chebnri ;
delete [] idom ;
double* beta = new double[ntv0*nrv0*(npm+2)] ;
p_coef = beta ;
double* tetlj = 0x0 ;
initialize_spectral_theta(mp, meudon.get_spectral_va().get_base(), tetlj) ;
p_func = tetlj ;
double* p_interm = alpha ;
//Second partial summation
//------------------------
for (int jtv=0; jtv<ntv0; jtv++) {
for (int irv=0; irv<nrv0; irv++) {
for (int mpm=0; mpm<npm+2; mpm++) {
*p_coef = 0 ;
for (int ltm=0; ltm<ntm; ltm++) {
*p_coef += (*p_interm) * (*p_func) ;
p_interm++ ;
p_func++ ;
}
p_coef++ ;
} // Loop on m
p_func -= (npm+2)*ntm ;
} //Loop on irv
p_interm = alpha ;
p_func += (npm+2)*ntm ;
} //Loop on jtv
delete [] alpha ;
delete [] tetlj ;
// Final summation
//----------------
p_interm = beta ;
double* expmk = 0x0 ;
initialize_spectral_phi(mp, meudon.get_spectral_va().get_base(), expmk) ;
p_func = expmk ;
p_coef = resu ;
for (int ip=0; ip<nfantome; ip++) {
for (int it=0; it<ntv2; it++) {
for (int ir=0; ir<nrv2; ir++){
*p_coef = 0. ;
p_coef++ ;
}
}
}
for (int ip=nfantome; ip<npv; ip++) {
for (int it=0; it<nfantome; it++) {
for (int ir=0; ir<nrv2; ir++) {
*p_coef = 0. ;
p_coef++ ;
}
}
for (int it=nfantome; it<ntv; it++) {
for (int ir=0; ir<nfantome; ir++) {
*p_coef = 0. ;
p_coef++ ;
}
for (int ir=nfantome; ir<nrv; ir++) {
*p_coef = 0. ;
for (int mpm=0; mpm<npm+2; mpm++) {
*p_coef += (*p_interm) * (*p_func) ;
p_interm++ ;
p_func++ ;
}
p_coef++ ;
p_func -= (npm+2) ;
}
for (int ir=nrv; ir<nrv2; ir++) {
*p_coef = 0. ;
p_coef++ ;
}
}
for (int it=ntv; it<ntv2; it++) {
for (int ir=0; ir<nrv2; ir++) {
*p_coef = 0. ;
p_coef++ ;
}
}
p_func += npm+2 ; //Next point in phi
p_interm = beta ;
}
for (int ip=npv; ip<npv2; ip++) {
for (int it=0; it<ntv2; it++) {
for (int ir=0; ir<nrv2; ir++){
*p_coef = 0. ;
p_coef++ ;
}
}
}
delete [] expmk ;
delete [] beta ;
}
void Gval_spher::initialize_spectral_r(const Map& mp, const Base_val& base,
int*& idom, double*& chebnri) const {
int nrv0 = dim.dim[0] ;
const Mg3d* mg = mp.get_mg() ;
int npm = mg->get_np(0) ;
int ntm = mg->get_nt(0) ;
assert (idom == 0x0) ;
idom = new int[nrv0] ;
double* xi = new double[nrv0] ;
int nrmax = 0 ;
for (int i=0; i<nrv0; i++) {
mp.val_lx(zr->t[i+nfantome], 0., 0., idom[i], xi[i]) ;
nrmax += mg->get_nr(idom[i]) ;
}
assert (chebnri == 0x0) ;
chebnri = new double[(npm+2)*ntm*nrmax] ;
double* p_out = chebnri ;
for (int irv=0; irv<nrv0; irv++) {
bool nucleus = (mg->get_type_r(idom[irv]) == RARE) ;
int nmax = (nucleus ? 2*mg->get_nr(idom[irv]) + 1
: mg->get_nr(idom[irv])) ;
double* cheb = new double[nmax] ;
cheb[0] = 1. ;
cheb[1] = xi[irv] ;
for (int ir=2; ir<nmax; ir++) {
cheb[ir] = 2*xi[irv]*cheb[ir-1] - cheb[ir-2] ;
}
int base_r = base.get_base_r(idom[irv]) ;
for (int ip=0; ip<npm+2; ip++) {
for (int it=0; it<ntm; it++) {
int fact = 1 ;
int par = 0 ;
if (nucleus) {
fact = 2 ;
switch (base_r) {
case R_CHEBP : {
break ;
}
case R_CHEBI : {
par = 1 ;
break ;
}
case R_CHEBPI_P : {
par = it % 2 ;
break ;
}
case R_CHEBPI_I : {
par = 1 - (it % 2) ;
break ;
}
case R_CHEBPIM_P : {
par = (ip/2) % 2 ;
break ;
}
case R_CHEBPIM_I : {
par = 1 - ((ip/2) % 2) ;
break ;
}
default : {
cout << "Gval_spher::initialize_spectral_r : " << '\n'
<< "Unexpected radial base !" << '\n'
<< "Base : " << base_r << endl ;
abort() ;
break ;
}
}
}
for (int ir=0; ir<mg->get_nr(idom[irv]); ir++) {
*p_out = cheb[fact*ir+par] ;
p_out++ ;
}
} // Loop on it
} // Loop on ip
delete [] cheb ;
}// Loop on irv
delete [] xi ;
}
void Gval_spher::initialize_spectral_theta(const Map& mp, const Base_val& base,
double*& tetlj) const {
int ntv0 = dim.dim[1] ;
const Mg3d* mg = mp.get_mg() ;
int npm = mg->get_np(0) ;
int ntm = mg->get_nt(0) ;
int base_t = base.get_base_t(0) ;
assert (tetlj == 0x0) ;
tetlj = new double[(npm+2)*ntv0*ntm] ;
double* p_out = tetlj ;
for (int jtv=0; jtv<ntv0; jtv++) {
double teta = tet->t[jtv+nfantome] ;
for (int mpm=0; mpm < npm+2; mpm++) {
for (int ltm=0; ltm<ntm; ltm++) {
switch (base_t) { //## One should use array of functions...
case T_COS : {
*p_out = cos(ltm*teta) ;
break ;
}
case T_SIN : {
*p_out = sin(ltm*teta) ;
break ;
}
case T_COS_P : {
*p_out = cos(2*ltm*teta) ;
break ;
}
case T_COS_I : {
*p_out = cos((2*ltm+1)*teta) ;
break ;
}
case T_SIN_P : {
*p_out = sin(2*ltm*teta) ;
break ;
}
case T_SIN_I : {
*p_out = sin((2*ltm+1)*teta) ;
break ;
}
case T_COSSIN_CP : {
*p_out = ( ((mpm/2) % 2 == 0) ? cos(2*ltm*teta)
: sin((2*ltm+1)*teta)) ;
break ;
}
case T_COSSIN_CI : {
*p_out = ( ((mpm/2) % 2 == 0) ? cos((2*ltm+1)*teta)
: sin(2*ltm*teta)) ;
break ;
}
case T_COSSIN_SP : {
*p_out = ( ((mpm/2) % 2 == 0) ? sin(2*ltm*teta)
: cos((2*ltm+1)*teta)) ;
break ;
}
case T_COSSIN_SI : {
*p_out = ( ((mpm/2) % 2 == 0) ? sin((2*ltm+1)*teta)
: cos(2*ltm*teta)) ;
break ;
}
case T_COSSIN_C : {
*p_out = ( ((mpm/2) % 2 == 0) ? cos(ltm*teta)
: sin(ltm*teta)) ;
break ;
}
case T_COSSIN_S : {
*p_out = ( ((mpm/2) % 2 == 0) ? sin(ltm*teta)
: cos(ltm*teta)) ;
break ;
}
default : {
cout << "Gval_spher::initialize_spectral_theta : " << '\n'
<< "Unexpected theta base !" << '\n'
<< "Base : " << base_t << endl ;
abort() ;
break ;
}
}
p_out++ ;
}
if ( (base_t == T_COS_I) || (base_t == T_SIN_P) || (base_t == T_SIN_I) )
{
p_out-- ;
*p_out = 0. ;
p_out++ ;
}
} //Loop on mpm
} //Loop on jtv
}
void Gval_spher::initialize_spectral_phi(const Map& mp, const Base_val& base,
double*& expmk) const {
int npv0 = dim.dim[2] ;
const Mg3d* mg = mp.get_mg() ;
int npm = mg->get_np(0) ;
int base_p = base.get_base_p(0) ;
assert (expmk == 0x0) ;
expmk = new double[(npm+2)*npv0] ;
double* p_out = expmk ;
for (int kpv=0; kpv<npv0; kpv++) {
double fi = phi->t[kpv+nfantome] ;
for (int mpm=0; mpm < npm+2; mpm++) {
switch (base_p) { //## One should use array of functions...
case P_COSSIN : {
int m = mpm / 2 ;
*p_out = ( (mpm%2 == 0) ? cos(m*fi) : sin(m*fi) ) ;
break ;
}
case P_COSSIN_P : {
int m = mpm / 2 ;
*p_out = ( (mpm%2 == 0) ? cos(2*m*fi) : sin(2*m*fi) ) ;
break ;
}
case P_COSSIN_I : {
int m = mpm / 2 ;
*p_out = ( (mpm%2 == 0) ? cos((2*m+1)*fi) : sin((2*m+1)*fi) ) ;
break ;
}
default : {
cout << "Gval_spher::initialize_spectral_phi : " << '\n'
<< "Unexpected phi base !" << '\n'
<< "Base : " << base_p << endl ;
abort() ;
break ;
}
}
p_out++ ;
} //Loop on mpm
} //Loop on kpv
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::get_base_r
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition base_val.h:400
Lorene::Dim_tbl::dim
int * dim
Array of dimensions (size: ndim).
Definition dim_tbl.h:102
Lorene::Grille_val::nfantome
int nfantome
The number of hidden cells (same on each side).
Definition grille_val.h:104
Lorene::Grille_val::somme_spectrale1
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...
Definition gval_from_spectral.C:90
Lorene::Grille_val::dim
Dim_tbl dim
The dimensions of the grid.
Definition grille_val.h:102
Lorene::Grille_val::zr
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes.
Definition grille_val.h:124
Lorene::Grille_val::zri
Tbl * zri
Arrays containing the values of coordinate z (or r) on the interfaces.
Definition grille_val.h:126
Lorene::Gval_cart::somme_spectrale2
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...
Definition gval_from_spectral.C:110
Lorene::Gval_cart::somme_spectrale3
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
Definition gval_from_spectral.C:158
Lorene::Gval_cart::x
Tbl * x
Arrays containing the values of coordinate x on the nodes.
Definition grille_val.h:343
Lorene::Gval_spher::somme_spectrale3
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const
Same as before but for the 3D case.
Definition gval_from_spectral.C:366
Lorene::Gval_spher::somme_spectrale2
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...
Definition gval_from_spectral.C:229
Lorene::Gval_spher::phi
Tbl * phi
Arrays containing the values of coordinate on the nodes.
Definition grille_val.h:591
Lorene::Gval_spher::teti
Tbl * teti
Arrays containing the values of coordinate on the interfaces.
Definition grille_val.h:589
Lorene::Gval_spher::somme_spectrale2ri
double * somme_spectrale2ri(const Scalar &meudon) const
Same as before but at radial grid interfaces.
Definition gval_from_spectral.C:278
Lorene::Gval_spher::somme_spectrale2ti
double * somme_spectrale2ti(const Scalar &meudon) const
Same as before but at angular grid interfaces.
Definition gval_from_spectral.C:322
Lorene::Gval_spher::tet
Tbl * tet
Arrays containing the values of coordinate on the nodes.
Definition grille_val.h:587
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::Valeur::val_point
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 ...
Definition valeur.C:882
Lorene::Valeur::val_point_jk
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...
Definition valeur.C:900
Lorene::Valeur::get_base
const Base_val & get_base() const
Return the bases for spectral expansions (member base ).
Definition valeur.h:480
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:169
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
P_COSSIN_P
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition type_parite.h:247
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition type_parite.h:210
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition type_parite.h:202
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition type_parite.h:194
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition type_parite.h:214
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition type_parite.h:200
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition type_parite.h:212
P_COSSIN
#define P_COSSIN
dev. standart
Definition type_parite.h:245
P_COSSIN_I
#define P_COSSIN_I
dev. sur Phi = 2*phi, freq. impaires
Definition type_parite.h:249
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition type_parite.h:176
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition type_parite.h:208
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition type_parite.h:206
T_COS
#define T_COS
dev. cos seulement
Definition type_parite.h:196
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
T_SIN
#define T_SIN
dev. sin seulement
Definition type_parite.h:198
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition type_parite.h:204
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition type_parite.h:192
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720