24char poisson_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson.C,v 1.9 2014/10/13 08:53:29 j_novak Exp $" ;
154#include "type_parite.h"
179Mtbl_cf sol_poisson(
const Map_af& mapping,
const Mtbl_cf& source,
int dzpuis,
186 assert ((source.get_mg()->get_type_r(0) == RARE) || (source.get_mg()->get_type_r(0) == FIN)) ;
187 assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
188 for (
int l=1 ; l<nz-1 ; l++)
189 assert(source.get_mg()->get_type_r(l) == FIN) ;
191 assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3) || (dzpuis==5)) ;
192 assert ((!match) || (dzpuis != 5)) ;
195 const Base_val& base = source.base ;
201 double alpha, beta, echelle ;
202 int l_quant, m_quant;
207 Tbl *sol_part = 0x0 ;
208 Matrice *operateur = 0x0 ;
209 Matrice *nondege = 0x0 ;
213 Mtbl_cf solution_part(source.get_mg(), base) ;
214 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
215 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
216 Mtbl_cf resultat(source.get_mg(), base) ;
218 solution_part.set_etat_qcq() ;
219 solution_hom_un.set_etat_qcq() ;
220 solution_hom_deux.set_etat_qcq() ;
221 resultat.annule_hard() ;
222 for (
int l=0 ; l<nz ; l++) {
223 solution_part.t[l]->set_etat_qcq() ;
224 solution_hom_un.t[l]->set_etat_qcq() ;
225 solution_hom_deux.t[l]->set_etat_qcq() ;
234 nr = source.get_mg()->get_nr(0) ;
235 nt = source.get_mg()->get_nt(0) ;
236 np = source.get_mg()->get_np(0) ;
241 alpha = mapping.get_alpha()[0] ;
242 beta = mapping.get_beta()[0] ;
244 for (
int k=0 ; k<np+1 ; k++)
245 for (
int j=0 ; j<nt ; j++)
246 if (nullite_plm(j, nt, k, np, base) == 1)
249 donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
250 assert( (source.get_mg()->get_type_r(0) == RARE) ||
254 operateur =
new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
255 (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
258 nondege =
new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
262 sol_hom =
new Tbl(solh(nr, l_quant, 0., base_r)) ;
268 for (
int i=0 ; i<nr ; i++)
269 so->set(i) = source(0, k, j, i) ;
271 sol_part =
new Tbl(solp(*operateur, *nondege, alpha, beta,
276 for (
int i=0 ; i<nr ; i++) {
277 solution_part.set(0, k, j, i) = (*sol_part)(i) ;
280 solution_hom_un.set(0, k, j, i) = (*sol_hom)(0, i) ;
281 solution_hom_deux.set(0, k, j, i) = (*sol_hom)(1, i) ;
284 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
285 solution_hom_deux.set(0, k, j, i) = 0. ;
295 if (match)
delete sol_hom ;
304 nr = source.get_mg()->get_nr(nz-1) ;
305 nt = source.get_mg()->get_nt(nz-1) ;
306 np = source.get_mg()->get_np(nz-1) ;
308 if (np > np_max) np_max = np ;
309 if (nt > nt_max) nt_max = nt ;
311 alpha = mapping.get_alpha()[nz-1] ;
312 beta = mapping.get_beta()[nz-1] ;
314 for (
int k=0 ; k<np+1 ; k++)
315 for (
int j=0 ; j<nt ; j++)
316 if (nullite_plm(j, nt, k, np, base) == 1)
320 donne_lm(nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
323 operateur =
new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis,
325 (*operateur) = combinaison(*operateur, l_quant, 0., dzpuis, base_r) ;
327 nondege =
new Matrice(prepa_nondege(*operateur, l_quant, 0.,
331 sol_hom =
new Tbl(solh(nr, l_quant, 0., base_r)) ;
336 for (
int i=0 ; i<nr ; i++)
337 so->set(i) = source(nz-1, k, j, i) ;
338 sol_part =
new Tbl(solp(*operateur, *nondege, alpha, beta,
339 *so, dzpuis, base_r)) ;
342 for (
int i=0 ; i<nr ; i++) {
343 solution_part.set(nz-1, k, j, i) = (*sol_part)(i) ;
345 solution_hom_un.set(nz-1, k, j, i) = (*sol_hom)(i) ;
346 solution_hom_deux.set(nz-1, k, j, i) = 0. ;
353 if (match)
delete sol_hom ;
361 for (
int zone=1 ; zone<nz-1 ; zone++) {
362 nr = source.get_mg()->get_nr(zone) ;
363 nt = source.get_mg()->get_nt(zone) ;
364 np = source.get_mg()->get_np(zone) ;
366 if (np > np_max) np_max = np ;
367 if (nt > nt_max) nt_max = nt ;
369 alpha = mapping.get_alpha()[zone] ;
370 beta = mapping.get_beta()[zone] ;
372 for (
int k=0 ; k<np+1 ; k++)
373 for (
int j=0 ; j<nt ; j++)
374 if (nullite_plm(j, nt, k, np, base) == 1)
377 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
380 operateur =
new Matrice(laplacien_mat
381 (nr, l_quant, beta/alpha, 0, base_r)) ;
383 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
387 nondege =
new Matrice(prepa_nondege(*operateur, l_quant,
388 beta/alpha, 0, base_r)) ;
391 sol_hom =
new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
397 for (
int i=0 ; i<nr ; i++)
398 so->set(i) = source(zone, k, j, i) ;
400 sol_part =
new Tbl (solp(*operateur, *nondege, alpha,
401 beta, *so, 0, base_r)) ;
405 for (
int i=0 ; i<nr ; i++) {
406 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
408 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
409 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
417 if (match)
delete sol_hom ;
430 int * indic =
new int [nz] ;
442 int sup_new, inf_new ;
445 for (
int k=0 ; k<np_max+1 ; k++)
446 for (
int j=0 ; j<nt_max ; j++)
447 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
454 for (
int zone=0 ; zone<nz ; zone++)
455 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
456 k, source.get_mg()->get_np(zone), base);
459 taille = indic[nz-1]+indic[0] ;
460 for (
int zone=1 ; zone<nz-1 ; zone++)
461 taille+=2*indic[zone] ;
468 sec_membre =
new Tbl(taille) ;
469 sec_membre->set_etat_qcq() ;
470 for (
int l=0 ; l<taille ; l++)
471 sec_membre->set(l) = 0 ;
473 systeme =
new Matrice(taille, taille) ;
474 systeme->set_etat_qcq() ;
475 for (
int l=0 ; l<taille ; l++)
476 for (
int c=0 ; c<taille ; c++)
477 systeme->set(l, c) = 0 ;
483 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
491 nr = source.get_mg()->get_nr(0) ;
493 alpha = mapping.get_alpha()[0] ;
495 systeme->set(ligne, colonne) = 1. ;
497 for (
int i=0 ; i<nr ; i++)
498 sec_membre->set(ligne) -= solution_part(0, k, j, i) ;
510 systeme->set(ligne, colonne) = 1./alpha*l_quant ;
515 for (
int i=0 ; i<nr ; i++)
516 sec_membre->set(ligne) -=
518 *solution_part(0, k, j, i) ;
521 for (
int i=0 ; i<nr ; i++)
522 sec_membre->set(ligne) -=
523 (2*i+1)*(2*i+1)/alpha
524 *solution_part(0, k, j, i) ;
536 for (
int zone=1 ; zone<nz-1 ; zone++)
537 if (indic[zone] == 1) {
539 nr = source.get_mg()->get_nr(zone) ;
540 alpha = mapping.get_alpha()[zone] ;
541 echelle = mapping.get_beta()[zone]/alpha ;
544 if (indic[zone-1] == 1) ligne -- ;
547 systeme->set(ligne, colonne) =
548 -
pow(echelle-1,
double(l_quant)) ;
551 systeme->set(ligne, colonne+1) =
552 -1/
pow(echelle-1,
double(l_quant+1)) ;
555 for (
int i=0 ; i<nr ; i++)
557 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
558 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
561 sup_new = colonne+1-ligne ;
562 if (sup_new > sup) sup = sup_new ;
569 if (indic[zone-1] == 1) {
571 systeme->set(ligne, colonne) =
572 -l_quant*
pow(echelle-1,
double(l_quant-1))/alpha ;
574 systeme->
set(ligne, colonne+1) =
575 (l_quant+1)/
pow(echelle-1,
double(l_quant+2))/alpha ;
580 for (
int i=0 ; i<nr ; i++)
582 sec_membre->
set(ligne) -=
583 i*i/alpha*solution_part(zone, k, j, i) ;
585 sec_membre->set(ligne) +=
586 i*i/alpha*solution_part(zone, k, j, i) ;
589 sup_new = colonne+1-ligne ;
590 if (sup_new > sup) sup = sup_new ;
597 systeme->set(ligne, colonne) =
598 pow(echelle+1,
double(l_quant)) ;
601 systeme->set(ligne, colonne+1) =
602 1/
pow(echelle+1,
double(l_quant+1)) ;
605 for (
int i=0 ; i<nr ; i++)
606 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
609 inf_new = ligne-colonne ;
610 if (inf_new > inf) inf = inf_new ;
616 if (indic[zone+1] == 1) {
619 systeme->set(ligne, colonne) =
620 l_quant*
pow(echelle+1,
double(l_quant-1))/alpha ;
623 systeme->
set(ligne, colonne+1) =
624 -(l_quant+1)/
pow(echelle+1,
double(l_quant+2))/alpha ;
627 for (
int i=0 ; i<nr ; i++)
628 sec_membre->
set(ligne) -=
629 i*i/alpha*solution_part(zone, k, j, i) ;
632 inf_new = ligne-colonne ;
633 if (inf_new > inf) inf = inf_new ;
644 if (indic[nz-1] == 1) {
646 nr = source.get_mg()->get_nr(nz-1) ;
649 alpha = mapping.get_alpha()[nz-1] ;
651 if (indic[nz-2] == 1) ligne -- ;
654 systeme->set(ligne, colonne) = -
pow(-2,
double(l_quant+1)) ;
656 for (
int i=0 ; i<nr ; i++)
658 sec_membre->set(ligne) += solution_part(nz-1, k, j, i) ;
659 else sec_membre->set(ligne) -= solution_part(nz-1, k, j, i) ;
662 if (indic[nz-2] == 1) {
665 systeme->set(ligne+1, colonne) =
666 alpha*(l_quant+1)*
pow(-2.,
double(l_quant+2)) ;
669 for (
int i=0 ; i<nr ; i++)
671 sec_membre->
set(ligne+1) -=
672 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
674 sec_membre->set(ligne+1) +=
675 -4*alpha*i*i*solution_part(nz-1, k, j, i) ;
678 if (sup == 0) sup = 1 ;
686 systeme->set_band(sup, inf) ;
689 Tbl facteurs(systeme->inverse(*sec_membre)) ;
696 nr = source.get_mg()->get_nr(0) ;
697 for (
int i=0 ; i<nr ; i++)
698 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
699 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
704 for (
int zone=1 ; zone<nz-1 ; zone++)
705 if (indic[zone] == 1) {
706 nr = source.get_mg()->get_nr(zone) ;
707 for (
int i=0 ; i<nr ; i++)
708 resultat.set(zone, k, j, i) =
709 solution_part(zone, k, j, i)
710 +facteurs(conte)*solution_hom_un(zone, k, j, i)
711 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
716 if (indic[nz-1] == 1) {
717 nr = source.get_mg()->get_nr(nz-1) ;
718 for (
int i=0 ; i<nr ; i++)
719 resultat.set(nz-1, k, j, i) =
720 solution_part(nz-1, k, j, i)
721 +facteurs(conte)*solution_hom_un(nz-1, k, j, i) ;
734 for (
int zone = 0; zone < nz; zone++)
735 for (
int k=0 ; k<source.get_mg()->get_np(zone)+1 ; k++)
736 for (
int j=0 ; j<source.get_mg()->get_nt(zone) ; j++)
737 if (nullite_plm(j,source.get_mg()->get_nt(zone) ,
738 k, source.get_mg()->get_np(zone), base) == 1)
739 for (
int i=0; i<source.get_mg()->get_nr(zone); i++)
740 resultat.set(zone, k, j, i) = solution_part(zone, k, j, i) ;
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 .
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEBP
base de Cheb. paire (rare) seulement