28 char map_af_poisson_regu_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson_regu.C,v 1.5 2014/10/13 08:53:03 j_novak Exp $" ;
110 double unsgam1,
Param& par,
Cmp& uu,
112 Cmp& source_regu,
Cmp& source_div)
const {
115 assert(source.
get_etat() != ETATNONDEF) ;
119 double aa = unsgam1 ;
126 assert(nr-k_div > 0) ;
136 assert(sourva.
get_etat() == ETATQCQ) ;
140 rho = *(sourva.
c_cf) ;
150 Tbl& ccf = *((rho.
c_cf)->t[0]) ;
152 Tbl nilm_cos(np/2+1, 2*nt, nr) ;
153 Tbl nilm_sin(np/2+1, 2*nt, nr) ;
158 for (
int k=0; k<=np; k+=2) {
162 for (
int j=0; j<nt; j++) {
171 for (
int i=0; i<nr; i++) {
173 nilm_cos.
set(m, l, i) = ccf(k, j, i) ;
174 nilm_sin.
set(m, l, i) = ccf(k+1, j, i) ;
186 Tbl cf_cil(2*nt, nr, k_div) ;
199 double* tmp1 =
new double[nr] ;
201 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
203 for (
int l=0; l<2*nt; l++) {
205 for (
int i=0; i<nr; i++) {
207 double xi = grid.
x[i] ;
209 tmp1[i] = (aa + k_deg + 1.) *
210 ( -(4. * l + 6.) *
pow(1. - xi * xi, aa + k_deg) *
pow(xi, l)
211 + 4. * (aa + k_deg) *
pow(1. - xi * xi, aa + k_deg - 1.) *
217 cfrchebp(deg, dim, tmp1, dim, tmp1) ;
220 cfrchebi(deg, dim, tmp1, dim, tmp1) ;
222 for (
int i=0; i<nr; i++) {
224 cf_cil.
set(l, i, k_deg-1) = tmp1[i] ;
233 Tbl alm_cos(np/2+1, 2*nt, k_div) ;
234 Tbl alm_sin(np/2+1, 2*nt, k_div) ;
248 for (
int k=0; k<=np; k+=2) {
252 for (
int j=0; j<nt; j++) {
261 for (
int i=0; i<k_div; i++) {
262 for (
int j2=0; j2<k_div; j2++) {
263 matrix.
set(i, j2) = cf_cil(l, nr-1-i, j2) ;
271 for (
int i=0; i<k_div; i++) {
272 rhs_cos.
set(i) = nilm_cos(m, l, nr-1-i) ;
273 rhs_sin.
set(i) = nilm_sin(m, l, nr-1-i) ;
279 for (
int i=0; i<k_div; i++) {
280 alm_cos.
set(m, l, i) = sol_cos(i) ;
281 alm_sin.
set(m, l, i) = sol_sin(i) ;
291 (source_div.
va).set_etat_cf_qcq() ;
292 (source_div.
va).c_cf->set_etat_qcq() ;
293 (source_div.
va).c_cf->t[0]->set_etat_qcq() ;
299 for (
int k=0; k<=np; k+=2) {
300 for (
int j=0; j<nt; j++) {
301 for (
int i=0; i<nr; i++) {
302 sdva_cf.
set(0, k, j, i) = 0 ;
303 sdva_cf.
set(0, k+1, j, i) = 0 ;
308 for (
int k=0; k<=np; k+=2) {
312 for (
int j=0; j<nt; j++) {
322 for (
int i=0; i<nr; i++) {
324 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
326 sdva_cf.
set(0, k, j, i) = sdva_cf(0, k, j, i)
327 + alm_cos(m, l, k_deg-1) * cf_cil(l, i, k_deg-1) ;
328 sdva_cf.
set(0, k+1, j, i) = sdva_cf(0, k+1, j, i)
329 + alm_sin(m, l, k_deg-1) * cf_cil(l, i, k_deg-1) ;
336 source_div.
annule(nzet, nzm1) ;
341 for (
int l=0; l<=nzm1; l++) {
342 int base_s_div_r = base_std.
b[l] &
MSQ_R ;
343 int base_s_div_p = base_std.
b[l] &
MSQ_P ;
345 base_s_div.
b[l] = base_s_div_r |
T_LEG_P | base_s_div_p ;
348 sdva_cf.
base = base_s_div ;
357 source_regu = source - source_div ;
367 (*this).poisson(source_regu, par, uu_regu) ;
373 Tbl cf_pil(2*nt, nr, k_div) ;
376 double* tmp2 =
new double[nr] ;
378 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
380 for (
int l=0; l<2*nt; l++) {
382 for (
int i=0; i<nr; i++) {
384 double xi = grid.
x[i] ;
385 tmp2[i] =
pow(xi, l) *
pow(1. - xi * xi, aa + 1. + k_deg) ;
390 cfrchebp(deg, dim, tmp2, dim, tmp2) ;
393 cfrchebi(deg, dim, tmp2, dim, tmp2) ;
395 for (
int i=0; i<nr; i++) {
397 cf_pil.
set(l, i, k_deg-1) = tmp2[i] ;
404 (uu_div.
va).set_etat_cf_qcq() ;
405 ((uu_div.
va).c_cf)->set_etat_qcq() ;
406 ((uu_div.
va).c_cf)->t[0]->set_etat_qcq() ;
412 for (
int k=0; k<=np; k+=2) {
413 for (
int j=0; j<nt; j++) {
414 for (
int i=0; i<nr; i++) {
415 udva_cf.
set(0, k, j, i) = 0 ;
416 udva_cf.
set(0, k+1, j, i) = 0 ;
421 for (
int k=0; k<=np; k+=2) {
425 for (
int j=0; j<nt; j++) {
435 for (
int i=0; i<nr; i++) {
437 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
439 udva_cf.
set(0, k, j, i) = udva_cf(0, k, j, i)
440 + alm_cos(m, l, k_deg-1) * cf_pil(l, i, k_deg-1) ;
441 udva_cf.
set(0, k+1, j, i) = udva_cf(0, k+1, j, i)
442 + alm_sin(m, l, k_deg-1) * cf_pil(l, i, k_deg-1) ;
449 uu_div.
annule(nzet, nzm1) ;
452 for (
int l=0; l<=nzm1; l++) {
453 int base_uu_r = base_std.
b[l] &
MSQ_R ;
454 int base_uu_p = base_std.
b[l] &
MSQ_P ;
456 base_uu_div.
b[l] = base_uu_r |
T_LEG_P | base_uu_p ;
459 udva_cf.
base = base_uu_div ;
473 uu = uu_regu + uu_div ;
482 (duu_div.
set(0).
va).set_etat_cf_qcq() ;
483 ((duu_div.
set(0).
va).c_cf)->set_etat_qcq() ;
484 ((duu_div.
set(0).
va).c_cf)->t[0]->set_etat_qcq() ;
487 (duu_div.
set(1).
va).set_etat_cf_qcq() ;
488 ((duu_div.
set(1).
va).c_cf)->set_etat_qcq() ;
489 ((duu_div.
set(1).
va).c_cf)->t[0]->set_etat_qcq() ;
492 (duu_div.
set(2).
va).set_etat_cf_qcq() ;
493 ((duu_div.
set(2).
va).c_cf)->set_etat_qcq() ;
494 ((duu_div.
set(2).
va).c_cf)->t[0]->set_etat_qcq() ;
508 Tbl cf_dril(2*nt, nr, k_div) ;
511 double* tmp3 =
new double[nr] ;
513 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
515 for (
int i=0; i<nr; i++) {
517 double xi = grid.
x[i] ;
518 tmp3[i] = -2. * (aa + 1. + k_deg) * xi
519 *
pow(1. - xi * xi, aa + k_deg) ;
523 cfrchebi(deg, dim, tmp3, dim, tmp3) ;
525 for (
int i=0; i<nr; i++) {
527 cf_dril.
set(0, i, k_deg-1) = tmp3[i] ;
531 for (
int l=1; l<2*nt; l++) {
533 for (
int i=0; i<nr; i++) {
535 double xi = grid.
x[i] ;
536 tmp3[i] = l *
pow(xi, l - 1.) *
pow(1. - xi * xi, aa + 1. + k_deg)
537 -2. * (aa + 1. + k_deg) *
pow(xi, l + 1.)
538 *
pow(1. - xi * xi, aa + k_deg) ;
543 cfrchebi(deg, dim, tmp3, dim, tmp3) ;
546 cfrchebp(deg, dim, tmp3, dim, tmp3) ;
548 for (
int i=0; i<nr; i++) {
550 cf_dril.
set(l, i, k_deg-1) = tmp3[i] ;
557 for (
int k=0; k<=np; k+=2) {
558 for (
int j=0; j<nt; j++) {
559 for (
int i=0; i<nr; i++) {
560 vr_cf.
set(0, k, j, i) = 0 ;
561 vr_cf.
set(0, k+1, j, i) = 0 ;
566 for (
int k=0; k<=np; k+=2) {
570 for (
int j=0; j<nt; j++) {
580 for (
int i=0; i<nr; i++) {
582 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
584 vr_cf.
set(0, k, j, i) = vr_cf(0, k, j, i)
585 + alm_cos(m, l, k_deg-1) * cf_dril(l, i, k_deg-1) ;
586 vr_cf.
set(0, k+1, j, i) = vr_cf(0, k+1, j, i)
587 + alm_sin(m, l, k_deg-1) * cf_dril(l, i, k_deg-1) ;
594 (duu_div.
set(0)).annule(nzet, nzm1) ;
600 for (
int l=0; l<=nzm1; l++) {
601 int base_duu_r_p = base_std.
b[l] &
MSQ_P ;
606 vr_cf.
base = base_duu_div_r ;
615 vr = duu_div(0).va *
alpha[0] *
alpha[0] * RR ;
616 vr.
base = sauve_base ;
622 Tbl cf_dpil(2*nt, nr, k_div) ;
625 double* tmp4 =
new double[nr] ;
627 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
629 for (
int i=0; i<nr; i++) {
633 cfrchebi(deg, dim, tmp4, dim, tmp4) ;
635 for (
int i=0; i<nr; i++) {
637 cf_dpil.
set(0, i, k_deg-1) = tmp4[i] ;
641 for (
int l=1; l<2*nt; l++) {
643 for (
int i=0; i<nr; i++) {
645 double xi = grid.
x[i] ;
646 tmp4[i] =
pow(xi, l - 1.) *
pow(1. - xi * xi, aa + 1. + k_deg) ;
651 cfrchebi(deg, dim, tmp4, dim, tmp4) ;
654 cfrchebp(deg, dim, tmp4, dim, tmp4) ;
656 for (
int i=0; i<nr; i++) {
658 cf_dpil.
set(l, i, k_deg-1) = tmp4[i] ;
665 for (
int k=0; k<=np; k+=2) {
666 for (
int j=0; j<nt; j++) {
667 for (
int i=0; i<nr; i++) {
668 vt_cf.
set(0, k, j, i) = 0 ;
669 vt_cf.
set(0, k+1, j, i) = 0 ;
670 vp_cf.
set(0, k, j, i) = 0 ;
671 vp_cf.
set(0, k+1, j, i) = 0 ;
676 for (
int k=0; k<=np; k+=2) {
680 for (
int j=0; j<nt; j++) {
690 for (
int i=0; i<nr; i++) {
692 for (
int k_deg=1; k_deg<=k_div; k_deg++) {
694 vt_cf.
set(0, k, j, i) = vt_cf(0, k, j, i)
695 + alm_cos(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
696 vt_cf.
set(0, k+1, j, i) = vt_cf(0, k+1, j, i)
697 + alm_sin(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
699 vp_cf.
set(0, k, j, i) = vp_cf(0, k, j, i)
700 + alm_cos(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
701 vp_cf.
set(0, k+1, j, i) = vp_cf(0, k+1, j, i)
702 + alm_sin(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
709 (duu_div.
set(1)).annule(nzet, nzm1) ;
710 (duu_div.
set(2)).annule(nzet, nzm1) ;
717 for (
int l=0; l<=nzm1; l++) {
718 int base_duu_p_p = base_std.
b[l] &
MSQ_P ;
723 vt_cf.
base = base_duu_div_p ;
731 for (
int l=0; l<=nzm1; l++) {
732 int base_duu_t_p = base_std.
b[l] &
MSQ_P ;
737 vp_cf.
base = base_duu_div_t ;
744 vt = (duu_div(1).va).
dsdt() ;
746 vp = (duu_div(2).va).
stdsdp() ;
751 vt = duu_div(1).va *
alpha[0] ;
753 vp = duu_div(2).va *
alpha[0] ;
758 duu_div.
set_triad( (*this).get_bvect_spher() ) ;
Bases of the spectral expansions.
int * b
Array (size: nzone ) of the spectral basis in each domain.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
int get_etat() const
Returns the logical state.
Valeur va
The numerical value of the Cmp
void annule(int l)
Sets the Cmp to zero in a given domain.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Map * get_mp() const
Returns the mapping.
void set_dzpuis(int)
Set a value to dzpuis.
Active physical coordinates and mapping derivatives.
3D grid class in one domain.
double * x
Array of values of at the nr collocation points.
virtual void poisson_regular(const Cmp &source, int k_div, int nzet, double unsgam1, Param &par, Cmp &uu, Cmp &uu_regu, Cmp &uu_div, Tenseur &duu_div, Cmp &source_regu, Cmp &source_div) const
Computes the solution of a scalar Poisson equation.
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
virtual void dsdt(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
virtual void stdsdp(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this 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.
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for 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.
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
Coefficients storage for the multi-domain spectral method.
Base_val base
Bases of the spectral expansions.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
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).
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Values and coefficients of a (real-value) function.
int get_etat() const
Returns the logical state.
void ylm()
Computes the coefficients of *this.
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
Cmp pow(const Cmp &, int)
Power .
#define MSQ_R
Extraction de l'info sur R.
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define T_LEG_P
fct. de Legendre associees paires
#define MSQ_P
Extraction de l'info sur Phi.