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 char 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 
87 namespace 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) +
112  hole1.gamt_point*hole1.decouple ) / (2.* hole1.nn) ;
113 
114  aa_auto_deux = ( hole2.beta_auto.ope_killing_conf(hole2.tgam) +
115  hole2.gamt_point*hole2.decouple ) / (2.* hole2.nn) ;
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 
451  aa_deux.change_triad(hole2.mp.get_bvect_spher()) ;
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 }
void decouple()
Calculates decouple which is used to obtain tkij_auto and tkij_comp.
Definition: binhor_kij.C:473
Coord xa
Absolute x coordinate.
Definition: map.h:730
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition: binhor_kij.C:88
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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:783
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Scalar n_comp
Lapse function .
Definition: isol_hor.h:918
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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:997
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
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:784
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
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:94
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 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
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:462
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
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:481
Vector beta_auto
Shift function .
Definition: isol_hor.h:944
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
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:348
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:890
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Coord ya
Absolute y coordinate.
Definition: map.h:731
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
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:791
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:654
Coord za
Absolute z coordinate.
Definition: map.h:732
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:69
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
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
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:223
Scalar nn
Lapse function .
Definition: isol_hor.h:921
Coord r
r coordinate centered on the grid
Definition: map.h:718