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 
90 namespace Lorene {
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
213 
214  }
215  else {
216 
217  // Computation of A^{ij}_auto
218 
219  aa_auto_un = ( hole1.beta_auto.ope_killing_conf(hole1.tgam) +
221  aa_auto_deux = ( hole2.beta_auto.ope_killing_conf(hole2.tgam) +
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 :
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 :
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 
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 
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
544  (xabs, yabs, zabs, air_un, theta, phi) ;
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 :
595  (xabs, yabs, zabs, air_un, theta, phi) ;
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 }
void decouple()
Calculates decouple which is used to obtain tkij_auto and tkij_comp.
Definition: binhor_kij.C:476
Coord xa
Absolute x coordinate.
Definition: map.h:748
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition: binhor_kij.C:91
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Multi-domain array.
Definition: mtbl.h:118
Scalar decouple
Function used to construct from the total .
Definition: isol_hor.h:1002
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Scalar n_comp
Lapse function .
Definition: isol_hor.h:918
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: scalar.C:1003
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
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
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Sym_tensor aa_auto
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
Definition: isol_hor.h:959
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
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
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
Sym_tensor aa
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
Definition: isol_hor.h:971
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
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...
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:611
Vector beta_auto
Shift function .
Definition: isol_hor.h:944
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Sym_tensor aa_comp
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
Definition: isol_hor.h:965
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Scalar n_auto
Lapse function .
Definition: isol_hor.h:915
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
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
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Coord ya
Absolute y coordinate.
Definition: map.h:749
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
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Coord za
Absolute z coordinate.
Definition: map.h:750
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition: isol_hor.h:986
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
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:305
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...
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
Scalar nn
Lapse function .
Definition: isol_hor.h:921
Coord r
r coordinate centered on the grid
Definition: map.h:736