23 char bhole_equations_bin_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $" ;
79 #include "utilitaires.h"
80 #include "graphique.h"
86 assert ((relax >0) && (relax<=1)) ;
88 cout <<
"-----------------------------------------------" << endl ;
89 cout <<
"Resolution LAPSE" << endl ;
95 Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
98 Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
116 dirichlet_binaire (source_un, source_deux, -1., -1.,
hole1.
n_auto.
set(),
136 assert ((relax>0) && (relax<=1)) ;
138 cout <<
"-----------------------------------------------" << endl ;
139 cout <<
"Resolution PSI" << endl ;
145 Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
148 Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
162 for (
int k=0 ; k<np_un ; k++)
163 for (
int j=0 ; j<nt_un ; j++)
171 for (
int k=0 ; k<np_deux ; k++)
172 for (
int j=0 ; j<nt_deux ; j++)
177 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
198 cout <<
"------------------------------------------------" << endl ;
199 cout <<
"Resolution shift : Omega = " <<
omega << endl ;
205 if (source_un.
get_etat() == ETATZERO) {
207 for (
int i=0 ; i<3 ; i++)
215 if (source_deux.
get_etat() == ETATZERO) {
217 for (
int i=0 ; i<3 ; i++)
223 for (
int i=0 ; i<3 ; i++) {
224 if (source_un(i).get_etat() != ETATZERO)
226 if (source_deux(i).get_etat() != ETATZERO)
232 assert ((orientation_un==0) || (orientation_un == M_PI)) ;
235 assert ((orientation_deux==0) || (orientation_deux == M_PI)) ;
237 int aligne_un = (orientation_un == 0) ? 1 : -1 ;
238 int aligne_deux = (orientation_deux == 0) ? 1 : -1 ;
243 for (
int i=0 ; i<3 ; i++) {
244 if (source_un(i).
get_etat() == ETATQCQ)
246 if (source_deux(i).
get_etat() == ETATQCQ)
264 xa_mtbl_un = source_un.
get_mp()->
xa ;
265 ya_mtbl_un = source_un.
get_mp()->
ya ;
272 x_mtbl_un = source_un.
get_mp()->
x ;
273 y_mtbl_un = source_un.
get_mp()->
y ;
282 for (
int k=0 ; k<np_un ; k++)
283 for (
int j=0 ; j<nt_un ; j++)
284 lim_x_un.
set(0, k, j, 0) = aligne_un*
omega*ya_mtbl_un(0, 0, 0, 0) + aligne_un*
hole1.
omega_local*y_mtbl_un(1,k,j,0) ;
285 lim_x_un.
base = *bases_un[0] ;
290 for (
int k=0 ; k<np_un ; k++)
291 for (
int j=0 ; j<nt_un ; j++)
292 lim_y_un.
set(0, k, j, 0) = -aligne_un*
omega*xa_mtbl_un(0, 0, 0, 0) - aligne_un*
hole1.
omega_local*x_mtbl_un(1,k,j,0) ;;
293 lim_y_un.
base = *bases_un[1] ;
297 for (
int k=0 ; k<np_un ; k++)
298 for (
int j=0 ; j<nt_un ; j++)
299 lim_z_un.
set(0, k, j, 0) = 0 ;
300 lim_z_un.
base = *bases_un[2] ;
311 xa_mtbl_deux = source_deux.
get_mp()->
xa ;
312 ya_mtbl_deux = source_deux.
get_mp()->
ya ;
319 x_mtbl_deux = source_deux.
get_mp()->
x ;
320 y_mtbl_deux = source_deux.
get_mp()->
y ;
325 for (
int k=0 ; k<np_deux ; k++)
326 for (
int j=0 ; j<nt_deux ; j++)
327 lim_x_deux.
set(0, k, j, 0) = aligne_deux*
omega*ya_mtbl_deux(0, 0, 0, 0) + aligne_deux*
hole2.
omega_local*y_mtbl_deux(1,k,j,0) ;
328 lim_x_deux.
base = *bases_deux[0] ;
333 for (
int k=0 ; k<np_deux ; k++)
334 for (
int j=0 ; j<nt_deux ; j++)
335 lim_y_deux.
set(0, k, j, 0) = -aligne_deux*
omega*xa_mtbl_deux(0, 0, 0, 0) - aligne_deux*
hole2.
omega_local*x_mtbl_deux(1,k,j,0) ;
336 lim_y_deux.
base = *bases_deux[1] ;
340 for (
int k=0 ; k<np_deux ; k++)
341 for (
int j=0 ; j<nt_deux ; j++)
342 lim_z_deux.
set(0, k, j, 0) = 0 ;
343 lim_z_deux.
base = *bases_deux[2] ;
345 for (
int i=0 ; i<3 ; i++) {
347 delete bases_deux[i] ;
350 delete [] bases_deux ;
356 poisson_vect_binaire (1./3., source_un, source_deux,
357 lim_x_un, lim_y_un, lim_z_un,
358 lim_x_deux, lim_y_deux, lim_z_deux,
361 for (
int i=0 ; i<3 ; i++) {
368 (1-relax)*shift_un_old ;
370 (1-relax)*shift_deux_old ;
377 cout <<
"Difference relatives : " << diff_un <<
" " << diff_deux << endl ;