103 source_r.
poisson(par_khi, khi) ;
153 for (
int i=0; i<3; i++)
154 assert(
cmp[i]->check_dzpuis(4)) ;
159 assert( mpaff != 0x0) ;
171 if (
cmp[0]->get_etat() == ETATZERO) {
180 int nz = mg.
get_nzone() ;
int nzm1 = nz - 1;
185 S_r.set_spectral_va().ylm() ;
202 double beta = mpaff->
get_beta()[0] ;
203 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
208 for (
int k=0 ; k<np+1 ; k++) {
209 for (
int j=0 ; j<nt ; j++) {
211 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
212 int dege1 = (l_q <3 ? 0 : 1) ;
214 int nr_eq1 = nr0 - dege1 ;
215 int nr_eq2 = nr0 - dege2 ;
216 int nrtot = nr_eq1 + nr_eq2 ;
225 for (
int lin=0; lin<nr_eq1; lin++) {
226 for (
int col=dege1; col<nr0; col++)
227 oper.
set(lin,col-dege1)
228 = md2(lin,col) + 6*mxd(lin,col) + 6*mid(lin,col) ;
229 for (
int col=dege2; col<nr0; col++)
230 oper.
set(lin,col-dege2+nr_eq1) = -mxd(lin,col) -2*mid(lin,col) ;
232 for (
int lin=0; lin<nr0; lin++) {
233 for (
int col=dege1; col<nr0; col++)
234 oper.
set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
235 for (
int col=dege2; col<nr0; col++)
236 oper.
set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) + 4*mid(lin,col) ;
242 for (
int i=0; i<nr_eq1; i++)
244 for (
int i=0; i<nr0; i++)
245 sec_membre.
set(i+nr_eq1) = 0 ;
255 for (
int i=0; i<dege1; i++)
257 for (
int i=dege1; i<nr0; i++)
258 res_eta.
set(i) = big_res(i-dege1) ;
259 res_eta.
set(nr0) = 0 ;
260 for (
int i=0; i<dege2; i++)
262 for (
int i=dege2; i<nr0; i++)
263 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
264 res_vr.
set(nr0) = 0 ;
268 multx2_1d(nr, &res_eta.
t, base_r) ;
269 multx2_1d(nr, &res_vr.
t, base_r) ;
272 Tbl sol_hom = solh(nr, l_q-1, 0., base_r) ;
273 for (
int i=0 ; i<nr ; i++) {
274 sol_part_eta.
set(0, k, j, i) = alpha*alpha*res_eta(i) ;
275 sol_part_vr.
set(0, k, j, i) = alpha*alpha*res_vr(i) ;
276 solution_hom_un.
set(0, k, j, i) = sol_hom(i) ;
277 solution_hom_deux.
set(0, k, j, i) = 0. ;
286 for (
int zone=1 ; zone<nzm1 ; zone++) {
289 assert(nt == mg.
get_nt(zone)) ;
290 assert(np == mg.
get_np(zone)) ;
293 double ech = beta / alpha ;
297 for (
int k=0 ; k<np+1 ; k++) {
298 for (
int j=0 ; j<nt ; j++) {
300 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
303 int nr_eq1 = nr - dege1 ;
304 int nr_eq2 = nr - dege2 ;
305 int nrtot = nr_eq1 + nr_eq2 + 1;
313 Diff_id id(base_r, nr+1) ;
const Matrice& mid =
id.get_matrice() ;
317 for (
int lin=0; lin<nr_eq1; lin++) {
318 for (
int col=dege1; col<nr; col++)
319 oper.
set(lin,col-dege1)
320 = mxd2(lin,col) + ech*md2(lin,col) + 2*md(lin,col) ;
321 for (
int col=dege2; col<nr+1; col++)
322 oper.
set(lin,col-dege2+nr_eq1) = -md(lin,col) ;
324 for (
int lin=0; lin<nr_eq2; lin++) {
325 for (
int col=dege1; col<nr; col++)
326 oper.
set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
327 for (
int col=dege2; col<nr+1; col++)
328 oper.
set(lin+nr_eq1, col-dege2+nr_eq1) =
329 mxd(lin,col) + ech*md(lin,col) + 2*mid(lin,col) ;
333 for (
int col=dege1; col<nr; col++)
334 oper.
set(nrtot-1, col-dege1) = 0 ;
335 for (
int col=dege2; col<nr+1; col++)
336 oper.
set(nrtot-1, col-dege2+nr_eq1) =
337 m2d2(0,col) + ech*(2*mxd2(0,col) + ech*md2(0,col))
338 +4*(mxd(0,col) +ech*md(0,col))
339 +(2 - l_q*(l_q+1))*mid(0,col) ;
346 for (
int i=0; i<5; i++) {
347 sr.
set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
350 for (
int i=5; i<nr; i++)
352 Tbl xsr= sr ;
Tbl x2sr= sr ;
354 multx2_1d(5, &x2sr.
t, base_r) ; multx_1d(5, &xsr.
t, base_r) ;
355 multx_1d(nr, &xseta.
t, base_r) ;
357 for (
int i=0; i<nr_eq1; i++)
358 sec_membre.
set(i) = alpha*(alpha*xseta(i) + beta*seta(i)) ;
359 for (
int i=0; i<nr_eq2; i++)
360 sec_membre.
set(i+nr_eq1) = 0 ;
361 sec_membre.
set(nr_eq1+nr_eq2) = alpha*alpha*x2sr(0) + 2*alpha*beta*xsr(0)
372 for (
int i=0; i<dege1; i++)
374 for (
int i=dege1; i<nr; i++)
375 res_eta.
set(i) = big_res(i-dege1) ;
376 for (
int i=0; i<dege2; i++)
378 for (
int i=dege2; i<nr; i++)
379 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
382 Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
383 Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
384 for (
int i=0 ; i<nr ; i++) {
385 sol_part_eta.
set(zone, k, j, i) = res_eta(i) ;
386 sol_part_vr.
set(zone, k, j, i) = res_vr(i) ;
387 solution_hom_un.
set(zone, k, j, i) = sol_hom1(0,i) ;
388 solution_hom_deux.
set(zone, k, j, i) = sol_hom2(1,i) ;
398 assert(nt == mg.
get_nt(nzm1)) ;
399 assert(np == mg.
get_np(nzm1)) ;
406 for (
int k=0 ; k<np+1 ; k++) {
407 for (
int j=0 ; j<nt ; j++) {
409 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
412 int nr_eq1 = nr0 - dege1 ;
413 int nr_eq2 = nr0 - dege2 ;
414 int nrtot = nr_eq1 + nr_eq2 ;
423 for (
int lin=0; lin<nr_eq1; lin++) {
424 for (
int col=dege1; col<nr0; col++)
425 oper.
set(lin,col-dege1)
426 = m2d2(lin,col) + 4*mxd(lin,col) + 2*mid(lin,col) ;
427 for (
int col=dege2; col<nr0; col++)
428 oper.
set(lin,col-dege2+nr_eq1) =
429 mxd(lin,col) + 2*mid(lin,col) ;
431 for (
int lin=0; lin<nr_eq2; lin++) {
432 for (
int col=dege1; col<nr0; col++)
433 oper.
set(lin+nr_eq1,col-dege1) = l_q*(l_q+1)*mid(lin,col) ;
434 for (
int col=dege2; col<nr0; col++)
435 oper.
set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) ;
441 for (
int i=0; i<nr_eq1; i++)
443 for (
int i=0; i<nr_eq2; i++)
444 sec_membre.
set(i+nr_eq1) = 0 ;
451 for (
int i=0; i<dege1; i++)
453 for (
int i=dege1; i<nr0; i++)
454 res_eta.
set(i) = big_res(i-dege1) ;
455 res_eta.
set(nr0) = 0 ;
456 res_eta.
set(nr0+1) = 0 ;
457 for (
int i=0; i<dege2; i++)
459 for (
int i=dege2; i<nr0; i++)
460 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
461 res_vr.
set(nr0) = 0 ;
462 res_vr.
set(nr0+1) = 0 ;
468 mult2_xm1_1d_cheb(nr, res_eta.
t, res_e2.
t) ;
469 mult2_xm1_1d_cheb(nr, res_vr.
t, res_v2.
t) ;
472 Tbl sol_hom = solh(nr, l_q+1, 0., base_r) ;
473 for (
int i=0 ; i<nr ; i++) {
474 sol_part_eta.
set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
475 sol_part_vr.
set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
476 solution_hom_un.
set(nzm1, k, j, i) = sol_hom(i) ;
477 solution_hom_deux.
set(nzm1, k, j, i) = 0. ;
497 int taille = 2*nzm1 ;
498 Tbl sec_membre(taille) ;
499 Matrice systeme(taille, taille) ;
501 int ligne ;
int colonne ;
505 for (
int k=0 ; k<np+1 ; k++)
506 for (
int j=0 ; j<nt ; j++) {
508 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
513 for (
int l=0; l<taille; l++)
514 for (
int c=0; c<taille; c++)
515 systeme.
set(l,c) = 0 ;
520 systeme.
set(ligne, colonne) = 1. ;
521 for (
int i=0 ; i<nr ; i++)
522 sec_membre.
set(ligne) -= sol_part_eta(0, k, j, i) ;
525 systeme.
set(ligne, colonne) = l_q;
526 for (
int i=0; i<nr; i++)
527 sec_membre.
set(ligne) -= sol_part_vr(0,k,j,i) ;
530 for (
int zone=1 ; zone<nzm1 ; zone++) {
533 double echelle = mpaff->
get_beta()[zone]/alpha ;
536 systeme.
set(ligne, colonne) = -
pow(echelle-1.,
double(l_q-1)) ;
538 systeme.
set(ligne, colonne+1) = -1/
pow(echelle-1.,
double(l_q+2)) ;
539 for (
int i=0 ; i<nr ; i++)
541 sec_membre.
set(ligne) += sol_part_eta(zone, k, j, i) ;
542 else sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
545 systeme.
set(ligne, colonne) = -l_q*
pow(echelle-1.,
double(l_q-1)) ;
546 systeme.
set(ligne, colonne+1) = (l_q+1)/
pow(echelle-1.,
double(l_q+2));
547 for (
int i=0 ; i<nr ; i++)
549 sec_membre.
set(ligne) += sol_part_vr(zone, k, j, i) ;
550 else sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i) ;
554 systeme.
set(ligne, colonne) =
pow(echelle+1.,
double(l_q-1)) ;
556 systeme.
set(ligne, colonne+1) = 1./
pow(echelle+1.,
double(l_q+2)) ;
557 for (
int i=0 ; i<nr ; i++)
558 sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
561 systeme.
set(ligne, colonne) = l_q*
pow(echelle+1.,
double(l_q-1)) ;
562 systeme.
set(ligne, colonne+1) = -double(l_q+1)
563 /
pow(echelle+1.,
double(l_q+2)) ;
564 for (
int i=0 ; i<nr ; i++)
565 sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i);
574 systeme.
set(ligne, colonne) = -
pow(-2,
double(l_q+2)) ;
575 for (
int i=0 ; i<nr ; i++)
576 if (i%2 == 0) sec_membre.
set(ligne) += sol_part_eta(nzm1, k, j, i) ;
577 else sec_membre.
set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
579 systeme.
set(ligne+1, colonne) = double(l_q+1)*
pow(-2,
double(l_q+2)) ;
580 for (
int i=0 ; i<nr ; i++)
581 if (i%2 == 0) sec_membre.
set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
582 else sec_membre.
set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
587 if (taille > 2) systeme.
set_band(2,2) ;
597 for (
int i=0 ; i<nr ; i++) {
598 cf_eta.
set(0, k, j, i) = sol_part_eta(0, k, j, i)
599 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
600 cf_vr.
set(0, k, j, i) = sol_part_vr(0, k, j, i)
601 +double(l_q)*facteurs(conte)*solution_hom_un(0, k, j, i) ;
604 for (
int zone=1 ; zone<nzm1 ; zone++) {
606 for (
int i=0 ; i<nr ; i++) {
607 cf_eta.
set(zone, k, j, i) =
608 sol_part_eta(zone, k, j, i)
609 +facteurs(conte)*solution_hom_un(zone, k, j, i)
610 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
611 cf_vr.
set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
612 +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
613 -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
618 for (
int i=0 ; i<nr ; i++) {
619 cf_eta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
620 +facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
621 cf_vr.
set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
622 -double(l_q+1)*facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Class for the elementary differential operator (see the base class Diff ).
const double * get_alpha() const
Returns the pointer on the array alpha.
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
void ylm_i()
Inverse of ylm()
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...
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void ylm()
Computes the coefficients of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Tensor field of valence 0 (or component of a tensorial field).
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
const Metric *const met_div
Metric with respect to which the divergence is defined.
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Class for the elementary differential operator (see the base class Diff ).
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Class for the elementary differential operator Identity (see the base class Diff ).
void div_r()
Division by r everywhere; dzpuis is not changed.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Base_val base
Bases on which the spectral expansion is performed.
Scalar ** cmp
Array of size n_comp of pointers onto the components.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
double * t
The array of double.
const double * get_beta() const
Returns the pointer on the array beta.
int get_nzone() const
Returns the number of domains.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Cmp pow(const Cmp &, int)
Power .
Vector_divfree poisson() const
Computes the solution of a vectorial Poisson equation with *this as a source: .
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
double & set(int j, int i)
Read/write of a particuliar element.
void set_band(int up, int low) const
Calculate the band storage of *std.
Spherical orthonormal vectorial bases (triads).
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Bases of the spectral expansions.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Class for the elementary differential operator (see the base class Diff ).
void set_vr_mu(const Scalar &vr_i, const Scalar &mu_i)
Sets the angular potentials (see member p_mu ), and the component of the vector.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Coefficients storage for the multi-domain spectral method.
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Class for the elementary differential operator (see the base class Diff ).
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Valeur & set_spectral_va()
Returns va (read/write version)
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
void annule_hard()
Sets the Tbl to zero in a hard way.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through , and .
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
const Valeur & get_spectral_va() const
Returns va (read only version)