LORENE
binhor_kij.C
1/*
2 * Copyright (c) 2004-2005 Francois Limousin
3 * Jose Luis Jaramillo
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24
25
26/*
27 * $Id: binhor_kij.C,v 1.12 2016/12/05 16:17:46 j_novak Exp $
28 * $Log: binhor_kij.C,v $
29 * Revision 1.12 2016/12/05 16:17:46 j_novak
30 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31 *
32 * Revision 1.11 2014/10/13 08:52:42 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.10 2014/10/06 15:13:01 j_novak
36 * Modified #include directives to use c++ syntax.
37 *
38 * Revision 1.9 2007/04/13 15:28:55 f_limousin
39 * Lots of improvements, generalisation to an arbitrary state of
40 * rotation, implementation of the spatial metric given by Samaya.
41 *
42 * Revision 1.8 2006/05/24 16:56:37 f_limousin
43 * Many small modifs.
44 *
45 * Revision 1.7 2005/09/13 18:33:15 f_limousin
46 * New function vv_bound_cart_bin(double) for computing binaries with
47 * berlin condition for the shift vector.
48 * Suppress all the symy and asymy in the importations.
49 *
50 * Revision 1.6 2005/04/29 14:02:44 f_limousin
51 * Important changes : manage the dependances between quantities (for
52 * instance psi and psi4). New function write_global(ost).
53 *
54 * Revision 1.5 2005/03/10 17:21:52 f_limousin
55 * Add the Berlin boundary condition for the shift.
56 * Some changes to avoid warnings.
57 *
58 * Revision 1.4 2005/03/03 13:49:35 f_limousin
59 * Add the spectral bases for both Scalars decouple.
60 *
61 * Revision 1.3 2005/02/07 10:48:00 f_limousin
62 * The extrinsic curvature can now be computed in the case N=0 on the
63 * horizon.
64 *
65 * Revision 1.2 2004/12/31 15:41:54 f_limousin
66 * Correction of an error
67 *
68 * Revision 1.1 2004/12/29 16:12:03 f_limousin
69 * First version
70 *
71 *
72 * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_kij.C,v 1.12 2016/12/05 16:17:46 j_novak Exp $
73 *
74 */
75
76
77//standard
78#include <cstdlib>
79#include <cmath>
80
81// Lorene
82#include "nbr_spx.h"
83#include "tenseur.h"
84#include "tensor.h"
85#include "isol_hor.h"
86#include "proto.h"
87#include "utilitaires.h"
88//#include "graphique.h"
89
90namespace Lorene {
91void Bin_hor::extrinsic_curvature () {
92
93
94 int nnt = hole1.mp.get_mg()->get_nt(1) ;
95 int nnp = hole1.mp.get_mg()->get_np(1) ;
96
97 int check ;
98 check = 0 ;
99 for (int k=0; k<nnp; k++)
100 for (int j=0; j<nnt; j++){
101 if ((hole1.n_auto+hole1.n_comp).val_grid_point(1, k, j , 0) < 1e-4){
102 check = 1 ;
103 break ;
104 }
105 }
106
107 Sym_tensor aa_auto_un (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
108 Sym_tensor aa_auto_deux (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
109
110 if (check == 0){
111
112 // Computation of A^{ij}_auto
113
114 aa_auto_un = ( hole1.beta_auto.ope_killing_conf(hole1.tgam) +
115 hole1.gamt_point*hole1.decouple ) / (2.* hole1.nn) ;
116
117 aa_auto_deux = ( hole2.beta_auto.ope_killing_conf(hole2.tgam) +
118 hole2.gamt_point*hole2.decouple ) / (2.* hole2.nn) ;
119
120
121 aa_auto_un.change_triad(hole1.mp.get_bvect_cart()) ;
122 aa_auto_deux.change_triad(hole2.mp.get_bvect_cart()) ;
123
124 for (int i=1 ; i<=3 ; i++)
125 for (int j=i ; j<=3 ; j++) {
126 if (aa_auto_un(i,j).get_etat() != ETATZERO)
127 aa_auto_un.set(i, j).raccord(3) ;
128 if (aa_auto_deux(i,j).get_etat() != ETATZERO)
129 aa_auto_deux.set(i, j).raccord(3) ;
130 }
131
132 aa_auto_un.change_triad(hole1.mp.get_bvect_spher()) ;
133 aa_auto_deux.change_triad(hole2.mp.get_bvect_spher()) ;
134
135 hole1.aa_auto = aa_auto_un ;
136 hole2.aa_auto = aa_auto_deux ;
137
138
139 // Computation of A^{ij}_comp
140
141 aa_auto_un.dec_dzpuis(2) ;
142 aa_auto_deux.dec_dzpuis(2) ;
143
144 Sym_tensor aa_comp_un (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
145 aa_comp_un.set_etat_qcq() ;
146 Sym_tensor aa_comp_deux (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
147 aa_comp_deux.set_etat_qcq() ;
148
149 aa_auto_deux.change_triad(hole2.mp.get_bvect_cart()) ;
150 aa_auto_deux.change_triad(hole1.mp.get_bvect_cart()) ;
151 assert(*(aa_auto_deux.get_triad()) == *(aa_comp_un.get_triad())) ;
152
153 // importations :
154 for (int i=1 ; i<=3 ; i++)
155 for (int j=i ; j<=3 ; j++) {
156 aa_comp_un.set(i, j).import(aa_auto_deux(i, j)) ;
157 aa_comp_un.set(i, j).set_spectral_va().set_base(aa_auto_deux(i, j).
158 get_spectral_va().get_base()) ;
159 }
160
161 aa_comp_un.inc_dzpuis(2) ;
162 aa_comp_un.change_triad(hole1.mp.get_bvect_spher()) ;
163
164 aa_auto_un.change_triad(hole1.mp.get_bvect_cart()) ;
165 aa_auto_un.change_triad(hole2.mp.get_bvect_cart()) ;
166 assert(*(aa_auto_un.get_triad()) == *(aa_comp_deux.get_triad())) ;
167 // importations :
168 for (int i=1 ; i<=3 ; i++)
169 for (int j=i ; j<=3 ; j++) {
170 aa_comp_deux.set(i, j).import(aa_auto_un(i, j)) ;
171 aa_comp_deux.set(i, j).set_spectral_va().set_base(aa_auto_un(i, j).
172 get_spectral_va().get_base()) ;
173 }
174
175 aa_comp_deux.inc_dzpuis(2) ;
176 aa_comp_deux.change_triad(hole2.mp.get_bvect_spher()) ;
177
178 /*
179 // Computation of A^{ij}_comp in the last domains
180 // -----------------------------------------------
181
182 int nz = hole1.mp.get_mg()->get_nzone() ;
183
184 Sym_tensor aa_comp_un_zec (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
185 aa_comp_un_zec.set_etat_qcq() ;
186 Sym_tensor aa_comp_deux_zec (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
187 aa_comp_deux_zec.set_etat_qcq() ;
188
189 aa_comp_un_zec = ( hole1.beta_comp().ope_killing_conf(hole1.tgam) +
190 hole1.gamt_point*(1.-hole1.decouple) ) / (2.* hole1.nn()) ;
191
192 aa_comp_deux_zec =( hole2.beta_comp().ope_killing_conf(hole2.tgam) +
193 hole2.gamt_point*(1.-hole2.decouple) ) / (2.* hole2.nn()) ;
194
195 for (int i=1 ; i<=3 ; i++)
196 for (int j=i ; j<=3 ; j++)
197 for (int l=nz-1 ; l<=nz-1 ; l++) {
198 if (aa_comp_un.set(i,j).get_etat() == ETATQCQ)
199 aa_comp_un.set(i,j).set_domain(l) =
200 aa_comp_un_zec(i,j).domain(l) ;
201 if (aa_comp_deux.set(i,j).get_etat() == ETATQCQ)
202 aa_comp_deux.set(i,j).set_domain(l)=
203 aa_comp_deux_zec(i,j).domain(l) ;
204 }
205 */
206
207 hole1.aa_comp = aa_comp_un ;
208 hole2.aa_comp = aa_comp_deux ;
209
210 // Computation of A^{ij}_ total
211 hole1.aa = hole1.aa_auto + hole1.aa_comp ;
212 hole2.aa = hole2.aa_auto + hole2.aa_comp ;
213
214 }
215 else {
216
217 // Computation of A^{ij}_auto
218
219 aa_auto_un = ( hole1.beta_auto.ope_killing_conf(hole1.tgam) +
220 hole1.gamt_point*hole1.decouple ) ;
221 aa_auto_deux = ( hole2.beta_auto.ope_killing_conf(hole2.tgam) +
222 hole2.gamt_point*hole2.decouple ) ;
223
224 aa_auto_un.change_triad(hole1.mp.get_bvect_cart()) ;
225 aa_auto_deux.change_triad(hole2.mp.get_bvect_cart()) ;
226
227 for (int i=1 ; i<=3 ; i++)
228 for (int j=1 ; j<=3 ; j++) {
229 if (aa_auto_un(i,j).get_etat() != ETATZERO)
230 aa_auto_un.set(i, j).raccord(3) ;
231 if (aa_auto_deux(i,j).get_etat() != ETATZERO)
232 aa_auto_deux.set(i, j).raccord(3) ;
233 }
234
235 // Computation of A^{ij}_comp
236
237 Sym_tensor aa_auto_1 (aa_auto_un) ;
238 Sym_tensor aa_auto_2 (aa_auto_deux) ;
239
240 aa_auto_1.dec_dzpuis(2) ;
241 aa_auto_2.dec_dzpuis(2) ;
242
243 Sym_tensor aa_comp_un (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
244 aa_comp_un.set_etat_qcq() ;
245 Sym_tensor aa_comp_deux (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
246 aa_comp_deux.set_etat_qcq() ;
247
248 aa_auto_2.change_triad(hole1.mp.get_bvect_cart()) ;
249 assert(*(aa_auto_2.get_triad()) == *(aa_comp_un.get_triad())) ;
250 // importations :
251 aa_comp_un.set(1, 1).import(aa_auto_2(1, 1)) ;
252 aa_comp_un.set(1, 2).import(aa_auto_2(1, 2)) ;
253 aa_comp_un.set(1, 3).import(aa_auto_2(1, 3)) ;
254 aa_comp_un.set(2, 2).import(aa_auto_2(2, 2)) ;
255 aa_comp_un.set(2, 3).import(aa_auto_2(2, 3)) ;
256 aa_comp_un.set(3, 3).import(aa_auto_2(3, 3)) ;
257
258 aa_comp_un.std_spectral_base() ;
259 aa_comp_un.inc_dzpuis(2) ;
260
261 aa_auto_1.change_triad(hole2.mp.get_bvect_cart()) ;
262 assert(*(aa_auto_1.get_triad()) == *(aa_comp_deux.get_triad())) ;
263 // importations :
264 aa_comp_deux.set(1, 1).import(aa_auto_1(1, 1)) ;
265 aa_comp_deux.set(1, 2).import(aa_auto_1(1, 2)) ;
266 aa_comp_deux.set(1, 3).import(aa_auto_1(1, 3)) ;
267 aa_comp_deux.set(2, 2).import(aa_auto_1(2, 2)) ;
268 aa_comp_deux.set(2, 3).import(aa_auto_1(2, 3)) ;
269 aa_comp_deux.set(3, 3).import(aa_auto_1(3, 3)) ;
270
271 aa_comp_deux.std_spectral_base() ;
272 aa_comp_deux.inc_dzpuis(2) ;
273
274 // Computation of A^{ij}_ total
275 Sym_tensor aa_tot_un (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
276 Sym_tensor aa_tot_deux (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
277 aa_tot_un = aa_auto_un + aa_comp_un ;
278 aa_tot_deux = aa_auto_deux + aa_comp_deux ;
279
280 Sym_tensor temp_aa_tot1 (aa_tot_un) ;
281 Sym_tensor temp_aa_tot2 (aa_tot_deux) ;
282
283 temp_aa_tot1.change_triad(hole1.mp.get_bvect_spher()) ;
284 temp_aa_tot2.change_triad(hole2.mp.get_bvect_spher()) ;
285
286 // Regularisation
287 // --------------
288
289 Sym_tensor aa_un (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
290 Sym_tensor aa_deux (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
291
292 int nz_un = hole1.mp.get_mg()->get_nzone() ;
293 int nz_deux = hole2.mp.get_mg()->get_nzone() ;
294
295 Scalar ntot_un (hole1.n_auto+hole1.n_comp) ;
296 ntot_un = division_xpun (ntot_un, 0) ;
297 ntot_un.raccord(1) ;
298
299 Scalar ntot_deux (hole2.n_auto+hole2.n_comp) ;
300 ntot_deux = division_xpun (ntot_deux, 0) ;
301 ntot_deux.raccord(1) ;
302
303 // THE TWO Aij are aligned of not !
304 double orientation_un = aa_auto_un.get_mp().get_rot_phi() ;
305 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
306 double orientation_deux = aa_auto_deux.get_mp().get_rot_phi() ;
307 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
308 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
309
310
311 // Loop on the composants :
312 for (int lig = 1 ; lig<=3 ; lig++)
313 for (int col = lig ; col<=3 ; col++) {
314
315 // The orientation
316 int ind = 1 ;
317 if (lig !=3)
318 ind *= -1 ;
319 if (col != 3)
320 ind *= -1 ;
321 if (same_orient == 1)
322 ind = 1 ;
323
324 // Close to hole one :
325 Scalar auxi_un (aa_tot_un(lig, col)/2.) ;
326 auxi_un.dec_dzpuis(2) ;
327 auxi_un = division_xpun (auxi_un, 0) ;
328 auxi_un = auxi_un / ntot_un ;
329 if (auxi_un.get_etat() != ETATZERO)
330 auxi_un.raccord(1) ;
331
332 // Close to hole two :
333 Scalar auxi_deux (aa_tot_deux(lig, col)/2.) ;
334 auxi_deux.dec_dzpuis(2) ;
335 auxi_deux = division_xpun (auxi_deux, 0) ;
336 auxi_deux = auxi_deux / ntot_deux ;
337 if (auxi_deux.get_etat() != ETATZERO)
338 auxi_deux.raccord(1) ;
339
340 // copy :
341 Scalar copie_un (aa_tot_un(lig, col)) ;
342 copie_un.dec_dzpuis(2) ;
343
344 Scalar copie_deux (aa_tot_deux(lig, col)) ;
345 copie_deux.dec_dzpuis(2) ;
346
347 double lim_un = hole1.mp.get_alpha()[1] + hole1.mp.get_beta()[1] ;
348 double lim_deux = hole2.mp.get_alpha()[1] + hole2.mp.get_beta()[1] ;
349
350 Mtbl xabs_un (hole1.mp.xa) ;
351 Mtbl yabs_un (hole1.mp.ya) ;
352 Mtbl zabs_un (hole1.mp.za) ;
353
354 Mtbl xabs_deux (hole2.mp.xa) ;
355 Mtbl yabs_deux (hole2.mp.ya) ;
356 Mtbl zabs_deux (hole2.mp.za) ;
357
358 double xabs, yabs, zabs, air, theta, phi ;
359
360 if (auxi_un.get_etat() != ETATZERO){
361 // Loop on the other zones :
362 for (int l=2 ; l<nz_un ; l++) {
363
364 int nr = hole1.mp.get_mg()->get_nr (l) ;
365
366 if (l==nz_un-1)
367 nr -- ;
368
369 int np = hole1.mp.get_mg()->get_np (l) ;
370 int nt = hole1.mp.get_mg()->get_nt (l) ;
371
372 for (int k=0 ; k<np ; k++)
373 for (int j=0 ; j<nt ; j++)
374 for (int i=0 ; i<nr ; i++) {
375
376 xabs = xabs_un (l, k, j, i) ;
377 yabs = yabs_un (l, k, j, i) ;
378 zabs = zabs_un (l, k, j, i) ;
379
380 // coordinates of the point in 2 :
381 hole2.mp.convert_absolute
382 (xabs, yabs, zabs, air, theta, phi) ;
383
384 if (air >= lim_deux)
385 // Far from hole two : no pb :
386 auxi_un.set_grid_point(l, k, j, i) =
387 copie_un.val_grid_point(l, k, j, i) /
388 ntot_un.val_grid_point(l, k, j, i)/2. ;
389 else
390 // close to hole two :
391 auxi_un.set_grid_point(l, k, j, i) =
392 ind * auxi_deux.val_point (air, theta, phi) ;
393
394 }
395
396 // Case infinity :
397 if (l==nz_un-1)
398 for (int k=0 ; k<np ; k++)
399 for (int j=0 ; j<nt ; j++)
400 auxi_un.set_grid_point(nz_un-1, k, j, nr) = 0 ;
401 }
402 }
403
404 if (auxi_deux.get_etat() != ETATZERO){
405 // The second hole :
406 for (int l=2 ; l<nz_deux ; l++) {
407
408 int nr = hole2.mp.get_mg()->get_nr (l) ;
409
410 if (l==nz_deux-1)
411 nr -- ;
412
413 int np = hole2.mp.get_mg()->get_np (l) ;
414 int nt = hole2.mp.get_mg()->get_nt (l) ;
415
416 for (int k=0 ; k<np ; k++)
417 for (int j=0 ; j<nt ; j++)
418 for (int i=0 ; i<nr ; i++) {
419
420 xabs = xabs_deux (l, k, j, i) ;
421 yabs = yabs_deux (l, k, j, i) ;
422 zabs = zabs_deux (l, k, j, i) ;
423
424 // coordinates of the point in 2 :
425 hole1.mp.convert_absolute
426 (xabs, yabs, zabs, air, theta, phi) ;
427
428 if (air >= lim_un)
429 // Far from hole one : no pb :
430 auxi_deux.set_grid_point(l, k, j, i) =
431 copie_deux.val_grid_point(l, k, j, i) /
432 ntot_deux.val_grid_point(l, k, j, i) /2.;
433 else
434 // close to hole one :
435 auxi_deux.set_grid_point(l, k, j, i) =
436 ind * auxi_un.val_point (air, theta, phi) ;
437 }
438 // Case infinity :
439 if (l==nz_deux-1)
440 for (int k=0 ; k<np ; k++)
441 for (int j=0 ; j<nt ; j++)
442 auxi_un.set_grid_point(nz_deux-1, k, j, nr) = 0 ;
443 }
444 }
445
446 auxi_un.inc_dzpuis(2) ;
447 auxi_deux.inc_dzpuis(2) ;
448
449 aa_un.set(lig, col) = auxi_un ;
450 aa_deux.set(lig, col) = auxi_deux ;
451 }
452
453 aa_un.change_triad(hole1.mp.get_bvect_spher()) ;
454 aa_deux.change_triad(hole2.mp.get_bvect_spher()) ;
455
456 hole1.aa = aa_un ;
457 hole2.aa = aa_deux ;
458
459 aa_auto_un.change_triad(hole1.mp.get_bvect_spher()) ;
460 aa_auto_deux.change_triad(hole2.mp.get_bvect_spher()) ;
461
462 for (int lig=1 ; lig<=3 ; lig++)
463 for (int col=lig ; col<=3 ; col++) {
464 aa_auto_un.set(lig, col) = aa_un(lig, col)*hole1.decouple ;
465 aa_auto_deux.set(lig, col) = aa_deux(lig, col)*hole2.decouple ;
466 }
467
468 hole1.aa_auto = aa_auto_un ;
469 hole2.aa_auto = aa_auto_deux ;
470
471 }
472
473}
474
475
476void Bin_hor::decouple () {
477
478 int nz_un = hole1.mp.get_mg()->get_nzone() ;
479 int nz_deux = hole2.mp.get_mg()->get_nzone() ;
480
481 // We determine R_limite :
482 double distance = hole1.mp.get_ori_x() - hole2.mp.get_ori_x() ;
483 double lim_un = distance/2. ;
484 double lim_deux = distance/2. ;
485 double int_un = distance/6. ;
486 double int_deux = distance/6. ;
487
488 // The functions used.
489 Scalar fonction_f_un (hole1.mp) ;
490 fonction_f_un = 0.5*pow(
491 cos((hole1.mp.r-int_un)*M_PI/2./(lim_un-int_un)), 2.)+0.5 ;
492 fonction_f_un.std_spectral_base();
493
494 Scalar fonction_g_un (hole1.mp) ;
495 fonction_g_un = 0.5*pow
496 (sin((hole1.mp.r-int_un)*M_PI/2./(lim_un-int_un)), 2.) ;
497 fonction_g_un.std_spectral_base();
498
499 Scalar fonction_f_deux (hole2.mp) ;
500 fonction_f_deux = 0.5*pow(
501 cos((hole2.mp.r-int_deux)*M_PI/2./(lim_deux-int_deux)), 2.)+0.5 ;
502 fonction_f_deux.std_spectral_base();
503
504 Scalar fonction_g_deux (hole2.mp) ;
505 fonction_g_deux = 0.5*pow(
506 sin((hole2.mp.r-int_deux)*M_PI/2./(lim_un-int_deux)), 2.) ;
507 fonction_g_deux.std_spectral_base();
508
509 // The functions total :
510 Scalar decouple_un (hole1.mp) ;
511 decouple_un.allocate_all() ;
512 Scalar decouple_deux (hole2.mp) ;
513 decouple_deux.allocate_all() ;
514
515 Mtbl xabs_un (hole1.mp.xa) ;
516 Mtbl yabs_un (hole1.mp.ya) ;
517 Mtbl zabs_un (hole1.mp.za) ;
518
519 Mtbl xabs_deux (hole2.mp.xa) ;
520 Mtbl yabs_deux (hole2.mp.ya) ;
521 Mtbl zabs_deux (hole2.mp.za) ;
522
523 double xabs, yabs, zabs, air_un, air_deux, theta, phi ;
524
525 for (int l=0 ; l<nz_un ; l++) {
526 int nr = hole1.mp.get_mg()->get_nr (l) ;
527
528 if (l==nz_un-1)
529 nr -- ;
530
531 int np = hole1.mp.get_mg()->get_np (l) ;
532 int nt = hole1.mp.get_mg()->get_nt (l) ;
533
534 for (int k=0 ; k<np ; k++)
535 for (int j=0 ; j<nt ; j++)
536 for (int i=0 ; i<nr ; i++) {
537
538 xabs = xabs_un (l, k, j, i) ;
539 yabs = yabs_un (l, k, j, i) ;
540 zabs = zabs_un (l, k, j, i) ;
541
542 // Coordinates of the point
543 hole1.mp.convert_absolute
544 (xabs, yabs, zabs, air_un, theta, phi) ;
545 hole2.mp.convert_absolute
546 (xabs, yabs, zabs, air_deux, theta, phi) ;
547
548 if (air_un <= lim_un)
549 if (air_un < int_un)
550 decouple_un.set_grid_point(l, k, j, i) = 1 ;
551 else
552 // Close to hole 1 :
553 decouple_un.set_grid_point(l, k, j, i) =
554 fonction_f_un.val_grid_point(l, k, j, i) ;
555 else
556 if (air_deux <= lim_deux)
557 if (air_deux < int_deux)
558 decouple_un.set_grid_point(l, k, j, i) = 0 ;
559 else
560 // Close to hole 2 :
561 decouple_un.set_grid_point(l, k, j, i) =
562 fonction_g_deux.val_point (air_deux, theta, phi) ;
563
564 else
565 // Far from each holes :
566 decouple_un.set_grid_point(l, k, j, i) = 0.5 ;
567 }
568
569 // Case infinity :
570 if (l==nz_un-1)
571 for (int k=0 ; k<np ; k++)
572 for (int j=0 ; j<nt ; j++)
573 decouple_un.set_grid_point(nz_un-1, k, j, nr)=0.5 ;
574 }
575
576 for (int l=0 ; l<nz_deux ; l++) {
577 int nr = hole2.mp.get_mg()->get_nr (l) ;
578
579 if (l==nz_deux-1)
580 nr -- ;
581
582 int np = hole2.mp.get_mg()->get_np (l) ;
583 int nt = hole2.mp.get_mg()->get_nt (l) ;
584
585 for (int k=0 ; k<np ; k++)
586 for (int j=0 ; j<nt ; j++)
587 for (int i=0 ; i<nr ; i++) {
588
589 xabs = xabs_deux (l, k, j, i) ;
590 yabs = yabs_deux (l, k, j, i) ;
591 zabs = zabs_deux (l, k, j, i) ;
592
593 // coordinates of the point :
594 hole1.mp.convert_absolute
595 (xabs, yabs, zabs, air_un, theta, phi) ;
596 hole2.mp.convert_absolute
597 (xabs, yabs, zabs, air_deux, theta, phi) ;
598
599 if (air_deux <= lim_deux)
600 if (air_deux < int_deux)
601 decouple_deux.set_grid_point(l, k, j, i) = 1 ;
602 else
603 // close to hole two :
604 decouple_deux.set_grid_point(l, k, j, i) =
605 fonction_f_deux.val_grid_point(l, k, j, i) ;
606 else
607 if (air_un <= lim_un)
608 if (air_un < int_un)
609 decouple_deux.set_grid_point(l, k, j, i) = 0 ;
610 else
611 // close to hole one :
612 decouple_deux.set_grid_point(l, k, j, i) =
613 fonction_g_un.val_point (air_un, theta, phi) ;
614
615 else
616 // Far from each hole :
617 decouple_deux.set_grid_point(l, k, j, i) = 0.5 ;
618 }
619
620 // Case infinity :
621 if (l==nz_deux-1)
622 for (int k=0 ; k<np ; k++)
623 for (int j=0 ; j<nt ; j++)
624 decouple_deux.set_grid_point(nz_un-1, k, j, nr)=0.5 ;
625 }
626
627 decouple_un.std_spectral_base() ;
628 decouple_deux.std_spectral_base() ;
629
630 hole1.decouple = decouple_un ;
631 hole2.decouple = decouple_deux ;
632
633}
634}
Lorene::Bin_hor::hole1
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Lorene::Bin_hor::hole2
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
Lorene::Bin_hor::decouple
void decouple()
Calculates decouple which is used to obtain tkij_auto and tkij_comp.
Definition binhor_kij.C:476
Lorene::Bin_hor::extrinsic_curvature
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition binhor_kij.C:91
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition scalar_raccord.C:65
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:790
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:643
Lorene::Scalar::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:473
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:610
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:560
Lorene::Scalar::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:371
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:690
Lorene::Scalar::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition scalar.C:896
Lorene::Scalar::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:421
Lorene::Scalar::import
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition scalar_import.C:71
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:226
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base ).
Definition valeur.C:813
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:72
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:97
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:874
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:817
Lorene::Tensor::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tensor_change_triad.C:80
Lorene::Tensor::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:825
Lorene::Tensor::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:879
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:490
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:663
Lorene::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:935
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:732