23char poisson_vect_frontiere_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_vect_frontiere.C,v 1.9 2014/10/13 08:53:30 j_novak Exp $" ;
85void poisson_vect_frontiere (
double lambda,
const Tenseur& source, Tenseur& shift,
86 const Valeur& lim_x,
const Valeur& lim_y,
const Valeur& lim_z,
87 int num_front,
double precision,
int itermax) {
92 int nt = lim_x.get_mg()->get_nt(num_front+1) ;
93 int np = lim_x.get_mg()->get_np(num_front+1) ;
94 int nz = lim_x.get_mg()->get_nzone() ;
96 if (shift.get_etat() == ETATZERO) {
97 shift.set_etat_qcq() ;
98 for (
int i=0 ; i<3 ; i++)
99 shift.set(i).annule_hard() ;
100 shift.set_std_base() ;
103 Tenseur so (source) ;
106 Tenseur cop_so (so) ;
107 cop_so.dec2_dzpuis() ;
108 cop_so.dec2_dzpuis() ;
110 Tenseur scal (*so.get_mp()) ;
111 scal.set_etat_qcq() ;
113 Cmp source_scal (
contract(cop_so.gradient(), 0, 1)()/(lambda+1)) ;
114 source_scal.inc2_dzpuis() ;
115 if (source_scal.get_etat()== ETATZERO) {
116 source_scal.annule_hard() ;
117 source_scal.std_base_scal() ;
118 source_scal.set_dzpuis(4) ;
121 Tenseur copie_so (so) ;
122 copie_so.dec_dzpuis() ;
124 Tenseur source_vect (*so.get_mp(), 1, CON, *source.get_triad()) ;
125 Tenseur auxi (*so.get_mp(), 1, COV, *source.get_triad()) ;
126 Cmp grad_shift (source_scal.get_mp()) ;
129 Valeur lim_scal (lim_x.get_mg()) ;
130 Tenseur shift_old (*shift.get_mp(), 1, CON, shift.get_mp()->get_bvect_cart()) ;
139 grad_shift =
contract(shift.gradient(), 0, 1)() ;
141 grad_shift.va.coef_i() ;
146 for (
int k=0 ; k<np ; k++)
147 for (
int j=0 ; j<nt ; j++)
148 lim_scal.set(num_front, k, j, 0) =
149 grad_shift.va (num_front+1, k, j, 0) ;
150 lim_scal.std_base_scal() ;
153 scal.set() = source_scal.poisson_dirichlet (lim_scal, num_front) ;
156 source_vect.set_etat_qcq() ;
157 auxi = scal.gradient() ;
159 for (
int i=0 ; i<3 ; i++)
160 source_vect.set(i) = copie_so(i) - lambda * auxi(i) ;
163 for (
int i=0 ; i<3 ; i++)
164 if (source_vect(i).get_etat()==ETATQCQ)
167 for (
int i=0 ; i<3 ; i++)
168 source_vect.set(i).annule_hard() ;
169 source_vect.set_std_base() ;
173 shift.set(0) = source_vect(0).poisson_dirichlet (lim_x, num_front) ;
174 shift.set(1) = source_vect(1).poisson_dirichlet (lim_y, num_front) ;
175 shift.set(2) = source_vect(2).poisson_dirichlet (lim_z, num_front) ;
178 for (
int i=0 ; i<3 ; i++)
179 if (
max(
norme(shift(i))) > precision) {
180 Tbl diff (
diffrelmax (shift(i), shift_old(i))) ;
181 for (
int j=num_front+1 ; j<nz ; j++)
186 cout <<
"Pas " << conte <<
" : Difference " << erreur << endl ;
189 if ((erreur <precision) || (conte > itermax))
197void poisson_vect_boundary (
double lambda,
const Vector& source,Vector& shift,
198 const Valeur& lim_x,
const Valeur& lim_y,
const Valeur& lim_z,
199 int num_front,
double precision,
int itermax) {
202 assert(source.get_mp().get_bvect_spher() == *(source.get_triad())) ;
203 assert(source.get_mp().get_bvect_spher() == *(shift.get_triad())) ;
207 int nt = lim_x.get_mg()->get_nt(num_front+1) ;
208 int np = lim_x.get_mg()->get_np(num_front+1) ;
209 int nz = lim_x.get_mg()->get_nzone() ;
211 Metric_flat ff(source.get_mp(), source.get_mp().get_bvect_spher()) ;
217 cop_so.dec_dzpuis(2) ;
218 cop_so.dec_dzpuis(2) ;
220 Scalar scal (so.get_mp()) ;
222 Scalar source_scal (
contract(cop_so.derive_cov(ff), 0, 1)/(lambda+1)) ;
223 source_scal.inc_dzpuis(2) ;
224 if (source_scal.get_etat()== ETATZERO) {
225 source_scal.annule_hard() ;
226 source_scal.std_spectral_base() ;
227 source_scal.set_dzpuis(4) ;
230 Vector copie_so (so) ;
231 copie_so.dec_dzpuis() ;
233 Vector source_vect (so.get_mp(), CON, *source.get_triad()) ;
234 Vector auxi (so.get_mp(), COV, *source.get_triad()) ;
235 Scalar grad_shift (source_scal.get_mp()) ;
238 Valeur lim_scal (lim_x.get_mg()) ;
239 Vector shift_old (shift.get_mp(), CON, shift.get_mp().get_bvect_cart()) ;
248 grad_shift =
contract(shift.derive_cov(ff), 0, 1) ;
249 grad_shift.dec_dzpuis(2) ;
250 grad_shift.set_spectral_va().coef_i() ;
254 for (
int k=0 ; k<np ; k++)
255 for (
int j=0 ; j<nt ; j++)
256 lim_scal.set(num_front, k, j, 0) =
257 grad_shift.get_spectral_va() (num_front+1, k, j, 0) ;
258 lim_scal.std_base_scal() ;
262 source_scal.filtre(4) ;
263 scal = source_scal.poisson_dirichlet (lim_scal, num_front) ;
266 source_vect.set_etat_qcq() ;
267 auxi = scal.derive_cov(ff) ;
269 for (
int i=1 ; i<=3 ; i++)
270 source_vect.set(i) = copie_so(i) - lambda * auxi(i) ;
273 for (
int i=1 ; i<=3 ; i++)
274 if (source_vect(i).get_etat()==ETATQCQ)
277 for (
int i=1 ; i<=3 ; i++)
278 source_vect.set(i).annule_hard() ;
279 source_vect.std_spectral_base() ;
282 shift.change_triad(source.get_mp().get_bvect_cart()) ;
283 source_vect.change_triad(source.get_mp().get_bvect_cart()) ;
285 for (
int i=1 ; i<=3 ; i++)
286 source_vect.set(i).filtre(4) ;
289 shift.set(1) = source_vect(1).poisson_dirichlet (lim_x, num_front) ;
290 shift.set(2) = source_vect(2).poisson_dirichlet (lim_y, num_front) ;
291 shift.set(3) = source_vect(3).poisson_dirichlet (lim_z, num_front) ;
293 shift.change_triad(source.get_mp().get_bvect_spher()) ;
294 source_vect.change_triad(source.get_mp().get_bvect_spher()) ;
297 for (
int i=1 ; i<=3 ; i++)
298 if (
max(
norme(shift(i))) > precision) {
299 Tbl diff (
diffrelmax (shift(i), shift_old(i))) ;
300 for (
int j=num_front+1 ; j<nz ; j++)
305 cout <<
"Pas " << conte <<
" : Difference " << erreur << endl ;
308 if ((erreur <precision) || (conte > itermax))
316void poisson_vect_binaire (
double lambda,
317 const Tenseur& source_un,
const Tenseur& source_deux,
318 const Valeur& bound_x_un,
const Valeur& bound_y_un,
319 const Valeur& bound_z_un,
const Valeur& bound_x_deux,
320 const Valeur& bound_y_deux,
const Valeur& bound_z_deux,
321 Tenseur& sol_un, Tenseur& sol_deux,
int num_front,
double precision) {
324 assert (sol_un.get_etat() != ETATNONDEF) ;
325 assert (sol_deux.get_etat() != ETATNONDEF) ;
329 assert (sol_un.get_mp() == source_un.get_mp()) ;
330 assert (sol_deux.get_mp() == source_deux.get_mp()) ;
332 double orientation_un = sol_un.get_mp()->get_rot_phi() ;
333 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
335 double orientation_deux = sol_deux.get_mp()->get_rot_phi() ;
336 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
338 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
341 if (sol_un.get_etat() == ETATZERO) {
342 sol_un.set_etat_qcq() ;
343 for (
int i=0 ; i<3 ; i++)
344 sol_un.set(i).annule_hard() ;
345 sol_un.set_std_base() ;
348 if (sol_deux.get_etat() == ETATZERO) {
349 sol_deux.set_etat_qcq() ;
350 for (
int i=0 ; i<3 ; i++)
351 sol_deux.set(i).annule_hard() ;
352 sol_deux.set_std_base() ;
355 Valeur limite_x_un (bound_x_un.get_mg()) ;
356 limite_x_un = bound_x_un ;
357 Valeur limite_y_un (bound_y_un.get_mg()) ;
358 limite_y_un = bound_y_un ;
359 Valeur limite_z_un (bound_z_un.get_mg()) ;
360 limite_z_un = bound_z_un ;
362 Valeur limite_x_deux (bound_x_deux.get_mg()) ;
363 limite_x_deux = bound_x_deux ;
364 Valeur limite_y_deux (bound_y_deux.get_mg()) ;
365 limite_y_deux = bound_y_deux ;
366 Valeur limite_z_deux (bound_z_deux.get_mg()) ;
367 limite_z_deux = bound_z_deux ;
369 Valeur limite_chi_un (bound_x_un.get_mg()) ;
371 limite_chi_un.std_base_scal() ;
373 Valeur limite_chi_deux (bound_x_deux.get_mg()) ;
374 limite_chi_deux = 0 ;
375 limite_chi_deux.std_base_scal() ;
377 Mtbl xa_mtbl_un (source_un.get_mp()->get_mg()) ;
378 xa_mtbl_un.set_etat_qcq() ;
379 Mtbl ya_mtbl_un (source_un.get_mp()->get_mg()) ;
380 ya_mtbl_un.set_etat_qcq() ;
381 Mtbl za_mtbl_un (source_un.get_mp()->get_mg()) ;
382 za_mtbl_un.set_etat_qcq() ;
383 Mtbl xa_mtbl_deux (source_deux.get_mp()->get_mg()) ;
384 xa_mtbl_deux.set_etat_qcq() ;
385 Mtbl ya_mtbl_deux (source_deux.get_mp()->get_mg()) ;
386 ya_mtbl_deux.set_etat_qcq() ;
387 Mtbl za_mtbl_deux (source_deux.get_mp()->get_mg()) ;
388 za_mtbl_deux.set_etat_qcq() ;
390 xa_mtbl_un = source_un.get_mp()->xa ;
391 ya_mtbl_un = source_un.get_mp()->ya ;
392 za_mtbl_un = source_un.get_mp()->za ;
394 xa_mtbl_deux = source_deux.get_mp()->xa ;
395 ya_mtbl_deux = source_deux.get_mp()->ya ;
396 za_mtbl_deux = source_deux.get_mp()->za ;
398 double xabs, yabs, zabs ;
399 double air, theta, phi ;
402 int nbrep_un = bound_x_un.get_mg()->get_np(num_front) ;
403 int nbret_un = bound_x_un.get_mg()->get_nt(num_front) ;
404 int nbrep_deux = bound_x_deux.get_mg()->get_np(num_front) ;
405 int nbret_deux = bound_x_deux.get_mg()->get_nt(num_front) ;
406 int nz_un = bound_x_un.get_mg()->get_nzone() ;
407 int nz_deux = bound_x_deux.get_mg()->get_nzone() ;
410 Tenseur cop_so_un (source_un) ;
411 cop_so_un.dec2_dzpuis() ;
412 cop_so_un.dec2_dzpuis() ;
414 Cmp source_scal_un (
contract (cop_so_un.gradient(), 0, 1)()/(lambda+1)) ;
415 if (source_scal_un.get_etat() == ETATZERO) {
416 source_scal_un.annule_hard() ;
417 source_scal_un.std_base_scal() ;
419 source_scal_un.inc2_dzpuis() ;
422 Tenseur cop_so_deux (source_deux) ;
423 cop_so_deux.dec2_dzpuis() ;
424 cop_so_deux.dec2_dzpuis() ;
426 Cmp source_scal_deux (
contract (cop_so_deux.gradient(), 0, 1)()/(lambda+1)) ;
427 if (source_scal_deux.get_etat() == ETATZERO) {
428 source_scal_deux.annule_hard() ;
429 source_scal_deux.std_base_scal() ;
431 source_scal_deux.inc2_dzpuis() ;
434 Tenseur copie_so_un (source_un) ;
435 copie_so_un.dec_dzpuis() ;
437 Tenseur copie_so_deux (source_deux) ;
438 copie_so_deux.dec_dzpuis() ;
441 Tenseur sol_un_old (sol_un) ;
442 Tenseur sol_deux_old (sol_deux) ;
450 Tenseur chi_un (source_scal_un.poisson_dirichlet (limite_chi_un, num_front)) ;
451 Tenseur chi_deux (source_scal_deux.poisson_dirichlet (limite_chi_deux, num_front)) ;
454 Tenseur source_vect_un (copie_so_un) ;
455 if (source_vect_un.get_etat() == ETATZERO) {
456 source_vect_un.set_etat_qcq() ;
457 for (
int i=0 ; i<3 ; i++) {
458 source_vect_un.set(i).annule_hard() ;
459 source_vect_un.set(i).set_dzpuis(3) ;
461 source_vect_un.set_std_base() ;
463 Tenseur grad_chi_un (chi_un.gradient()) ;
464 grad_chi_un.inc_dzpuis() ;
465 for (
int i=0 ; i<3 ; i++)
466 source_vect_un.set(i) = source_vect_un(i)-lambda*grad_chi_un(i) ;
468 Tenseur source_vect_deux (copie_so_deux) ;
469 if (source_vect_deux.get_etat() == ETATZERO) {
470 source_vect_deux.set_etat_qcq() ;
471 for (
int i=0 ; i<3 ; i++) {
472 source_vect_deux.set(i).annule_hard() ;
473 source_vect_deux.set(i).set_dzpuis(3) ;
475 source_vect_deux.set_std_base() ;
477 Tenseur grad_chi_deux (chi_deux.gradient()) ;
478 grad_chi_deux.inc_dzpuis() ;
479 for (
int i=0 ; i<3 ; i++)
480 source_vect_deux.set(i) = source_vect_deux(i)-lambda*grad_chi_deux(i) ;
483 sol_un_old = sol_un ;
484 sol_deux_old = sol_deux ;
487 sol_un.set(0) = source_vect_un(0).poisson_dirichlet (limite_x_un, num_front) ;
488 sol_un.set(1) = source_vect_un(1).poisson_dirichlet (limite_y_un, num_front) ;
489 sol_un.set(2) = source_vect_un(2).poisson_dirichlet (limite_z_un, num_front) ;
490 sol_deux.set(0) = source_vect_deux(0).poisson_dirichlet (limite_x_deux, num_front) ;
491 sol_deux.set(1) = source_vect_deux(1).poisson_dirichlet (limite_y_deux, num_front) ;
492 sol_deux.set(2) = source_vect_deux(2).poisson_dirichlet (limite_z_deux, num_front) ;
496 Cmp div_shift_un (
contract(sol_un.gradient(), 0, 1)()) ;
497 div_shift_un.dec2_dzpuis() ;
498 div_shift_un.va.coef_i() ;
501 for (
int k=0 ; k<nbrep_un ; k++)
502 for (
int j=0 ; j<nbret_un ; j++)
503 limite_chi_un.set(num_front, k, j, 0) =
504 div_shift_un.va (num_front+1, k, j, 0) ;
505 limite_chi_un.std_base_scal() ;
507 Cmp div_shift_deux (
contract(sol_deux.gradient(), 0, 1)()) ;
508 div_shift_deux.dec2_dzpuis() ;
509 div_shift_deux.va.coef_i() ;
511 limite_chi_deux = 1 ;
512 for (
int k=0 ; k<nbrep_deux ; k++)
513 for (
int j=0 ; j<nbret_deux ; j++)
514 limite_chi_deux.set(num_front, k, j, 0) =
515 div_shift_deux.va (num_front+1, k, j, 0) ;
516 limite_chi_deux.std_base_scal() ;
520 for (
int k=0 ; k<nbrep_un ; k++)
521 for (
int j=0 ; j<nbret_un ; j++) {
522 xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
523 yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
524 zabs = za_mtbl_un (num_front+1, k, j, 0) ;
526 source_deux.get_mp()->convert_absolute
527 (xabs, yabs, zabs, air, theta, phi) ;
529 valeur = sol_deux(0).val_point(air, theta, phi) ;
530 limite_x_un.set(num_front, k, j, 0) =
531 bound_x_un(num_front, k, j, 0) - same_orient*valeur ;
533 valeur = sol_deux(1).val_point(air, theta, phi) ;
534 limite_y_un.set(num_front, k, j, 0) =
535 bound_y_un(num_front, k, j, 0) - same_orient*valeur ;
537 valeur = sol_deux(2).val_point(air, theta, phi) ;
538 limite_z_un.set(num_front, k, j, 0) =
539 bound_z_un(num_front, k, j, 0) - valeur ;
543 for (
int k=0 ; k<nbrep_deux ; k++)
544 for (
int j=0 ; j<nbret_deux ; j++) {
545 xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
546 yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
547 zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
549 source_un.get_mp()->convert_absolute
550 (xabs, yabs, zabs, air, theta, phi) ;
552 valeur = sol_un(0).val_point(air, theta, phi) ;
553 limite_x_deux.set(num_front, k, j, 0) =
554 bound_x_deux(num_front, k, j, 0) - same_orient*valeur ;
556 valeur = sol_un(1).val_point(air, theta, phi) ;
557 limite_y_deux.set(num_front, k, j, 0) =
558 bound_y_deux(num_front, k, j, 0) - same_orient*valeur ;
560 valeur = sol_un(2).val_point(air, theta, phi) ;
561 limite_z_deux.set(num_front, k, j, 0) =
562 bound_z_deux(num_front, k, j, 0) - valeur ;
567 for (
int i=0 ; i<3 ; i++) {
568 Tbl diff_un (
diffrelmax (sol_un_old(i), sol_un(i))) ;
569 for (
int j=num_front+1 ; j<nz_un ; j++)
570 if (erreur<diff_un(j))
571 erreur = diff_un(j) ;
574 for (
int i=0 ; i<3 ; i++) {
575 Tbl diff_deux (
diffrelmax (sol_deux_old(i), sol_deux(i))) ;
576 for (
int j=num_front+1 ; j<nz_deux ; j++)
577 if (erreur<diff_deux(j))
578 erreur = diff_deux(j) ;
581 cout <<
"Pas " << conte <<
" : Difference " << erreur << endl ;
583 if (erreur < precision)
588void poisson_vect_binaire (
double lambda,
589 const Vector& source_un,
const Vector& source_deux,
590 const Valeur& bound_x_un,
const Valeur& bound_y_un,
591 const Valeur& bound_z_un,
const Valeur& bound_x_deux,
592 const Valeur& bound_y_deux,
const Valeur& bound_z_deux,
593 Vector& sol_un, Vector& sol_deux,
int num_front,
double precision) {
595 sol_un.change_triad(source_un.get_mp().get_bvect_cart()) ;
596 sol_deux.change_triad(source_deux.get_mp().get_bvect_cart()) ;
600 assert (sol_un.get_mp() == source_un.get_mp()) ;
601 assert (sol_deux.get_mp() == source_deux.get_mp()) ;
603 double orientation_un = sol_un.get_mp().get_rot_phi() ;
604 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
606 double orientation_deux = sol_deux.get_mp().get_rot_phi() ;
607 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
609 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
611 Valeur limite_x_un (bound_x_un.get_mg()) ;
612 limite_x_un = bound_x_un ;
613 Valeur limite_y_un (bound_y_un.get_mg()) ;
614 limite_y_un = bound_y_un ;
615 Valeur limite_z_un (bound_z_un.get_mg()) ;
616 limite_z_un = bound_z_un ;
618 Valeur limite_x_deux (bound_x_deux.get_mg()) ;
619 limite_x_deux = bound_x_deux ;
620 Valeur limite_y_deux (bound_y_deux.get_mg()) ;
621 limite_y_deux = bound_y_deux ;
622 Valeur limite_z_deux (bound_z_deux.get_mg()) ;
623 limite_z_deux = bound_z_deux ;
625 Valeur limite_chi_un (bound_x_un.get_mg()) ;
627 limite_chi_un.std_base_scal() ;
629 Valeur limite_chi_deux (bound_x_deux.get_mg()) ;
630 limite_chi_deux = 0 ;
631 limite_chi_deux.std_base_scal() ;
633 Mtbl xa_mtbl_un (source_un.get_mp().get_mg()) ;
634 xa_mtbl_un.set_etat_qcq() ;
635 Mtbl ya_mtbl_un (source_un.get_mp().get_mg()) ;
636 ya_mtbl_un.set_etat_qcq() ;
637 Mtbl za_mtbl_un (source_un.get_mp().get_mg()) ;
638 za_mtbl_un.set_etat_qcq() ;
639 Mtbl xa_mtbl_deux (source_deux.get_mp().get_mg()) ;
640 xa_mtbl_deux.set_etat_qcq() ;
641 Mtbl ya_mtbl_deux (source_deux.get_mp().get_mg()) ;
642 ya_mtbl_deux.set_etat_qcq() ;
643 Mtbl za_mtbl_deux (source_deux.get_mp().get_mg()) ;
644 za_mtbl_deux.set_etat_qcq() ;
646 xa_mtbl_un = source_un.get_mp().xa ;
647 ya_mtbl_un = source_un.get_mp().ya ;
648 za_mtbl_un = source_un.get_mp().za ;
650 xa_mtbl_deux = source_deux.get_mp().xa ;
651 ya_mtbl_deux = source_deux.get_mp().ya ;
652 za_mtbl_deux = source_deux.get_mp().za ;
654 double xabs, yabs, zabs ;
655 double air, theta, phi ;
658 int nbrep_un = bound_x_un.get_mg()->get_np(num_front) ;
659 int nbret_un = bound_x_un.get_mg()->get_nt(num_front) ;
660 int nbrep_deux = bound_x_deux.get_mg()->get_np(num_front) ;
661 int nbret_deux = bound_x_deux.get_mg()->get_nt(num_front) ;
662 int nz_un = bound_x_un.get_mg()->get_nzone() ;
663 int nz_deux = bound_x_deux.get_mg()->get_nzone() ;
665 const Metric_flat& ff_un (source_un.get_mp().flat_met_cart()) ;
666 const Metric_flat& ff_deux (source_deux.get_mp().flat_met_cart()) ;
669 Vector cop_so_un (source_un) ;
670 cop_so_un.dec_dzpuis(2) ;
671 cop_so_un.dec_dzpuis(2) ;
672 cop_so_un.change_triad(source_un.get_mp().get_bvect_cart()) ;
675 Scalar source_scal_un (
contract (cop_so_un.derive_cov(ff_un), 0, 1)/(lambda+1)) ;
676 if (source_scal_un.get_etat() == ETATZERO) {
677 source_scal_un.annule_hard() ;
678 source_scal_un.std_spectral_base() ;
680 source_scal_un.inc_dzpuis(2) ;
683 Vector cop_so_deux (source_deux) ;
684 cop_so_deux.dec_dzpuis(2) ;
685 cop_so_deux.dec_dzpuis(2) ;
686 cop_so_deux.change_triad(source_deux.get_mp().get_bvect_cart()) ;
688 Scalar source_scal_deux (
contract (cop_so_deux.derive_cov(ff_deux), 0, 1)/(lambda+1)) ;
689 if (source_scal_deux.get_etat() == ETATZERO) {
690 source_scal_deux.annule_hard() ;
691 source_scal_deux.std_spectral_base() ;
693 source_scal_deux.inc_dzpuis(2) ;
696 Vector copie_so_un (source_un) ;
697 copie_so_un.dec_dzpuis() ;
698 copie_so_un.change_triad(source_un.get_mp().get_bvect_cart()) ;
700 Vector copie_so_deux (source_deux) ;
701 copie_so_deux.dec_dzpuis() ;
702 copie_so_deux.change_triad(source_deux.get_mp().get_bvect_cart()) ;
705 Vector sol_un_old (sol_un) ;
706 Vector sol_deux_old (sol_deux) ;
714 Scalar chi_un (source_scal_un.poisson_dirichlet (limite_chi_un, num_front)) ;
715 Scalar chi_deux (source_scal_deux.poisson_dirichlet (limite_chi_deux, num_front)) ;
719 Vector source_vect_un (copie_so_un) ;
720 Vector grad_chi_un (chi_un.derive_con(ff_un)) ;
721 grad_chi_un.inc_dzpuis() ;
722 source_vect_un = source_vect_un - lambda*grad_chi_un ;
725 Vector source_vect_deux (copie_so_deux) ;
726 Vector grad_chi_deux (chi_deux.derive_con(ff_deux)) ;
727 grad_chi_deux.inc_dzpuis() ;
728 source_vect_deux = source_vect_deux - lambda*grad_chi_deux ;
730 sol_un_old = sol_un ;
731 sol_deux_old = sol_deux ;
734 sol_un.set(1) = source_vect_un(1).poisson_dirichlet (limite_x_un, num_front) ;
735 sol_un.set(2) = source_vect_un(2).poisson_dirichlet (limite_y_un, num_front) ;
736 sol_un.set(3) = source_vect_un(3).poisson_dirichlet (limite_z_un, num_front) ;
737 sol_deux.set(1) = source_vect_deux(1).poisson_dirichlet (limite_x_deux, num_front) ;
738 sol_deux.set(2) = source_vect_deux(2).poisson_dirichlet (limite_y_deux, num_front) ;
739 sol_deux.set(3) = source_vect_deux(3).poisson_dirichlet (limite_z_deux, num_front) ;
743 Scalar div_shift_un (
contract(sol_un.derive_cov(ff_un), 0, 1)) ;
744 div_shift_un.dec_dzpuis(2) ;
745 div_shift_un.get_spectral_va().coef_i() ;
748 for (
int k=0 ; k<nbrep_un ; k++)
749 for (
int j=0 ; j<nbret_un ; j++)
750 limite_chi_un.set(num_front, k, j, 0) =
751 div_shift_un.get_spectral_va() (num_front+1, k, j, 0) ;
752 limite_chi_un.std_base_scal() ;
754 Scalar div_shift_deux (
contract(sol_deux.derive_cov(ff_deux), 0, 1)) ;
755 div_shift_deux.dec_dzpuis(2) ;
756 div_shift_deux.get_spectral_va().coef_i() ;
758 limite_chi_deux = 1 ;
759 for (
int k=0 ; k<nbrep_deux ; k++)
760 for (
int j=0 ; j<nbret_deux ; j++)
761 limite_chi_deux.set(num_front, k, j, 0) =
762 div_shift_deux.get_spectral_va() (num_front+1, k, j, 0) ;
763 limite_chi_deux.std_base_scal() ;
767 for (
int k=0 ; k<nbrep_un ; k++)
768 for (
int j=0 ; j<nbret_un ; j++) {
769 xabs = xa_mtbl_un (num_front+1, k, j, 0) ;
770 yabs = ya_mtbl_un (num_front+1, k, j, 0) ;
771 zabs = za_mtbl_un (num_front+1, k, j, 0) ;
773 source_deux.get_mp().convert_absolute
774 (xabs, yabs, zabs, air, theta, phi) ;
776 valeur = sol_deux(1).val_point(air, theta, phi) ;
777 limite_x_un.set(num_front, k, j, 0) =
778 bound_x_un(num_front, k, j, 0) - same_orient*valeur ;
780 valeur = sol_deux(2).val_point(air, theta, phi) ;
781 limite_y_un.set(num_front, k, j, 0) =
782 bound_y_un(num_front, k, j, 0) - same_orient*valeur ;
784 valeur = sol_deux(3).val_point(air, theta, phi) ;
785 limite_z_un.set(num_front, k, j, 0) =
786 bound_z_un(num_front, k, j, 0) - valeur ;
790 for (
int k=0 ; k<nbrep_deux ; k++)
791 for (
int j=0 ; j<nbret_deux ; j++) {
792 xabs = xa_mtbl_deux (num_front+1, k, j, 0) ;
793 yabs = ya_mtbl_deux (num_front+1, k, j, 0) ;
794 zabs = za_mtbl_deux (num_front+1, k, j, 0) ;
796 source_un.get_mp().convert_absolute
797 (xabs, yabs, zabs, air, theta, phi) ;
799 valeur = sol_un(1).val_point(air, theta, phi) ;
800 limite_x_deux.set(num_front, k, j, 0) =
801 bound_x_deux(num_front, k, j, 0) - same_orient*valeur ;
803 valeur = sol_un(2).val_point(air, theta, phi) ;
804 limite_y_deux.set(num_front, k, j, 0) =
805 bound_y_deux(num_front, k, j, 0) - same_orient*valeur ;
807 valeur = sol_un(3).val_point(air, theta, phi) ;
808 limite_z_deux.set(num_front, k, j, 0) =
809 bound_z_deux(num_front, k, j, 0) - valeur ;
814 for (
int i=1 ; i<=3 ; i++) {
815 Tbl diff_un (
diffrelmax (sol_un_old(i), sol_un(i))) ;
816 for (
int j=num_front+1 ; j<nz_un ; j++)
817 if (erreur<diff_un(j))
818 erreur = diff_un(j) ;
821 for (
int i=1 ; i<=3 ; i++) {
822 Tbl diff_deux (
diffrelmax (sol_deux_old(i), sol_deux(i))) ;
823 for (
int j=num_front+1 ; j<nz_deux ; j++)
824 if (erreur<diff_deux(j))
825 erreur = diff_deux(j) ;
828 cout <<
"Pas " << conte <<
" : Difference " << erreur << endl ;
839 if (erreur < precision)
void dec2_dzpuis()
dzpuis -= 2 ;
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .