31 char vector_poisson_boundary_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary.C,v 1.3 2014/10/13 08:53:45 j_novak Exp $" ;
60 #include "param_elliptic.h" 62 #include "utilitaires.h" 67 int num_front,
double fact_dir,
double fact_neu,
72 for (
int i=0; i<3; i++)
73 assert(
cmp[i]->check_dzpuis(4)) ;
78 assert( mpaff != 0x0) ;
81 if (fabs(lam + 1.) < 1.e-8) {
82 cout <<
"Not implemented yet !!" << endl ;
96 int nz = mg.
get_nzone() ;
int nzm1 = nz - 1;
102 S_r.set_spectral_va().ylm() ;
117 for (
int l=0; l<nz; l++)
132 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
133 double alpha, beta, ech ;
135 assert(num_front+1 < nzm1) ;
136 for (
int zone=num_front+1 ; zone<nzm1 ; zone++) {
142 assert(nt == mg.
get_nt(zone)) ;
143 assert(np == mg.
get_np(zone)) ;
147 for (
int k=0 ; k<np+1 ; k++) {
148 for (
int j=0 ; j<nt ; j++) {
150 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
153 int nr_eq1 = nr - dege1 ;
154 int nr_eq2 = nr - dege2 ;
155 int nrtot = nr_eq1 + nr_eq2 ;
167 for (
int lin=0; lin<nr_eq1; lin++) {
168 for (
int col=dege1; col<nr; col++)
169 oper.
set(lin,col-dege1)
170 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
171 + 2*(mxd(lin,col) + ech*md(lin,col))
172 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
173 for (
int col=dege2; col<nr; col++)
174 oper.
set(lin,col-dege2+nr_eq1)
175 = lam*(mxd(lin,col) + ech*md(lin,col)) + 2*(1+lam)*mid(lin,col) ;
177 for (
int lin=0; lin<nr_eq2; lin++) {
178 for (
int col=dege1; col<nr; col++)
179 oper.
set(lin+nr_eq1,col-dege1)
180 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
181 - (lam+2)*mid(lin,col)) ;
182 for (
int col=dege2; col<nr; col++)
183 oper.
set(lin+nr_eq1, col-dege2+nr_eq1)
184 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
185 + ech*ech*md2(lin,col)
186 + 2*(mxd(lin,col) + ech*md(lin,col)))
187 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
195 for (
int i=0; i<nr; i++) {
196 sr.
set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
199 Tbl xsr= sr ;
Tbl x2sr= sr ;
200 Tbl xseta= seta ;
Tbl x2seta = seta ;
201 multx2_1d(nr, &x2sr.
t, base_r) ; multx_1d(nr, &xsr.
t, base_r) ;
202 multx2_1d(nr, &x2seta.
t, base_r) ; multx_1d(nr, &xseta.
t, base_r) ;
204 for (
int i=0; i<nr_eq1; i++)
205 sec_membre.
set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
207 for (
int i=0; i<nr_eq2; i++)
208 sec_membre.
set(i+nr_eq1) = beta*beta*sr(i)
209 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
219 for (
int i=0; i<dege1; i++)
221 for (
int i=dege1; i<nr; i++)
222 res_eta.
set(i) = big_res(i-dege1) ;
223 for (
int i=0; i<dege2; i++)
225 for (
int i=dege2; i<nr; i++)
226 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
229 Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
230 Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
231 for (
int i=0 ; i<nr ; i++) {
232 sol_part_eta.
set(zone, k, j, i) = res_eta(i) ;
233 sol_part_vr.
set(zone, k, j, i) = res_vr(i) ;
234 solution_hom_un.
set(zone, k, j, i) = sol_hom1(0,i) ;
235 solution_hom_deux.
set(zone, k, j, i) = sol_hom2(1,i) ;
236 solution_hom_trois.
set(zone, k, j, i) = sol_hom2(0,i) ;
237 solution_hom_quatre.
set(zone, k, j, i) = sol_hom1(1,i) ;
247 assert(nt == mg.
get_nt(nzm1)) ;
248 assert(np == mg.
get_np(nzm1)) ;
255 for (
int k=0 ; k<np+1 ; k++) {
256 for (
int j=0 ; j<nt ; j++) {
258 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
261 int nr_eq1 = nr0 - dege1 ;
262 int nr_eq2 = nr0 - dege2 ;
263 int nrtot = nr_eq1 + nr_eq2 ;
272 for (
int lin=0; lin<nr_eq1; lin++) {
273 for (
int col=dege1; col<nr0; col++)
274 oper.
set(lin,col-dege1)
275 = m2d2(lin,col) + 4*mxd(lin,col)
276 + (2-(lam+1)*l_q*(l_q+1))*mid(lin,col) ;
277 for (
int col=dege2; col<nr0; col++)
278 oper.
set(lin,col-dege2+nr_eq1) =
279 -lam*mxd(lin,col) + 2*mid(lin,col) ;
281 for (
int lin=0; lin<nr_eq2; lin++) {
282 for (
int col=dege1; col<nr0; col++)
283 oper.
set(lin+nr_eq1,col-dege1)
284 = l_q*(l_q+1)*(lam*mxd(lin,col) + (3*lam+2)*mid(lin,col)) ;
285 for (
int col=dege2; col<nr0; col++)
286 oper.
set(lin+nr_eq1, col-dege2+nr_eq1)
287 = (lam+1)*(m2d2(lin,col) + 4*mxd(lin,col))
288 - l_q*(l_q+1)*mid(lin,col) ;
294 for (
int i=0; i<nr_eq1; i++)
296 for (
int i=0; i<nr_eq2; i++)
297 sec_membre.
set(i+nr_eq1) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
304 for (
int i=0; i<dege1; i++)
306 for (
int i=dege1; i<nr0; i++)
307 res_eta.
set(i) = big_res(i-dege1) ;
308 res_eta.
set(nr0) = 0 ;
309 res_eta.
set(nr0+1) = 0 ;
310 for (
int i=0; i<dege2; i++)
312 for (
int i=dege2; i<nr0; i++)
313 res_vr.
set(i) = big_res(i-dege2+nr_eq1) ;
314 res_vr.
set(nr0) = 0 ;
315 res_vr.
set(nr0+1) = 0 ;
321 mult2_xm1_1d_cheb(nr, res_eta.
t, res_e2.
t) ;
322 mult2_xm1_1d_cheb(nr, res_vr.
t, res_v2.
t) ;
325 Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
326 Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
327 for (
int i=0 ; i<nr ; i++) {
328 sol_part_eta.
set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
329 sol_part_vr.
set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
330 solution_hom_un.
set(nzm1, k, j, i) = 0. ;
331 solution_hom_deux.
set(nzm1, k, j, i) = sol_hom2(i) ;
332 solution_hom_trois.
set(nzm1, k, j, i) = 0. ;
333 solution_hom_quatre.
set(nzm1, k, j, i) = sol_hom1(i) ;
353 int taille = 4*(nzm1-num_front)-2 ;
354 Tbl sec_membre(taille) ;
355 Matrice systeme(taille, taille) ;
357 int ligne ;
int colonne ;
361 for (
int k=0 ; k<np+1 ; k++)
362 for (
int j=0 ; j<nt ; j++) {
364 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
366 double f3_eta = lam*l_q + 3*lam + 2 ;
367 double f4_eta = 2 + 2*lam - lam*l_q ;
368 double f3_vr = (l_q+1)*(lam*l_q - 2) ;
369 double f4_vr = l_q*(lam*l_q + lam + 2) ;
373 for (
int l=0; l<taille; l++)
374 for (
int c=0; c<taille; c++)
375 systeme.
set(l,c) = 0 ;
378 nr = mg.
get_nr(num_front+1) ;
379 alpha = mpaff->
get_alpha()[num_front+1] ;
380 double echelle = mpaff->
get_beta()[num_front+1]/alpha ;
383 systeme.
set(ligne, colonne) =
pow(echelle-1.,
double(l_q-1)) ;
386 systeme.
set(ligne, colonne+1) = 1/
pow(echelle-1.,
double(l_q+2)) ;
389 systeme.
set(ligne, colonne+2) = f3_eta*
pow(echelle-1.,
double(l_q+1)) ;
391 systeme.
set(ligne, colonne+3) = f4_eta/
pow(echelle-1.,
double(l_q)) ;
392 for (
int i=0 ; i<nr ; i++)
394 sec_membre.
set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
395 else sec_membre.
set(ligne) += sol_part_eta(num_front+1, k, j, i) ;
396 sec_membre.
set(ligne) += bound_eta(num_front+1, k, j, 0) ;
400 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 ;
401 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 ;
402 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 ;
403 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 ;
404 for (
int i=0 ; i<nr ; i++)
406 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) ;
407 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) ;
408 sec_membre.
set(ligne) += bound_vr(num_front+1, k, j, 0) ;
416 systeme.
set(ligne, colonne) =
pow(echelle+1.,
double(l_q-1)) ;
418 systeme.
set(ligne, colonne+1) = 1./
pow(echelle+1.,
double(l_q+2)) ;
420 systeme.
set(ligne, colonne+2) = f3_eta*
pow(echelle+1.,
double(l_q+1));
422 systeme.
set(ligne, colonne+3) = f4_eta/
pow(echelle+1.,
double(l_q)) ;
423 for (
int i=0 ; i<nr ; i++)
424 sec_membre.
set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
427 systeme.
set(ligne, colonne) = l_q*
pow(echelle+1.,
double(l_q-1)) ;
428 systeme.
set(ligne, colonne+1)
429 = -double(l_q+1) /
pow(echelle+1.,
double(l_q+2));
430 systeme.
set(ligne, colonne+2) = f3_vr*
pow(echelle+1.,
double(l_q+1)) ;
431 systeme.
set(ligne, colonne+3) = f4_vr/
pow(echelle+1.,
double(l_q));
432 for (
int i=0 ; i<nr ; i++)
433 sec_membre.
set(ligne) -= sol_part_vr(num_front+1, k, j, i) ;
439 systeme.
set(ligne, colonne)
440 = (l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
442 systeme.
set(ligne, colonne+1)
443 = -(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
445 systeme.
set(ligne, colonne+2)
446 = f3_eta*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha;
448 systeme.
set(ligne, colonne+3)
449 = -f4_eta*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
450 for (
int i=0 ; i<nr ; i++)
451 sec_membre.
set(ligne) -= i*i/alpha*sol_part_eta(num_front+1, k, j, i) ;
454 systeme.
set(ligne, colonne)
455 = l_q*(l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
456 systeme.
set(ligne, colonne+1)
457 = (l_q+1)*(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
458 systeme.
set(ligne, colonne+2)
459 = f3_vr*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha ;
460 systeme.
set(ligne, colonne+3)
461 = -f4_vr*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
462 for (
int i=0 ; i<nr ; i++)
463 sec_membre.
set(ligne) -= i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
469 if (num_front+2<nzm1){
470 for (
int zone=num_front+2 ; zone<nzm1 ; zone++) {
473 echelle = mpaff->
get_beta()[zone]/alpha ;
476 systeme.
set(ligne, colonne) = -
pow(echelle-1.,
double(l_q-1)) ;
478 systeme.
set(ligne, colonne+1) = -1/
pow(echelle-1.,
double(l_q+2)) ;
480 systeme.
set(ligne, colonne+2) = -f3_eta*
pow(echelle-1.,
double(l_q+1));
482 systeme.
set(ligne, colonne+3) = -f4_eta/
pow(echelle-1.,
double(l_q)) ;
483 for (
int i=0 ; i<nr ; i++)
485 sec_membre.
set(ligne) += sol_part_eta(zone, k, j, i) ;
486 else sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
489 systeme.
set(ligne, colonne) = -l_q*
pow(echelle-1.,
double(l_q-1)) ;
490 systeme.
set(ligne, colonne+1) = (l_q+1)/
pow(echelle-1.,
double(l_q+2));
491 systeme.
set(ligne, colonne+2) = -f3_vr*
pow(echelle-1.,
double(l_q+1)) ;
492 systeme.
set(ligne, colonne+3) = -f4_vr/
pow(echelle-1.,
double(l_q));
493 for (
int i=0 ; i<nr ; i++)
495 sec_membre.
set(ligne) += sol_part_vr(zone, k, j, i) ;
496 else sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i) ;
500 systeme.
set(ligne, colonne)
501 = -(l_q-1)*
pow(echelle-1.,
double(l_q-2))/alpha ;
503 systeme.
set(ligne, colonne+1)
504 = (l_q+2)/
pow(echelle-1.,
double(l_q+3))/alpha ;
506 systeme.
set(ligne, colonne+2)
507 = -f3_eta*(l_q+1)*
pow(echelle-1.,
double(l_q))/alpha;
509 systeme.
set(ligne, colonne+3)
510 = (f4_eta*l_q/
pow(echelle-1.,
double(l_q+1)))/alpha ;
511 for (
int i=0 ; i<nr ; i++)
512 if (i%2 == 0) sec_membre.
set(ligne)
513 -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
514 else sec_membre.
set(ligne) +=
515 i*i/alpha*sol_part_eta(zone, k, j, i) ;
518 systeme.
set(ligne, colonne)
519 = -l_q*(l_q-1)*
pow(echelle-1.,
double(l_q-2))/alpha ;
520 systeme.
set(ligne, colonne+1)
521 = -(l_q+1)*(l_q+2)/
pow(echelle-1.,
double(l_q+3))/alpha ;
522 systeme.
set(ligne, colonne+2)
523 = -f3_vr*(l_q+1)*
pow(echelle-1.,
double(l_q))/alpha ;
524 systeme.
set(ligne, colonne+3)
525 = (f4_vr*l_q/
pow(echelle-1.,
double(l_q+1)))/alpha ;
526 for (
int i=0 ; i<nr ; i++)
527 if (i%2 == 0) sec_membre.
set(ligne)
528 -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
529 else sec_membre.
set(ligne) +=
530 i*i/alpha*sol_part_vr(zone, k, j, i) ;
534 systeme.
set(ligne, colonne) =
pow(echelle+1.,
double(l_q-1)) ;
536 systeme.
set(ligne, colonne+1) = 1./
pow(echelle+1.,
double(l_q+2)) ;
538 systeme.
set(ligne, colonne+2) = f3_eta*
pow(echelle+1.,
double(l_q+1));
540 systeme.
set(ligne, colonne+3) = f4_eta/
pow(echelle+1.,
double(l_q)) ;
541 for (
int i=0 ; i<nr ; i++)
542 sec_membre.
set(ligne) -= sol_part_eta(zone, k, j, i) ;
545 systeme.
set(ligne, colonne) = l_q*
pow(echelle+1.,
double(l_q-1)) ;
546 systeme.
set(ligne, colonne+1)
547 = -double(l_q+1) /
pow(echelle+1.,
double(l_q+2));
548 systeme.
set(ligne, colonne+2) = f3_vr*
pow(echelle+1.,
double(l_q+1)) ;
549 systeme.
set(ligne, colonne+3) = f4_vr/
pow(echelle+1.,
double(l_q));
550 for (
int i=0 ; i<nr ; i++)
551 sec_membre.
set(ligne) -= sol_part_vr(zone, k, j, i) ;
555 systeme.
set(ligne, colonne)
556 = (l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
558 systeme.
set(ligne, colonne+1)
559 = -(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
561 systeme.
set(ligne, colonne+2)
562 = f3_eta*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha;
564 systeme.
set(ligne, colonne+3)
565 = -f4_eta*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
566 for (
int i=0 ; i<nr ; i++)
567 sec_membre.
set(ligne) -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
570 systeme.
set(ligne, colonne)
571 = l_q*(l_q-1) *
pow(echelle+1.,
double(l_q-2))/alpha ;
572 systeme.
set(ligne, colonne+1)
573 = (l_q+1)*(l_q+2) /
pow(echelle+1.,
double(l_q+3))/alpha ;
574 systeme.
set(ligne, colonne+2)
575 = f3_vr*(l_q+1) *
pow(echelle+1.,
double(l_q))/alpha ;
576 systeme.
set(ligne, colonne+3)
577 = -f4_vr*l_q /
pow(echelle+1.,
double(l_q+1))/alpha ;
578 for (
int i=0 ; i<nr ; i++)
579 sec_membre.
set(ligne) -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
590 systeme.
set(ligne, colonne) = -
pow(-2,
double(l_q+2)) ;
592 systeme.
set(ligne, colonne+1) = -f4_eta*
pow(-2,
double(l_q)) ;
593 for (
int i=0 ; i<nr ; i++)
594 if (i%2 == 0) sec_membre.
set(ligne) += sol_part_eta(nzm1, k, j, i) ;
595 else sec_membre.
set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
597 systeme.
set(ligne+1, colonne) = double(l_q+1)*
pow(-2,
double(l_q+2)) ;
598 systeme.
set(ligne+1, colonne+1) = -f4_vr*
pow(-2,
double(l_q)) ;
599 for (
int i=0 ; i<nr ; i++)
600 if (i%2 == 0) sec_membre.
set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
601 else sec_membre.
set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
605 systeme.
set(ligne, colonne) = alpha*(l_q+2)*
pow(-2,
double(l_q+3)) ;
607 systeme.
set(ligne, colonne+1) = alpha*l_q*f4_eta*
pow(-2,
double(l_q+1)) ;
608 for (
int i=0 ; i<nr ; i++)
609 if (i%2 == 0) sec_membre.
set(ligne)
610 -= -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
611 else sec_membre.
set(ligne)
612 += -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
614 systeme.
set(ligne+1, colonne)
615 = -alpha*double((l_q+1)*(l_q+2))*
pow(-2,
double(l_q+3)) ;
616 systeme.
set(ligne+1, colonne+1)
617 = alpha*double(l_q)*f4_vr*
pow(-2,
double(l_q+1)) ;
618 for (
int i=0 ; i<nr ; i++)
619 if (i%2 == 0) sec_membre.
set(ligne+1)
620 -= -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
621 else sec_membre.
set(ligne+1)
622 += -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
627 if (taille > 2) systeme.
set_band(5,5) ;
637 for (
int zone=1 ; zone<nzm1 ; zone++) {
639 for (
int i=0 ; i<nr ; i++) {
640 cf_eta.
set(zone, k, j, i) =
641 sol_part_eta(zone, k, j, i)
642 +facteurs(conte)*solution_hom_un(zone, k, j, i)
643 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i)
644 +facteurs(conte+2)*f3_eta*solution_hom_trois(zone, k, j, i)
645 +facteurs(conte+3)*f4_eta*solution_hom_quatre(zone, k, j, i) ;
646 cf_vr.
set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
647 +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
648 -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i)
649 +f3_vr*facteurs(conte+2)*solution_hom_trois(zone, k, j, i)
650 +f4_vr*facteurs(conte+3)*solution_hom_quatre(zone, k, j, i) ;
655 for (
int i=0 ; i<nr ; i++) {
656 cf_eta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
657 +facteurs(conte)*solution_hom_deux(nzm1, k, j, i)
658 +f4_eta*facteurs(conte+1)*solution_hom_quatre(nzm1, k, j, i) ;
659 cf_vr.
set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
660 -double(l_q+1)*facteurs(conte)*solution_hom_deux(nzm1, k, j, i)
661 +f4_vr*facteurs(conte+1)*solution_hom_quatre(nzm1, k, j, i) ;
672 const Valeur& limit_mu (temp_mu) ;
674 resu.
set_vr_eta_mu(vr, 0*het,
mu().poisson_dirichlet(limit_mu, num_front)) ;
680 Vector Vector::poisson_dirichlet(
double lam,
const Valeur& bound_vr,
682 int num_front)
const {
685 resu =
poisson_robin(lam, bound_vr, bound_vt, bound_vp, 1., 0., num_front) ;
693 int num_front)
const {
696 resu =
poisson_robin(lam, bound_vr, bound_vt, bound_vp, 0., 1., num_front) ;
704 double fact_dir,
double fact_neu,
705 int num_front)
const {
709 Valeur limit_vr (bound_vr) ;
722 for (
int l=0; l<nz; l++)
731 cout <<
"temp_vp" << endl <<
norme(temp_vp) << endl ;
735 Scalar vtstant (temp_vt) ;
737 source_eta = temp_vt.
dsdt() + vtstant + temp_vp.
stdsdp() ;
741 Scalar vpstant (temp_vp) ;
743 source_mu = temp_vp.
dsdt() + vpstant - temp_vt.
stdsdp() ;
749 Valeur limit_mu ((*mp).get_mg()->get_angu() ) ;
751 int nnt = (*mp).get_mg()->get_nt(1) ;
753 for (
int k=0 ; k<nnp ; k++)
754 for (
int j=0 ; j<nnt ; j++)
755 limit_mu.set(1, k, j, 0) = temp_mu.val_grid_point(1, k, j, 0) ;
756 limit_mu.set_base(temp_mu.get_spectral_va().get_base()) ;
760 Mtbl_cf lim_mu (*(limit_mu.c_cf)) ;
763 Valeur limit_eta ((*mp).get_mg()->get_angu() ) ;
765 for (
int k=0 ; k<nnp ; k++)
766 for (
int j=0 ; j<nnt ; j++)
767 limit_eta.
set(1, k, j, 0) = temp_eta.val_grid_point(1, k, j, 0) ;
768 limit_eta.set_base(temp_eta.get_spectral_va().get_base()) ;
772 Mtbl_cf lim_eta (*(limit_eta.c_cf)) ;
776 bool nullite = true ;
777 for (
int i=0; i<3; i++) {
778 assert(
cmp[i]->check_dzpuis(4)) ;
779 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).