26 char scalar_raccord_zec_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_zec.C,v 1.8 2014/10/13 08:53:47 j_novak Exp $" ;
70 #include "utilitaires.h"
80 assert (
etat != ETATNONDEF) ;
81 if ((
etat == ETATZERO) || (
etat == ETATUN))
87 cout <<
"Le mapping doit etre affine" << endl ;
97 double r_cont = -1./2./alpha ;
100 Tbl coef (nbre+2*lmax, nr) ;
103 int* deg =
new int[3] ;
104 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
105 double* auxi =
new double[nr] ;
107 for (
int conte=0 ; conte<nbre+2*lmax ; conte++) {
108 for (
int i=0 ; i<nr ; i++)
109 auxi[i] =
pow(-1-
cos(M_PI*i/(nr-1)), (conte+puis)) ;
111 cfrcheb(deg, deg, auxi, deg, auxi) ;
112 for (
int i=0 ; i<nr ; i++)
113 coef.
set(conte, i) = auxi[i]*
pow (alpha, conte+puis) ;
118 Tbl valeurs (nbre, nt, np+1) ;
122 double* res_val =
new double[1] ;
124 for (
int conte=0 ; conte<nbre ; conte++) {
131 for (
int k=0 ; k<np+1 ; k++)
132 for (
int j=0 ; j<nt ; j++)
133 if (nullite_plm(j, nt, k, np, courant.
va.
base) == 1) {
135 for (
int i=0 ; i<nr ; i++)
136 auxi[i] = (*courant.
va.
c_cf)(nz-2, k, j, i) ;
140 som_r_cheb (auxi, nr, 1, 1, 1, res_val) ;
143 cout <<
"Cas non prevu dans raccord_zec" << endl ;
147 valeurs.
set(conte, k, j) = res_val[0] ;
151 courant = copie.
dsdr() ;
165 int base_r, l_quant, m_quant ;
166 for (
int k=0 ; k<np+1 ; k++)
167 for (
int j=0 ; j<nt ; j++)
168 if (nullite_plm(j, nt, k, np,
va.
base) == 1) {
170 donne_lm (nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
177 for (
int col=0 ; col<nbre ; col++)
178 for (
int lig=0 ; lig<nbre ; lig++) {
180 int facteur = (lig%2==0) ? 1 : -1 ;
181 for (
int conte=0 ; conte<lig ; conte++)
182 facteur *= puis+col+conte+2*l_quant ;
183 systeme.
set(lig, col) = facteur/
pow(r_cont, puis+col+lig+2*l_quant) ;
189 Tbl sec_membre (nbre) ;
191 for (
int conte=0 ; conte<nbre ; conte++)
192 sec_membre.
set(conte) = valeurs(conte, k, j) ;
196 for (
int conte=0 ; conte<nbre ; conte++)
197 for (
int i=0 ; i<nr ; i++)
199 inv(conte)*coef(conte+2*l_quant, i) ;
201 else for (
int i=0 ; i<nr ; i++)
217 if (
etat == ETATZERO) return ;
224 int np = mgrid.
get_np(nzm1) ;
225 int nt = mgrid.
get_nt(nzm1) ;
226 int nr_ced = mgrid.
get_nr(nzm1) ;
227 int nr_shell = mgrid.
get_nr(nzm1-1) ;
229 assert(mgrid.
get_np(nzm1-1) == np) ;
230 assert(mgrid.
get_nt(nzm1-1) == nt) ;
237 cout <<
"Scalar::smooth_decay: present version supports only \n"
238 <<
" affine mappings !" << endl ;
250 assert(
va.
c_cf->
t[nzm1] != 0x0) ;
256 int nbr1[] = {nr_shell, nr_ced} ;
257 int nbt1[] = {1, 1} ;
258 int nbp1[] = {1, 1} ;
259 int typr1[] = {mgrid.
get_type_r(nzm1-1), UNSURR} ;
260 Mg3d grid1d(2, nbr1, typr1, nbt1, SYM, nbp1, SYM) ;
264 double rbound = mapaf->
get_alpha()[nzm1-1]
266 double rmin = - mapaf->
get_alpha()[nzm1-1]
268 double bound1[] = {rmin, rbound, __infinity} ;
270 Map_af map1d(grid1d, bound1) ;
284 for (
int k=0; k<=np; k++) {
285 for (
int j=0; j<nt; j++) {
287 if (nullite_plm(j, nt, k, np, base) != 1) {
288 for (
int i=0 ; i<nr_ced ; i++) {
289 coefresu.
set(k, j, i) = 0 ;
293 int base_r_ced, base_r_shell , l_quant, m_quant ;
294 donne_lm(nz, nzm1, j, k, base, m_quant, l_quant, base_r_ced) ;
295 donne_lm(nz, nzm1-1, j, k, base, m_quant, l_quant, base_r_shell) ;
297 int nn0 = l_quant + nn ;
304 for (
int i=0; i<nr_shell; i++) {
306 va.
c_cf->operator()(nzm1-1, k, j, i) ;
323 for (
int p=1; p<=kk; p++) {
335 for (
int im=0; im<=kk; im++) {
337 double fact = (im%2 == 0) ? 1. : -1. ;
338 fact /=
pow(rbound, nn0 + im) ;
340 for (
int jm=0; jm<=kk; jm++) {
343 for (
int ir=0; ir<im; ir++) {
344 prod *= nn0 + jm + ir ;
347 mat.
set(im, jm) = fact * prod /
pow(rbound, jm) ;
362 const Coord& r = map1d.
r ;
363 for (
int p=0; p<=kk; p++) {
364 tmp = aa(p) /
pow(r, nn0 + p) ;
375 for (
int i=0; i<nr_ced; i++) {
376 coefresu.
set(k,j,i) = 0 ;
380 assert( vv.
va.
c_cf != 0x0 ) ;
381 for (
int i=0; i<nr_ced; i++) {
382 coefresu.
set(k,j,i) = vv.
va.
c_cf->operator()(1,0,0,i) ;
404 for (
int p=0; p<=kk; p++) {
407 for (
int k=0; k<np; k++) {
408 for (
int j=0; j<nt; j++) {
411 if (diff > err) err = diff ;
415 "Scalar::smooth_decay: Max error matching of " << p
416 <<
"-th derivative: " << err << 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.
Active physical coordinates and mapping derivatives.
const double * get_beta() const
Returns the pointer on the array beta.
const double * get_alpha() const
Returns the pointer on the array alpha.
Coord r
r coordinate centered on the grid
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int j, int i)
Read/write of a particuliar element.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void set_band(int up, int low) const
Calculate the band storage of *std.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
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.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Tensor field of valence 0 (or component of a tensorial field).
friend Scalar cos(const Scalar &)
Cosine.
int get_dzpuis() const
Returns dzpuis.
void smooth_decay(int k, int n)
Performs a matching of the last non-compactified shell with a decaying function where is the spher...
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
const Scalar & dsdr() const
Returns of *this .
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
void set_dzpuis(int)
Modifies the dzpuis flag.
void raccord_c1_zec(int puis, int nbre, int lmax)
Performs the matching of the external domain with respect to the last shell using function like wit...
friend Scalar pow(const Scalar &, int)
Power .
Valeur va
The numerical value of the Scalar
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
void annule_hard()
Sets the Tbl to zero in a hard way.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
int get_etat() const
Returns the logical state.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
void ylm()
Computes the coefficients of *this.
Mtbl * c
Values of the function at the points of the multi-grid
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
void ylm_i()
Inverse of ylm()
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Base_val base
Bases on which the spectral expansion is performed.
#define R_CHEB
base de Chebychev ordinaire (fin)
const Map *const mp
Mapping on which the numerical values at the grid points are defined.