63 #include "param_elliptic.h" 65 #include "utilitaires.h" 70 int num_front,
double fact_dir,
double fact_neu,
75 for (
int i=0; i<3; i++)
76 assert(
cmp[i]->check_dzpuis(4)) ;
81 assert( mpaff != 0x0) ;
84 if (fabs(lam + 1.) < 1.e-8) {
85 cout <<
"Not implemented yet !!" << endl ;
99 int nz = mg.
get_nzone() ;
int nzm1 = nz - 1;
102 if (S_r.get_etat() == ETATZERO) S_r.
annule_hard() ;
105 S_r.set_spectral_va().ylm() ;
120 for (
int l=0; l<nz; l++)
135 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
136 double alpha, beta, ech ;
138 assert(num_front+1 < nzm1) ;
139 for (
int zone=num_front+1 ; zone<nzm1 ; zone++) {
145 assert(nt == mg.
get_nt(zone)) ;
146 assert(np == mg.
get_np(zone)) ;
150 for (
int k=0 ; k<np+1 ; k++) {
151 for (
int j=0 ; j<nt ; j++) {
153 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
156 int nr_eq1 = nr - dege1 ;
157 int nr_eq2 = nr - dege2 ;
158 int nrtot = nr_eq1 + nr_eq2 ;
170 for (
int lin=0; lin<nr_eq1; lin++) {
171 for (
int col=dege1; col<nr; col++)
172 oper.
set(lin,col-dege1)
173 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
174 + 2*(mxd(lin,col) + ech*md(lin,col))
175 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
176 for (
int col=dege2; col<nr; col++)
177 oper.
set(lin,col-dege2+nr_eq1)
178 = lam*(mxd(lin,col) + ech*md(lin,col)) + 2*(1+lam)*mid(lin,col) ;
180 for (
int lin=0; lin<nr_eq2; lin++) {
181 for (
int col=dege1; col<nr; col++)
182 oper.
set(lin+nr_eq1,col-dege1)
183 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
184 - (lam+2)*mid(lin,col)) ;
185 for (
int col=dege2; col<nr; col++)
186 oper.
set(lin+nr_eq1, col-dege2+nr_eq1)
187 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
188 + ech*ech*md2(lin,col)
189 + 2*(mxd(lin,col) + ech*md(lin,col)))
190 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
198 for (
int i=0; i<nr; i++) {
199 sr.
set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
202 Tbl xsr= sr ;
Tbl x2sr= sr ;
203 Tbl xseta= seta ;
Tbl x2seta = seta ;
204 multx2_1d(nr, &x2sr.
t, base_r) ; multx_1d(nr, &xsr.
t, base_r) ;
205 multx2_1d(nr, &x2seta.
t, base_r) ; multx_1d(nr, &xseta.
t, base_r) ;
207 for (
int i=0; i<nr_eq1; i++)
208 sec_membre.
set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
210 for (
int i=0; i<nr_eq2; i++)
211 sec_membre.
set(i+nr_eq1) = beta*beta*sr(i)
212 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
222 for (
int i=0; i<dege1; i++)
224 for (
int i=dege1; i<nr; i++)
225 res_eta.
set(i) = big_res(i-dege1) ;
226 for (
int i=0; i<dege2; i++)
228 for (
int i=dege2; i<nr; i++)
229 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
232 Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
233 Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
234 for (
int i=0 ; i<nr ; i++) {
235 sol_part_eta.
set(zone, k, j, i) = res_eta(i) ;
236 sol_part_vr.
set(zone, k, j, i) = res_vr(i) ;
237 solution_hom_un.
set(zone, k, j, i) = sol_hom1(0,i) ;
238 solution_hom_deux.
set(zone, k, j, i) = sol_hom2(1,i) ;
239 solution_hom_trois.
set(zone, k, j, i) = sol_hom2(0,i) ;
240 solution_hom_quatre.
set(zone, k, j, i) = sol_hom1(1,i) ;
250 assert(nt == mg.
get_nt(nzm1)) ;
251 assert(np == mg.
get_np(nzm1)) ;
258 for (
int k=0 ; k<np+1 ; k++) {
259 for (
int j=0 ; j<nt ; j++) {
261 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
264 int nr_eq1 = nr0 - dege1 ;
265 int nr_eq2 = nr0 - dege2 ;
266 int nrtot = nr_eq1 + nr_eq2 ;
275 for (
int lin=0; lin<nr_eq1; lin++) {
276 for (
int col=dege1; col<nr0; col++)
277 oper.
set(lin,col-dege1)
278 = m2d2(lin,col) + 4*mxd(lin,col)
279 + (2-(lam+1)*l_q*(l_q+1))*mid(lin,col) ;
280 for (
int col=dege2; col<nr0; col++)
281 oper.
set(lin,col-dege2+nr_eq1) =
282 -lam*mxd(lin,col) + 2*mid(lin,col) ;
284 for (
int lin=0; lin<nr_eq2; lin++) {
285 for (
int col=dege1; col<nr0; col++)
286 oper.
set(lin+nr_eq1,col-dege1)
287 = l_q*(l_q+1)*(lam*mxd(lin,col) + (3*lam+2)*mid(lin,col)) ;
288 for (
int col=dege2; col<nr0; col++)
289 oper.
set(lin+nr_eq1, col-dege2+nr_eq1)
290 = (lam+1)*(m2d2(lin,col) + 4*mxd(lin,col))
291 - l_q*(l_q+1)*mid(lin,col) ;
297 for (
int i=0; i<nr_eq1; i++)
299 for (
int i=0; i<nr_eq2; i++)
300 sec_membre.
set(i+nr_eq1) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
307 for (
int i=0; i<dege1; i++)
309 for (
int i=dege1; i<nr0; i++)
310 res_eta.
set(i) = big_res(i-dege1) ;
311 res_eta.
set(nr0) = 0 ;
312 res_eta.
set(nr0+1) = 0 ;
313 for (
int i=0; i<dege2; i++)
315 for (
int i=dege2; i<nr0; i++)
316 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
317 res_vr.
set(nr0) = 0 ;
318 res_vr.
set(nr0+1) = 0 ;
324 mult2_xm1_1d_cheb(nr, res_eta.
t, res_e2.
t) ;
325 mult2_xm1_1d_cheb(nr, res_vr.
t, res_v2.
t) ;
328 Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
329 Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
330 for (
int i=0 ; i<nr ; i++) {
331 sol_part_eta.
set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
332 sol_part_vr.
set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
333 solution_hom_un.
set(nzm1, k, j, i) = 0. ;
334 solution_hom_deux.
set(nzm1, k, j, i) = sol_hom2(i) ;
335 solution_hom_trois.
set(nzm1, k, j, i) = 0. ;
336 solution_hom_quatre.
set(nzm1, k, j, i) = sol_hom1(i) ;
356 int taille = 4*(nzm1-num_front)-2 ;
357 Tbl sec_membre(taille) ;
358 Matrice systeme(taille, taille) ;
360 int ligne ;
int colonne ;
364 for (
int k=0 ; k<np+1 ; k++)
365 for (
int j=0 ; j<nt ; j++) {
367 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
369 double f3_eta = lam*l_q + 3*lam + 2 ;
370 double f4_eta = 2 + 2*lam - lam*l_q ;
371 double f3_vr = (l_q+1)*(lam*l_q - 2) ;
372 double f4_vr = l_q*(lam*l_q + lam + 2) ;
376 for (
int l=0; l<taille; l++)
377 for (
int c=0; c<taille; c++)
378 systeme.
set(l,c) = 0 ;
381 nr = mg.
get_nr(num_front+1) ;
382 alpha = mpaff->
get_alpha()[num_front+1] ;
383 double echelle = mpaff->
get_beta()[num_front+1]/alpha ;
386 systeme.
set(ligne, colonne) =
pow(echelle-1.,
double(l_q-1)) ;
389 systeme.
set(ligne, colonne+1) = 1/
pow(echelle-1.,
double(l_q+2)) ;
392 systeme.
set(ligne, colonne+2) = f3_eta*
pow(echelle-1.,
double(l_q+1)) ;
394 systeme.
set(ligne, colonne+3) = f4_eta/
pow(echelle-1.,
double(l_q)) ;
395 for (
int i=0 ; i<nr ; i++)
397 sec_membre.
set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
398 else sec_membre.
set(ligne) += sol_part_eta(num_front+1, k, j, i) ;
399 sec_membre.
set(ligne) += bound_eta(num_front+1, k, j, 0) ;
403 systeme.
set(ligne, colonne) = fact_dir*l_q*
pow(echelle-1.,
double(l_q-1)) + fact_neu*l_q*(l_q-1)*
pow(echelle-1.,
double(l_q-2))/alpha ;
404 systeme.
set(ligne, colonne+1) = -fact_dir*(l_q+1)/
pow(echelle-1.,
double(l_q+2)) + fact_neu*(l_q+1)*(l_q+2)/
pow(echelle-1.,
double(l_q+3))/alpha ;
405 systeme.
set(ligne, colonne+2) = fact_dir*f3_vr*
pow(echelle-1.,
double(l_q+1)) + fact_neu*f3_vr*(l_q+1)*
pow(echelle-1.,
double(l_q))/alpha ;
406 systeme.
set(ligne, colonne+3) = fact_dir*f4_vr/
pow(echelle-1.,
double(l_q)) - fact_neu*(f4_vr*l_q/
pow(echelle-1.,
double(l_q+1)))/alpha ;
407 for (
int i=0 ; i<nr ; i++)
409 sec_membre.
set(ligne) -= fact_dir*sol_part_vr(num_front+1, k, j, i) - fact_neu*i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
410 else sec_membre.
set(ligne) += fact_dir*sol_part_vr(num_front+1, k, j, i) - fact_neu*i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
411 sec_membre.
set(ligne) += bound_vr(num_front+1, k, j, 0) ;
419 systeme.
set(ligne, colonne) =
pow(echelle+1.,
double(l_q-1)) ;
421 systeme.
set(ligne, colonne+1) = 1./
pow(echelle+1.,
double(l_q+2)) ;
423 systeme.
set(ligne, colonne+2) = f3_eta*
pow(echelle+1.,
double(l_q+1));
425 systeme.
set(ligne, colonne+3) = f4_eta/
pow(echelle+1.,
double(l_q)) ;
426 for (
int i=0 ; i<nr ; i++)
427 sec_membre.
set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
430 systeme.
set(ligne, colonne) = l_q*
pow(echelle+1.,
double(l_q-1)) ;
431 systeme.
set(ligne, colonne+1)
432 = -double(l_q+1) /
pow(echelle+1.,
double(l_q+2));
433 systeme.
set(ligne, colonne+2) = f3_vr*
pow(echelle+1.,
double(l_q+1)) ;
434 systeme.
set(ligne, colonne+3) = f4_vr/
pow(echelle+1.,
double(l_q));
435 for (
int i=0 ; i<nr ; i++)
436 sec_membre.
set(ligne) -= sol_part_vr(num_front+1, k, j, i) ;
442 systeme.
set(ligne, colonne)
443 = (l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
445 systeme.
set(ligne, colonne+1)
446 = -(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
448 systeme.
set(ligne, colonne+2)
449 = f3_eta*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha;
451 systeme.
set(ligne, colonne+3)
452 = -f4_eta*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
453 for (
int i=0 ; i<nr ; i++)
454 sec_membre.
set(ligne) -= i*i/alpha*sol_part_eta(num_front+1, k, j, i) ;
457 systeme.
set(ligne, colonne)
458 = l_q*(l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
459 systeme.
set(ligne, colonne+1)
460 = (l_q+1)*(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
461 systeme.
set(ligne, colonne+2)
462 = f3_vr*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha ;
463 systeme.
set(ligne, colonne+3)
464 = -f4_vr*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
465 for (
int i=0 ; i<nr ; i++)
466 sec_membre.
set(ligne) -= i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
472 if (num_front+2<nzm1){
473 for (
int zone=num_front+2 ; zone<nzm1 ; zone++) {
476 echelle = mpaff->
get_beta()[zone]/alpha ;
479 systeme.
set(ligne, colonne) = -
pow(echelle-1.,
double(l_q-1)) ;
481 systeme.
set(ligne, colonne+1) = -1/
pow(echelle-1.,
double(l_q+2)) ;
483 systeme.
set(ligne, colonne+2) = -f3_eta*
pow(echelle-1.,
double(l_q+1));
485 systeme.
set(ligne, colonne+3) = -f4_eta/
pow(echelle-1.,
double(l_q)) ;
486 for (
int i=0 ; i<nr ; i++)
488 sec_membre.
set(ligne) += sol_part_eta(zone, k, j, i) ;
489 else sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
492 systeme.
set(ligne, colonne) = -l_q*
pow(echelle-1.,
double(l_q-1)) ;
493 systeme.
set(ligne, colonne+1) = (l_q+1)/
pow(echelle-1.,
double(l_q+2));
494 systeme.
set(ligne, colonne+2) = -f3_vr*
pow(echelle-1.,
double(l_q+1)) ;
495 systeme.
set(ligne, colonne+3) = -f4_vr/
pow(echelle-1.,
double(l_q));
496 for (
int i=0 ; i<nr ; i++)
498 sec_membre.
set(ligne) += sol_part_vr(zone, k, j, i) ;
499 else sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i) ;
503 systeme.
set(ligne, colonne)
504 = -(l_q-1)*
pow(echelle-1.,
double(l_q-2))/alpha ;
506 systeme.
set(ligne, colonne+1)
507 = (l_q+2)/
pow(echelle-1.,
double(l_q+3))/alpha ;
509 systeme.
set(ligne, colonne+2)
510 = -f3_eta*(l_q+1)*
pow(echelle-1.,
double(l_q))/alpha;
512 systeme.
set(ligne, colonne+3)
513 = (f4_eta*l_q/
pow(echelle-1.,
double(l_q+1)))/alpha ;
514 for (
int i=0 ; i<nr ; i++)
515 if (i%2 == 0) sec_membre.
set(ligne)
516 -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
517 else sec_membre.
set(ligne) +=
518 i*i/alpha*sol_part_eta(zone, k, j, i) ;
521 systeme.
set(ligne, colonne)
522 = -l_q*(l_q-1)*
pow(echelle-1.,
double(l_q-2))/alpha ;
523 systeme.
set(ligne, colonne+1)
524 = -(l_q+1)*(l_q+2)/
pow(echelle-1.,
double(l_q+3))/alpha ;
525 systeme.
set(ligne, colonne+2)
526 = -f3_vr*(l_q+1)*
pow(echelle-1.,
double(l_q))/alpha ;
527 systeme.
set(ligne, colonne+3)
528 = (f4_vr*l_q/
pow(echelle-1.,
double(l_q+1)))/alpha ;
529 for (
int i=0 ; i<nr ; i++)
530 if (i%2 == 0) sec_membre.
set(ligne)
531 -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
532 else sec_membre.
set(ligne) +=
533 i*i/alpha*sol_part_vr(zone, k, j, i) ;
537 systeme.
set(ligne, colonne) =
pow(echelle+1.,
double(l_q-1)) ;
539 systeme.
set(ligne, colonne+1) = 1./
pow(echelle+1.,
double(l_q+2)) ;
541 systeme.
set(ligne, colonne+2) = f3_eta*
pow(echelle+1.,
double(l_q+1));
543 systeme.
set(ligne, colonne+3) = f4_eta/
pow(echelle+1.,
double(l_q)) ;
544 for (
int i=0 ; i<nr ; i++)
545 sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
548 systeme.
set(ligne, colonne) = l_q*
pow(echelle+1.,
double(l_q-1)) ;
549 systeme.
set(ligne, colonne+1)
550 = -double(l_q+1) /
pow(echelle+1.,
double(l_q+2));
551 systeme.
set(ligne, colonne+2) = f3_vr*
pow(echelle+1.,
double(l_q+1)) ;
552 systeme.
set(ligne, colonne+3) = f4_vr/
pow(echelle+1.,
double(l_q));
553 for (
int i=0 ; i<nr ; i++)
554 sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i) ;
558 systeme.
set(ligne, colonne)
559 = (l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
561 systeme.
set(ligne, colonne+1)
562 = -(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
564 systeme.
set(ligne, colonne+2)
565 = f3_eta*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha;
567 systeme.
set(ligne, colonne+3)
568 = -f4_eta*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
569 for (
int i=0 ; i<nr ; i++)
570 sec_membre.
set(ligne) -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
573 systeme.
set(ligne, colonne)
574 = l_q*(l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
575 systeme.
set(ligne, colonne+1)
576 = (l_q+1)*(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
577 systeme.
set(ligne, colonne+2)
578 = f3_vr*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha ;
579 systeme.
set(ligne, colonne+3)
580 = -f4_vr*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
581 for (
int i=0 ; i<nr ; i++)
582 sec_membre.
set(ligne) -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
593 systeme.
set(ligne, colonne) = -
pow(-2,
double(l_q+2)) ;
595 systeme.
set(ligne, colonne+1) = -f4_eta*
pow(-2,
double(l_q)) ;
596 for (
int i=0 ; i<nr ; i++)
597 if (i%2 == 0) sec_membre.
set(ligne) += sol_part_eta(nzm1, k, j, i) ;
598 else sec_membre.
set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
600 systeme.
set(ligne+1, colonne) = double(l_q+1)*
pow(-2,
double(l_q+2)) ;
601 systeme.
set(ligne+1, colonne+1) = -f4_vr*
pow(-2,
double(l_q)) ;
602 for (
int i=0 ; i<nr ; i++)
603 if (i%2 == 0) sec_membre.
set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
604 else sec_membre.
set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
608 systeme.
set(ligne, colonne) = alpha*(l_q+2)*
pow(-2,
double(l_q+3)) ;
610 systeme.
set(ligne, colonne+1) = alpha*l_q*f4_eta*
pow(-2,
double(l_q+1)) ;
611 for (
int i=0 ; i<nr ; i++)
612 if (i%2 == 0) sec_membre.
set(ligne)
613 -= -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
614 else sec_membre.
set(ligne)
615 += -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
617 systeme.
set(ligne+1, colonne)
618 = -alpha*double((l_q+1)*(l_q+2))*
pow(-2,
double(l_q+3)) ;
619 systeme.
set(ligne+1, colonne+1)
620 = alpha*double(l_q)*f4_vr*
pow(-2,
double(l_q+1)) ;
621 for (
int i=0 ; i<nr ; i++)
622 if (i%2 == 0) sec_membre.
set(ligne+1)
623 -= -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
624 else sec_membre.
set(ligne+1)
625 += -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
630 if (taille > 2) systeme.
set_band(5,5) ;
640 for (
int zone=1 ; zone<nzm1 ; zone++) {
642 for (
int i=0 ; i<nr ; i++) {
643 cf_eta.
set(zone, k, j, i) =
644 sol_part_eta(zone, k, j, i)
645 +facteurs(conte)*solution_hom_un(zone, k, j, i)
646 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i)
647 +facteurs(conte+2)*f3_eta*solution_hom_trois(zone, k, j, i)
648 +facteurs(conte+3)*f4_eta*solution_hom_quatre(zone, k, j, i) ;
649 cf_vr.
set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
650 +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
651 -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i)
652 +f3_vr*facteurs(conte+2)*solution_hom_trois(zone, k, j, i)
653 +f4_vr*facteurs(conte+3)*solution_hom_quatre(zone, k, j, i) ;
658 for (
int i=0 ; i<nr ; i++) {
659 cf_eta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
660 +facteurs(conte)*solution_hom_deux(nzm1, k, j, i)
661 +f4_eta*facteurs(conte+1)*solution_hom_quatre(nzm1, k, j, i) ;
662 cf_vr.
set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
663 -double(l_q+1)*facteurs(conte)*solution_hom_deux(nzm1, k, j, i)
664 +f4_vr*facteurs(conte+1)*solution_hom_quatre(nzm1, k, j, i) ;
675 const Valeur& limit_mu (temp_mu) ;
677 resu.
set_vr_eta_mu(vr, 0*het,
mu().poisson_dirichlet(limit_mu, num_front)) ;
683 Vector Vector::poisson_dirichlet(
double lam,
const Valeur& bound_vr,
685 int num_front)
const {
688 resu =
poisson_robin(lam, bound_vr, bound_vt, bound_vp, 1., 0., num_front) ;
696 int num_front)
const {
699 resu =
poisson_robin(lam, bound_vr, bound_vt, bound_vp, 0., 1., num_front) ;
707 double fact_dir,
double fact_neu,
708 int num_front)
const {
712 Valeur limit_vr (bound_vr) ;
725 for (
int l=0; l<nz; l++)
734 cout <<
"temp_vp" << endl <<
norme(temp_vp) << endl ;
738 Scalar vtstant (temp_vt) ;
740 source_eta = temp_vt.
dsdt() + vtstant + temp_vp.
stdsdp() ;
744 Scalar vpstant (temp_vp) ;
746 source_mu = temp_vp.
dsdt() + vpstant - temp_vt.
stdsdp() ;
752 Valeur limit_mu ((*mp).get_mg()->get_angu() ) ;
754 int nnt = (*mp).get_mg()->get_nt(1) ;
756 for (
int k=0 ; k<nnp ; k++)
757 for (
int j=0 ; j<nnt ; j++)
758 limit_mu.set(1, k, j, 0) = temp_mu.val_grid_point(1, k, j, 0) ;
759 limit_mu.set_base(temp_mu.get_spectral_va().get_base()) ;
763 Mtbl_cf lim_mu (*(limit_mu.c_cf)) ;
766 Valeur limit_eta ((*mp).get_mg()->get_angu() ) ;
768 for (
int k=0 ; k<nnp ; k++)
769 for (
int j=0 ; j<nnt ; j++)
770 limit_eta.
set(1, k, j, 0) = temp_eta.val_grid_point(1, k, j, 0) ;
771 limit_eta.set_base(temp_eta.get_spectral_va().get_base()) ;
775 Mtbl_cf lim_eta (*(limit_eta.c_cf)) ;
779 bool nullite = true ;
780 for (
int i=0; i<3; i++) {
781 assert(
cmp[i]->check_dzpuis(4)) ;
782 if (
cmp[i]->get_etat() != ETATZERO || bound_vr.
get_etat() != ETATZERO ||
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.
Class for the elementary differential operator (see the base class Diff ).
const double * get_alpha() const
Returns the pointer on the array alpha.
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...
void coef() const
Computes the coeffcients of *this.
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.
void poisson_boundary(double lambda, const Mtbl_cf &limit_vr, const Mtbl_cf &limit_eta, const Mtbl_cf &limit_mu, int num_front, double fact_dir, double fact_neu, Vector &resu) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void ylm()
Computes the coefficients of *this.
const Scalar & dsdt() const
Returns 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).
Vector poisson_robin(double lambda, const Valeur &limit_vr, const Valeur &limit_vt, const Valeur &limit_vp, double fact_dir, double fact_neu, int num_front) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Values and coefficients of a (real-value) function.
Class for the elementary differential operator (see the base class Diff ).
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tensor field of valence 1.
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
void annule_hard()
Sets the Scalar to zero in a hard way.
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Class for the elementary differential operator Identity (see the base class Diff ).
int get_etat() const
Returns the logical state.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Base_val base
Bases on which the spectral expansion is performed.
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through potentials and (see members p_eta and p_mu ), as well as the compon...
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.
Scalar sol_elliptic_boundary(Param_elliptic ¶ms, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
Resolution of a general elliptic equation, putting zero at infinity and with inner boundary condition...
int get_nzone() const
Returns the number of domains.
const Scalar & stdsdp() const
Returns of *this .
Cmp pow(const Cmp &, int)
Power .
This class contains the parameters needed to call the general elliptic solver.
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.
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Class for the elementary differential operator (see the base class Diff ).
Vector poisson_neumann(double lambda, const Valeur &limit_vr, const Valeur &limit_vt, const Valeur &limit_vp, int num_front) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
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 ).
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
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)
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 div_tant()
Division by .
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.
const Valeur & get_spectral_va() const
Returns va (read only version)
void set_poisson_vect_r(int zone, bool only_l_zero=false)
Sets the operator to in all domains, for ; and to in all domains otherwise.
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).