78 void Vector::poisson_block(
double lam, Vector& resu)
const {
80 const Map_af* mpaff =
dynamic_cast<const Map_af*
>(
mp) ;
82 for (
int i=0; i<3; i++)
83 assert(
cmp[i]->check_dzpuis(4)) ;
85 const Base_vect_spher* bvs =
dynamic_cast<const Base_vect_spher*
>(
triad) ;
88 assert( mpaff != 0x0) ;
91 if (fabs(lam + 1.) < 1.e-8) {
93 Vector_divfree sou(*
mp, *
triad, mets) ;
94 for (
int i=1; i<=3; i++) sou.set(i) = *
cmp[i-1] ;
101 const Mg3d& mg = *mpaff->get_mg() ;
102 int nz = mg.get_nzone() ;
int nzm1 = nz - 1;
103 int nt = mg.get_nt(0) ;
104 int np = mg.get_np(0) ;
105 assert(mg.get_type_r(0) == RARE) ;
106 assert(mg.get_type_r(nzm1) == UNSURR) ;
107 Scalar S_r = *
cmp[0] ;
108 Scalar S_eta =
eta() ;
111 bool all_zero = false ;
112 if (S_r.get_etat() == ETATZERO) {
113 if (S_eta.get_etat() == ETATZERO) {
115 het.set_etat_zero() ;
120 S_r.set_spectral_base(S_eta.get_spectral_base()) ;
123 if ((S_eta.get_etat() == ETATZERO)&&(!all_zero)) {
124 S_eta.annule_hard() ;
125 S_eta.set_spectral_base(S_r.get_spectral_base()) ;
128 S_r.set_spectral_va().ylm() ;
129 S_eta.set_spectral_va().ylm() ;
130 const Base_val& base = S_eta.get_spectral_va().base ;
131 Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
132 Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
133 Mtbl_cf sol_hom_un_eta(mg, base) ; sol_hom_un_eta.annule_hard() ;
134 Mtbl_cf sol_hom_deux_eta(mg, base) ; sol_hom_deux_eta.annule_hard() ;
135 Mtbl_cf sol_hom_trois_eta(mg, base) ; sol_hom_trois_eta.annule_hard() ;
136 Mtbl_cf sol_hom_quatre_eta(mg, base) ; sol_hom_quatre_eta.annule_hard() ;
137 Mtbl_cf sol_hom_un_vr(mg, base) ; sol_hom_un_vr.annule_hard() ;
138 Mtbl_cf sol_hom_deux_vr(mg, base) ; sol_hom_deux_vr.annule_hard() ;
139 Mtbl_cf sol_hom_trois_vr(mg, base) ; sol_hom_trois_vr.annule_hard() ;
140 Mtbl_cf sol_hom_quatre_vr(mg, base) ; sol_hom_quatre_vr.annule_hard() ;
144 Scalar sou_l0 = (*
cmp[0]) / (1. + lam) ;
145 sou_l0.set_spectral_va().ylm() ;
146 Scalar vrl0 = pois_vect_r0(sou_l0) ;
153 int nr = mg.get_nr(0) ;
154 double alpha = mpaff->get_alpha()[0] ;
double alp2 = alpha*alpha ;
155 double beta = mpaff->get_beta()[0] ;
156 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
160 for (
int k=0 ; k<np+1 ; k++) {
161 for (
int j=0 ; j<nt ; j++) {
162 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
163 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0) ) {
164 int aa = 0 ;
int bb = 0 ;
int nr0 = 0 ;
176 Matrice oper(2*nr, 2*nr) ; oper.set_etat_qcq() ;
177 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
178 Diff_dsdx2 d2(base_r, nr) ;
const Matrice& md2 = d2.get_matrice() ;
179 Diff_sxdsdx xd(base_r, nr) ;
const Matrice& mxd = xd.get_matrice() ;
180 Diff_sx2 s2(base_r, nr) ;
const Matrice& ms2 = s2.get_matrice() ;
184 for (
int lin=0; lin<nr; lin++) {
185 for (
int col=0; col<nr; col++)
187 (md2(lin,col) + 2*mxd(lin,col)
188 -(lam+1)*l_q*(l_q+1)*ms2(lin,col)) / alp2 ;
189 for (
int col=0; col<nr; col++)
190 oper.set(lin,col+nr) =
191 (lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
193 for (
int lin=0; lin<nr; lin++) {
194 for (
int col=0; col<nr; col++)
195 oper.set(lin+nr,col) =
196 (-lam*l_q*(l_q+1)*mxd(lin,col)
197 +(lam+2)*l_q*(l_q+1)*ms2(lin,col)) / alp2 ;
198 for (
int col=0; col<nr; col++)
199 oper.set(lin+nr, col+nr) =
200 ((lam+1)*(md2(lin,col) + 2*mxd(lin,col))
201 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
203 bool pb_eta = ( ( fabs( lam*
double(l_q+3) + 2 ) < 0.01) && (l_q <=2) ) ;
205 for (
int col=0; col<nr; col++)
206 oper.set(nr0-1, col) = 1 ;
207 for (
int col=0; col<nr; col++)
208 oper.set(nr0-1,nr+col) = 0 ;
210 if ((l_q > 2)||pb_eta) {
211 for (
int col=0; col<nr; col++)
212 oper.set(nr+nr0-1, col) = 0 ;
213 for (
int col=0; col<nr; col++)
214 oper.set(nr+nr0-1,nr+col) = 1 ;
217 Matrice op2(nrtot, nrtot) ;
219 for (
int i=0; i<nr0; i++) {
220 for (
int col=0; col<nr0;col++)
221 op2.set(i,col) = (aa*col+bb)*oper(i,col+1)
222 + (aa*(col+1)+bb)*oper(i,col) ;
223 for (
int col=0; col<nr0;col++)
224 op2.set(i,col+nr0) = (aa*col+bb)*oper(i,col+nr+1)
225 + (aa*(col+1)+bb)*oper(i,col+nr) ;
228 for (
int i=nr0; i<nrtot; i++) {
229 for (
int col=0; col<nr0;col++)
230 op2.set(i,col) = (aa*col+bb)*oper(i+d0,col+1)
231 + (aa*(col+1)+bb)*oper(i+d0,col) ;
232 for (
int col=0; col<nr0;col++)
233 op2.set(i,col+nr0) = (aa*col+bb)*oper(i+d0,col+nr+1)
234 + (aa*(col+1)+bb)*oper(i+d0,col+nr) ;
240 for (
int i=0; i<nr0; i++)
241 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(0, k, j, i) ;
242 if (!pb_eta) sec_membre.set(nr0-1) = 0 ;
243 for (
int i=0; i<nr0; i++)
244 sec_membre.set(i+nr0)
245 = (*S_r.get_spectral_va().c_cf)(0, k, j, i) ;
246 if ((l_q > 2)||pb_eta) sec_membre.set(nrtot-1) = 0 ;
250 Tbl big_res = op2.inverse(sec_membre) ;
254 sol_part_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
255 for (
int i=1; i<nr0; i++)
256 sol_part_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
257 + (aa*(i+1)+bb)*big_res(i);
258 sol_part_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1) + bb) ;
259 sol_part_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
260 for (
int i=1; i<nr0; i++)
261 sol_part_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
262 + (aa*(i+1)+bb)*big_res(nr0+i);
263 sol_part_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
265 sol_part_eta.set(0, k, j, nr-1) = 0 ;
266 sol_part_vr.set(0, k, j, nr-1) = 0 ;
271 double fac_eta = lam*double(l_q+3) + 2. ;
272 double fac_vr = double(l_q+1)*(lam*l_q - 2.) ;
273 Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
274 Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
275 for (
int i=0 ; i<nr ; i++) {
276 sol_hom_un_eta.set(0, k, j, i) = sol_hom1(i) ;
277 sol_hom_un_vr.set(0, k, j, i) = l_q*sol_hom1(i) ;
278 sol_hom_trois_eta.set(0, k, j, i) = fac_eta*sol_hom2(i) ;
279 sol_hom_trois_vr.set(0, k, j, i) = fac_vr*sol_hom2(i) ;
283 for (
int i=0; i<nrtot; i++) sec_membre.set(i) = 0 ;
285 sec_membre.set(nr0-1) = 1 ;
286 big_res = op2.inverse(sec_membre) ;
288 sol_hom_un_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
289 for (
int i=1; i<nr0; i++)
290 sol_hom_un_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
291 + (aa*(i+1)+bb)*big_res(i);
292 sol_hom_un_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1) + bb) ;
293 sol_hom_un_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
294 for (
int i=1; i<nr0; i++)
295 sol_hom_un_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
296 + (aa*(i+1)+bb)*big_res(nr0+i);
297 sol_hom_un_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
299 sol_hom_un_eta.set(0, k, j, nr-1) = 0 ;
300 sol_hom_un_vr.set(0, k, j, nr-1) = 0 ;
303 sec_membre.set(nr0-1) = 0 ;
304 sec_membre.set(nrtot-1) = 1 ;
305 big_res = op2.inverse(sec_membre) ;
307 sol_hom_trois_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
308 for (
int i=1; i<nr0; i++)
309 sol_hom_trois_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
310 + (aa*(i+1)+bb)*big_res(i);
311 sol_hom_trois_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1)+bb) ;
312 sol_hom_trois_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
313 for (
int i=1; i<nr0; i++)
314 sol_hom_trois_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
315 + (aa*(i+1)+bb)*big_res(nr0+i);
316 sol_hom_trois_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
318 sol_hom_trois_eta.set(0, k, j, nr-1) = 0 ;
319 sol_hom_trois_vr.set(0, k, j, nr-1) = 0 ;
330 for (
int zone=1 ; zone<nzm1 ; zone++) {
331 nr = mg.get_nr(zone) ;
333 assert(nt == mg.get_nt(zone)) ;
334 assert(np == mg.get_np(zone)) ;
335 alpha = mpaff->get_alpha()[zone] ;
336 beta = mpaff->get_beta()[zone] ;
337 double ech = beta / alpha ;
341 for (
int k=0 ; k<np+1 ; k++) {
342 for (
int j=0 ; j<nt ; j++) {
343 base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
344 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
347 int nr_eq1 = nr - dege1 ;
348 int nr_eq2 = nr - dege2 ;
350 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
351 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
352 Diff_x2dsdx2 x2d2(base_r, nr);
const Matrice& m2d2 = x2d2.get_matrice() ;
353 Diff_xdsdx2 xd2(base_r, nr) ;
const Matrice& mxd2 = xd2.get_matrice() ;
354 Diff_dsdx2 d2(base_r, nr) ;
const Matrice& md2 = d2.get_matrice() ;
355 Diff_xdsdx xd(base_r, nr) ;
const Matrice& mxd = xd.get_matrice() ;
356 Diff_dsdx d1(base_r, nr) ;
const Matrice& md = d1.get_matrice() ;
357 Diff_id id(base_r, nr) ;
const Matrice& mid =
id.get_matrice() ;
361 for (
int lin=0; lin<nr_eq1; lin++) {
362 for (
int col=0; col<nr; col++)
364 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
365 + 2*(mxd(lin,col) + ech*md(lin,col))
366 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
367 for (
int col=0; col<nr; col++)
369 = lam*(mxd(lin,col) + ech*md(lin,col))
370 + 2*(1+lam)*mid(lin,col) ;
372 oper.set(nr_eq1, 0) = 1 ;
373 for (
int col=1; col<2*nr; col++)
374 oper.set(nr_eq1, col) = 0 ;
375 oper.set(nr_eq1+1, 0) = 0 ;
376 oper.set(nr_eq1+1, 1) = 1 ;
377 for (
int col=2; col<2*nr; col++)
378 oper.set(nr_eq1+1, col) = 0 ;
379 for (
int lin=0; lin<nr_eq2; lin++) {
380 for (
int col=0; col<nr; col++)
382 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
383 - (lam+2)*mid(lin,col)) ;
384 for (
int col=0; col<nr; col++)
385 oper.set(lin+nr, col+nr)
386 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
387 + ech*ech*md2(lin,col)
388 + 2*(mxd(lin,col) + ech*md(lin,col)))
389 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
391 for (
int col=0; col<nr; col++)
392 oper.set(nr+nr_eq2, col) = 0 ;
393 oper.set(nr+nr_eq2, nr) = 1 ;
394 for (
int col=nr+1; col<2*nr; col++)
395 oper.set(nr+nr_eq2, col) = 0 ;
396 for (
int col=0; col<nr+1; col++)
397 oper.set(nr+nr_eq2+1, col) = 0 ;
398 oper.set(nr+nr_eq2+1, nr+1) = 1 ;
399 for (
int col=nr+2; col<2*nr; col++)
400 oper.set(nr+nr_eq2+1, col) = 0 ;
406 Tbl sr(nr) ; sr.set_etat_qcq() ;
407 Tbl seta(nr) ; seta.set_etat_qcq() ;
408 for (
int i=0; i<nr; i++) {
409 sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
410 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
412 Tbl xsr= sr ; Tbl x2sr= sr ;
413 Tbl xseta= seta ; Tbl x2seta = seta ;
414 multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
415 multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
417 for (
int i=0; i<nr_eq1; i++)
418 sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
420 sec_membre.set(nr_eq1) = 0 ;
421 sec_membre.set(nr_eq1+1) = 0 ;
422 for (
int i=0; i<nr_eq2; i++)
423 sec_membre.set(i+nr) = beta*beta*sr(i)
424 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
425 sec_membre.set(nr+nr_eq2) = 0 ;
426 sec_membre.set(nr+nr_eq2+1) = 0 ;
430 Tbl big_res = oper.inverse(sec_membre) ;
431 for (
int i=0; i<nr; i++) {
432 sol_part_eta.set(zone, k, j, i) = big_res(i) ;
433 sol_part_vr.set(zone, k, j, i) = big_res(nr+i) ;
438 sec_membre.annule_hard() ;
439 sec_membre.set(nr_eq1) = 1 ;
440 big_res = oper.inverse(sec_membre) ;
441 for (
int i=0 ; i<nr ; i++) {
442 sol_hom_un_eta.set(zone, k, j, i) = big_res(i) ;
443 sol_hom_un_vr.set(zone, k, j, i) = big_res(nr+i) ;
445 sec_membre.set(nr_eq1) = 0 ;
446 sec_membre.set(nr_eq1+1) = 1 ;
447 big_res = oper.inverse(sec_membre) ;
448 for (
int i=0 ; i<nr ; i++) {
449 sol_hom_deux_eta.set(zone, k, j, i) = big_res(i) ;
450 sol_hom_deux_vr.set(zone, k, j, i) = big_res(nr+i) ;
452 sec_membre.set(nr_eq1+1) = 0 ;
453 sec_membre.set(nr+nr_eq2) = 1 ;
454 big_res = oper.inverse(sec_membre) ;
455 for (
int i=0 ; i<nr ; i++) {
456 sol_hom_trois_eta.set(zone, k, j, i) = big_res(i) ;
457 sol_hom_trois_vr.set(zone, k, j, i) = big_res(nr+i) ;
459 sec_membre.set(nr+nr_eq2) = 0 ;
460 sec_membre.set(nr+nr_eq2+1) = 1 ;
461 big_res = oper.inverse(sec_membre) ;
462 for (
int i=0 ; i<nr ; i++) {
463 sol_hom_quatre_eta.set(zone, k, j, i) = big_res(i) ;
464 sol_hom_quatre_vr.set(zone, k, j, i) = big_res(nr+i) ;
474 nr = mg.get_nr(nzm1) ;
475 assert(nt == mg.get_nt(nzm1)) ;
476 assert(np == mg.get_np(nzm1)) ;
477 alpha = mpaff->get_alpha()[nzm1] ; alp2 = alpha*alpha ;
482 for (
int k=0 ; k<np+1 ; k++) {
483 for (
int j=0 ; j<nt ; j++) {
484 base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
485 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
487 int dege2 = (l_q == 1 ? 2 : 3);
488 int nr_eq1 = nr - dege1 ;
489 int nr_eq2 = nr - dege2 ;
491 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
492 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
493 Diff_dsdx2 d2(base_r, nr) ;
const Matrice& md2 = d2.get_matrice() ;
494 Diff_sxdsdx sxd(base_r, nr) ;
const Matrice& mxd = sxd.get_matrice() ;
495 Diff_sx2 sx2(base_r, nr) ;
const Matrice& ms2 = sx2.get_matrice() ;
499 for (
int lin=0; lin<nr_eq1; lin++) {
500 for (
int col=0; col<nr; col++)
502 = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
503 for (
int col=0; col<nr; col++)
504 oper.set(lin,col+nr) =
505 (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
507 for (
int col=0; col<nr; col++)
508 oper.set(nr_eq1, col) = 1 ;
509 for (
int col=nr; col<nrtot; col++)
510 oper.set(nr_eq1, col) = 0 ;
512 for (
int col=0; col<nr; col++) {
513 oper.set(nr_eq1+1, col) = pari*col*col ;
516 for (
int col=nr; col<nrtot; col++)
517 oper.set(nr_eq1+1, col) = 0 ;
518 oper.set(nr_eq1+2, 0) = 1 ;
519 for (
int col=1; col<nrtot; col++)
520 oper.set(nr_eq1+2, col) = 0 ;
521 for (
int lin=0; lin<nr_eq2; lin++) {
522 for (
int col=0; col<nr; col++)
523 oper.set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
524 + (lam+2)*ms2(lin,col))) / alp2 ;
525 for (
int col=0; col<nr; col++)
526 oper.set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
527 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
529 for (
int col=0; col<nr; col++)
530 oper.set(nr+nr_eq2, col) = 0 ;
531 for (
int col=nr; col<nrtot; col++)
532 oper.set(nr+nr_eq2, col) = 1 ;
533 for (
int col=0; col<nr; col++)
534 oper.set(nr+nr_eq2+1, col) = 0 ;
536 for (
int col=0; col<nr; col++) {
537 oper.set(nr+nr_eq2+1, nr+col) = pari*col*col ;
541 for (
int col=0; col<nr; col++)
542 oper.set(nr+nr_eq2+2, col) = 0 ;
543 oper.set(nr+nr_eq2+2, nr) = 1 ;
544 for (
int col=nr+1; col<nrtot; col++)
545 oper.set(nr+nr_eq2+2, col) = 0 ;
551 for (
int i=0; i<nr_eq1; i++)
552 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
553 for (
int i=nr_eq1; i<nr; i++)
554 sec_membre.set(i) = 0 ;
555 for (
int i=0; i<nr_eq2; i++)
556 sec_membre.set(i+nr) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
557 for (
int i=nr_eq2; i<nr; i++)
558 sec_membre.set(nr+i) = 0 ;
559 Tbl big_res = oper.inverse(sec_membre) ;
563 for (
int i=0; i<nr; i++) {
564 sol_part_eta.set(nzm1, k, j, i) = big_res(i) ;
565 sol_part_vr.set(nzm1, k, j, i) = big_res(i+nr) ;
571 Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
572 Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
573 double fac_eta = lam + 2. ;
574 double fac_vr = 2*lam + 2. ;
575 for (
int i=0 ; i<nr ; i++) {
576 sol_hom_deux_eta.set(nzm1, k, j, i) = sol_hom2(i) ;
577 sol_hom_quatre_eta.set(nzm1, k, j, i) = fac_eta*sol_hom1(i) ;
578 sol_hom_deux_vr.set(nzm1, k, j, i) = -2*sol_hom2(i) ;
579 sol_hom_quatre_vr.set(nzm1, k, j, i) = fac_vr*sol_hom1(i) ;
583 sec_membre.annule_hard() ;
584 sec_membre.set(nr-1) = 1 ;
585 big_res = oper.inverse(sec_membre) ;
587 for (
int i=0; i<nr; i++) {
588 sol_hom_deux_eta.set(nzm1, k, j, i) = big_res(i) ;
589 sol_hom_deux_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
591 sec_membre.set(nr-1) = 0 ;
592 sec_membre.set(2*nr-1) = 1 ;
593 big_res = oper.inverse(sec_membre) ;
594 for (
int i=0; i<nr; i++) {
595 sol_hom_quatre_eta.set(nzm1, k, j, i) = big_res(i) ;
596 sol_hom_quatre_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
608 vr.set_spectral_base(base) ;
609 vr.set_spectral_va().set_etat_cf_qcq() ;
610 Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
611 cf_vr.annule_hard() ;
613 het.set_spectral_base(base) ;
614 het.set_spectral_va().set_etat_cf_qcq() ;
615 Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
616 cf_eta.annule_hard() ;
617 int taille = 4*nzm1 ;
618 Tbl sec_membre(taille) ;
619 Matrice systeme(taille, taille) ;
620 systeme.set_etat_qcq() ;
621 int ligne ;
int colonne ;
625 for (
int k=0 ; k<np+1 ; k++)
626 for (
int j=0 ; j<nt ; j++) {
627 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
628 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
632 systeme.annule_hard() ;
633 sec_membre.annule_hard() ;
637 alpha = mpaff->get_alpha()[0] ;
639 systeme.set(ligne, colonne) = sol_hom_un_eta.val_out_bound_jk(0, j, k) ;
640 systeme.set(ligne, colonne+1) = sol_hom_trois_eta.val_out_bound_jk(0, j, k) ;
641 sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
644 systeme.set(ligne, colonne) = sol_hom_un_vr.val_out_bound_jk(0, j, k) ;
645 systeme.set(ligne, colonne+1) = sol_hom_trois_vr.val_out_bound_jk(0, j, k) ;
646 sec_membre.set(ligne) = -sol_part_vr.val_out_bound_jk(0,j,k) ;
650 int pari = (base_r ==
R_CHEBP ? 0 : 1) ;
651 for (
int i=0; i<nr; i++) {
652 systeme.set(ligne, colonne)
653 += (2*i+pari)*(2*i+pari)*sol_hom_un_eta(0, k, j, i)/alpha ;
654 systeme.set(ligne, colonne+1)
655 += (2*i+pari)*(2*i+pari)*sol_hom_trois_eta(0, k, j, i)/alpha ;
656 sec_membre.set(ligne)
657 -= (2*i+pari)*(2*i+pari)* sol_part_eta(0, k, j, i)/alpha ;
661 for (
int i=0; i<nr; i++) {
662 systeme.set(ligne, colonne)
663 += (2*i+pari)*(2*i+pari)*sol_hom_un_vr(0, k, j, i)/alpha ;
664 systeme.set(ligne, colonne+1)
665 += (2*i+pari)*(2*i+pari)*sol_hom_trois_vr(0, k, j, i)/alpha ;
666 sec_membre.set(ligne)
667 -= (2*i+pari)*(2*i+pari)* sol_part_vr(0, k, j, i)/alpha ;
672 for (
int zone=1 ; zone<nzm1 ; zone++) {
673 nr = mg.get_nr(zone) ;
674 alpha = mpaff->get_alpha()[zone] ;
676 systeme.set(ligne, colonne)
677 = -sol_hom_un_eta.val_in_bound_jk(zone, j, k) ;
678 systeme.set(ligne, colonne+1)
679 = -sol_hom_deux_eta.val_in_bound_jk(zone, j, k) ;
680 systeme.set(ligne, colonne+2)
681 = -sol_hom_trois_eta.val_in_bound_jk(zone, j, k) ;
682 systeme.set(ligne, colonne+3)
683 = -sol_hom_quatre_eta.val_in_bound_jk(zone, j, k) ;
684 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
687 systeme.set(ligne, colonne)
688 = -sol_hom_un_vr.val_in_bound_jk(zone, j, k) ;
689 systeme.set(ligne, colonne+1)
690 = -sol_hom_deux_vr.val_in_bound_jk(zone, j, k) ;
691 systeme.set(ligne, colonne+2)
692 = -sol_hom_trois_vr.val_in_bound_jk(zone, j, k) ;
693 systeme.set(ligne, colonne+3)
694 = -sol_hom_quatre_vr.val_in_bound_jk(zone, j, k) ;
695 sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
700 for (
int i=0; i<nr; i++) {
701 systeme.set(ligne, colonne)
702 -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
703 systeme.set(ligne, colonne+1)
704 -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
705 systeme.set(ligne, colonne+2)
706 -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
707 systeme.set(ligne, colonne+3)
708 -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
709 sec_membre.set(ligne)
710 += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
716 for (
int i=0; i<nr; i++) {
717 systeme.set(ligne, colonne)
718 -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
719 systeme.set(ligne, colonne+1)
720 -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
721 systeme.set(ligne, colonne+2)
722 -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
723 systeme.set(ligne, colonne+3)
724 -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
725 sec_membre.set(ligne)
726 += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
731 systeme.set(ligne, colonne)
732 += sol_hom_un_eta.val_out_bound_jk(zone, j, k) ;
733 systeme.set(ligne, colonne+1)
734 += sol_hom_deux_eta.val_out_bound_jk(zone, j, k) ;
735 systeme.set(ligne, colonne+2)
736 += sol_hom_trois_eta.val_out_bound_jk(zone, j, k) ;
737 systeme.set(ligne, colonne+3)
738 += sol_hom_quatre_eta.val_out_bound_jk(zone, j, k) ;
739 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
742 systeme.set(ligne, colonne)
743 += sol_hom_un_vr.val_out_bound_jk(zone, j, k) ;
744 systeme.set(ligne, colonne+1)
745 += sol_hom_deux_vr.val_out_bound_jk(zone, j, k) ;
746 systeme.set(ligne, colonne+2)
747 += sol_hom_trois_vr.val_out_bound_jk(zone, j, k) ;
748 systeme.set(ligne, colonne+3)
749 += sol_hom_quatre_vr.val_out_bound_jk(zone, j, k) ;
750 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
755 for (
int i=0; i<nr; i++) {
756 systeme.set(ligne, colonne)
757 += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
758 systeme.set(ligne, colonne+1)
759 += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
760 systeme.set(ligne, colonne+2)
761 += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
762 systeme.set(ligne, colonne+3)
763 += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
764 sec_membre.set(ligne)
765 -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
771 for (
int i=0; i<nr; i++) {
772 systeme.set(ligne, colonne)
773 += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
774 systeme.set(ligne, colonne+1)
775 += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
776 systeme.set(ligne, colonne+2)
777 += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
778 systeme.set(ligne, colonne+3)
779 += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
780 sec_membre.set(ligne)
781 -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
787 nr = mg.get_nr(nzm1) ;
789 alpha = mpaff->get_alpha()[nzm1] ;
792 systeme.set(ligne, colonne)
793 -= sol_hom_deux_eta.val_in_bound_jk(nzm1, j, k) ;
794 systeme.set(ligne, colonne+1)
795 -= sol_hom_quatre_eta.val_in_bound_jk(nzm1, j, k) ;
796 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nzm1, j, k) ;
798 systeme.set(ligne+1, colonne)
799 -= sol_hom_deux_vr.val_in_bound_jk(nzm1, j, k) ;
800 systeme.set(ligne+1, colonne+1)
801 -= sol_hom_quatre_vr.val_in_bound_jk(nzm1, j, k) ;
802 sec_membre.set(ligne+1) += sol_part_vr.val_in_bound_jk(nzm1, j, k) ;
807 for (
int i=0; i<nr; i++) {
808 systeme.set(ligne, colonne)
809 -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
810 systeme.set(ligne, colonne+1)
811 -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
812 sec_membre.set(ligne)
813 += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
819 for (
int i=0; i<nr; i++) {
820 systeme.set(ligne, colonne)
821 -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
822 systeme.set(ligne, colonne+1)
823 -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
824 sec_membre.set(ligne)
825 += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
833 Tbl facteurs(systeme.inverse(sec_membre)) ;
840 for (
int i=0 ; i<nr ; i++) {
841 cf_eta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
842 +facteurs(conte)*sol_hom_un_eta(0, k, j, i)
843 +facteurs(conte+1)*sol_hom_trois_eta(0, k, j, i) ;
844 cf_vr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
845 +facteurs(conte)*sol_hom_un_vr(0, k, j, i)
846 +facteurs(conte+1)*sol_hom_trois_vr(0, k, j, i) ;
849 for (
int zone=1 ; zone<nzm1 ; zone++) {
850 nr = mg.get_nr(zone) ;
851 for (
int i=0 ; i<nr ; i++) {
852 cf_eta.set(zone, k, j, i) =
853 sol_part_eta(zone, k, j, i)
854 +facteurs(conte)*sol_hom_un_eta(zone, k, j, i)
855 +facteurs(conte+1)*sol_hom_deux_eta(zone, k, j, i)
856 +facteurs(conte+2)*sol_hom_trois_eta(zone, k, j, i)
857 +facteurs(conte+3)*sol_hom_quatre_eta(zone, k, j, i) ;
858 cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
859 +facteurs(conte)*sol_hom_un_vr(zone, k, j, i)
860 +facteurs(conte+1)*sol_hom_deux_vr(zone, k, j, i)
861 +facteurs(conte+2)*sol_hom_trois_vr(zone, k, j, i)
862 +facteurs(conte+3)*sol_hom_quatre_vr(zone, k, j, i) ;
866 nr = mg.get_nr(nz-1) ;
867 for (
int i=0 ; i<nr ; i++) {
868 cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
869 +facteurs(conte)*sol_hom_deux_eta(nzm1, k, j, i)
870 +facteurs(conte+1)*sol_hom_quatre_eta(nzm1, k, j, i) ;
871 cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
872 +facteurs(conte)*sol_hom_deux_vr(nzm1, k, j, i)
873 +facteurs(conte+1)*sol_hom_quatre_vr(nzm1, k, j, i) ;
878 vr.set_spectral_va().ylm_i() ;
880 het.set_spectral_va().ylm_i() ;
882 resu.set_vr_eta_mu(vr, het,
mu().
poisson()) ;
883 if ((nt==1)&&(np==1)) {
884 resu.set(2).set_etat_zero() ;
885 resu.set(3).set_etat_zero() ;
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
Scalar ** cmp
Array of size n_comp of pointers onto the components.
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...