26char 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 & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
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).
Tensor field of valence 0 (or component of a tensorial field).
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...
friend Scalar cos(const Scalar &)
Cosine.
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...
Valeur va
The numerical value of the Scalar
friend Scalar pow(const Scalar &, int)
Power .
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).
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)
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
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()
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.