25char poisson_compact_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_compact.C,v 1.6 2014/10/13 08:53:29 j_novak Exp $" ;
72#include "type_parite.h"
75#include "utilitaires.h"
90Mtbl_cf sol_poisson_compact(
const Mtbl_cf& source,
double a,
double b,
94 assert (source.get_etat() != ETATNONDEF) ;
100 const int nmax = 200 ;
101 static Matrice* tab_op[nmax] ;
102 static int nb_deja_fait = 0 ;
103 static int l_deja_fait[nmax] ;
104 static int n_deja_fait[nmax] ;
107 for (
int i=0 ; i<nb_deja_fait ; i++)
112 int nz = source.get_mg()->get_nzone() ;
115 int nr = source.get_mg()->get_nr(0) ;
116 int nt = source.get_mg()->get_nt(0) ;
117 int np = source.get_mg()->get_np(0) ;
124 Mtbl_cf solution(source.get_mg(), source.base) ;
125 solution.set_etat_qcq() ;
126 solution.t[0]->set_etat_qcq() ;
128 for (
int k=0 ; k<np+1 ; k++)
129 for (
int j=0 ; j<nt ; j++)
130 if (nullite_plm(j, nt, k, np, source.base) == 1)
133 donne_lm(nz, 0, j, k, source.base, m_quant, l_quant, base_r) ;
138 for (
int i=0 ; i<nr ; i++)
139 solution.set(0, k, j, i) = 0 ;
148 int taille = (l_quant == 1) ? nr : nr-1 ;
150 Matrice operateur (taille, taille) ;
151 for (
int conte=0 ; conte<nb_deja_fait ; conte++)
152 if ((l_deja_fait[conte]== l_quant) && (n_deja_fait[conte] == nr))
156 if (nb_deja_fait >= nmax) {
157 cout <<
"sol_poisson_compact : trop de matrices ..." << endl;
161 operateur = a*lap_cpt_mat (nr, l_quant, base_r)
162 + b*xdsdx_mat(nr, l_quant, base_r) ;
163 operateur = combinaison_cpt (operateur, l_quant, base_r) ;
165 l_deja_fait[nb_deja_fait] = l_quant ;
166 n_deja_fait[nb_deja_fait] = nr ;
167 tab_op[nb_deja_fait] =
new Matrice(operateur) ;
173 operateur = *tab_op[indice] ;
179 for (
int i=0 ; i<taille ; i++)
180 so.set(i) = source(0, k, j, i) ;
181 so = combinaison_cpt (so, base_r) ;
183 Tbl part (operateur.inverse(so)) ;
186 for (
int i=0 ; i<nr ; i++)
187 solution.set(0, k, j, i) = part(i) ;
189 solution.set(0, k, j, nr-1) = 0 ;
190 for (
int i=nr-2 ; i>=0 ; i--)
192 solution.set(0, k, j, i) = part(i) ;
193 solution.set(0, k, j, i+1) += part(i) ;
196 solution.set(0, k, j, i) = part(i)*(2*i+3) ;
197 solution.set(0, k, j, i+1) += part(i)*(2*i+1) ;
203 for (
int i=0 ; i<nr ; i++)
204 solution.set(0, k, j, i) = 0 ;
207 for (
int j=0; j<nt; j++)
208 for(
int i=0 ; i<nr ; i++)
209 solution.set(0, np+1, j, i) = 0 ;
211 for (
int zone=1 ; zone<nz ; zone++)
212 solution.t[zone]->set_etat_zero() ;
224Mtbl_cf sol_poisson_compact(
const Map_af& mapping,
const Mtbl_cf& source,
225 const Tbl& ac,
const Tbl& bc,
bool ) {
228 int nzet = ac.get_dim(1) ;
229 assert(nzet<=source.get_mg()->get_nzone()) ;
232 assert (source.get_mg()->get_type_r(0) == RARE) ;
233 for (
int l=1 ; l<nzet ; l++)
234 assert(source.get_mg()->get_type_r(l) == FIN) ;
237 const Base_val& base = source.base ;
240 Mtbl_cf resultat(source.get_mg(), base) ;
241 resultat.annule_hard() ;
246 double alpha, beta, echelle ;
247 int l_quant, m_quant;
252 for (
int l=0 ; l<nzet ; l++) {
253 nr = mapping.get_mg()->get_nr(l) ;
255 if (nr > max_nr) max_nr = nr ;
258 Matrice systeme (size, size) ;
259 systeme.set_etat_qcq() ;
260 Tbl sec_membre (size) ;
262 soluce.set_etat_qcq() ;
264 np = mapping.get_mg()->get_np(0) ;
265 nt = mapping.get_mg()->get_nt(0) ;
269 for (
int k=0 ; k<np+1 ; k++)
270 for (
int j=0 ; j<nt ; j++)
271 if (nullite_plm(j, nt, k, np, base) == 1) {
273 systeme.annule_hard() ;
274 sec_membre.annule_hard() ;
276 int column_courant = 0 ;
277 int ligne_courant = 0 ;
283 nr = mapping.get_mg()->get_nr(0) ;
284 alpha = mapping.get_alpha()[0] ;
285 base.give_quant_numbers (0, k, j, m_quant, l_quant, base_r) ;
288 for (
int i=0 ; i<size ; i++)
293 Diff_dsdx2 d2_n(base_r, nr) ;
294 Diff_sxdsdx sxd_n(base_r, nr) ;
295 Diff_sx2 sx2_n(base_r, nr) ;
296 Diff_x2dsdx2 x2d2_n(base_r,nr) ;
297 Diff_xdsdx xd_n(base_r,nr) ;
298 Diff_id id_n(base_r,nr) ;
300 work =
new Matrice( ac(0,0) * ( d2_n + 2.*sxd_n -l_quant*(l_quant+1)*sx2_n )
301 + ac(0,2) * ( x2d2_n + 2.*xd_n -l_quant*(l_quant+1)*id_n )
302 + alpha * bc(0,1) * xd_n ) ;
313 for (
int col=0 ; col<nr ; col++)
315 systeme.set(ligne_courant, col+column_courant) = 1 ;
317 systeme.set(ligne_courant, col+column_courant) = -1 ;
321 for (
int col=0 ; col<nr ; col++)
323 systeme.set(ligne_courant, col+column_courant) = 2*col+1 ;
325 systeme.set(ligne_courant, col+column_courant) = -(2*col+1) ;
331 for (
int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
332 for (
int col=0 ; col<nr ; col++)
333 systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
334 sec_membre.set(lig+ligne_courant) = alpha*alpha*source(0, k, j, lig) ;
338 ligne_courant += nr-1-nbr_cl ;
341 for (
int col=0 ; col<nr ; col++) {
343 systeme.set(ligne_courant, col+column_courant) = 1 ;
346 systeme.set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
348 systeme.set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
350 column_courant += nr ;
356 for (
int l=1 ; l<nzet ; l++) {
358 nr = mapping.get_mg()->get_nr(l) ;
359 alpha = mapping.get_alpha()[l] ;
360 beta = mapping.get_beta()[l] ;
361 echelle = beta/alpha ;
362 double bsa = echelle ;
363 double bsa2 = bsa*bsa ;
365 base.give_quant_numbers (l, k, j, m_quant, l_quant, base_r) ;
367 Diff_dsdx dx(base_r, nr) ;
368 Diff_xdsdx xdx(base_r, nr) ;
369 Diff_x2dsdx x2dx(base_r, nr) ;
370 Diff_x3dsdx x3dx(base_r, nr) ;
371 Diff_dsdx2 dx2(base_r, nr) ;
372 Diff_xdsdx2 xdx2(base_r, nr) ;
373 Diff_x2dsdx2 x2dx2(base_r, nr) ;
374 Diff_x3dsdx2 x3dx2(base_r, nr) ;
375 Diff_id id(base_r,nr) ;
376 Diff_mx mx(base_r,nr) ;
378 work =
new Matrice ( ac(l,0) * ( x2dx2 + 2.*bsa*xdx2 + bsa2*dx2 + 2.*xdx
379 + 2.*bsa*dx - l_quant*(l_quant+1)*
id )
380 + ac(l,1) * ( x3dx2 + 2.*bsa*x2dx2 + bsa2*xdx2 + 2.*x2dx
381 + 2.*bsa*xdx - l_quant*(l_quant+1) *mx )
382 + alpha * ( bc(l,0) * ( x2dx + 2.*bsa*xdx + bsa2*dx )
383 + bc(l,1) * ( x3dx + 2.*bsa*x2dx + bsa2*xdx ) ) ) ;
386 for (
int col=0 ; col<nr ; col++) {
389 systeme.set(ligne_courant, col+column_courant) = -1 ;
391 systeme.set(ligne_courant, col+column_courant) = 1 ;
394 systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
396 systeme.set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
402 source_aux.set_etat_qcq() ;
403 for (
int i=0 ; i<nr ; i++)
404 source_aux.set(i) = source(l,k,j,i)*alpha*alpha ;
405 Tbl xso(source_aux) ;
406 Tbl xxso(source_aux) ;
407 multx_1d(nr, &xso.t,
R_CHEB) ;
408 multx_1d(nr, &xxso.t,
R_CHEB) ;
409 multx_1d(nr, &xxso.t,
R_CHEB) ;
410 source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
414 for (
int lig=0 ; lig<nr-1 ; lig++) {
415 for (
int col=0 ; col<nr ; col++)
416 systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
417 sec_membre.set(lig+ligne_courant) = source_aux(lig) ;
424 ligne_courant += nr-2 ;
428 for (
int col=0 ; col<nr ; col++) {
430 systeme.set(ligne_courant, col+column_courant) = 1 ;
432 systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
436 column_courant += nr ;
445 systeme.set_band (size, size) ;
451 soluce = systeme.inverse(sec_membre) ;
460 for (
int l=0 ; l<nzet ; l++) {
461 nr = mapping.get_mg()->get_nr(l) ;
462 for (
int i=0 ; i<nr ; i++) {
463 resultat.set(l,k,j,i) = soluce(conte) ;
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement