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;
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++) {
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++)
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++)
723 for (
int j=0; j<
mp->
get_mg()->get_nt(l); j++)
724 for (
int k=0; k<
mp->
get_mg()->get_np(l); k++) {
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() ) ;
750 int nnp = (*mp).get_mg()->get_np(1) ;
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++)
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++)
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 ||
Bases of the spectral expansions.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Spherical orthonormal vectorial bases (triads).
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator Identity (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
const double * get_beta() const
Returns the pointer on the array beta.
const double * get_alpha() const
Returns the pointer on the array alpha.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int j, int i)
Read/write of a particuliar element.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void set_band(int up, int low) const
Calculate the band storage of *std.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Coefficients storage for the multi-domain spectral method.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
This class contains the parameters needed to call the general elliptic solver.
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.
Tensor field of valence 0 (or component of a tensorial field).
const Scalar & dsdt() const
Returns of *this .
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
const Scalar & stdsdp() const
Returns of *this .
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
void div_tant()
Division by .
void annule_hard()
Sets the Scalar to zero in a hard way.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
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...
const Valeur & get_spectral_va() const
Returns va (read only version)
Valeur & set_spectral_va()
Returns va (read/write version)
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
void annule_hard()
Sets the Tbl to zero in a hard way.
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double * t
The array of double.
Values and coefficients of a (real-value) function.
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
int get_etat() const
Returns the logical state.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
void ylm()
Computes the coefficients of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
void ylm_i()
Inverse of ylm()
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Base_val base
Bases on which the spectral expansion is performed.
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Tensor field of valence 1.
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...
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
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...
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
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...
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...
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
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.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.