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;
182 S_r.set_spectral_va().ylm() ;
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++) {
344 sr.
set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, 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) ;
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)