23char poisson_tau_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_tau.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $" ;
67#include "type_parite.h"
91Mtbl_cf sol_poisson_tau(
const Map_af& mapping,
const Mtbl_cf& source,
int dzpuis)
97 assert ((source.get_mg()->get_type_r(0) == RARE) || (source.get_mg()->get_type_r(0) == FIN)) ;
98 assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
99 for (
int l=1 ; l<nz-1 ; l++)
100 assert(source.get_mg()->get_type_r(l) == FIN) ;
102 assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3)) ;
105 const Base_val& base = source.base ;
108 Mtbl_cf resultat(source.get_mg(), base) ;
109 resultat.annule_hard() ;
114 double alpha, beta, echelle ;
115 int l_quant, m_quant;
120 for (
int l=0 ; l<nz ; l++) {
121 nr = mapping.get_mg()->get_nr(l) ;
127 Matrice systeme (size, size) ;
128 systeme.set_etat_qcq() ;
129 Tbl sec_membre (size) ;
131 np = mapping.get_mg()->get_np(0) ;
132 nt = mapping.get_mg()->get_nt(0) ;
136 for (
int k=0 ; k<np+1 ; k++)
137 for (
int j=0 ; j<nt ; j++)
138 if (nullite_plm(j, nt, k, np, base) == 1) {
143 systeme.annule_hard() ;
144 sec_membre.annule_hard() ;
146 int column_courant = 0 ;
147 int ligne_courant = 0 ;
153 nr = mapping.get_mg()->get_nr(0) ;
154 alpha = mapping.get_alpha()[0] ;
155 base.give_quant_numbers (0, k, j, m_quant, l_quant, base_r) ;
156 work =
new Matrice (laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
160 if (source.get_mg()->get_type_r(0) == RARE) {
166 for (
int col=0 ; col<nr ; col++)
168 systeme.set(ligne_courant, col+column_courant) = 1 ;
170 systeme.set(ligne_courant, col+column_courant) = -1 ;
174 for (
int col=0 ; col<nr ; col++)
176 systeme.set(ligne_courant, col+column_courant) = 2*col+1 ;
178 systeme.set(ligne_courant, col+column_courant) = -(2*col+1) ;
189 for (
int col=0 ; col<nr ; col++) {
190 systeme.set(ligne_courant, col+column_courant) = col*(col+1)*(col+2)*(col+3)/
double(12)*(2*(col%2)-1);
193 else if (l_quant == 1) {
195 for (
int col=0 ; col<nr ; col++) {
196 systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
201 for (
int col=0 ; col<nr ; col++) {
202 systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
203 systeme.set(ligne_courant+1, col+column_courant) = col*(col+1)*(col+2)*(col+3)/
double(12)*(2*(col%2)-1) ;
207 ligne_courant += nbr_cl ;
210 for (
int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
211 for (
int col=0 ; col<nr ; col++)
212 systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
213 sec_membre.set(lig+ligne_courant) = alpha*alpha*source(0, k, j, lig) ;
217 ligne_courant += nr-1-nbr_cl ;
220 for (
int col=0 ; col<nr ; col++) {
221 if (source.get_mg()->get_type_r(0) == RARE) {
223 systeme.set(ligne_courant, col+column_courant) = 1 ;
226 systeme.set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
229 systeme.set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
234 systeme.set(ligne_courant, col+column_courant) = 1 ;
236 systeme.set(ligne_courant+1, col+column_courant) = col*(col+3)/
double(2)/alpha ;
240 column_courant += nr ;
248 for (
int l=1 ; l<nz-1 ; l++) {
250 nr = mapping.get_mg()->get_nr(l) ;
251 alpha = mapping.get_alpha()[l] ;
252 beta = mapping.get_beta()[l] ;
253 echelle = beta/alpha ;
255 base.give_quant_numbers (l, k, j, m_quant, l_quant, base_r) ;
256 work =
new Matrice (laplacien_mat(nr, l_quant, echelle, 0, base_r)) ;
259 for (
int col=0 ; col<nr ; col++) {
262 systeme.set(ligne_courant, col+column_courant) = -1 ;
264 systeme.set(ligne_courant, col+column_courant) = 1 ;
267 systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
269 systeme.set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
277 source_aux.set_etat_qcq() ;
278 for (
int i=0 ; i<nr ; i++)
279 source_aux.set(i) = source(l,k,j,i)*alpha*alpha ;
280 Tbl xso(source_aux) ;
281 Tbl xxso(source_aux) ;
282 multx_1d(nr, &xso.t,
R_CHEB) ;
283 multx_1d(nr, &xxso.t,
R_CHEB) ;
284 multx_1d(nr, &xxso.t,
R_CHEB) ;
285 source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
287 for (
int lig=0 ; lig<nr-2 ; lig++) {
288 for (
int col=0 ; col<nr ; col++)
289 systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
290 sec_membre.set(lig+ligne_courant) = source_aux(lig) ;
294 ligne_courant += nr-2 ;
296 for (
int col=0 ; col<nr ; col++) {
298 systeme.set(ligne_courant, col+column_courant) = 1 ;
300 systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
303 column_courant += nr ;
310 nr = mapping.get_mg()->get_nr(nz-1) ;
311 alpha = mapping.get_alpha()[nz-1] ;
312 beta = mapping.get_beta()[nz-1] ;
314 base.give_quant_numbers (nz-1, k, j, m_quant, l_quant, base_r) ;
315 work =
new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis, base_r)) ;
318 for (
int col=0 ; col<nr ; col++) {
321 systeme.set(ligne_courant, col+column_courant) = -1 ;
323 systeme.set(ligne_courant, col+column_courant) = 1 ;
326 systeme.set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
328 systeme.set(ligne_courant+1, col+column_courant) = 4*alpha*col*col ;
339 for (
int col=0 ; col<nr ; col++)
340 systeme.set(ligne_courant, col+column_courant) = 1 ;
345 for (
int col=0 ; col<nr ; col++)
346 systeme.set(ligne_courant, col+column_courant) = 1 ;
348 for (
int col=0 ; col<nr ; col++)
349 systeme.set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
356 for (
int col=0 ; col<nr ; col++)
357 systeme.set(ligne_courant, col+column_courant) = 1 ;
364 for (
int col=0 ; col<nr ; col++)
365 systeme.set(ligne_courant, col+column_courant) = 1 ;
369 cout <<
"Unknown dzpuis in sol_poisson_tau ..." << endl ;
373 ligne_courant += nbr_cl ;
379 indic = alpha*alpha ;
389 for (
int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
390 for (
int col=0 ; col<nr ; col++)
391 systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
392 sec_membre.set(lig+ligne_courant) = indic*source(nz-1, k, j, lig) ;
397 systeme.set_band (max_nr, max_nr) ;
399 Tbl soluce (systeme.inverse(sec_membre)) ;
403 for (
int l=0 ; l<nz ; l++) {
404 nr = mapping.get_mg()->get_nr(l) ;
405 for (
int i=0 ; i<nr ; i++) {
406 resultat.set(l,k,j,i) = soluce(conte) ;
int get_nzone() const
Returns the number of domains.
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEB
base de Chebychev ordinaire (fin)