23char dalembert_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/dalembert.C,v 1.13 2014/10/13 08:53:28 j_novak Exp $" ;
113Mtbl_cf sol_dalembert(Param& par,
const Map_af& mapping,
const Mtbl_cf& source)
118 bool ced = (source.get_mg()->get_type_r(nz-1) == UNSURR ) ;
119 int nz0 = (ced ? nz - 1 : nz ) ;
120 assert ((source.get_mg()->get_type_r(0) == RARE)||(source.get_mg()->get_type_r(0) == FIN)) ;
121 for (
int l=1 ; l<nz0 ; l++) {
122 assert(source.get_mg()->get_type_r(l) == FIN) ;
123 assert(source.get_mg()->get_nt(l) == source.get_mg()->get_nt(0)) ;
124 assert(source.get_mg()->get_np(l) == source.get_mg()->get_np(0)) ;
126 assert (par.get_n_double() > 0) ;
127 assert (par.get_n_tbl_mod() > 1) ;
132 if (par.get_n_int() > 1) {
134 l_min = par.get_int(1) ;
138 const Base_val& base = source.base ;
142 int base_r, type_dal ;
144 int l_quant, m_quant;
145 nt = source.get_mg()->get_nt(0) ;
146 np = source.get_mg()->get_np(0) ;
155 Mtbl_cf solution_part(source.get_mg(), base) ;
156 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
157 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
158 Mtbl_cf resultat(source.get_mg(), base) ;
160 solution_part.set_etat_qcq() ;
161 solution_hom_un.set_etat_qcq() ;
162 solution_hom_deux.set_etat_qcq() ;
163 resultat.annule_hard() ;
166 double* bc1 = &par.get_double_mod(1) ;
167 double* bc2 = &par.get_double_mod(2) ;
168 Tbl* tbc3 = &par.get_tbl_mod(1) ;
170 for (
int l=0 ; l<nz ; l++) {
171 solution_part.t[l]->annule_hard() ;
172 solution_hom_un.t[l]->annule_hard() ;
173 solution_hom_deux.t[l]->annule_hard() ;
180 nr = source.get_mg()->get_nr(lz) ;
183 alpha = mapping.get_alpha()[lz] ;
185 for (
int k=0 ; k<np+1 ; k++) {
186 for (
int j=0 ; j<nt ; j++) {
188 base.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
189 assert( (source.get_mg()->get_type_r(0) == RARE) ||
192 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_quant >=l_min) )
195 par.get_tbl_mod().set(4,lz) = 2*par.get_tbl_mod()(2,lz) ;
196 par.get_tbl_mod().set(5,lz) = 2*par.get_tbl_mod()(3,lz) ;
197 par.get_tbl_mod().set(6,lz) = 2*par.get_tbl_mod()(1,lz) ;
198 par.get_tbl_mod().set(7,lz) =
199 -l_quant*(l_quant+1)*par.get_tbl_mod()(3,lz) ;
200 par.get_tbl_mod().set(8,lz) =
201 -l_quant*(l_quant+1)*par.get_tbl_mod()(2,lz) ;
202 par.get_tbl_mod().set(9,lz) =
203 -l_quant*(l_quant+1)*par.get_tbl_mod()(1,lz) ;
205 Matrice operateur(nr,nr) ;
207 get_operateur_dal(par, lz, base_r, type_dal, operateur) ;
211 for (
int i=0 ; i<nr ; i++)
212 so->set(i) = source(lz, k, j, i) ;
216 sol_part =
new Tbl(dal_inverse(base_r, type_dal, operateur,
220 sol_hom =
new Tbl(dal_inverse(base_r, type_dal, operateur,
224 for (
int i=0 ; i<nr ; i++) {
225 solution_part.set(lz, k, j, i) = (*sol_part)(i) ;
226 solution_hom_un.set(lz, k, j, i) = (*sol_hom)(i) ;
227 solution_hom_deux.set(lz, k, j, i) = 0. ;
234 double part, dpart, hom, dhom;
235 Tbl der_part(3,1,nr) ;
236 der_part.set_etat_qcq() ;
237 for (
int i=0; i<nr; i++)
238 der_part.set(0,0,i) = (*sol_part)(i) ;
239 Tbl der_hom(3,1,nr) ;
240 der_hom.set_etat_qcq() ;
241 for (
int i=0; i<nr; i++)
242 der_hom.set(0,0,i) = (*sol_hom)(i) ;
245 som_r_chebp(sol_part->t, nr, 1, 1, 1., &part) ;
246 _dsdx_r_chebp(&der_part, base_pipo) ;
247 som_r_chebi(der_part.t, nr, 1, 1, 1., &dpart) ;
248 som_r_chebp(sol_hom->t, nr, 1, 1, 1., &hom) ;
249 _dsdx_r_chebp(&der_hom, base_pipo) ;
250 som_r_chebi(der_hom.t, nr, 1, 1, 1., &dhom) ;
253 som_r_chebi(sol_part->t, nr, 1, 1, 1., &part) ;
254 _dsdx_r_chebi(&der_part, base_pipo) ;
255 som_r_chebp(der_part.t, nr, 1, 1, 1., &dpart) ;
256 som_r_chebi(sol_hom->t, nr, 1, 1, 1., &hom) ;
257 _dsdx_r_chebi(&der_hom, base_pipo) ;
258 som_r_chebp(der_hom.t, nr, 1, 1, 1., &dhom) ;
261 part = part*(*bc1) + dpart*(*bc2)/alpha ;
262 hom = hom*(*bc1) + dhom*(*bc2)/alpha ;
263 double lambda = ((*tbc3)(k,j) - part) / hom ;
264 for (
int i=0 ; i<nr ; i++)
265 resultat.set(lz, k, j, i) =
266 solution_part(lz, k, j, i)
267 +lambda*solution_hom_un(lz, k, j, i) ;
280 for (lz=1 ; lz<nz0 ; lz++) {
281 nr = source.get_mg()->get_nr(lz) ;
283 alpha = mapping.get_alpha()[lz] ;
284 beta = mapping.get_beta()[lz] ;
286 for (
int k=0 ; k<np+1 ; k++)
287 for (
int j=0 ; j<nt ; j++) {
289 base.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
291 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_quant >=l_min) )
294 par.get_tbl_mod().set(4,lz) = 2*par.get_tbl_mod()(2,lz) ;
295 par.get_tbl_mod().set(5,lz) = 2*par.get_tbl_mod()(3,lz) ;
296 par.get_tbl_mod().set(6,lz) = 2*par.get_tbl_mod()(1,lz) ;
297 par.get_tbl_mod().set(7,lz) =
298 -l_quant*(l_quant+1)*par.get_tbl_mod()(3,lz) ;
299 par.get_tbl_mod().set(8,lz) =
300 -l_quant*(l_quant+1)*par.get_tbl_mod()(2,lz) ;
301 par.get_tbl_mod().set(9,lz) =
302 -l_quant*(l_quant+1)*par.get_tbl_mod()(1,lz) ;
304 Matrice operateur(nr,nr) ;
306 get_operateur_dal(par, lz, base_r, type_dal, operateur) ;
310 for (
int i=0; i<nr; i++) so->set(i) = 0. ;
312 sol_hom =
new Tbl(dal_inverse(base_r, type_dal, operateur, *so,
316 sol_hom2 =
new Tbl(dal_inverse(base_r, type_dal, operateur, *so,
320 double *tmp =
new double[nr] ;
321 for (
int i=0 ; i<nr ; i++)
322 tmp[i] = source(lz, k, j, i) ;
324 for (
int i=0; i<nr; i++) so->set(i) = beta*tmp[i] ;
325 multx_1d(nr, &tmp,
R_CHEB) ;
326 for (
int i=0; i<nr; i++) so->set(i) += alpha*tmp[i] ;
329 for (
int i=0; i<nr; i++) so->set(i) = beta*beta*tmp[i] ;
330 multx_1d(nr, &tmp,
R_CHEB) ;
331 for (
int i=0; i<nr; i++) so->set(i) += 2*alpha*beta*tmp[i] ;
332 multx_1d(nr, &tmp,
R_CHEB) ;
333 for (
int i=0; i<nr; i++) so->set(i) += alpha*alpha*tmp[i] ;
338 sol_part =
new Tbl (dal_inverse(base_r, type_dal, operateur,
341 for (
int i=0 ; i<nr ; i++) {
342 solution_part.set(lz, k, j, i) = (*sol_part)(i) ;
343 solution_hom_un.set(lz, k, j, i) = (*sol_hom)(i) ;
344 solution_hom_deux.set(lz, k, j, i) = (*sol_hom2)(i) ;
362 int taille = 2*nz0 - 1 ;
364 deuz.set_etat_qcq() ;
365 Matrice systeme(taille,taille) ;
366 systeme.set_etat_qcq() ;
368 int inf = (nz0>2) ? 2 : 1 ;
369 for (
int k=0; k<np+1; k++) {
370 for (
int j=0; j<nt; j++) {
372 base.give_quant_numbers(0, k, j, m_quant, l_quant, base_r) ;
373 if ( (nullite_plm(j, nt, k, np, base)) && (l_quant + dl >= l_min) ) {
375 int parite = (base_r ==
R_CHEBP) ? 0 : 1 ;
378 for (l=0; l<taille; l++)
379 for (c=0; c<taille; c++) systeme.set(l,c) = xx ;
380 for (l=0; l<taille; l++) deuz.set(l) = xx ;
385 nr = source.get_mg()->get_nr(0) ;
386 alpha = mapping.get_alpha()[0] ;
388 for (
int i=0; i<nr; i++)
389 systeme.set(l,c) += solution_hom_un(0, k, j, i) ;
390 for (
int i=0; i<nr; i++) deuz.set(l) -= solution_part(0, k, j, i) ;
394 for (
int i=0; i<nr; i++)
395 xx +=(2*i+parite)*(2*i+parite)
396 *solution_hom_un(0, k, j, i) ;
397 systeme.set(l,c) += xx/alpha ;
399 for (
int i=0; i<nr; i++) xx -= (2*i+parite)*
400 (2*i+parite)*solution_part(0, k, j, i) ;
401 deuz.set(l) += xx/alpha ;
406 for (lz=1; lz<nz0; lz++) {
407 nr = source.get_mg()->get_nr(lz) ;
408 alpha = mapping.get_alpha()[lz] ;
411 for (
int i=0; i<nr; i++)
413 systeme.set(l,c) -= solution_hom_un(lz, k, j, i) ;
415 systeme.set(l,c) += solution_hom_un(lz, k, j, i) ;
417 for (
int i=0; i<nr; i++)
419 systeme.set(l,c) -= solution_hom_deux(lz, k, j, i) ;
421 systeme.set(l,c) += solution_hom_deux(lz, k, j, i) ;
422 for (
int i=0; i<nr; i++)
423 if (i%2 == 0) deuz.set(l) += solution_part(lz, k, j, i) ;
424 else deuz.set(l) -= solution_part(lz, k, j, i) ;
428 for (
int i=0; i<nr; i++)
430 xx += i*i*solution_hom_un(lz, k, j, i) ;
432 xx -= i*i*solution_hom_un(lz, k, j, i) ;
433 systeme.set(l,c) += xx/alpha ;
436 for (
int i=0; i<nr; i++)
438 xx += i*i*solution_hom_deux(lz, k, j, i) ;
440 xx -= i*i*solution_hom_deux(lz, k, j, i) ;
441 systeme.set(l,c) += xx/alpha ;
443 for (
int i=0; i<nr; i++)
444 if (i%2 == 0) xx -= i*i*solution_part(lz, k, j, i) ;
445 else xx += i*i*solution_part(lz, k, j, i) ;
446 deuz.set(l) += xx/alpha ;
450 for (
int i=0; i<nr; i++)
452 ((*bc1)+(*bc2)*i*i/alpha)*solution_hom_un(lz, k, j, i) ;
454 for (
int i=0; i<nr; i++)
456 ((*bc1)+(*bc2)*i*i/alpha)*solution_hom_deux(lz, k, j, i) ;
457 for (
int i=0; i<nr; i++)
459 ((*bc1)+(*bc2)*i*i/alpha)*solution_part(lz, k, j, i) ;
460 deuz.set(l) += (*tbc3)(k,j) ;
463 for (
int i=0; i<nr; i++)
464 systeme.set(l,c) += solution_hom_un(lz, k, j, i) ;
466 for (
int i=0; i<nr; i++)
467 systeme.set(l,c) += solution_hom_deux(lz, k, j, i) ;
468 for (
int i=0; i<nr; i++)
469 deuz.set(l) -= solution_part(lz, k, j, i) ;
472 for (
int i=0; i<nr; i++) xx += i*i*solution_hom_un(lz, k, j, i) ;
473 systeme.set(l,c) += xx/alpha ;
476 for (
int i=0; i<nr; i++)
477 xx += i*i*solution_hom_deux(lz, k, j, i) ;
478 systeme.set(l,c) += xx/alpha ;
480 for (
int i=0; i<nr; i++)
481 xx -= i*i*solution_part(lz, k, j, i) ;
482 deuz.set(l) += xx/alpha ;
490 systeme.set_band(sup, inf) ;
492 Tbl facteur(systeme.inverse(deuz)) ;
495 nr = source.get_mg()->get_nr(0) ;
496 for (
int i=0; i<nr; i++)
497 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
498 + facteur(0)*solution_hom_un(0, k, j, i) ;
501 for (lz=1; lz<nz0; lz++) {
502 nr = source.get_mg()->get_nr(lz) ;
503 for (
int i=0; i<nr; i++)
504 resultat.set(lz, k, j, i) = solution_part(lz, k, j, i)
505 + facteur(2*lz-1)*solution_hom_un(lz, k, j, i)
506 + facteur(2*lz)*solution_hom_deux(lz, k, j, i) ;
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 O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define ORDRE1_LARGE
Operateur du premier ordre .
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement