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