23char poisson_frontiere_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_frontiere.C,v 1.5 2014/10/13 08:53:29 j_novak Exp $" ;
76#include "type_parite.h"
77#include "utilitaires.h"
108Mtbl_cf sol_poisson_frontiere(
const Map_af& mapping,
const Mtbl_cf& source,
109 const Mtbl_cf& limite,
int type_raccord,
int num_front,
int dzpuis,
double fact_dir,
double fact_neu)
116 assert ((num_front>=0) && (num_front<nz-2)) ;
117 assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
118 for (
int l=num_front+1 ; l<nz-1 ; l++)
119 assert(source.get_mg()->get_type_r(l) == FIN) ;
121 assert (source.get_etat() != ETATNONDEF) ;
122 assert (limite.get_etat() != ETATNONDEF) ;
124 assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3)) ;
125 assert ((type_raccord == 1) || (type_raccord == 2)|| (type_raccord == 3)) ;
128 const Base_val& base = source.base ;
133 double alpha, beta, echelle ;
134 int l_quant, m_quant;
145 Mtbl_cf solution_part(source.get_mg(), base) ;
146 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
147 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
148 Mtbl_cf resultat(source.get_mg(), base) ;
150 solution_part.set_etat_qcq() ;
151 solution_hom_un.set_etat_qcq() ;
152 solution_hom_deux.set_etat_qcq() ;
153 resultat.set_etat_qcq() ;
155 for (
int l=0 ; l<nz ; l++) {
156 solution_part.t[l]->set_etat_qcq() ;
157 solution_hom_un.t[l]->set_etat_qcq() ;
158 solution_hom_deux.t[l]->set_etat_qcq() ;
159 resultat.t[l]->set_etat_qcq() ;
160 for (
int k=0 ; k<source.get_mg()->get_np(l)+2 ; k++)
161 for (
int j=0 ; j<source.get_mg()->get_nt(l) ; j++)
162 for (
int i=0 ; i<source.get_mg()->get_nr(l) ; i++)
163 resultat.set(l, k, j, i) = 0 ;
175 nr = source.get_mg()->get_nr(nz-1) ;
176 nt = source.get_mg()->get_nt(nz-1) ;
177 np = source.get_mg()->get_np(nz-1) ;
179 if (np > np_max) np_max = np ;
180 if (nt > nt_max) nt_max = nt ;
182 alpha = mapping.get_alpha()[nz-1] ;
183 beta = mapping.get_beta()[nz-1] ;
185 for (
int k=0 ; k<np+1 ; k++)
186 for (
int j=0 ; j<nt ; j++)
187 if (nullite_plm(j, nt, k, np, base) == 1)
191 donne_lm(nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
194 operateur =
new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis,
196 (*operateur) = combinaison(*operateur, l_quant, 0., dzpuis, base_r) ;
199 nondege =
new Matrice(prepa_nondege(*operateur, l_quant, 0.,
203 sol_hom =
new Tbl(solh(nr, l_quant, 0., base_r)) ;
208 for (
int i=0 ; i<nr ; i++)
209 so->set(i) = source(nz-1, k, j, i) ;
210 sol_part =
new Tbl(solp(*operateur, *nondege, alpha, beta,
211 *so, dzpuis, base_r)) ;
214 for (
int i=0 ; i<nr ; i++) {
215 solution_part.set(nz-1, k, j, i) = (*sol_part)(i) ;
216 solution_hom_un.set(nz-1, k, j, i) = (*sol_hom)(i) ;
217 solution_hom_deux.set(nz-1, k, j, i) = 0. ;
231 for (
int zone=num_front+1 ; zone<nz-1 ; zone++) {
232 nr = source.get_mg()->get_nr(zone) ;
233 nt = source.get_mg()->get_nt(zone) ;
234 np = source.get_mg()->get_np(zone) ;
236 if (np > np_max) np_max = np ;
237 if (nt > nt_max) nt_max = nt ;
239 alpha = mapping.get_alpha()[zone] ;
240 beta = mapping.get_beta()[zone] ;
242 for (
int k=0 ; k<np+1 ; k++)
243 for (
int j=0 ; j<nt ; j++)
244 if (nullite_plm(j, nt, k, np, base) == 1)
247 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
250 operateur =
new Matrice(laplacien_mat
251 (nr, l_quant, beta/alpha, 0, base_r)) ;
253 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
257 nondege =
new Matrice(prepa_nondege(*operateur, l_quant,
258 beta/alpha, 0, base_r)) ;
261 sol_hom =
new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
266 for (
int i=0 ; i<nr ; i++)
267 so->set(i) = source(zone, k, j, i) ;
269 sol_part =
new Tbl (solp(*operateur, *nondege, alpha,
270 beta, *so, 0, base_r)) ;
274 for (
int i=0 ; i<nr ; i++) {
275 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
276 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
277 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
293 nr = source.get_mg()->get_nr(num_front+1) ;
294 nt = source.get_mg()->get_nt(num_front+1) ;
295 np = source.get_mg()->get_np(num_front+1) ;
299 alpha = mapping.get_alpha()[num_front+1] ;
300 beta = mapping.get_beta()[num_front+1] ;
301 echelle = beta/alpha ;
303 for (
int k=0 ; k<np+1 ; k++)
304 for (
int j=0 ; j<nt ; j++)
305 if (nullite_plm(j, nt, k, np, base) == 1)
308 donne_lm(nz, num_front+1, j, k, base, m_quant, l_quant, base_r) ;
310 switch (type_raccord) {
315 for (
int i=0 ; i<nr ; i++)
317 somme += solution_part(num_front+1, k, j, i) ;
319 somme -= solution_part(num_front+1, k, j, i) ;
320 facteur = (limite(num_front, k, j, 0)-somme)
321 *
pow(echelle-1, l_quant+1) ;
323 for (
int i=0 ; i<nr ; i++)
324 solution_part.
set(num_front+1, k, j, i) +=
325 facteur*solution_hom_deux(num_front+1, k, j, i) ;
328 facteur = -
pow(echelle-1, 2*l_quant+1) ;
329 for (
int i=0 ; i<nr ; i++)
330 solution_hom_un.set(num_front+1, k, j, i) +=
331 facteur*solution_hom_deux(num_front+1, k, j, i) ;
339 for (
int i=0 ; i<nr ; i++)
341 somme -= i*i/alpha*solution_part(num_front+1, k, j, i) ;
343 somme += i*i/alpha*solution_part(num_front+1, k, j, i) ;
344 facteur = (somme-limite(num_front, k, j, 0))
345 * alpha*
pow(echelle-1, l_quant+2)/(l_quant+1) ;
346 for (
int i=0 ; i<nr ; i++)
347 solution_part.
set(num_front+1, k, j, i) +=
348 facteur*solution_hom_deux(num_front+1, k, j, i) ;
351 facteur =
pow(echelle-1, 2*l_quant+1)*l_quant/(l_quant+1) ;
352 for (
int i=0 ; i<nr ; i++)
353 solution_hom_un.
set(num_front+1, k, j, i) +=
354 facteur*solution_hom_deux(num_front+1, k, j, i) ;
360 for (
int i=0 ; i<nr ; i++)
362 somme += solution_part(num_front+1, k, j, i) *
363 fact_dir - fact_neu *
364 i*i/alpha*solution_part(num_front+1, k, j, i) ;
366 somme += - solution_part(num_front+1, k, j, i) *
367 fact_dir + fact_neu *
368 i*i/alpha*solution_part(num_front+1, k, j, i) ;
371 somme2 = fact_dir *
pow(echelle-1, -l_quant-1) -
372 fact_neu/alpha*
pow(echelle-1, -l_quant-2)*(l_quant+1) ;
374 facteur = (limite(num_front, k, j, 0)- somme) / somme2 ;
376 for (
int i=0 ; i<nr ; i++)
377 solution_part.
set(num_front+1, k, j, i) +=
378 facteur*solution_hom_deux(num_front+1, k, j, i) ;
382 somme1 = fact_dir *
pow(echelle-1, l_quant) +
383 fact_neu / alpha *
pow(echelle-1, l_quant-1) *
385 facteur = - somme1 / somme2 ;
386 for (
int i=0 ; i<nr ; i++)
387 solution_hom_un.
set(num_front+1, k, j, i) +=
388 facteur*solution_hom_deux(num_front+1, k, j, i) ;
393 cout <<
"Diantres nous ne devrions pas etre ici ! " << endl ;
399 for (
int i=0 ; i<nr ; i++)
400 solution_hom_deux.set(num_front+1, k, j, i) = 0 ;
411 int * indic =
new int [nz] ;
423 int sup_new, inf_new ;
426 for (
int k=0 ; k<np_max+1 ; k++)
427 for (
int j=0 ; j<nt_max ; j++)
428 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
435 for (
int zone=num_front+1 ; zone<nz ; zone++)
436 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
437 k, source.get_mg()->get_np(zone), base);
440 taille = indic[nz-1]+indic[num_front+1] ;
441 for (
int zone=num_front+2 ; zone<nz-1 ; zone++)
442 taille+=2*indic[zone] ;
449 sec_membre =
new Tbl(taille) ;
450 sec_membre->set_etat_qcq() ;
451 for (
int l=0 ; l<taille ; l++)
452 sec_membre->set(l) = 0 ;
454 systeme =
new Matrice(taille, taille) ;
455 systeme->set_etat_qcq() ;
456 for (
int l=0 ; l<taille ; l++)
457 for (
int c=0 ; c<taille ; c++)
458 systeme->set(l, c) = 0 ;
464 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
471 if (indic[num_front+1] == 1) {
472 nr = source.get_mg()->get_nr(num_front+1) ;
474 alpha = mapping.get_alpha()[num_front+1] ;
475 beta = mapping.get_beta()[num_front+1] ;
476 echelle = beta/alpha ;
480 for (
int i=0 ; i<nr ; i++)
481 somme += solution_hom_un (num_front+1, k, j, i) ;
482 systeme->set(ligne, colonne) = somme ;
485 for (
int i=0 ; i<nr ; i++)
486 sec_membre->set(ligne) -= solution_part(num_front+1, k, j, i) ;
492 if (indic[num_front+1] == 1) {
496 for (
int i=0 ; i<nr ; i++)
498 solution_hom_un(num_front+1, k, j, i) ;
499 systeme->set(ligne, colonne) = somme ;
502 for (
int i=0 ; i<nr ; i++)
503 sec_membre->set(ligne) -= i*i/alpha
504 *solution_part(num_front+1, k, j, i) ;
516 for (
int zone=num_front+2 ; zone<nz-1 ; zone++)
517 if (indic[zone] == 1) {
519 nr = source.get_mg()->get_nr(zone) ;
520 alpha = mapping.get_alpha()[zone] ;
521 echelle = mapping.get_beta()[zone]/alpha ;
524 if (indic[zone-1] == 1) ligne -- ;
527 systeme->set(ligne, colonne) =
528 -
pow(echelle-1,
double(l_quant)) ;
531 systeme->set(ligne, colonne+1) =
532 -1/
pow(echelle-1,
double(l_quant+1)) ;
535 for (
int i=0 ; i<nr ; i++)
537 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
538 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
541 sup_new = colonne+1-ligne ;
542 if (sup_new > sup) sup = sup_new ;
549 if (indic[zone-1] == 1) {
551 systeme->set(ligne, colonne) =
552 -l_quant*
pow(echelle-1,
double(l_quant-1))/alpha ;
554 systeme->
set(ligne, colonne+1) =
555 (l_quant+1)/
pow(echelle-1,
double(l_quant+2))/alpha ;
560 for (
int i=0 ; i<nr ; i++)
562 sec_membre->
set(ligne) -=
563 i*i/alpha*solution_part(zone, k, j, i) ;
565 sec_membre->set(ligne) +=
566 i*i/alpha*solution_part(zone, k, j, i) ;
569 sup_new = colonne+1-ligne ;
570 if (sup_new > sup) sup = sup_new ;
577 systeme->set(ligne, colonne) =
578 pow(echelle+1,
double(l_quant)) ;
581 systeme->set(ligne, colonne+1) =
582 1/
pow(echelle+1,
double(l_quant+1)) ;
585 for (
int i=0 ; i<nr ; i++)
586 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
589 inf_new = ligne-colonne ;
590 if (inf_new > inf) inf = inf_new ;
596 if (indic[zone+1] == 1) {
599 systeme->set(ligne, colonne) =
600 l_quant*
pow(echelle+1,
double(l_quant-1))/alpha ;
603 systeme->
set(ligne, colonne+1) =
604 -(l_quant+1)/
pow(echelle+1,
double(l_quant+2))/alpha ;
607 for (
int i=0 ; i<nr ; i++)
608 sec_membre->
set(ligne) -=
609 i*i/alpha*solution_part(zone, k, j, i) ;
612 inf_new = ligne-colonne ;
613 if (inf_new > inf) inf = inf_new ;
624 if (indic[nz-1] == 1) {
626 nr = source.get_mg()->get_nr(nz-1) ;
629 alpha = mapping.get_alpha()[nz-1] ;
631 if (indic[nz-2] == 1) ligne -- ;
634 systeme->set(ligne, colonne) = -
pow(-2,
double(l_quant+1)) ;
636 for (
int i=0 ; i<nr ; i++)
638 sec_membre->set(ligne) += solution_part(nz-1, k, j, i) ;
639 else sec_membre->set(ligne) -= solution_part(nz-1, k, j, i) ;
642 if (indic[nz-2] == 1) {
645 systeme->set(ligne+1, colonne) =
646 alpha*(l_quant+1)*
pow(-2.,
double(l_quant+2)) ;
649 for (
int i=0 ; i<nr ; i++)
651 sec_membre->
set(ligne+1) -=
652 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
654 sec_membre->set(ligne+1) +=
655 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
658 if (sup == 0) sup = 1 ;
666 systeme->set_band(sup, inf) ;
669 Tbl facteurs(systeme->inverse(*sec_membre)) ;
675 if (indic[num_front+1] == 1) {
676 nr = source.get_mg()->get_nr(num_front+1) ;
677 for (
int i=0 ; i<nr ; i++)
678 resultat.set(num_front+1, k, j, i) =
679 solution_part(num_front+1, k, j, i)
680 +facteurs(conte)*solution_hom_un(num_front+1, k, j, i) ;
685 for (
int zone=num_front+2 ; zone<nz-1 ; zone++)
686 if (indic[zone] == 1) {
687 nr = source.get_mg()->get_nr(zone) ;
688 for (
int i=0 ; i<nr ; i++)
689 resultat.set(zone, k, j, i) =
690 solution_part(zone, k, j, i)
691 +facteurs(conte)*solution_hom_un(zone, k, j, i)
692 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
697 if (indic[nz-1] == 1) {
698 nr = source.get_mg()->get_nr(nz-1) ;
699 for (
int i=0 ; i<nr ; i++)
700 resultat.set(nz-1, k, j, i) =
701 solution_part(nz-1, k, j, i)
702 +facteurs(conte)*solution_hom_un(nz-1, k, j, i) ;
714 for (
int l=0 ; l<num_front+1 ; l++)
715 resultat.t[l]->set_etat_zero() ;
Tbl & set(int l)
Read/write of the value in a given domain.
int get_nzone() const
Returns the number of domains.
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Cmp pow(const Cmp &, int)
Power .