28 char vector_df_poisson_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $" ;
100 source_r.
poisson(par_khi, khi) ;
150 for (
int i=0; i<3; i++)
151 assert(
cmp[i]->check_dzpuis(4)) ;
156 assert( mpaff != 0x0) ;
168 if (
cmp[0]->get_etat() == ETATZERO) {
177 int nz = mg.
get_nzone() ;
int nzm1 = nz - 1;
199 double beta = mpaff->
get_beta()[0] ;
200 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
205 for (
int k=0 ; k<np+1 ; k++) {
206 for (
int j=0 ; j<nt ; j++) {
208 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
209 int dege1 = (l_q <3 ? 0 : 1) ;
211 int nr_eq1 = nr0 - dege1 ;
212 int nr_eq2 = nr0 - dege2 ;
213 int nrtot = nr_eq1 + nr_eq2 ;
222 for (
int lin=0; lin<nr_eq1; lin++) {
223 for (
int col=dege1; col<nr0; col++)
224 oper.
set(lin,col-dege1)
225 = md2(lin,col) + 6*mxd(lin,col) + 6*mid(lin,col) ;
226 for (
int col=dege2; col<nr0; col++)
227 oper.
set(lin,col-dege2+nr_eq1) = -mxd(lin,col) -2*mid(lin,col) ;
229 for (
int lin=0; lin<nr0; lin++) {
230 for (
int col=dege1; col<nr0; col++)
231 oper.
set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
232 for (
int col=dege2; col<nr0; col++)
233 oper.
set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) + 4*mid(lin,col) ;
239 for (
int i=0; i<nr_eq1; i++)
241 for (
int i=0; i<nr0; i++)
242 sec_membre.
set(i+nr_eq1) = 0 ;
252 for (
int i=0; i<dege1; i++)
254 for (
int i=dege1; i<nr0; i++)
255 res_eta.
set(i) = big_res(i-dege1) ;
256 res_eta.
set(nr0) = 0 ;
257 for (
int i=0; i<dege2; i++)
259 for (
int i=dege2; i<nr0; i++)
260 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
261 res_vr.
set(nr0) = 0 ;
265 multx2_1d(nr, &res_eta.
t, base_r) ;
266 multx2_1d(nr, &res_vr.
t, base_r) ;
269 Tbl sol_hom = solh(nr, l_q-1, 0., base_r) ;
270 for (
int i=0 ; i<nr ; i++) {
271 sol_part_eta.
set(0, k, j, i) = alpha*alpha*res_eta(i) ;
272 sol_part_vr.
set(0, k, j, i) = alpha*alpha*res_vr(i) ;
273 solution_hom_un.
set(0, k, j, i) = sol_hom(i) ;
274 solution_hom_deux.
set(0, k, j, i) = 0. ;
283 for (
int zone=1 ; zone<nzm1 ; zone++) {
286 assert(nt == mg.
get_nt(zone)) ;
287 assert(np == mg.
get_np(zone)) ;
290 double ech = beta / alpha ;
294 for (
int k=0 ; k<np+1 ; k++) {
295 for (
int j=0 ; j<nt ; j++) {
297 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
300 int nr_eq1 = nr - dege1 ;
301 int nr_eq2 = nr - dege2 ;
302 int nrtot = nr_eq1 + nr_eq2 + 1;
310 Diff_id id(base_r, nr+1) ;
const Matrice& mid =
id.get_matrice() ;
314 for (
int lin=0; lin<nr_eq1; lin++) {
315 for (
int col=dege1; col<nr; col++)
316 oper.
set(lin,col-dege1)
317 = mxd2(lin,col) + ech*md2(lin,col) + 2*md(lin,col) ;
318 for (
int col=dege2; col<nr+1; col++)
319 oper.
set(lin,col-dege2+nr_eq1) = -md(lin,col) ;
321 for (
int lin=0; lin<nr_eq2; lin++) {
322 for (
int col=dege1; col<nr; col++)
323 oper.
set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
324 for (
int col=dege2; col<nr+1; col++)
325 oper.
set(lin+nr_eq1, col-dege2+nr_eq1) =
326 mxd(lin,col) + ech*md(lin,col) + 2*mid(lin,col) ;
330 for (
int col=dege1; col<nr; col++)
331 oper.
set(nrtot-1, col-dege1) = 0 ;
332 for (
int col=dege2; col<nr+1; col++)
333 oper.
set(nrtot-1, col-dege2+nr_eq1) =
334 m2d2(0,col) + ech*(2*mxd2(0,col) + ech*md2(0,col))
335 +4*(mxd(0,col) +ech*md(0,col))
336 +(2 - l_q*(l_q+1))*mid(0,col) ;
343 for (
int i=0; i<5; i++) {
347 for (
int i=5; i<nr; i++)
349 Tbl xsr= sr ;
Tbl x2sr= sr ;
351 multx2_1d(5, &x2sr.
t, base_r) ; multx_1d(5, &xsr.
t, base_r) ;
352 multx_1d(nr, &xseta.
t, base_r) ;
354 for (
int i=0; i<nr_eq1; i++)
355 sec_membre.
set(i) = alpha*(alpha*xseta(i) + beta*seta(i)) ;
356 for (
int i=0; i<nr_eq2; i++)
357 sec_membre.
set(i+nr_eq1) = 0 ;
358 sec_membre.
set(nr_eq1+nr_eq2) = alpha*alpha*x2sr(0) + 2*alpha*beta*xsr(0)
369 for (
int i=0; i<dege1; i++)
371 for (
int i=dege1; i<nr; i++)
372 res_eta.
set(i) = big_res(i-dege1) ;
373 for (
int i=0; i<dege2; i++)
375 for (
int i=dege2; i<nr; i++)
376 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
379 Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
380 Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
381 for (
int i=0 ; i<nr ; i++) {
382 sol_part_eta.
set(zone, k, j, i) = res_eta(i) ;
383 sol_part_vr.
set(zone, k, j, i) = res_vr(i) ;
384 solution_hom_un.
set(zone, k, j, i) = sol_hom1(0,i) ;
385 solution_hom_deux.
set(zone, k, j, i) = sol_hom2(1,i) ;
395 assert(nt == mg.
get_nt(nzm1)) ;
396 assert(np == mg.
get_np(nzm1)) ;
403 for (
int k=0 ; k<np+1 ; k++) {
404 for (
int j=0 ; j<nt ; j++) {
406 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
409 int nr_eq1 = nr0 - dege1 ;
410 int nr_eq2 = nr0 - dege2 ;
411 int nrtot = nr_eq1 + nr_eq2 ;
420 for (
int lin=0; lin<nr_eq1; lin++) {
421 for (
int col=dege1; col<nr0; col++)
422 oper.
set(lin,col-dege1)
423 = m2d2(lin,col) + 4*mxd(lin,col) + 2*mid(lin,col) ;
424 for (
int col=dege2; col<nr0; col++)
425 oper.
set(lin,col-dege2+nr_eq1) =
426 mxd(lin,col) + 2*mid(lin,col) ;
428 for (
int lin=0; lin<nr_eq2; lin++) {
429 for (
int col=dege1; col<nr0; col++)
430 oper.
set(lin+nr_eq1,col-dege1) = l_q*(l_q+1)*mid(lin,col) ;
431 for (
int col=dege2; col<nr0; col++)
432 oper.
set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) ;
438 for (
int i=0; i<nr_eq1; i++)
440 for (
int i=0; i<nr_eq2; i++)
441 sec_membre.
set(i+nr_eq1) = 0 ;
448 for (
int i=0; i<dege1; i++)
450 for (
int i=dege1; i<nr0; i++)
451 res_eta.
set(i) = big_res(i-dege1) ;
452 res_eta.
set(nr0) = 0 ;
453 res_eta.
set(nr0+1) = 0 ;
454 for (
int i=0; i<dege2; i++)
456 for (
int i=dege2; i<nr0; i++)
457 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
458 res_vr.
set(nr0) = 0 ;
459 res_vr.
set(nr0+1) = 0 ;
465 mult2_xm1_1d_cheb(nr, res_eta.
t, res_e2.
t) ;
466 mult2_xm1_1d_cheb(nr, res_vr.
t, res_v2.
t) ;
469 Tbl sol_hom = solh(nr, l_q+1, 0., base_r) ;
470 for (
int i=0 ; i<nr ; i++) {
471 sol_part_eta.
set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
472 sol_part_vr.
set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
473 solution_hom_un.
set(nzm1, k, j, i) = sol_hom(i) ;
474 solution_hom_deux.
set(nzm1, k, j, i) = 0. ;
494 int taille = 2*nzm1 ;
495 Tbl sec_membre(taille) ;
496 Matrice systeme(taille, taille) ;
498 int ligne ;
int colonne ;
502 for (
int k=0 ; k<np+1 ; k++)
503 for (
int j=0 ; j<nt ; j++) {
505 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
510 for (
int l=0; l<taille; l++)
511 for (
int c=0; c<taille; c++)
512 systeme.
set(l,c) = 0 ;
517 systeme.
set(ligne, colonne) = 1. ;
518 for (
int i=0 ; i<nr ; i++)
519 sec_membre.
set(ligne) -= sol_part_eta(0, k, j, i) ;
522 systeme.
set(ligne, colonne) = l_q;
523 for (
int i=0; i<nr; i++)
524 sec_membre.
set(ligne) -= sol_part_vr(0,k,j,i) ;
527 for (
int zone=1 ; zone<nzm1 ; zone++) {
530 double echelle = mpaff->
get_beta()[zone]/alpha ;
533 systeme.
set(ligne, colonne) = -
pow(echelle-1.,
double(l_q-1)) ;
535 systeme.
set(ligne, colonne+1) = -1/
pow(echelle-1.,
double(l_q+2)) ;
536 for (
int i=0 ; i<nr ; i++)
538 sec_membre.
set(ligne) += sol_part_eta(zone, k, j, i) ;
539 else sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
542 systeme.
set(ligne, colonne) = -l_q*
pow(echelle-1.,
double(l_q-1)) ;
543 systeme.
set(ligne, colonne+1) = (l_q+1)/
pow(echelle-1.,
double(l_q+2));
544 for (
int i=0 ; i<nr ; i++)
546 sec_membre.
set(ligne) += sol_part_vr(zone, k, j, i) ;
547 else sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i) ;
551 systeme.
set(ligne, colonne) =
pow(echelle+1.,
double(l_q-1)) ;
553 systeme.
set(ligne, colonne+1) = 1./
pow(echelle+1.,
double(l_q+2)) ;
554 for (
int i=0 ; i<nr ; i++)
555 sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
558 systeme.
set(ligne, colonne) = l_q*
pow(echelle+1.,
double(l_q-1)) ;
559 systeme.
set(ligne, colonne+1) = -double(l_q+1)
560 /
pow(echelle+1.,
double(l_q+2)) ;
561 for (
int i=0 ; i<nr ; i++)
562 sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i);
571 systeme.
set(ligne, colonne) = -
pow(-2,
double(l_q+2)) ;
572 for (
int i=0 ; i<nr ; i++)
573 if (i%2 == 0) sec_membre.
set(ligne) += sol_part_eta(nzm1, k, j, i) ;
574 else sec_membre.
set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
576 systeme.
set(ligne+1, colonne) = double(l_q+1)*
pow(-2,
double(l_q+2)) ;
577 for (
int i=0 ; i<nr ; i++)
578 if (i%2 == 0) sec_membre.
set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
579 else sec_membre.
set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
584 if (taille > 2) systeme.
set_band(2,2) ;
594 for (
int i=0 ; i<nr ; i++) {
595 cf_eta.
set(0, k, j, i) = sol_part_eta(0, k, j, i)
596 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
597 cf_vr.
set(0, k, j, i) = sol_part_vr(0, k, j, i)
598 +double(l_q)*facteurs(conte)*solution_hom_un(0, k, j, i) ;
601 for (
int zone=1 ; zone<nzm1 ; zone++) {
603 for (
int i=0 ; i<nr ; i++) {
604 cf_eta.
set(zone, k, j, i) =
605 sol_part_eta(zone, k, j, i)
606 +facteurs(conte)*solution_hom_un(zone, k, j, i)
607 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
608 cf_vr.
set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
609 +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
610 -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
615 for (
int i=0 ; i<nr ; i++) {
616 cf_eta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
617 +facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
618 cf_vr.
set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
619 -double(l_q+1)*facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
Bases of the spectral expansions.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Spherical orthonormal vectorial bases (triads).
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator Identity (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
const double * get_beta() const
Returns the pointer on the array beta.
const double * get_alpha() const
Returns the pointer on the array alpha.
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.
Coefficients storage for the multi-domain spectral method.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Tensor field of valence 0 (or component of a tensorial field).
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void div_r()
Division by r everywhere; dzpuis is not changed.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
const Valeur & get_spectral_va() const
Returns va (read only version)
Valeur & set_spectral_va()
Returns va (read/write version)
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
void annule_hard()
Sets the Tbl to zero in a hard way.
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).
double * t
The array of double.
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
void ylm()
Computes the coefficients of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
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.
Vector_divfree poisson() const
Computes the solution of a vectorial Poisson equation with *this as a source:
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through , and .
const Metric *const met_div
Metric with respect to which the divergence is defined.
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
Cmp pow(const Cmp &, int)
Power .
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Scalar ** cmp
Array of size n_comp of pointers onto the components.
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.