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 ;
#define R_CHEB
base de Chebychev ordinaire (fin)
double & set(int i)
Read/write of a particular element (index i) (1D case)
Bases of the spectral expansions.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Valeur va
The numerical value of the Scalar
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...
int get_etat() const
Returns the logical state.
const Scalar & dsdr() const
Returns of *this .
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
void annule_hard()
Sets the Tbl to zero in a hard way.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Coord r
r coordinate centered on the grid
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Tensor field of valence 0 (or component of a tensorial field).
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Mtbl * c
Values of the function at the points of the multi-grid
void coef() const
Computes the coeffcients of *this.
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...
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Active physical coordinates and mapping derivatives.
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
void set_band(int up, int low) const
Calculate the band storage of *std.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
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.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void set_dzpuis(int)
Modifies the dzpuis flag.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const double * get_beta() const
Returns the pointer on the array beta.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Base_val base
Bases on which the spectral expansion is performed.
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
int get_dzpuis() const
Returns dzpuis.
void ylm()
Computes the coefficients of *this.
double & set(int j, int i)
Read/write of a particuliar element.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
friend Scalar pow(const Scalar &, int)
Power .
const double * get_alpha() const
Returns the pointer on the array alpha.
void smooth_decay(int k, int n)
Performs a matching of the last non-compactified shell with a decaying function where is the spher...
friend Scalar cos(const Scalar &)
Cosine.
void ylm_i()
Inverse of ylm()
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.