28 char vector_poisson_boundary2_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $" ;
71 #include "param_elliptic.h"
73 #include "graphique.h"
75 #include "utilitaires.h"
80 void Vector::poisson_boundary2(
double lam,
Vector& resu,
Scalar boundvr,
Scalar boundeta,
Scalar boundmu,
double dir_vr,
double neum_vr,
double dir_eta,
double neum_eta,
double dir_mu,
double neum_mu)
const {
84 for (
int i=0; i<3; i++)
85 assert(
cmp[i]->check_dzpuis(4)) ;
90 assert( mpaff != 0x0) ;
93 if (fabs(lam + 1.) < 1.e-8) {
95 cout <<
" lambda = -1 is not supported by this code!" << endl;
102 Scalar bound_vr2 = boundvr;
104 Scalar bound_eta2 = boundeta;
118 int nz = mg.
get_nzone() ;
int nzm1 = nz - 1;
164 for (
int l=0; l<nz; l++){
182 double alpha = mpaff->
get_alpha()[1] ;
double alp2 = alpha*alpha ;
183 double beta = mpaff->
get_beta()[1] ;
185 int l_q = 0 ;
int m_q = 0;
int base_r = 0 ;
186 for (
int zone=1 ; zone<nzm1 ; zone++) {
189 assert(nt == mg.
get_nt(zone)) ;
190 assert(np == mg.
get_np(zone)) ;
193 double ech = beta / alpha ;
197 for (
int k=0 ; k<np+1 ; k++) {
198 for (
int j=0 ; j<nt ; j++) {
200 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
203 int nr_eq1 = nr - dege1 ;
204 int nr_eq2 = nr - dege2 ;
217 for (
int lin=0; lin<nr_eq1; lin++) {
218 for (
int col=0; col<nr; col++)
220 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
221 + 2*(mxd(lin,col) + ech*md(lin,col))
222 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
223 for (
int col=0; col<nr; col++)
225 = lam*(mxd(lin,col) + ech*md(lin,col))
226 + 2*(1+lam)*mid(lin,col) ;
228 oper.
set(nr_eq1, 0) = 1 ;
229 for (
int col=1; col<2*nr; col++)
230 oper.
set(nr_eq1, col) = 0 ;
231 oper.
set(nr_eq1+1, 0) = 0 ;
232 oper.
set(nr_eq1+1, 1) = 1 ;
233 for (
int col=2; col<2*nr; col++)
234 oper.
set(nr_eq1+1, col) = 0 ;
235 for (
int lin=0; lin<nr_eq2; lin++) {
236 for (
int col=0; col<nr; col++)
238 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
239 - (lam+2)*mid(lin,col)) ;
240 for (
int col=0; col<nr; col++)
241 oper.
set(lin+nr, col+nr)
242 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
243 + ech*ech*md2(lin,col)
244 + 2*(mxd(lin,col) + ech*md(lin,col)))
245 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
247 for (
int col=0; col<nr; col++)
248 oper.
set(nr+nr_eq2, col) = 0 ;
249 oper.
set(nr+nr_eq2, nr) = 1 ;
250 for (
int col=nr+1; col<2*nr; col++)
251 oper.
set(nr+nr_eq2, col) = 0 ;
252 for (
int col=0; col<nr+1; col++)
253 oper.
set(nr+nr_eq2+1, col) = 0 ;
254 oper.
set(nr+nr_eq2+1, nr+1) = 1 ;
255 for (
int col=nr+2; col<2*nr; col++)
256 oper.
set(nr+nr_eq2+1, col) = 0 ;
264 for (
int i=0; i<nr; i++) {
268 Tbl xsr= sr ;
Tbl x2sr= sr ;
269 Tbl xseta= seta ;
Tbl x2seta = seta ;
270 multx2_1d(nr, &x2sr.
t, base_r) ; multx_1d(nr, &xsr.
t, base_r) ;
271 multx2_1d(nr, &x2seta.
t, base_r) ; multx_1d(nr, &xseta.
t, base_r) ;
273 for (
int i=0; i<nr_eq1; i++)
274 sec_membre.
set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
276 sec_membre.
set(nr_eq1) = 0 ;
277 sec_membre.
set(nr_eq1+1) = 0 ;
278 for (
int i=0; i<nr_eq2; i++)
279 sec_membre.
set(i+nr) = beta*beta*sr(i)
280 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
281 sec_membre.
set(nr+nr_eq2) = 0 ;
282 sec_membre.
set(nr+nr_eq2+1) = 0 ;
287 for (
int i=0; i<nr; i++) {
288 sol_part_eta.
set(zone, k, j, i) = big_res(i) ;
289 sol_part_vr.
set(zone, k, j, i) = big_res(nr+i) ;
295 sec_membre.
set(nr_eq1) = 1 ;
296 big_res = oper.
inverse(sec_membre) ;
297 for (
int i=0 ; i<nr ; i++) {
298 sol_hom_un_eta.
set(zone, k, j, i) = big_res(i) ;
299 sol_hom_un_vr.
set(zone, k, j, i) = big_res(nr+i) ;
301 sec_membre.
set(nr_eq1) = 0 ;
302 sec_membre.
set(nr_eq1+1) = 1 ;
303 big_res = oper.
inverse(sec_membre) ;
304 for (
int i=0 ; i<nr ; i++) {
305 sol_hom_deux_eta.
set(zone, k, j, i) = big_res(i) ;
306 sol_hom_deux_vr.
set(zone, k, j, i) = big_res(nr+i) ;
308 sec_membre.
set(nr_eq1+1) = 0 ;
309 sec_membre.
set(nr+nr_eq2) = 1 ;
310 big_res = oper.
inverse(sec_membre) ;
311 for (
int i=0 ; i<nr ; i++) {
312 sol_hom_trois_eta.
set(zone, k, j, i) = big_res(i) ;
313 sol_hom_trois_vr.
set(zone, k, j, i) = big_res(nr+i) ;
315 sec_membre.
set(nr+nr_eq2) = 0 ;
316 sec_membre.
set(nr+nr_eq2+1) = 1 ;
317 big_res = oper.
inverse(sec_membre) ;
318 for (
int i=0 ; i<nr ; i++) {
319 sol_hom_quatre_eta.
set(zone, k, j, i) = big_res(i) ;
320 sol_hom_quatre_vr.
set(zone, k, j, i) = big_res(nr+i) ;
331 assert(nt == mg.
get_nt(nzm1)) ;
332 assert(np == mg.
get_np(nzm1)) ;
333 alpha = mpaff->
get_alpha()[nzm1] ; alp2 = alpha*alpha ;
338 for (
int k=0 ; k<np+1 ; k++) {
339 for (
int j=0 ; j<nt ; j++) {
341 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
343 int dege2 = (l_q == 1 ? 2 : 3);
344 int nr_eq1 = nr - dege1 ;
345 int nr_eq2 = nr - dege2 ;
355 for (
int lin=0; lin<nr_eq1; lin++) {
356 for (
int col=0; col<nr; col++)
358 = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
359 for (
int col=0; col<nr; col++)
360 oper.
set(lin,col+nr) =
361 (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
363 for (
int col=0; col<nr; col++)
364 oper.
set(nr_eq1, col) = 1 ;
365 for (
int col=nr; col<nrtot; col++)
366 oper.
set(nr_eq1, col) = 0 ;
368 for (
int col=0; col<nr; col++) {
369 oper.
set(nr_eq1+1, col) = pari*col*col ;
372 for (
int col=nr; col<nrtot; col++)
373 oper.
set(nr_eq1+1, col) = 0 ;
374 oper.
set(nr_eq1+2, 0) = 1 ;
375 for (
int col=1; col<nrtot; col++)
376 oper.
set(nr_eq1+2, col) = 0 ;
377 for (
int lin=0; lin<nr_eq2; lin++) {
378 for (
int col=0; col<nr; col++)
379 oper.
set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
380 + (lam+2)*ms2(lin,col))) / alp2 ;
381 for (
int col=0; col<nr; col++)
382 oper.
set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
383 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
385 for (
int col=0; col<nr; col++)
386 oper.
set(nr+nr_eq2, col) = 0 ;
387 for (
int col=nr; col<nrtot; col++)
388 oper.
set(nr+nr_eq2, col) = 1 ;
389 for (
int col=0; col<nr; col++)
390 oper.
set(nr+nr_eq2+1, col) = 0 ;
392 for (
int col=0; col<nr; col++) {
393 oper.
set(nr+nr_eq2+1, nr+col) = pari*col*col ;
397 for (
int col=0; col<nr; col++)
398 oper.
set(nr+nr_eq2+2, col) = 0 ;
399 oper.
set(nr+nr_eq2+2, nr) = 1 ;
400 for (
int col=nr+1; col<nrtot; col++)
401 oper.
set(nr+nr_eq2+2, col) = 0 ;
407 for (
int i=0; i<nr_eq1; i++)
409 for (
int i=nr_eq1; i<nr; i++)
410 sec_membre.
set(i) = 0 ;
411 for (
int i=0; i<nr_eq2; i++)
413 for (
int i=nr_eq2; i<nr; i++)
414 sec_membre.
set(nr+i) = 0 ;
419 for (
int i=0; i<nr; i++) {
420 sol_part_eta.
set(nzm1, k, j, i) = big_res(i) ;
421 sol_part_vr.
set(nzm1, k, j, i) = big_res(i+nr) ;
427 Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
428 Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
429 double fac_eta = lam + 2. ;
430 double fac_vr = 2*lam + 2. ;
431 for (
int i=0 ; i<nr ; i++) {
432 sol_hom_deux_eta.
set(nzm1, k, j, i) = sol_hom2(i) ;
433 sol_hom_quatre_eta.
set(nzm1, k, j, i) = fac_eta*sol_hom1(i) ;
434 sol_hom_deux_vr.
set(nzm1, k, j, i) = -2*sol_hom2(i) ;
435 sol_hom_quatre_vr.
set(nzm1, k, j, i) = fac_vr*sol_hom1(i) ;
440 sec_membre.
set(nr-1) = 1 ;
441 big_res = oper.
inverse(sec_membre) ;
443 for (
int i=0; i<nr; i++) {
444 sol_hom_deux_eta.
set(nzm1, k, j, i) = big_res(i) ;
445 sol_hom_deux_vr.
set(nzm1, k, j, i) = big_res(nr+i) ;
447 sec_membre.
set(nr-1) = 0 ;
448 sec_membre.
set(2*nr-1) = 1 ;
449 big_res = oper.
inverse(sec_membre) ;
450 for (
int i=0; i<nr; i++) {
451 sol_hom_quatre_eta.
set(nzm1, k, j, i) = big_res(i) ;
452 sol_hom_quatre_vr.
set(nzm1, k, j, i) = big_res(nr+i) ;
473 int taille = 4*nzm1 -2 ;
474 Tbl sec_membre(taille) ;
475 Matrice systeme(taille, taille) ;
477 int ligne ;
int colonne ;
481 for (
int k=0 ; k<np+1 ; k++)
482 for (
int j=0 ; j<nt ; j++) {
484 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
499 systeme.
set(ligne, colonne)
501 systeme.
set(ligne, colonne+1)
503 systeme.
set(ligne, colonne+2)
505 systeme.
set(ligne, colonne+3)
516 for (
int i=0; i<nr; i++) {
517 systeme.
set(ligne, colonne)
518 -= neum_eta*pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
519 systeme.
set(ligne, colonne+1)
520 -= neum_eta*pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
521 systeme.
set(ligne, colonne+2)
522 -= neum_eta*pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
523 systeme.
set(ligne, colonne+3)
524 -= neum_eta*pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
525 sec_membre.
set(ligne)
526 += neum_eta*pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
530 sec_membre.
set(ligne) -= (*bound_eta).val_in_bound_jk(1,j,k);
536 systeme.
set(ligne, colonne)
538 systeme.
set(ligne, colonne+1)
540 systeme.
set(ligne, colonne+2)
542 systeme.
set(ligne, colonne+3)
550 for (
int i=0; i<nr; i++) {
551 systeme.
set(ligne, colonne)
552 -= neum_vr*pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
553 systeme.
set(ligne, colonne+1)
554 -= neum_vr*pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
555 systeme.
set(ligne, colonne+2)
556 -= neum_vr*pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
557 systeme.
set(ligne, colonne+3)
558 -= neum_vr*pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
559 sec_membre.
set(ligne)
560 += neum_vr*pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
564 sec_membre.
set(ligne) -= (*bound_vr).val_in_bound_jk(1,j,k);
572 systeme.
set(ligne, colonne)
574 systeme.
set(ligne, colonne+1)
576 systeme.
set(ligne, colonne+2)
578 systeme.
set(ligne, colonne+3)
583 systeme.
set(ligne, colonne)
585 systeme.
set(ligne, colonne+1)
587 systeme.
set(ligne, colonne+2)
589 systeme.
set(ligne, colonne+3)
596 for (
int i=0; i<nr; i++) {
597 systeme.
set(ligne, colonne)
598 += pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
599 systeme.
set(ligne, colonne+1)
600 += pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
601 systeme.
set(ligne, colonne+2)
602 += pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
603 systeme.
set(ligne, colonne+3)
604 += pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
605 sec_membre.
set(ligne)
606 -= pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
612 for (
int i=0; i<nr; i++) {
613 systeme.
set(ligne, colonne)
614 += pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
615 systeme.
set(ligne, colonne+1)
616 += pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
617 systeme.
set(ligne, colonne+2)
618 += pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
619 systeme.
set(ligne, colonne+3)
620 += pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
621 sec_membre.
set(ligne)
622 -= pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
629 cout <<
"WARNING!! Mapping must contain at least 2 shells!!" << endl;}
630 for (
int zone=2 ; zone<nzm1 ; zone++) {
635 systeme.
set(ligne, colonne)
637 systeme.
set(ligne, colonne+1)
639 systeme.
set(ligne, colonne+2)
641 systeme.
set(ligne, colonne+3)
646 systeme.
set(ligne, colonne)
648 systeme.
set(ligne, colonne+1)
650 systeme.
set(ligne, colonne+2)
652 systeme.
set(ligne, colonne+3)
659 for (
int i=0; i<nr; i++) {
660 systeme.
set(ligne, colonne)
661 -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
662 systeme.
set(ligne, colonne+1)
663 -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
664 systeme.
set(ligne, colonne+2)
665 -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
666 systeme.
set(ligne, colonne+3)
667 -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
668 sec_membre.
set(ligne)
669 += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
675 for (
int i=0; i<nr; i++) {
676 systeme.
set(ligne, colonne)
677 -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
678 systeme.
set(ligne, colonne+1)
679 -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
680 systeme.
set(ligne, colonne+2)
681 -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
682 systeme.
set(ligne, colonne+3)
683 -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
684 sec_membre.
set(ligne)
685 += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
690 systeme.
set(ligne, colonne)
692 systeme.
set(ligne, colonne+1)
694 systeme.
set(ligne, colonne+2)
696 systeme.
set(ligne, colonne+3)
701 systeme.
set(ligne, colonne)
703 systeme.
set(ligne, colonne+1)
705 systeme.
set(ligne, colonne+2)
707 systeme.
set(ligne, colonne+3)
714 for (
int i=0; i<nr; i++) {
715 systeme.
set(ligne, colonne)
716 += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
717 systeme.
set(ligne, colonne+1)
718 += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
719 systeme.
set(ligne, colonne+2)
720 += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
721 systeme.
set(ligne, colonne+3)
722 += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
723 sec_membre.
set(ligne)
724 -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
730 for (
int i=0; i<nr; i++) {
731 systeme.
set(ligne, colonne)
732 += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
733 systeme.
set(ligne, colonne+1)
734 += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
735 systeme.
set(ligne, colonne+2)
736 += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
737 systeme.
set(ligne, colonne+3)
738 += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
739 sec_membre.
set(ligne)
740 -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
752 systeme.
set(ligne, colonne)
754 systeme.
set(ligne, colonne+1)
758 systeme.
set(ligne+1, colonne)
760 systeme.
set(ligne+1, colonne+1)
767 for (
int i=0; i<nr; i++) {
768 systeme.
set(ligne, colonne)
769 -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
770 systeme.
set(ligne, colonne+1)
771 -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
772 sec_membre.
set(ligne)
773 += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
779 for (
int i=0; i<nr; i++) {
780 systeme.
set(ligne, colonne)
781 -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
782 systeme.
set(ligne, colonne+1)
783 -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
784 sec_membre.
set(ligne)
785 += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
800 for (
int i=0 ; i<nr ; i++) {
801 cf_eta.
set(0, k, j, i) = 0.;
802 cf_vr.
set(0, k, j, i) = 0.;
805 for (
int zone=1 ; zone<nzm1 ; zone++) {
807 for (
int i=0 ; i<nr ; i++) {
808 cf_eta.
set(zone, k, j, i) =
809 sol_part_eta(zone, k, j, i)
810 +facteurs(conte)*sol_hom_un_eta(zone, k, j, i)
811 +facteurs(conte+1)*sol_hom_deux_eta(zone, k, j, i)
812 +facteurs(conte+2)*sol_hom_trois_eta(zone, k, j, i)
813 +facteurs(conte+3)*sol_hom_quatre_eta(zone, k, j, i) ;
814 cf_vr.
set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
815 +facteurs(conte)*sol_hom_un_vr(zone, k, j, i)
816 +facteurs(conte+1)*sol_hom_deux_vr(zone, k, j, i)
817 +facteurs(conte+2)*sol_hom_trois_vr(zone, k, j, i)
818 +facteurs(conte+3)*sol_hom_quatre_vr(zone, k, j, i) ;
823 for (
int i=0 ; i<nr ; i++) {
824 cf_eta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
825 +facteurs(conte)*sol_hom_deux_eta(nzm1, k, j, i)
826 +facteurs(conte+1)*sol_hom_quatre_eta(nzm1, k, j, i) ;
827 cf_vr.
set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
828 +facteurs(conte)*sol_hom_deux_vr(nzm1, k, j, i)
829 +facteurs(conte+1)*sol_hom_quatre_vr(nzm1, k, j, i) ;
849 resu.
set_vr_eta_mu(vr, het,
mu().sol_elliptic_boundary(param_mu, (*bound_mu), dir_mu , neum_mu)) ;