27char poisson_falloff_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_falloff.C,v 1.3 2014/10/13 08:53:29 j_novak Exp $" ;
57#include "type_parite.h"
82Mtbl_cf sol_poisson_falloff(
const Map_af& mapping,
const Mtbl_cf& source,
const int k_falloff)
88 assert (source.get_mg()->get_type_r(0) == RARE) ;
90 for (
int l=1 ; l<nz ; l++)
91 assert(source.get_mg()->get_type_r(l) == FIN) ;
97 const Base_val& base = source.base ;
103 double alpha, beta, echelle ;
104 int l_quant, m_quant;
115 Mtbl_cf solution_part(source.get_mg(), base) ;
116 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
117 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
118 Mtbl_cf resultat(source.get_mg(), base) ;
120 solution_part.set_etat_qcq() ;
121 solution_hom_un.set_etat_qcq() ;
122 solution_hom_deux.set_etat_qcq() ;
123 resultat.set_etat_qcq() ;
125 for (
int l=0 ; l<nz ; l++) {
126 solution_part.t[l]->set_etat_qcq() ;
127 solution_hom_un.t[l]->set_etat_qcq() ;
128 solution_hom_deux.t[l]->set_etat_qcq() ;
129 resultat.t[l]->set_etat_qcq() ;
130 for (
int k=0 ; k<source.get_mg()->get_np(l)+2 ; k++)
131 for (
int j=0 ; j<source.get_mg()->get_nt(l) ; j++)
132 for (
int i=0 ; i<source.get_mg()->get_nr(l) ; i++)
133 resultat.set(l, k, j, i) = 0 ;
143 nr = source.get_mg()->get_nr(0) ;
144 nt = source.get_mg()->get_nt(0) ;
145 np = source.get_mg()->get_np(0) ;
150 alpha = mapping.get_alpha()[0] ;
151 beta = mapping.get_beta()[0] ;
153 for (
int k=0 ; k<np+1 ; k++)
154 for (
int j=0 ; j<nt ; j++)
155 if (nullite_plm(j, nt, k, np, base) == 1)
158 donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
161 operateur =
new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
162 (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
165 nondege =
new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
168 sol_hom =
new Tbl(solh(nr, l_quant, 0., base_r)) ;
173 for (
int i=0 ; i<nr ; i++)
174 so->set(i) = source(0, k, j, i) ;
176 sol_part =
new Tbl(solp(*operateur, *nondege, alpha, beta,
181 for (
int i=0 ; i<nr ; i++) {
182 solution_part.set(0, k, j, i) = (*sol_part)(i) ;
183 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
184 solution_hom_deux.set(0, k, j, i) = 0. ;
202 for (
int zone=1 ; zone<nz ; zone++) {
203 nr = source.get_mg()->get_nr(zone) ;
204 nt = source.get_mg()->get_nt(zone) ;
205 np = source.get_mg()->get_np(zone) ;
207 if (np > np_max) np_max = np ;
208 if (nt > nt_max) nt_max = nt ;
210 alpha = mapping.get_alpha()[zone] ;
211 beta = mapping.get_beta()[zone] ;
213 for (
int k=0 ; k<np+1 ; k++)
214 for (
int j=0 ; j<nt ; j++)
215 if (nullite_plm(j, nt, k, np, base) == 1)
218 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
221 operateur =
new Matrice(laplacien_mat
222 (nr, l_quant, beta/alpha, 0, base_r)) ;
224 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
228 nondege =
new Matrice(prepa_nondege(*operateur, l_quant,
229 beta/alpha, 0, base_r)) ;
232 sol_hom =
new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
237 for (
int i=0 ; i<nr ; i++)
238 so->set(i) = source(zone, k, j, i) ;
240 sol_part =
new Tbl (solp(*operateur, *nondege, alpha,
241 beta, *so, 0, base_r)) ;
245 for (
int i=0 ; i<nr ; i++) {
246 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
247 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
248 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
265 int * indic =
new int [nz] ;
277 int sup_new, inf_new ;
280 for (
int k=0 ; k<np_max+1 ; k++)
281 for (
int j=0 ; j<nt_max ; j++)
282 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
289 for (
int zone=0 ; zone<nz ; zone++)
290 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
291 k, source.get_mg()->get_np(zone), base);
295 for (
int zone=1 ; zone<nz ; zone++)
296 taille+=2*indic[zone] ;
303 sec_membre =
new Tbl(taille) ;
304 sec_membre->set_etat_qcq() ;
305 for (
int l=0 ; l<taille ; l++)
306 sec_membre->set(l) = 0 ;
308 systeme =
new Matrice(taille, taille) ;
309 systeme->set_etat_qcq() ;
310 for (
int l=0 ; l<taille ; l++)
311 for (
int c=0 ; c<taille ; c++)
312 systeme->set(l, c) = 0 ;
318 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
326 nr = source.get_mg()->get_nr(0) ;
328 alpha = mapping.get_alpha()[0] ;
330 systeme->set(ligne, colonne) = 1. ;
332 for (
int i=0 ; i<nr ; i++)
333 sec_membre->set(ligne) -= solution_part(0, k, j, i) ;
341 systeme->set(ligne, colonne) = 1./alpha*l_quant ;
346 for (
int i=0 ; i<nr ; i++)
347 sec_membre->set(ligne) -=
349 *solution_part(0, k, j, i) ;
352 for (
int i=0 ; i<nr ; i++)
353 sec_membre->set(ligne) -=
354 (2*i+1)*(2*i+1)/alpha
355 *solution_part(0, k, j, i) ;
367 for (
int zone=1 ; zone<nz ; zone++)
368 if (indic[zone] == 1) {
370 nr = source.get_mg()->get_nr(zone) ;
371 alpha = mapping.get_alpha()[zone] ;
372 echelle = mapping.get_beta()[zone]/alpha ;
375 if (indic[zone-1] == 1) ligne -- ;
378 systeme->set(ligne, colonne) =
379 -
pow(echelle-1,
double(l_quant)) ;
382 systeme->set(ligne, colonne+1) =
383 -1/
pow(echelle-1,
double(l_quant+1)) ;
386 for (
int i=0 ; i<nr ; i++)
388 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
389 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
392 sup_new = colonne+1-ligne ;
393 if (sup_new > sup) sup = sup_new ;
400 if (indic[zone-1] == 1) {
402 systeme->set(ligne, colonne) =
403 -l_quant*
pow(echelle-1,
double(l_quant-1))/alpha ;
405 systeme->
set(ligne, colonne+1) =
406 (l_quant+1)/
pow(echelle-1,
double(l_quant+2))/alpha ;
411 for (
int i=0 ; i<nr ; i++)
413 sec_membre->
set(ligne) -=
414 i*i/alpha*solution_part(zone, k, j, i) ;
416 sec_membre->set(ligne) +=
417 i*i/alpha*solution_part(zone, k, j, i) ;
420 sup_new = colonne+1-ligne ;
421 if (sup_new > sup) sup = sup_new ;
431 systeme->set(ligne, colonne) =
432 pow(echelle+1,
double(l_quant)) ;
435 systeme->set(ligne, colonne+1) =
436 1/
pow(echelle+1,
double(l_quant+1)) ;
439 for (
int i=0 ; i<nr ; i++)
440 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ;
443 inf_new = ligne-colonne ;
444 if (inf_new > inf) inf = inf_new ;
450 if (indic[zone+1] == 1) {
453 systeme->set(ligne, colonne) =
454 l_quant*
pow(echelle+1,
double(l_quant-1))/alpha ;
457 systeme->
set(ligne, colonne+1) =
458 -(l_quant+1)/
pow(echelle+1,
double(l_quant+2))/alpha ;
461 for (
int i=0 ; i<nr ; i++)
462 sec_membre->
set(ligne) -=
463 i*i/alpha*solution_part(zone, k, j, i) ;
466 inf_new = ligne-colonne ;
467 if (inf_new > inf) inf = inf_new ;
482 systeme->set(ligne, colonne) =
483 double(k_falloff+l_quant)*
pow(echelle+1,
double(l_quant)) ;
486 systeme->set(ligne, colonne+1) =
487 double(k_falloff-l_quant-1)/
pow(echelle+1,
double(l_quant+1)) ;
491 for (
int i=0 ; i<nr ; i++)
492 sec_membre->set(ligne) -=
493 (k_falloff+(echelle+1)*i*i)*solution_part(zone, k, j, i) ;
496 inf_new = ligne-colonne ;
497 if (inf_new > inf) inf = inf_new ;
507 systeme->set_band(sup, inf) ;
510 Tbl facteurs(systeme->inverse(*sec_membre)) ;
517 nr = source.get_mg()->get_nr(0) ;
518 for (
int i=0 ; i<nr ; i++)
519 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
520 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
525 for (
int zone=1 ; zone<nz ; zone++)
526 if (indic[zone] == 1) {
527 nr = source.get_mg()->get_nr(zone) ;
528 for (
int i=0 ; i<nr ; i++)
529 resultat.set(zone, k, j, i) =
530 solution_part(zone, k, j, i)
531 +facteurs(conte)*solution_hom_un(zone, k, j, i)
532 +facteurs(conte+1)*solution_hom_deux(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_CHEBP
base de Cheb. paire (rare) seulement