28char vector_poisson_block_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_block.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $" ;
75void Vector::poisson_block(
double lam, Vector& resu)
const {
77 const Map_af* mpaff =
dynamic_cast<const Map_af*
>(
mp) ;
79 for (
int i=0; i<3; i++)
80 assert(
cmp[i]->check_dzpuis(4)) ;
82 const Base_vect_spher* bvs =
dynamic_cast<const Base_vect_spher*
>(
triad) ;
85 assert( mpaff != 0x0) ;
88 if (fabs(lam + 1.) < 1.e-8) {
89 const Metric_flat& mets =
mp->flat_met_spher() ;
90 Vector_divfree sou(*
mp, *
triad, mets) ;
91 for (
int i=1; i<=3; i++) sou.set(i) = *
cmp[i-1] ;
98 const Mg3d& mg = *mpaff->get_mg() ;
99 int nz = mg.get_nzone() ;
int nzm1 = nz - 1;
100 int nt = mg.get_nt(0) ;
101 int np = mg.get_np(0) ;
102 assert(mg.get_type_r(0) == RARE) ;
103 assert(mg.get_type_r(nzm1) == UNSURR) ;
104 Scalar S_r = *
cmp[0] ;
105 Scalar S_eta =
eta() ;
108 bool all_zero = false ;
109 if (S_r.get_etat() == ETATZERO) {
110 if (S_eta.get_etat() == ETATZERO) {
112 het.set_etat_zero() ;
117 S_r.set_spectral_base(S_eta.get_spectral_base()) ;
120 if ((S_eta.get_etat() == ETATZERO)&&(!all_zero)) {
121 S_eta.annule_hard() ;
122 S_eta.set_spectral_base(S_r.get_spectral_base()) ;
125 S_r.set_spectral_va().ylm() ;
126 S_eta.set_spectral_va().ylm() ;
127 const Base_val& base = S_eta.get_spectral_va().base ;
128 Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
129 Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
130 Mtbl_cf sol_hom_un_eta(mg, base) ; sol_hom_un_eta.annule_hard() ;
131 Mtbl_cf sol_hom_deux_eta(mg, base) ; sol_hom_deux_eta.annule_hard() ;
132 Mtbl_cf sol_hom_trois_eta(mg, base) ; sol_hom_trois_eta.annule_hard() ;
133 Mtbl_cf sol_hom_quatre_eta(mg, base) ; sol_hom_quatre_eta.annule_hard() ;
134 Mtbl_cf sol_hom_un_vr(mg, base) ; sol_hom_un_vr.annule_hard() ;
135 Mtbl_cf sol_hom_deux_vr(mg, base) ; sol_hom_deux_vr.annule_hard() ;
136 Mtbl_cf sol_hom_trois_vr(mg, base) ; sol_hom_trois_vr.annule_hard() ;
137 Mtbl_cf sol_hom_quatre_vr(mg, base) ; sol_hom_quatre_vr.annule_hard() ;
141 Scalar sou_l0 = (*
cmp[0]) / (1. + lam) ;
142 sou_l0.set_spectral_va().ylm() ;
143 Scalar vrl0 = pois_vect_r0(sou_l0) ;
150 int nr = mg.get_nr(0) ;
151 double alpha = mpaff->get_alpha()[0] ;
double alp2 = alpha*alpha ;
152 double beta = mpaff->get_beta()[0] ;
153 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
157 for (
int k=0 ; k<np+1 ; k++) {
158 for (
int j=0 ; j<nt ; j++) {
159 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
160 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0) ) {
161 int aa = 0 ;
int bb = 0 ;
int nr0 = 0 ;
173 Matrice oper(2*nr, 2*nr) ; oper.set_etat_qcq() ;
174 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
175 Diff_dsdx2 d2(base_r, nr) ;
const Matrice& md2 = d2.get_matrice() ;
176 Diff_sxdsdx xd(base_r, nr) ;
const Matrice& mxd = xd.get_matrice() ;
177 Diff_sx2 s2(base_r, nr) ;
const Matrice& ms2 = s2.get_matrice() ;
181 for (
int lin=0; lin<nr; lin++) {
182 for (
int col=0; col<nr; col++)
184 (md2(lin,col) + 2*mxd(lin,col)
185 -(lam+1)*l_q*(l_q+1)*ms2(lin,col)) / alp2 ;
186 for (
int col=0; col<nr; col++)
187 oper.set(lin,col+nr) =
188 (lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
190 for (
int lin=0; lin<nr; lin++) {
191 for (
int col=0; col<nr; col++)
192 oper.set(lin+nr,col) =
193 (-lam*l_q*(l_q+1)*mxd(lin,col)
194 +(lam+2)*l_q*(l_q+1)*ms2(lin,col)) / alp2 ;
195 for (
int col=0; col<nr; col++)
196 oper.set(lin+nr, col+nr) =
197 ((lam+1)*(md2(lin,col) + 2*mxd(lin,col))
198 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
200 bool pb_eta = ( ( fabs( lam*
double(l_q+3) + 2 ) < 0.01) && (l_q <=2) ) ;
202 for (
int col=0; col<nr; col++)
203 oper.set(nr0-1, col) = 1 ;
204 for (
int col=0; col<nr; col++)
205 oper.set(nr0-1,nr+col) = 0 ;
207 if ((l_q > 2)||pb_eta) {
208 for (
int col=0; col<nr; col++)
209 oper.set(nr+nr0-1, col) = 0 ;
210 for (
int col=0; col<nr; col++)
211 oper.set(nr+nr0-1,nr+col) = 1 ;
214 Matrice op2(nrtot, nrtot) ;
216 for (
int i=0; i<nr0; i++) {
217 for (
int col=0; col<nr0;col++)
218 op2.set(i,col) = (aa*col+bb)*oper(i,col+1)
219 + (aa*(col+1)+bb)*oper(i,col) ;
220 for (
int col=0; col<nr0;col++)
221 op2.set(i,col+nr0) = (aa*col+bb)*oper(i,col+nr+1)
222 + (aa*(col+1)+bb)*oper(i,col+nr) ;
225 for (
int i=nr0; i<nrtot; i++) {
226 for (
int col=0; col<nr0;col++)
227 op2.set(i,col) = (aa*col+bb)*oper(i+d0,col+1)
228 + (aa*(col+1)+bb)*oper(i+d0,col) ;
229 for (
int col=0; col<nr0;col++)
230 op2.set(i,col+nr0) = (aa*col+bb)*oper(i+d0,col+nr+1)
231 + (aa*(col+1)+bb)*oper(i+d0,col+nr) ;
237 for (
int i=0; i<nr0; i++)
238 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(0, k, j, i) ;
239 if (!pb_eta) sec_membre.set(nr0-1) = 0 ;
240 for (
int i=0; i<nr0; i++)
241 sec_membre.set(i+nr0)
242 = (*S_r.get_spectral_va().c_cf)(0, k, j, i) ;
243 if ((l_q > 2)||pb_eta) sec_membre.set(nrtot-1) = 0 ;
247 Tbl big_res = op2.inverse(sec_membre) ;
251 sol_part_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
252 for (
int i=1; i<nr0; i++)
253 sol_part_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
254 + (aa*(i+1)+bb)*big_res(i);
255 sol_part_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1) + bb) ;
256 sol_part_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
257 for (
int i=1; i<nr0; i++)
258 sol_part_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
259 + (aa*(i+1)+bb)*big_res(nr0+i);
260 sol_part_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
262 sol_part_eta.set(0, k, j, nr-1) = 0 ;
263 sol_part_vr.set(0, k, j, nr-1) = 0 ;
268 double fac_eta = lam*double(l_q+3) + 2. ;
269 double fac_vr = double(l_q+1)*(lam*l_q - 2.) ;
270 Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
271 Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
272 for (
int i=0 ; i<nr ; i++) {
273 sol_hom_un_eta.set(0, k, j, i) = sol_hom1(i) ;
274 sol_hom_un_vr.set(0, k, j, i) = l_q*sol_hom1(i) ;
275 sol_hom_trois_eta.set(0, k, j, i) = fac_eta*sol_hom2(i) ;
276 sol_hom_trois_vr.set(0, k, j, i) = fac_vr*sol_hom2(i) ;
280 for (
int i=0; i<nrtot; i++) sec_membre.set(i) = 0 ;
282 sec_membre.set(nr0-1) = 1 ;
283 big_res = op2.inverse(sec_membre) ;
285 sol_hom_un_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
286 for (
int i=1; i<nr0; i++)
287 sol_hom_un_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
288 + (aa*(i+1)+bb)*big_res(i);
289 sol_hom_un_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1) + bb) ;
290 sol_hom_un_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
291 for (
int i=1; i<nr0; i++)
292 sol_hom_un_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
293 + (aa*(i+1)+bb)*big_res(nr0+i);
294 sol_hom_un_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
296 sol_hom_un_eta.set(0, k, j, nr-1) = 0 ;
297 sol_hom_un_vr.set(0, k, j, nr-1) = 0 ;
300 sec_membre.set(nr0-1) = 0 ;
301 sec_membre.set(nrtot-1) = 1 ;
302 big_res = op2.inverse(sec_membre) ;
304 sol_hom_trois_eta.set(0, k, j, 0) = (aa+bb)*big_res(0) ;
305 for (
int i=1; i<nr0; i++)
306 sol_hom_trois_eta.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(i-1)
307 + (aa*(i+1)+bb)*big_res(i);
308 sol_hom_trois_eta.set(0, k, j, nr0) = big_res(nr0-1)*(aa*(nr0-1)+bb) ;
309 sol_hom_trois_vr.set(0, k, j, 0) = (aa+bb)*big_res(nr0) ;
310 for (
int i=1; i<nr0; i++)
311 sol_hom_trois_vr.set(0, k, j, i) = (aa*(i-1)+bb)*big_res(nr0+i-1)
312 + (aa*(i+1)+bb)*big_res(nr0+i);
313 sol_hom_trois_vr.set(0, k, j, nr0) = big_res(nrtot-1)*(aa*(nr0-1)+bb) ;
315 sol_hom_trois_eta.set(0, k, j, nr-1) = 0 ;
316 sol_hom_trois_vr.set(0, k, j, nr-1) = 0 ;
327 for (
int zone=1 ; zone<nzm1 ; zone++) {
328 nr = mg.get_nr(zone) ;
330 assert(nt == mg.get_nt(zone)) ;
331 assert(np == mg.get_np(zone)) ;
332 alpha = mpaff->get_alpha()[zone] ;
333 beta = mpaff->get_beta()[zone] ;
334 double ech = beta / alpha ;
338 for (
int k=0 ; k<np+1 ; k++) {
339 for (
int j=0 ; j<nt ; j++) {
340 base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
341 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
344 int nr_eq1 = nr - dege1 ;
345 int nr_eq2 = nr - dege2 ;
347 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
348 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
349 Diff_x2dsdx2 x2d2(base_r, nr);
const Matrice& m2d2 = x2d2.get_matrice() ;
350 Diff_xdsdx2 xd2(base_r, nr) ;
const Matrice& mxd2 = xd2.get_matrice() ;
351 Diff_dsdx2 d2(base_r, nr) ;
const Matrice& md2 = d2.get_matrice() ;
352 Diff_xdsdx xd(base_r, nr) ;
const Matrice& mxd = xd.get_matrice() ;
353 Diff_dsdx d1(base_r, nr) ;
const Matrice& md = d1.get_matrice() ;
354 Diff_id id(base_r, nr) ;
const Matrice& mid =
id.get_matrice() ;
358 for (
int lin=0; lin<nr_eq1; lin++) {
359 for (
int col=0; col<nr; col++)
361 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
362 + 2*(mxd(lin,col) + ech*md(lin,col))
363 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
364 for (
int col=0; col<nr; col++)
366 = lam*(mxd(lin,col) + ech*md(lin,col))
367 + 2*(1+lam)*mid(lin,col) ;
369 oper.set(nr_eq1, 0) = 1 ;
370 for (
int col=1; col<2*nr; col++)
371 oper.set(nr_eq1, col) = 0 ;
372 oper.set(nr_eq1+1, 0) = 0 ;
373 oper.set(nr_eq1+1, 1) = 1 ;
374 for (
int col=2; col<2*nr; col++)
375 oper.set(nr_eq1+1, col) = 0 ;
376 for (
int lin=0; lin<nr_eq2; lin++) {
377 for (
int col=0; col<nr; col++)
379 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
380 - (lam+2)*mid(lin,col)) ;
381 for (
int col=0; col<nr; col++)
382 oper.set(lin+nr, col+nr)
383 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
384 + ech*ech*md2(lin,col)
385 + 2*(mxd(lin,col) + ech*md(lin,col)))
386 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
388 for (
int col=0; col<nr; col++)
389 oper.set(nr+nr_eq2, col) = 0 ;
390 oper.set(nr+nr_eq2, nr) = 1 ;
391 for (
int col=nr+1; col<2*nr; col++)
392 oper.set(nr+nr_eq2, col) = 0 ;
393 for (
int col=0; col<nr+1; col++)
394 oper.set(nr+nr_eq2+1, col) = 0 ;
395 oper.set(nr+nr_eq2+1, nr+1) = 1 ;
396 for (
int col=nr+2; col<2*nr; col++)
397 oper.set(nr+nr_eq2+1, col) = 0 ;
403 Tbl sr(nr) ; sr.set_etat_qcq() ;
404 Tbl seta(nr) ; seta.set_etat_qcq() ;
405 for (
int i=0; i<nr; i++) {
406 sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
407 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
409 Tbl xsr= sr ; Tbl x2sr= sr ;
410 Tbl xseta= seta ; Tbl x2seta = seta ;
411 multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
412 multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
414 for (
int i=0; i<nr_eq1; i++)
415 sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
417 sec_membre.set(nr_eq1) = 0 ;
418 sec_membre.set(nr_eq1+1) = 0 ;
419 for (
int i=0; i<nr_eq2; i++)
420 sec_membre.set(i+nr) = beta*beta*sr(i)
421 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
422 sec_membre.set(nr+nr_eq2) = 0 ;
423 sec_membre.set(nr+nr_eq2+1) = 0 ;
427 Tbl big_res = oper.inverse(sec_membre) ;
428 for (
int i=0; i<nr; i++) {
429 sol_part_eta.set(zone, k, j, i) = big_res(i) ;
430 sol_part_vr.set(zone, k, j, i) = big_res(nr+i) ;
435 sec_membre.annule_hard() ;
436 sec_membre.set(nr_eq1) = 1 ;
437 big_res = oper.inverse(sec_membre) ;
438 for (
int i=0 ; i<nr ; i++) {
439 sol_hom_un_eta.set(zone, k, j, i) = big_res(i) ;
440 sol_hom_un_vr.set(zone, k, j, i) = big_res(nr+i) ;
442 sec_membre.set(nr_eq1) = 0 ;
443 sec_membre.set(nr_eq1+1) = 1 ;
444 big_res = oper.inverse(sec_membre) ;
445 for (
int i=0 ; i<nr ; i++) {
446 sol_hom_deux_eta.set(zone, k, j, i) = big_res(i) ;
447 sol_hom_deux_vr.set(zone, k, j, i) = big_res(nr+i) ;
449 sec_membre.set(nr_eq1+1) = 0 ;
450 sec_membre.set(nr+nr_eq2) = 1 ;
451 big_res = oper.inverse(sec_membre) ;
452 for (
int i=0 ; i<nr ; i++) {
453 sol_hom_trois_eta.set(zone, k, j, i) = big_res(i) ;
454 sol_hom_trois_vr.set(zone, k, j, i) = big_res(nr+i) ;
456 sec_membre.set(nr+nr_eq2) = 0 ;
457 sec_membre.set(nr+nr_eq2+1) = 1 ;
458 big_res = oper.inverse(sec_membre) ;
459 for (
int i=0 ; i<nr ; i++) {
460 sol_hom_quatre_eta.set(zone, k, j, i) = big_res(i) ;
461 sol_hom_quatre_vr.set(zone, k, j, i) = big_res(nr+i) ;
471 nr = mg.get_nr(nzm1) ;
472 assert(nt == mg.get_nt(nzm1)) ;
473 assert(np == mg.get_np(nzm1)) ;
474 alpha = mpaff->get_alpha()[nzm1] ; alp2 = alpha*alpha ;
479 for (
int k=0 ; k<np+1 ; k++) {
480 for (
int j=0 ; j<nt ; j++) {
481 base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
482 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
484 int dege2 = (l_q == 1 ? 2 : 3);
485 int nr_eq1 = nr - dege1 ;
486 int nr_eq2 = nr - dege2 ;
488 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
489 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
490 Diff_dsdx2 d2(base_r, nr) ;
const Matrice& md2 = d2.get_matrice() ;
491 Diff_sxdsdx sxd(base_r, nr) ;
const Matrice& mxd = sxd.get_matrice() ;
492 Diff_sx2 sx2(base_r, nr) ;
const Matrice& ms2 = sx2.get_matrice() ;
496 for (
int lin=0; lin<nr_eq1; lin++) {
497 for (
int col=0; col<nr; col++)
499 = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
500 for (
int col=0; col<nr; col++)
501 oper.set(lin,col+nr) =
502 (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
504 for (
int col=0; col<nr; col++)
505 oper.set(nr_eq1, col) = 1 ;
506 for (
int col=nr; col<nrtot; col++)
507 oper.set(nr_eq1, col) = 0 ;
509 for (
int col=0; col<nr; col++) {
510 oper.set(nr_eq1+1, col) = pari*col*col ;
513 for (
int col=nr; col<nrtot; col++)
514 oper.set(nr_eq1+1, col) = 0 ;
515 oper.set(nr_eq1+2, 0) = 1 ;
516 for (
int col=1; col<nrtot; col++)
517 oper.set(nr_eq1+2, col) = 0 ;
518 for (
int lin=0; lin<nr_eq2; lin++) {
519 for (
int col=0; col<nr; col++)
520 oper.set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
521 + (lam+2)*ms2(lin,col))) / alp2 ;
522 for (
int col=0; col<nr; col++)
523 oper.set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
524 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
526 for (
int col=0; col<nr; col++)
527 oper.set(nr+nr_eq2, col) = 0 ;
528 for (
int col=nr; col<nrtot; col++)
529 oper.set(nr+nr_eq2, col) = 1 ;
530 for (
int col=0; col<nr; col++)
531 oper.set(nr+nr_eq2+1, col) = 0 ;
533 for (
int col=0; col<nr; col++) {
534 oper.set(nr+nr_eq2+1, nr+col) = pari*col*col ;
538 for (
int col=0; col<nr; col++)
539 oper.set(nr+nr_eq2+2, col) = 0 ;
540 oper.set(nr+nr_eq2+2, nr) = 1 ;
541 for (
int col=nr+1; col<nrtot; col++)
542 oper.set(nr+nr_eq2+2, col) = 0 ;
548 for (
int i=0; i<nr_eq1; i++)
549 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
550 for (
int i=nr_eq1; i<nr; i++)
551 sec_membre.set(i) = 0 ;
552 for (
int i=0; i<nr_eq2; i++)
553 sec_membre.set(i+nr) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
554 for (
int i=nr_eq2; i<nr; i++)
555 sec_membre.set(nr+i) = 0 ;
556 Tbl big_res = oper.inverse(sec_membre) ;
560 for (
int i=0; i<nr; i++) {
561 sol_part_eta.set(nzm1, k, j, i) = big_res(i) ;
562 sol_part_vr.set(nzm1, k, j, i) = big_res(i+nr) ;
568 Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
569 Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
570 double fac_eta = lam + 2. ;
571 double fac_vr = 2*lam + 2. ;
572 for (
int i=0 ; i<nr ; i++) {
573 sol_hom_deux_eta.set(nzm1, k, j, i) = sol_hom2(i) ;
574 sol_hom_quatre_eta.set(nzm1, k, j, i) = fac_eta*sol_hom1(i) ;
575 sol_hom_deux_vr.set(nzm1, k, j, i) = -2*sol_hom2(i) ;
576 sol_hom_quatre_vr.set(nzm1, k, j, i) = fac_vr*sol_hom1(i) ;
580 sec_membre.annule_hard() ;
581 sec_membre.set(nr-1) = 1 ;
582 big_res = oper.inverse(sec_membre) ;
584 for (
int i=0; i<nr; i++) {
585 sol_hom_deux_eta.set(nzm1, k, j, i) = big_res(i) ;
586 sol_hom_deux_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
588 sec_membre.set(nr-1) = 0 ;
589 sec_membre.set(2*nr-1) = 1 ;
590 big_res = oper.inverse(sec_membre) ;
591 for (
int i=0; i<nr; i++) {
592 sol_hom_quatre_eta.set(nzm1, k, j, i) = big_res(i) ;
593 sol_hom_quatre_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
605 vr.set_spectral_base(base) ;
606 vr.set_spectral_va().set_etat_cf_qcq() ;
607 Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
608 cf_vr.annule_hard() ;
610 het.set_spectral_base(base) ;
611 het.set_spectral_va().set_etat_cf_qcq() ;
612 Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
613 cf_eta.annule_hard() ;
614 int taille = 4*nzm1 ;
615 Tbl sec_membre(taille) ;
616 Matrice systeme(taille, taille) ;
617 systeme.set_etat_qcq() ;
618 int ligne ;
int colonne ;
622 for (
int k=0 ; k<np+1 ; k++)
623 for (
int j=0 ; j<nt ; j++) {
624 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
625 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
629 systeme.annule_hard() ;
630 sec_membre.annule_hard() ;
634 alpha = mpaff->get_alpha()[0] ;
636 systeme.set(ligne, colonne) = sol_hom_un_eta.val_out_bound_jk(0, j, k) ;
637 systeme.set(ligne, colonne+1) = sol_hom_trois_eta.val_out_bound_jk(0, j, k) ;
638 sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
641 systeme.set(ligne, colonne) = sol_hom_un_vr.val_out_bound_jk(0, j, k) ;
642 systeme.set(ligne, colonne+1) = sol_hom_trois_vr.val_out_bound_jk(0, j, k) ;
643 sec_membre.set(ligne) = -sol_part_vr.val_out_bound_jk(0,j,k) ;
647 int pari = (base_r ==
R_CHEBP ? 0 : 1) ;
648 for (
int i=0; i<nr; i++) {
649 systeme.set(ligne, colonne)
650 += (2*i+pari)*(2*i+pari)*sol_hom_un_eta(0, k, j, i)/alpha ;
651 systeme.set(ligne, colonne+1)
652 += (2*i+pari)*(2*i+pari)*sol_hom_trois_eta(0, k, j, i)/alpha ;
653 sec_membre.set(ligne)
654 -= (2*i+pari)*(2*i+pari)* sol_part_eta(0, k, j, i)/alpha ;
658 for (
int i=0; i<nr; i++) {
659 systeme.set(ligne, colonne)
660 += (2*i+pari)*(2*i+pari)*sol_hom_un_vr(0, k, j, i)/alpha ;
661 systeme.set(ligne, colonne+1)
662 += (2*i+pari)*(2*i+pari)*sol_hom_trois_vr(0, k, j, i)/alpha ;
663 sec_membre.set(ligne)
664 -= (2*i+pari)*(2*i+pari)* sol_part_vr(0, k, j, i)/alpha ;
669 for (
int zone=1 ; zone<nzm1 ; zone++) {
670 nr = mg.get_nr(zone) ;
671 alpha = mpaff->get_alpha()[zone] ;
673 systeme.set(ligne, colonne)
674 = -sol_hom_un_eta.val_in_bound_jk(zone, j, k) ;
675 systeme.set(ligne, colonne+1)
676 = -sol_hom_deux_eta.val_in_bound_jk(zone, j, k) ;
677 systeme.set(ligne, colonne+2)
678 = -sol_hom_trois_eta.val_in_bound_jk(zone, j, k) ;
679 systeme.set(ligne, colonne+3)
680 = -sol_hom_quatre_eta.val_in_bound_jk(zone, j, k) ;
681 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
684 systeme.set(ligne, colonne)
685 = -sol_hom_un_vr.val_in_bound_jk(zone, j, k) ;
686 systeme.set(ligne, colonne+1)
687 = -sol_hom_deux_vr.val_in_bound_jk(zone, j, k) ;
688 systeme.set(ligne, colonne+2)
689 = -sol_hom_trois_vr.val_in_bound_jk(zone, j, k) ;
690 systeme.set(ligne, colonne+3)
691 = -sol_hom_quatre_vr.val_in_bound_jk(zone, j, k) ;
692 sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
697 for (
int i=0; i<nr; i++) {
698 systeme.set(ligne, colonne)
699 -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
700 systeme.set(ligne, colonne+1)
701 -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
702 systeme.set(ligne, colonne+2)
703 -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
704 systeme.set(ligne, colonne+3)
705 -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
706 sec_membre.set(ligne)
707 += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
713 for (
int i=0; i<nr; i++) {
714 systeme.set(ligne, colonne)
715 -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
716 systeme.set(ligne, colonne+1)
717 -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
718 systeme.set(ligne, colonne+2)
719 -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
720 systeme.set(ligne, colonne+3)
721 -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
722 sec_membre.set(ligne)
723 += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
728 systeme.set(ligne, colonne)
729 += sol_hom_un_eta.val_out_bound_jk(zone, j, k) ;
730 systeme.set(ligne, colonne+1)
731 += sol_hom_deux_eta.val_out_bound_jk(zone, j, k) ;
732 systeme.set(ligne, colonne+2)
733 += sol_hom_trois_eta.val_out_bound_jk(zone, j, k) ;
734 systeme.set(ligne, colonne+3)
735 += sol_hom_quatre_eta.val_out_bound_jk(zone, j, k) ;
736 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
739 systeme.set(ligne, colonne)
740 += sol_hom_un_vr.val_out_bound_jk(zone, j, k) ;
741 systeme.set(ligne, colonne+1)
742 += sol_hom_deux_vr.val_out_bound_jk(zone, j, k) ;
743 systeme.set(ligne, colonne+2)
744 += sol_hom_trois_vr.val_out_bound_jk(zone, j, k) ;
745 systeme.set(ligne, colonne+3)
746 += sol_hom_quatre_vr.val_out_bound_jk(zone, j, k) ;
747 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
752 for (
int i=0; i<nr; i++) {
753 systeme.set(ligne, colonne)
754 += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
755 systeme.set(ligne, colonne+1)
756 += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
757 systeme.set(ligne, colonne+2)
758 += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
759 systeme.set(ligne, colonne+3)
760 += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
761 sec_membre.set(ligne)
762 -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
768 for (
int i=0; i<nr; i++) {
769 systeme.set(ligne, colonne)
770 += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
771 systeme.set(ligne, colonne+1)
772 += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
773 systeme.set(ligne, colonne+2)
774 += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
775 systeme.set(ligne, colonne+3)
776 += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
777 sec_membre.set(ligne)
778 -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
784 nr = mg.get_nr(nzm1) ;
786 alpha = mpaff->get_alpha()[nzm1] ;
789 systeme.set(ligne, colonne)
790 -= sol_hom_deux_eta.val_in_bound_jk(nzm1, j, k) ;
791 systeme.set(ligne, colonne+1)
792 -= sol_hom_quatre_eta.val_in_bound_jk(nzm1, j, k) ;
793 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nzm1, j, k) ;
795 systeme.set(ligne+1, colonne)
796 -= sol_hom_deux_vr.val_in_bound_jk(nzm1, j, k) ;
797 systeme.set(ligne+1, colonne+1)
798 -= sol_hom_quatre_vr.val_in_bound_jk(nzm1, j, k) ;
799 sec_membre.set(ligne+1) += sol_part_vr.val_in_bound_jk(nzm1, j, k) ;
804 for (
int i=0; i<nr; i++) {
805 systeme.set(ligne, colonne)
806 -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
807 systeme.set(ligne, colonne+1)
808 -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
809 sec_membre.set(ligne)
810 += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
816 for (
int i=0; i<nr; i++) {
817 systeme.set(ligne, colonne)
818 -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
819 systeme.set(ligne, colonne+1)
820 -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
821 sec_membre.set(ligne)
822 += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
830 Tbl facteurs(systeme.inverse(sec_membre)) ;
837 for (
int i=0 ; i<nr ; i++) {
838 cf_eta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
839 +facteurs(conte)*sol_hom_un_eta(0, k, j, i)
840 +facteurs(conte+1)*sol_hom_trois_eta(0, k, j, i) ;
841 cf_vr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
842 +facteurs(conte)*sol_hom_un_vr(0, k, j, i)
843 +facteurs(conte+1)*sol_hom_trois_vr(0, k, j, i) ;
846 for (
int zone=1 ; zone<nzm1 ; zone++) {
847 nr = mg.get_nr(zone) ;
848 for (
int i=0 ; i<nr ; i++) {
849 cf_eta.set(zone, k, j, i) =
850 sol_part_eta(zone, k, j, i)
851 +facteurs(conte)*sol_hom_un_eta(zone, k, j, i)
852 +facteurs(conte+1)*sol_hom_deux_eta(zone, k, j, i)
853 +facteurs(conte+2)*sol_hom_trois_eta(zone, k, j, i)
854 +facteurs(conte+3)*sol_hom_quatre_eta(zone, k, j, i) ;
855 cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
856 +facteurs(conte)*sol_hom_un_vr(zone, k, j, i)
857 +facteurs(conte+1)*sol_hom_deux_vr(zone, k, j, i)
858 +facteurs(conte+2)*sol_hom_trois_vr(zone, k, j, i)
859 +facteurs(conte+3)*sol_hom_quatre_vr(zone, k, j, i) ;
863 nr = mg.get_nr(nz-1) ;
864 for (
int i=0 ; i<nr ; i++) {
865 cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
866 +facteurs(conte)*sol_hom_deux_eta(nzm1, k, j, i)
867 +facteurs(conte+1)*sol_hom_quatre_eta(nzm1, k, j, i) ;
868 cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
869 +facteurs(conte)*sol_hom_deux_vr(nzm1, k, j, i)
870 +facteurs(conte+1)*sol_hom_quatre_vr(nzm1, k, j, i) ;
875 vr.set_spectral_va().ylm_i() ;
877 het.set_spectral_va().ylm_i() ;
879 resu.set_vr_eta_mu(vr, het,
mu().
poisson()) ;
880 if ((nt==1)&&(np==1)) {
881 resu.set(2).set_etat_zero() ;
882 resu.set(3).set_etat_zero() ;
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
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.