LORENE
bound_hor.C
1 /*
2  * Method of class Isol_hor to compute boundary conditions
3  *
4  * (see file isol_hor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Jose Luis Jaramillo
10  * Francois Limousin
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 char bound_hor_C[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/bound_hor.C,v 1.36 2014/10/13 08:53:00 j_novak Exp $" ;
30 
31 /*
32  * $Id: bound_hor.C,v 1.36 2014/10/13 08:53:00 j_novak Exp $
33  * $Log: bound_hor.C,v $
34  * Revision 1.36 2014/10/13 08:53:00 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.35 2014/10/06 15:13:10 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.34 2008/08/27 11:22:25 j_novak
41  * Minor modifications
42  *
43  * Revision 1.33 2008/08/19 06:42:00 j_novak
44  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
45  * cast-type operations, and constant strings that must be defined as const char*
46  *
47  * Revision 1.32 2007/04/13 15:28:35 f_limousin
48  * Lots of improvements, generalisation to an arbitrary state of
49  * rotation, implementation of the spatial metric given by Samaya.
50  *
51  * Revision 1.31 2006/08/01 14:35:48 f_limousin
52  * Many small modifs
53  *
54  * Revision 1.30 2006/02/22 17:02:04 f_limousin
55  * Removal of warnings
56  *
57  * Revision 1.29 2006/02/22 16:29:55 jl_jaramillo
58  * corrections on the relaxation and boundary conditions
59  *
60  * Revision 1.27 2005/10/24 16:44:40 jl_jaramillo
61  * Cook boundary condition ans minot bound of kss
62  *
63  * Revision 1.26 2005/10/21 16:20:55 jl_jaramillo
64  * Version for the paper JaramL05
65  *
66  * Revision 1.25 2005/09/13 18:33:17 f_limousin
67  * New function vv_bound_cart_bin(double) for computing binaries with
68  * berlin condition for the shift vector.
69  * Suppress all the symy and asymy in the importations.
70  *
71  * Revision 1.24 2005/09/12 12:33:54 f_limousin
72  * Compilation Warning - Change of convention for the angular velocity
73  * Add Berlin boundary condition in the case of binary horizons.
74  *
75  * Revision 1.23 2005/07/08 13:15:23 f_limousin
76  * Improvements of boundary_vv_cart(), boundary_nn_lapl().
77  * Add a fonction to compute the departure of axisymmetry.
78  *
79  * Revision 1.22 2005/06/09 08:05:32 f_limousin
80  * Implement a new function sol_elliptic_boundary() and
81  * Vector::poisson_boundary(...) which solve the vectorial poisson
82  * equation (method 6) with an inner boundary condition.
83  *
84  * Revision 1.21 2005/05/12 14:48:07 f_limousin
85  * New boundary condition for the lapse : boundary_nn_lapl().
86  *
87  * Revision 1.20 2005/04/29 14:04:17 f_limousin
88  * Implementation of boundary_vv_x (y,z) for binary black holes.
89  *
90  * Revision 1.19 2005/04/19 16:40:51 jl_jaramillo
91  * change of sign of ang_vel in vv_bound_cart. Convention of Phys. Rep.
92  *
93  * Revision 1.18 2005/04/08 12:16:52 f_limousin
94  * Function set_psi(). And dependance in phi.
95  *
96  * Revision 1.17 2005/04/02 15:49:21 f_limousin
97  * New choice (Lichnerowicz) for aaquad. New member data nz.
98  *
99  * Revision 1.16 2005/03/22 13:25:36 f_limousin
100  * Small changes. The angular velocity and A^{ij} are computed
101  * with a differnet sign.
102  *
103  * Revision 1.15 2005/03/10 10:19:42 f_limousin
104  * Add the regularisation of the shift in the case of a single black hole
105  * and lapse zero on the horizon.
106  *
107  * Revision 1.14 2005/03/03 10:00:55 f_limousin
108  * The funtions beta_boost_x() and beta_boost_z() have been added.
109  *
110  * Revision 1.13 2005/02/24 17:21:04 f_limousin
111  * Suppression of the function beta_bound_cart() and direct computation
112  * of boundary_beta_x, y and z.
113  *
114  * Revision 1.12 2004/12/31 15:34:37 f_limousin
115  * Modifications to avoid warnings
116  *
117  * Revision 1.11 2004/12/22 18:15:16 f_limousin
118  * Many different changes.
119  *
120  * Revision 1.10 2004/11/24 19:30:58 jl_jaramillo
121  * Berlin boundary conditions vv_bound_cart
122  *
123  * Revision 1.9 2004/11/18 09:49:44 jl_jaramillo
124  * Some new conditions for the shift (Neumann + Dirichlet)
125  *
126  * Revision 1.8 2004/11/05 10:52:26 f_limousin
127  * Replace double aa by double cc in argument of boundary_beta_x
128  * boundary_beta_y and boundary_beta_z to avoid warnings.
129  *
130  * Revision 1.7 2004/10/29 15:42:14 jl_jaramillo
131  * Static shift boundary conbdition
132  *
133  * Revision 1.6 2004/10/01 16:46:51 f_limousin
134  * Added a pure Dirichlet boundary condition
135  *
136  * Revision 1.5 2004/09/28 16:06:41 f_limousin
137  * Correction of an error when taking the bases of the boundary
138  * condition for the shift.
139  *
140  * Revision 1.4 2004/09/17 13:36:23 f_limousin
141  * Add some new boundary conditions
142  *
143  * Revision 1.2 2004/09/09 16:53:49 f_limousin
144  * Add the two lines $Id: bound_hor.C,v 1.36 2014/10/13 08:53:00 j_novak Exp $Log: for CVS.
145  *
146  *
147  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/bound_hor.C,v 1.36 2014/10/13 08:53:00 j_novak Exp $
148  *
149  */
150 
151 // C++ headers
152 #include "headcpp.h"
153 
154 // C headers
155 #include <cstdlib>
156 #include <cassert>
157 
158 // Lorene headers
159 #include "time_slice.h"
160 #include "isol_hor.h"
161 #include "metric.h"
162 #include "evolution.h"
163 #include "unites.h"
164 #include "graphique.h"
165 #include "utilitaires.h"
166 #include "param.h"
167 
168 
169 // Dirichlet boundary condition for Psi
170 //-------------------------------------
171 // ONE HAS TO GUARANTEE THAT BETA IS NOT ZERO, BUT IT IS PROPORTIONAL TO THE RADIAL VECTOR
172 
173 namespace Lorene {
175 
176  Scalar tmp = - 6 * contract(beta(), 0, psi().derive_cov(ff), 0) ;
177  tmp = tmp / (contract(beta().derive_cov(ff), 0, 1) - nn() * trk() ) - 1 ;
178 
179  // We have substracted 1, since we solve for zero condition at infinity
180  //and then we add 1 to the solution
181 
182  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
183 
184  int nnp = mp.get_mg()->get_np(1) ;
185  int nnt = mp.get_mg()->get_nt(1) ;
186 
187  psi_bound = 1 ;
188 
189  for (int k=0 ; k<nnp ; k++)
190  for (int j=0 ; j<nnt ; j++)
191  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
192 
193  psi_bound.std_base_scal() ;
194 
195  return psi_bound ;
196 
197 }
198 
199 
200 // Neumann boundary condition for Psi
201 //-------------------------------------
202 
204 
205  // Introduce 2-trace gamma tilde dot
206  Scalar tmp = - 1./ 6. * psi() * (beta().divergence(ff) - nn() * trk() )
207  - beta()(2)* psi().derive_cov(ff)(2) - beta()(3)* psi().derive_cov(ff)(3) ;
208 
209  tmp = tmp / beta()(1) ;
210 
211  // in this case you don't have to substract any value
212 
213  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
214 
215  int nnp = mp.get_mg()->get_np(1) ;
216  int nnt = mp.get_mg()->get_nt(1) ;
217 
218  psi_bound = 1 ;
219 
220  for (int k=0 ; k<nnp ; k++)
221  for (int j=0 ; j<nnt ; j++)
222  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
223 
224  psi_bound.std_base_scal() ;
225 
226  return psi_bound ;
227 
228 }
229 
230 
232 
233  Scalar tmp = psi() * psi() * psi() * trk()
234  - contract(k_dd(), 0, 1, tradial_vect_hor() * tradial_vect_hor(), 0, 1)
235  / psi()
236  - 4.* contract(tradial_vect_hor(), 0, psi().derive_cov(ff), 0) ;
237 
238  // rho = 1 is the standart case.
239  double rho = 1. ;
240  tmp += (tradial_vect_hor().divergence(ff)) * (rho - 1.)*psi() ;
241 
242  tmp = tmp / (rho * tradial_vect_hor().divergence(ff)) - 1. ;
243 
244  tmp.std_spectral_base() ;
245  tmp.inc_dzpuis(2) ;
246 
247  // We have substracted 1, since we solve for zero condition at infinity
248  //and then we add 1 to the solution
249 
250  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
251 
252  int nnp = mp.get_mg()->get_np(1) ;
253  int nnt = mp.get_mg()->get_nt(1) ;
254 
255  psi_bound = 1 ;
256 
257  for (int k=0 ; k<nnp ; k++)
258  for (int j=0 ; j<nnt ; j++)
259  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
260 
261  psi_bound.std_base_scal() ;
262 
263  return psi_bound ;
264 
265 }
266 
268 
269  Scalar tmp = psi() * psi() * psi() * trk()
270  - contract(k_dd(), 0, 1, tradial_vect_hor() * tradial_vect_hor(), 0, 1)
271  / psi()
273  - 4 * ( tradial_vect_hor()(2) * psi().derive_cov(ff)(2)
274  + tradial_vect_hor()(3) * psi().derive_cov(ff)(3) ) ;
275 
276  tmp = tmp / (4 * tradial_vect_hor()(1)) ;
277 
278  // in this case you don't have to substract any value
279 
280  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
281 
282  int nnp = mp.get_mg()->get_np(1) ;
283  int nnt = mp.get_mg()->get_nt(1) ;
284 
285  psi_bound = 1 ;
286 
287  for (int k=0 ; k<nnp ; k++)
288  for (int j=0 ; j<nnt ; j++)
289  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
290 
291  psi_bound.std_base_scal() ;
292 
293  return psi_bound ;
294 
295 }
296 
298 
299  int nnp = mp.get_mg()->get_np(1) ;
300  int nnt = mp.get_mg()->get_nt(1) ;
301 
302  Valeur psi_bound (mp.get_mg()->get_angu()) ;
303 
304  psi_bound = 1 ; // Juste pour affecter dans espace des configs ;
305 
306 // if (psi_comp_evol.is_known(jtime)) {
307 // for (int k=0 ; k<nnp ; k++)
308 // for (int j=0 ; j<nnt ; j++)
309 // psi_bound.set(0, k, j, 0) = - 0.5/radius*(psi_auto().val_grid_point(1, k, j, 0) + psi_comp().val_grid_point(1, k, j, 0)) ;
310 // }
311 // else {
312  for (int k=0 ; k<nnp ; k++)
313  for (int j=0 ; j<nnt ; j++)
314  psi_bound.set(0, k, j, 0) = - 0.5/radius*psi().val_grid_point(1, k, j, 0) ;
315 // }
316 
317 
318  psi_bound.std_base_scal() ;
319 
320  return psi_bound ;
321 
322 }
323 
325 
326  Scalar rho (mp) ;
327  rho = 5. ;
328  rho.std_spectral_base() ;
329 
330 
331  Scalar tmp(mp) ;
332  // tmp = nn()+1. - 1 ;
333 
334 
335  //Scalar b_tilde_tmp = b_tilde() ;
336  //b_tilde_tmp.set_domain(0) = 1. ;
337  //tmp = pow(nn()/b_tilde_tmp, 0.5) ;
338 
339 
340  tmp = 1.5 ;
341  tmp.std_spectral_base() ;
342 
343  //tmp = 1/ (2* nn()) ;
344 
345  tmp = (tmp + rho * psi())/(1 + rho) ;
346 
347  tmp = tmp - 1. ;
348 
349 
350 
351  Valeur psi_bound (mp.get_mg()->get_angu()) ;
352 
353  psi_bound = 1 ; // Juste pour affecter dans espace des configs ;
354 
355  int nnp = mp.get_mg()->get_np(1) ;
356  int nnt = mp.get_mg()->get_nt(1) ;
357 
358  for (int k=0 ; k<nnp ; k++)
359  for (int j=0 ; j<nnt ; j++)
360  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
361 
362  psi_bound.std_base_scal() ;
363 
364  return psi_bound ;
365 
366 }
367 
368 // Dirichlet boundary condition on nn using the extrinsic curvature
369 // (No time evolution taken into account! Make this)
370 //--------------------------------------------------------------------------
372 
373  Scalar tmp(mp) ;
374  double rho = 0. ;
375 
376  // Scalar kk_rr = 0.8*psi().derive_cov(mp.flat_met_spher())(1) / psi() ;
377  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
378  , k_dd(), 0, 1 ) ;
379 
380  Scalar k_kerr (mp) ;
381  k_kerr = 0.1 ;//1.*kappa_hor() ;
382  k_kerr.std_spectral_base() ;
383  k_kerr.inc_dzpuis(2) ;
384 
385  Scalar temp (rho*nn()) ;
386  temp.inc_dzpuis(2) ;
387 
388  tmp = contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0) + temp
389  - k_kerr ;
390 
391  tmp = tmp / (kk_rr + rho) - 1;
392 
393  Scalar diN (contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0)) ;
394  cout << "k_kerr = " << k_kerr.val_grid_point(1, 0, 0, 0) << endl ;
395  cout << "D^i N = " << diN.val_grid_point(1, 0, 0, 0) << endl ;
396  cout << "kss = " << kk_rr.val_grid_point(1, 0, 0, 0) << endl ;
397 
398  // We have substracted 1, since we solve for zero condition at infinity
399  //and then we add 1 to the solution
400 
401  int nnp = mp.get_mg()->get_np(1) ;
402  int nnt = mp.get_mg()->get_nt(1) ;
403 
404  Valeur nn_bound (mp.get_mg()->get_angu()) ;
405 
406  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
407 
408 
409  for (int k=0 ; k<nnp ; k++)
410  for (int j=0 ; j<nnt ; j++)
411  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
412 
413  nn_bound.std_base_scal() ;
414 
415  return nn_bound ;
416 
417 }
418 
419 
420 const Valeur Isol_hor::boundary_nn_Neu_kk(int step) const {
421 
422  const Vector& dnnn = nn().derive_cov(ff) ;
423  double rho = 5. ;
424 
425  step = 100 ; // Truc bidon pour eviter warning
426 
427  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
428  , k_dd(), 0, 1 ) ;
429 
430  Scalar k_kerr (mp) ;
431  k_kerr = kappa_hor() ;
432  //k_kerr = 0.06 ;
433  k_kerr.std_spectral_base() ;
434  k_kerr.inc_dzpuis(2) ;
435 
436  Scalar k_hor (mp) ;
437  k_hor = kappa_hor() ;
438  k_hor.std_spectral_base() ;
439 
440  // Vector beta_tilde_d (beta().down(0, met_gamt)) ;
441  //Scalar naass = 1./2. * contract( met_gamt.radial_vect() * met_gamt.radial_vect(), 0, 1, beta_tilde_d.ope_killing_conf(met_gamt) , 0, 1 ) ;
442 
444  0, 1, aa_auto().up_down(met_gamt), 0, 1) ;
445 
446  Scalar traceK = 1./3. * nn() * trk() *
448  , met_gamt.cov() , 0, 1 ) ;
449 
450  Scalar sdN (contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0)) ;
451 
452 
453 
454  Scalar tmp = psi() * psi() * ( k_kerr + naass + 1.* traceK)
455  - met_gamt.radial_vect()(2) * dnnn(2)
456  - met_gamt.radial_vect()(3) * dnnn(3)
457  + ( rho - 1 ) * met_gamt.radial_vect()(1) * dnnn(1) ;
458 
459 
460  tmp = tmp / (rho * met_gamt.radial_vect()(1)) ;
461 
462 
463 
464 
465  Scalar rhs ( sdN - nn() * kk_rr) ;
466  cout << "kappa_pres = " << k_kerr.val_grid_point(1, 0, 0, 0) << endl ;
467  cout << "kappa_hor = " << k_hor.val_grid_point(1, 0, 0, 0) << endl ;
468  cout << "sDN = " << sdN.val_grid_point(1, 0, 0, 0) << endl ;
469 
470  // in this case you don't have to substract any value
471 
472  int nnp = mp.get_mg()->get_np(1) ;
473  int nnt = mp.get_mg()->get_nt(1) ;
474 
475  Valeur nn_bound (mp.get_mg()->get_angu()) ;
476 
477  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
478 
479  for (int k=0 ; k<nnp ; k++)
480  for (int j=0 ; j<nnt ; j++)
481  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
482 
483  nn_bound.std_base_scal() ;
484 
485  return nn_bound ;
486 
487 }
488 
490 
491  const Vector& dnnn = nn().derive_cov(ff) ;
492  double rho = 1. ;
493 
494 
495  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
496  , k_dd(), 0, 1 ) ;
497 
498  Sym_tensor qq_uu = gam_uu() - gam().radial_vect() * gam().radial_vect() ;
499 
500  // Free function L_theta = h
501  //--------------------------
502  Scalar L_theta (mp) ;
503  L_theta = .0;
504  L_theta.std_spectral_base() ;
505  L_theta.inc_dzpuis(4) ;
506 
507  //Divergence of Omega
508  //-------------------
509  Vector hom1 = nn().derive_cov(met_gamt)/nn() ;
510  hom1 = contract(qq_uu.down(1, gam()), 0, hom1, 0) ;
511 
512  Vector hom2 = -contract( k_dd(), 1, gam().radial_vect() , 0) ;
513  hom2 = contract(qq_uu.down(1, gam()), 0, hom2, 0) ;
514 
515  Vector hom = hom1 + hom2 ;
516 
517  Scalar div_omega = 1.*contract( qq_uu, 0, 1, (1.*hom1 + 1.*hom2).derive_cov(gam()), 0, 1) ;
518  div_omega.inc_dzpuis() ;
519 
520  //---------------------
521 
522 
523  // Two-dimensional Ricci scalar
524  //-----------------------------
525 
526  Scalar rr (mp) ;
527  rr = mp.r ;
528 
529  Scalar Ricci_conf(mp) ;
530  Ricci_conf = 2 / rr / rr ;
531  Ricci_conf.std_spectral_base() ;
532 
533  Scalar Ricci(mp) ;
534  Scalar log_psi (log(psi())) ;
535  log_psi.std_spectral_base() ;
536  Ricci = pow(psi(), -4) * (Ricci_conf - 4*log_psi.lapang())/rr /rr ;
537  Ricci.std_spectral_base() ;
538  Ricci.inc_dzpuis(4) ;
539  //-------------------------------
540 
541  Scalar theta_k = -1/(2*nn()) * (gam().radial_vect().divergence(gam()) -
542  kk_rr + trk() ) ;
543 
544  int nnp = mp.get_mg()->get_np(1) ;
545  int nnt = mp.get_mg()->get_nt(1) ;
546  /*
547  for (int k=0 ; k<nnp ; k++)
548  for (int j=0 ; j<nnt ; j++){
549  cout << theta_k.val_grid_point(1, k, j, 0) << endl ;
550  cout << nn().val_grid_point(1, k, j, 0) << endl ;
551  }
552  */
553 
554 
555 
556  //Source
557  //------
558  Scalar source = div_omega + contract( qq_uu, 0, 1, hom * hom , 0, 1) - 0.5 * Ricci - L_theta;
559  source = source / theta_k ;
560 
561  Scalar tmp = ( source + nn() * kk_rr + rho * contract(gam().radial_vect(), 0,
562  nn().derive_cov(gam()), 0))/(1+rho)
563  - gam().radial_vect()(2) * dnnn(2) - gam().radial_vect()(3) * dnnn(3) ;
564 
565  tmp = tmp / gam().radial_vect()(1) ;
566 
567  // in this case you don't have to substract any value
568 
569  Valeur nn_bound (mp.get_mg()->get_angu()) ;
570 
571  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
572 
573 
574  for (int k=0 ; k<nnp ; k++)
575  for (int j=0 ; j<nnt ; j++)
576  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
577 
578  nn_bound.std_base_scal() ;
579 
580  return nn_bound ;
581 
582 }
583 
584 
585 
586 const Valeur Isol_hor::boundary_nn_Dir_eff(double cc)const {
587 
588  Scalar tmp(mp) ;
589 
590  tmp = - nn().derive_cov(ff)(1) ;
591 
592  // rho = 1 is the standart case
593  double rho = 1. ;
594  tmp.dec_dzpuis(2) ;
595  tmp += cc * (rho -1)*nn() ;
596  tmp = tmp / (rho*cc) ;
597 
598  tmp = tmp - 1. ;
599 
600  // We have substracted 1, since we solve for zero condition at infinity
601  //and then we add 1 to the solution
602 
603  int nnp = mp.get_mg()->get_np(1) ;
604  int nnt = mp.get_mg()->get_nt(1) ;
605 
606  Valeur nn_bound (mp.get_mg()->get_angu()) ;
607 
608  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
609 
610  for (int k=0 ; k<nnp ; k++)
611  for (int j=0 ; j<nnt ; j++)
612  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
613 
614  nn_bound.std_base_scal() ;
615 
616  return nn_bound ;
617 
618 }
619 
620 
621 
622 const Valeur Isol_hor::boundary_nn_Neu_eff(double cc)const {
623 
624  Scalar tmp = - cc * nn() ;
625  // Scalar tmp = - nn()/psi()*psi().dsdr() ;
626 
627  // in this case you don't have to substract any value
628 
629  int nnp = mp.get_mg()->get_np(1) ;
630  int nnt = mp.get_mg()->get_nt(1) ;
631 
632  Valeur nn_bound (mp.get_mg()->get_angu()) ;
633 
634  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
635 
636  for (int k=0 ; k<nnp ; k++)
637  for (int j=0 ; j<nnt ; j++)
638  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
639 
640  nn_bound.std_base_scal() ;
641 
642  return nn_bound ;
643 
644 }
645 
646 
647 const Valeur Isol_hor::boundary_nn_Dir(double cc)const {
648 
649  Scalar rho(mp);
650  rho = 0. ; // 0 is the standard case
651 
652  Scalar tmp(mp) ;
653  tmp = cc ;
654 
655  /*
656  if (b_tilde().val_grid_point(1, 0, 0, 0) < 0.08)
657  tmp = 0.25 ;
658  else {
659  cout << "OK" << endl ;
660  // des_profile(nn(), 0, 20, M_PI/2, M_PI) ;
661  rho = 5. ;
662  tmp = b_tilde()*psi()*psi() ;
663  }
664  */
665 
666  //tmp = 1./(2*psi()) ;
667  // tmp = - psi() * nn().dsdr() / (psi().dsdr()) ;
668 
669  // We have substracted 1, since we solve for zero condition at infinity
670  //and then we add 1 to the solution
671 
672  tmp = (tmp + rho * nn())/(1 + rho) ;
673 
674  tmp = tmp - 1 ;
675 
676  int nnp = mp.get_mg()->get_np(1) ;
677  int nnt = mp.get_mg()->get_nt(1) ;
678 
679  Valeur nn_bound (mp.get_mg()->get_angu()) ;
680 
681  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
682 
683  for (int k=0 ; k<nnp ; k++)
684  for (int j=0 ; j<nnt ; j++)
685  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
686 
687  nn_bound.std_base_scal() ;
688 
689  return nn_bound ;
690 
691 }
692 
694 
695  int nnp = mp.get_mg()->get_np(1) ;
696  int nnt = mp.get_mg()->get_nt(1) ;
697 
698 
699  //Preliminaries
700  //-------------
701 
702  //Metrics on S^2
703  //--------------
704 
705  Sym_tensor q_uu = gam_uu() - gam().radial_vect() * gam().radial_vect() ;
706  Sym_tensor q_dd = q_uu.up_down(gam()) ;
707 
708  Scalar det_q = q_dd(2,2) * q_dd(3,3) - q_dd(2,3) * q_dd(3,2) ;
709  Scalar square_q (pow(det_q, 0.5)) ;
710  square_q.std_spectral_base() ;
711 
712  Sym_tensor qhat_uu = square_q * q_uu ;
713 
714  Sym_tensor ffr_uu = ff.con() - ff.radial_vect() * ff.radial_vect() ;
715  Sym_tensor ffr_dd = ff.cov() - ff.radial_vect().down(0, ff) * ff.radial_vect().down(0,ff) ;
716 
717  Sym_tensor h_uu = qhat_uu - ffr_uu ;
718 
719  //------------------
720 
721 
722  // Source of the "round" laplacian:
723  //---------------------------------
724 
725  //Source 1: "Divergence" term
726  //---------------------------
727  Vector kqs = contract(k_dd(), 1, gam().radial_vect(), 0 ) ;
728  kqs = contract( q_uu.down(0, gam()), 1, kqs, 0) ;
729  Scalar div_kqs = contract( q_uu, 0, 1, kqs.derive_cov(gam()), 0, 1) ;
730 
731  Scalar source1 = div_kqs ;
732  source1 *= square_q ;
733 
734 
735  // Source 2: correction term
736  //--------------------------
737 
738  Vector corr = - contract(h_uu, 1, nn().derive_cov(ff), 0) / nn() ;
739  Scalar source2 = contract( ffr_dd, 0, 1, corr.derive_con(ff), 0, 1 ) ;
740 
741 
742  // Source 3: (physical) divergence of Omega
743  //-----------------------------------------
744 
745  Scalar div_omega(mp) ;
746 
747  //Source from (L_l\theta_k=0)
748 
749  Scalar L_theta (mp) ;
750  L_theta = 0. ;
751  L_theta.std_spectral_base() ;
752 
753  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
754  , k_dd(), 0, 1 ) ;
755 
756  // Scalar kappa (mp) ;
757  // kappa = kappa_hor() ;
758  // kappa.std_spectral_base() ;
759  // kappa.set_dzpuis(2) ;
760 
761 
762  Scalar kappa = contract(gam().radial_vect(), 0, nn().derive_cov(gam()), 0) - nn() * kk_rr ;
763  Scalar theta_k = -1/(2*nn()) * (gam().radial_vect().divergence(gam()) -
764  kk_rr + trk() ) ;
765 
766  Sym_tensor qqq = gam_uu() - gam().radial_vect() * gam().radial_vect() ;
767 
768  Vector hom = nn().derive_cov(met_gamt) - contract( k_dd(), 1, gam().radial_vect() , 0) ;
769  hom = contract(qqq.down(1, gam()), 0, hom, 0) ;
770 
771  Scalar rr(mp) ;
772  rr = mp.r ;
773 
774  Scalar Ricci_conf = 2 / rr /rr ;
775  Ricci_conf.std_spectral_base() ;
776 
777  Scalar log_psi (log(psi())) ;
778  log_psi.std_spectral_base() ;
779  Scalar Ricci = pow(psi(), -4) * (Ricci_conf - 4*log_psi.lapang())/rr /rr ;
780  Ricci.std_spectral_base() ;
781  Ricci.inc_dzpuis(4) ;
782 
783 
784  div_omega = L_theta - contract(qqq, 0, 1, hom * hom, 0, 1) + 0.5 * Ricci
785  + theta_k * kappa ;
786  div_omega.dec_dzpuis() ;
787 
788  //End of Source from (L_l\theta_k=0)
789  //----------------------------------
790 
791  div_omega = 0. ;
792  // div_omega = 0.01*log(1/(2*psi())) ;
793  div_omega.std_spectral_base() ;
794  div_omega.inc_dzpuis(3) ;
795 
796 
797  //Construction of source3 (correction of div_omega by the square root of the determinant)
798  Scalar source3 = div_omega ;
799  source3 *= square_q ;
800 
801 
802  //"Correction" to the div_omega term (no Y_00 component must be present)
803 
804  double corr_const = mp.integrale_surface(source3, radius + 1e-15)/(4. * M_PI) ;
805  cout << "Constant part of div_omega = " << corr_const <<endl ;
806 
807  Scalar corr_div_omega (mp) ;
808  corr_div_omega = corr_const ;
809  corr_div_omega.std_spectral_base() ;
810  corr_div_omega.set_dzpuis(3) ;
811 
812  source3 -= corr_div_omega ;
813 
814 
815 
816  //Source
817  //------
818 
819  Scalar source = source1 + source2 + source3 ;
820 
821 
822 
823  //Resolution of round angular laplacian
824  //-------------------------------------
825 
826  Scalar source_tmp = source ;
827 
828  Scalar logn (mp) ;
829  logn = source.poisson_angu() ;
830 
831  double cc = 0.2 ; // Integration constant
832 
833  Scalar lapse (mp) ;
834  lapse = (exp(logn)*cc) ;
835  lapse.std_spectral_base() ;
836 
837 
838  //Test of Divergence of Omega
839  //---------------------------
840 
841  ofstream err ;
842  err.open ("err_laplBC.d", ofstream::app ) ;
843 
844  hom = nn().derive_cov(gam())/nn()
845  - contract( k_dd(), 1, gam().radial_vect() , 0) ;
846  hom = contract(q_uu.down(1, gam()), 0, hom, 0) ;
847 
848  Scalar div_omega_test = contract( q_uu, 0, 1, hom.derive_cov(gam()), 0, 1) ;
849 
850 
851  Scalar err_lapl = div_omega_test - div_omega ;
852 
853  double max_err = err_lapl.val_grid_point(1, 0, 0, 0) ;
854  double min_err = err_lapl.val_grid_point(1, 0, 0, 0) ;
855 
856  for (int k=0 ; k<nnp ; k++)
857  for (int j=0 ; j<nnt ; j++){
858  if (err_lapl.val_grid_point(1, k, j, 0) > max_err)
859  max_err = err_lapl.val_grid_point(1, k, j, 0) ;
860  if (err_lapl.val_grid_point(1, k, j, 0) < min_err)
861  min_err = err_lapl.val_grid_point(1, k, j, 0) ;
862  }
863 
864  err << mer << " " << max_err << " " << min_err << endl ;
865 
866  err.close() ;
867 
868 
869 
870  //Construction of the Valeur
871  //--------------------------
872 
873  lapse = lapse - 1. ;
874 
875  // int nnp = mp.get_mg()->get_np(1) ;
876  // int nnt = mp.get_mg()->get_nt(1) ;
877 
878  Valeur nn_bound (mp.get_mg()->get_angu()) ;
879 
880  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
881 
882  for (int k=0 ; k<nnp ; k++)
883  for (int j=0 ; j<nnt ; j++)
884  nn_bound.set(0, k, j, 0) = lapse.val_grid_point(1, k, j, 0) ;
885 
886  nn_bound.std_base_scal() ;
887 
888  return nn_bound ;
889 
890 }
891 
892 
893 // Component r of boundary value of beta (using expression in terms
894 // of radial vector)
895 // --------------------------------------
896 
898 
899  Scalar tmp (mp) ;
900 
901  tmp = nn() * radial_vect_hor()(1) ;
902 
903  Base_val base_tmp (radial_vect_hor()(1).get_spectral_va().base) ;
904 
905  int nnp = mp.get_mg()->get_np(1) ;
906  int nnt = mp.get_mg()->get_nt(1) ;
907 
908  Valeur bnd_beta_r (mp.get_mg()->get_angu()) ;
909 
910  bnd_beta_r = 1. ; // Juste pour affecter dans espace des configs ;
911 
912  for (int k=0 ; k<nnp ; k++)
913  for (int j=0 ; j<nnt ; j++)
914  bnd_beta_r.set(1, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
915 
916  bnd_beta_r.set_base_r(0, base_tmp.get_base_r(0)) ;
917  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
918  bnd_beta_r.set_base_r(l, base_tmp.get_base_r(l)) ;
919  bnd_beta_r.set_base_r((*mp.get_mg()).get_nzone()-1, base_tmp.get_base_r((*mp.get_mg()).get_nzone()-1)) ;
920  bnd_beta_r.set_base_t(tmp.get_spectral_va().get_base().get_base_t(1)) ;
921  bnd_beta_r.set_base_p(tmp.get_spectral_va().get_base().get_base_p(1)) ;
922 
923 // bnd_beta_r.set_base(*(mp.get_mg()->std_base_vect_spher()[0])) ;
924 
925  cout << "norme de lim_vr" << endl << norme(bnd_beta_r) << endl ;
926  cout << "bases" << endl << bnd_beta_r.base << endl ;
927 
928  return bnd_beta_r ;
929 
930 
931 }
932 
933 
934 // Component theta of boundary value of beta (using expression in terms
935 // of radial vector)
936 // ------------------------------------------
937 
939 
940  Scalar tmp(mp) ;
941 
942  tmp = nn() * radial_vect_hor()(2) ;
943  Base_val base_tmp (radial_vect_hor()(2).get_spectral_va().base) ;
944 
945 
946  int nnp = mp.get_mg()->get_np(1) ;
947  int nnt = mp.get_mg()->get_nt(1) ;
948 
949  Valeur bnd_beta_theta (mp.get_mg()->get_angu()) ;
950 
951  bnd_beta_theta = 1. ; // Juste pour affecter dans espace des configs ;
952 
953  for (int k=0 ; k<nnp ; k++)
954  for (int j=0 ; j<nnt ; j++)
955  bnd_beta_theta.set(1, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
956 
957 // bnd_beta_theta.set_base(*(mp.get_mg()->std_base_vect_spher()[1])) ;
958 
959  bnd_beta_theta.set_base_r(0, base_tmp.get_base_r(0)) ;
960  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
961  bnd_beta_theta.set_base_r(l, base_tmp.get_base_r(l)) ;
962  bnd_beta_theta.set_base_r((*mp.get_mg()).get_nzone()-1, base_tmp.get_base_r((*mp.get_mg()).get_nzone()-1)) ;
963  bnd_beta_theta.set_base_t(base_tmp.get_base_t(1)) ;
964  bnd_beta_theta.set_base_p(base_tmp.get_base_p(1)) ;
965 
966  cout << "bases" << endl << bnd_beta_theta.base << endl ;
967 
968 
969  return bnd_beta_theta ;
970 
971 
972 }
973 
974 // Component phi of boundary value of beta (using expression in terms
975 // of radial vector)
976 // -----------------------------------------------------------
977 
978 const Valeur Isol_hor::boundary_beta_phi(double om)const {
979 
980  Scalar tmp (mp) ;
981 
982  Scalar ang_vel(mp) ;
983  ang_vel = om ;
984  ang_vel.std_spectral_base() ;
985  ang_vel.mult_rsint() ;
986 
987  tmp = nn() * radial_vect_hor()(3) - ang_vel ;
988  Base_val base_tmp (ang_vel.get_spectral_va().base) ;
989 
990  int nnp = mp.get_mg()->get_np(1) ;
991  int nnt = mp.get_mg()->get_nt(1) ;
992 
993  Valeur bnd_beta_phi (mp.get_mg()->get_angu()) ;
994 
995  bnd_beta_phi = 1. ; // Juste pour affecter dans espace des configs ;
996 
997  for (int k=0 ; k<nnp ; k++)
998  for (int j=0 ; j<nnt ; j++)
999  bnd_beta_phi.set(1, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1000 
1001  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
1002 
1003 // bnd_beta_phi.set_base(*(mp.get_mg()->std_base_vect_spher()[2])) ;
1004 
1005  bnd_beta_phi.set_base_r(0, base_tmp.get_base_r(0)) ;
1006  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
1007  bnd_beta_phi.set_base_r(l, base_tmp.get_base_r(l)) ;
1008  bnd_beta_phi.set_base_r((*mp.get_mg()).get_nzone()-1, base_tmp.get_base_r((*mp.get_mg()).get_nzone()-1)) ;
1009  bnd_beta_phi.set_base_t(base_tmp.get_base_t(1)) ;
1010  bnd_beta_phi.set_base_p(base_tmp.get_base_p(1)) ;
1011 
1012 // cout << "bound beta_phi" << endl << bnd_beta_phi << endl ;
1013 
1014  return bnd_beta_phi ;
1015 
1016 }
1017 
1018 // Component x of boundary value of beta
1019 //--------------------------------------
1020 
1021 const Valeur Isol_hor:: boundary_beta_x(double om)const {
1022 
1023  // Les alignemenents pour le signe des CL.
1024  double orientation = mp.get_rot_phi() ;
1025  assert ((orientation == 0) || (orientation == M_PI)) ;
1026  int aligne = (orientation == 0) ? 1 : -1 ;
1027 
1028  int nnp = mp.get_mg()->get_np(1) ;
1029  int nnt = mp.get_mg()->get_nt(1) ;
1030 
1031  Vector tmp_vect = nn() * gam().radial_vect() ;
1032  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1033 
1034  //Isol_hor boundary conditions
1035 
1036  Valeur lim_x (mp.get_mg()->get_angu()) ;
1037 
1038  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
1039 
1040  Mtbl ya_mtbl (mp.get_mg()) ;
1041  ya_mtbl.set_etat_qcq() ;
1042  ya_mtbl = mp.ya ;
1043 
1044  Scalar cosp (mp) ;
1045  cosp = mp.cosp ;
1046  Scalar cost (mp) ;
1047  cost = mp.cost ;
1048  Scalar sinp (mp) ;
1049  sinp = mp.sinp ;
1050  Scalar sint (mp) ;
1051  sint = mp.sint ;
1052 
1053  Scalar dep_phi (mp) ;
1054  dep_phi = 0.0*sint*cosp ;
1055 
1056  for (int k=0 ; k<nnp ; k++)
1057  for (int j=0 ; j<nnt ; j++)
1058  lim_x.set(0, k, j, 0) = aligne * om * ya_mtbl(1, k, j, 0) * (1 +
1059  dep_phi.val_grid_point(1, k, j, 0))
1060  + tmp_vect(1).val_grid_point(1, k, j, 0) ;
1061 
1062  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1063 
1064  return lim_x ;
1065 }
1066 
1067 
1068 // Component y of boundary value of beta
1069 //--------------------------------------
1070 
1071 const Valeur Isol_hor:: boundary_beta_y(double om)const {
1072 
1073  // Les alignemenents pour le signe des CL.
1074  double orientation = mp.get_rot_phi() ;
1075  assert ((orientation == 0) || (orientation == M_PI)) ;
1076  int aligne = (orientation == 0) ? 1 : -1 ;
1077 
1078 
1079  int nnp = mp.get_mg()->get_np(1) ;
1080  int nnt = mp.get_mg()->get_nt(1) ;
1081 
1082  Vector tmp_vect = nn() * gam().radial_vect() ;
1083  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1084 
1085  //Isol_hor boundary conditions
1086 
1087  Valeur lim_y (mp.get_mg()->get_angu()) ;
1088 
1089  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
1090 
1091  Mtbl xa_mtbl (mp.get_mg()) ;
1092  xa_mtbl.set_etat_qcq() ;
1093  xa_mtbl = mp.xa ;
1094 
1095  Scalar cosp (mp) ;
1096  cosp = mp.cosp ;
1097  Scalar cost (mp) ;
1098  cost = mp.cost ;
1099  Scalar sinp (mp) ;
1100  sinp = mp.sinp ;
1101  Scalar sint (mp) ;
1102  sint = mp.sint ;
1103 
1104  Scalar dep_phi (mp) ;
1105  dep_phi = 0.0*sint*cosp ;
1106 
1107  for (int k=0 ; k<nnp ; k++)
1108  for (int j=0 ; j<nnt ; j++)
1109  lim_y.set(0, k, j, 0) = - aligne * om * xa_mtbl(1, k, j, 0) *(1 +
1110  dep_phi.val_grid_point(1, k, j, 0))
1111  + tmp_vect(2).val_grid_point(1, k, j, 0) ;
1112 
1113  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1114 
1115  return lim_y ;
1116 }
1117 
1118 // Component z of boundary value of beta
1119 //--------------------------------------
1120 
1122 
1123  int nnp = mp.get_mg()->get_np(1) ;
1124  int nnt = mp.get_mg()->get_nt(1) ;
1125 
1126  Vector tmp_vect = nn() * gam().radial_vect() ;
1127  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1128 
1129  //Isol_hor boundary conditions
1130 
1131  Valeur lim_z (mp.get_mg()->get_angu()) ;
1132 
1133  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
1134 
1135  for (int k=0 ; k<nnp ; k++)
1136  for (int j=0 ; j<nnt ; j++)
1137  lim_z.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
1138 
1139  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1140 
1141  return lim_z ;
1142 }
1143 
1145 
1146  Valeur lim_x (mp.get_mg()->get_angu()) ;
1147  lim_x = boost_x ; // Juste pour affecter dans espace des configs ;
1148 
1149  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1150 
1151  return lim_x ;
1152 }
1153 
1154 
1156 
1157  Valeur lim_z (mp.get_mg()->get_angu()) ;
1158  lim_z = boost_z ; // Juste pour affecter dans espace des configs ;
1159 
1160  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1161 
1162  return lim_z ;
1163 }
1164 
1165 
1166 // Neumann boundary condition for b_tilde
1167 //---------------------------------------
1168 
1170 
1171  // Introduce 2-trace gamma tilde dot
1172 
1173  Vector s_tilde = met_gamt.radial_vect() ;
1174 
1175  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1176 
1177  //des_profile(hh_tilde, 1.00001, 10, M_PI/2., 0., "H_tilde") ;
1178 
1179  Scalar tmp (mp) ;
1180 
1181  tmp = + b_tilde() * hh_tilde - 2 * ( s_tilde(2) * b_tilde().derive_cov(ff)(2)
1182  + s_tilde(3) * b_tilde().derive_cov(ff)(3) ) ;
1183 
1184  Scalar constant (mp) ;
1185  constant = 0. ;
1186  constant.std_spectral_base() ;
1187  constant.inc_dzpuis(2) ;
1188 
1189  tmp += constant ;
1190  tmp = tmp / (2 * s_tilde(1) ) ;
1191 
1192  // in this case you don't have to substract any value
1193 
1194  Valeur b_tilde_bound (mp.get_mg()->get_angu() ) ;
1195 
1196  int nnp = mp.get_mg()->get_np(1) ;
1197  int nnt = mp.get_mg()->get_nt(1) ;
1198 
1199  b_tilde_bound = 1 ;
1200 
1201  for (int k=0 ; k<nnp ; k++)
1202  for (int j=0 ; j<nnt ; j++)
1203  b_tilde_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1204 
1205  b_tilde_bound.std_base_scal() ;
1206 
1207  return b_tilde_bound ;
1208 
1209 }
1210 
1211 
1213 
1214  Vector s_tilde = met_gamt.radial_vect() ;
1215 
1216  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1217 
1218  Scalar tmp = 2 * contract (s_tilde, 0, b_tilde().derive_cov(ff) , 0) ;
1219 
1220  Scalar constant (mp) ;
1221  constant = -1. ;
1222  constant.std_spectral_base() ;
1223  constant.inc_dzpuis(2) ;
1224 
1225  tmp -= constant ;
1226 
1227  tmp = tmp / hh_tilde ;
1228 
1229  // des_profile(tmp, 1.00001, 10, M_PI/2., 0., "tmp") ;
1230 
1231  // We have substracted 1, since we solve for zero condition at infinity
1232  //and then we add 1 to the solution
1233 
1234  Valeur b_tilde_bound (mp.get_mg()->get_angu() ) ;
1235 
1236  int nnp = mp.get_mg()->get_np(1) ;
1237  int nnt = mp.get_mg()->get_nt(1) ;
1238 
1239  b_tilde_bound = 1 ;
1240 
1241  for (int k=0 ; k<nnp ; k++)
1242  for (int j=0 ; j<nnt ; j++)
1243  b_tilde_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1244 
1245  b_tilde_bound.std_base_scal() ;
1246 
1247  return b_tilde_bound ;
1248 
1249 }
1250 
1251 const Vector Isol_hor::vv_bound_cart(double om) const{
1252 
1253  // Preliminaries
1254  //--------------
1255 
1256  // HH_tilde
1257  Vector s_tilde = met_gamt.radial_vect() ;
1258 
1259  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1260  hh_tilde.dec_dzpuis(2) ;
1261 
1262  // Tangential part of the shift
1263  Vector tmp_vect = b_tilde() * s_tilde ;
1264 
1265  Vector VV = b_tilde() * s_tilde - beta() ;
1266 
1267  // "Acceleration" term V^a \tilde{D}_a ln M
1268  Scalar accel = 2 * contract( VV, 0, contract( s_tilde, 0, s_tilde.down(0,
1269  met_gamt).derive_cov(met_gamt), 1), 0) ;
1270 
1271 
1272  // Divergence of V^a
1273  Sym_tensor qq_spher = met_gamt.con() - s_tilde * s_tilde ;
1274  Scalar div_VV = contract( qq_spher.down(0, met_gamt), 0, 1, VV.derive_cov(met_gamt), 0, 1) ;
1275 
1276 
1277  Scalar cosp (mp) ;
1278  cosp = mp.cosp ;
1279  Scalar cost (mp) ;
1280  cost = mp.cost ;
1281  Scalar sinp (mp) ;
1282  sinp = mp.sinp ;
1283  Scalar sint (mp) ;
1284  sint = mp.sint ;
1285 
1286  Scalar dep_phi (mp) ;
1287  dep_phi = 0.0*sint*cosp ;
1288 
1289  // Les alignemenents pour le signe des CL.
1290  double orientation = mp.get_rot_phi() ;
1291  assert ((orientation == 0) || (orientation == M_PI)) ;
1292  int aligne = (orientation == 0) ? 1 : -1 ;
1293 
1294  Vector angular (mp, CON, mp.get_bvect_cart()) ;
1295  Scalar yya (mp) ;
1296  yya = mp.ya ;
1297  Scalar xxa (mp) ;
1298  xxa = mp.xa ;
1299 
1300  angular.set(1) = - aligne * om * yya * (1 + dep_phi) ;
1301  angular.set(2) = aligne * om * xxa * (1 + dep_phi) ;
1302  angular.set(3).annule_hard() ;
1303 
1304 
1305  angular.set(1).set_spectral_va()
1306  .set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1307  angular.set(2).set_spectral_va()
1308  .set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1309  angular.set(3).set_spectral_va()
1310  .set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1311 
1312  angular.change_triad(mp.get_bvect_spher()) ;
1313 
1314 /*
1315  Scalar ang_vel (mp) ;
1316  ang_vel = om * (1 + dep_phi) ;
1317  ang_vel.std_spectral_base() ;
1318  ang_vel.mult_rsint() ;
1319 */
1320 
1321  Scalar kss (mp) ;
1322  kss = - 0.2 ;
1323  // kss.std_spectral_base() ;
1324  // kss.inc_dzpuis(2) ;
1325 
1326 
1327  //kss from from L_lN=0 condition
1328  //------------------------------
1329  kss = - 0.15 / nn() ;
1330  kss.inc_dzpuis(2) ;
1331  kss += contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0) / nn() ;
1332 
1333 
1334  cout << "kappa_hor = " << kappa_hor() << endl ;
1335 
1336  /*
1337  // Apparent horizon condition
1338  //---------------------------
1339  kss = trk() ;
1340  kss -= (4* contract(psi().derive_cov(met_gamt), 0, met_gamt.radial_vect(),
1341  0)/psi() +
1342  contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) /
1343  (psi() * psi()) ;
1344  */
1345 
1346 
1347  // beta^r component
1348  //-----------------
1349  double rho = 5. ; // rho>1 ; rho=1 "pure Dirichlet" version
1350  Scalar beta_r_corr = (rho - 1) * b_tilde() * hh_tilde;
1351  beta_r_corr.inc_dzpuis(2) ;
1352  Scalar nn_KK = nn() * trk() ;
1353  nn_KK.inc_dzpuis(2) ;
1354 
1355  beta_r_corr.set_dzpuis(2) ;
1356  nn_KK.set_dzpuis(2) ;
1357  accel.set_dzpuis(2) ;
1358  div_VV.set_dzpuis(2) ;
1359 
1360  Scalar b_tilde_new (mp) ;
1361  b_tilde_new = 2 * contract(s_tilde, 0, b_tilde().derive_cov(ff), 0)
1362  + beta_r_corr
1363  - 3 * nn() * kss + nn_KK + accel + div_VV ;
1364 
1365  b_tilde_new = b_tilde_new / (hh_tilde * rho) ;
1366 
1367  b_tilde_new.dec_dzpuis(2) ;
1368 
1369  tmp_vect.set(1) = b_tilde_new * s_tilde(1) + angular(1) ;
1370  tmp_vect.set(2) = b_tilde_new * s_tilde(2) + angular(2) ;
1371  tmp_vect.set(3) = b_tilde_new * s_tilde(3) + angular(3) ;
1372 
1373  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1374 
1375  return tmp_vect ;
1376 }
1377 
1378 
1379 const Vector Isol_hor::vv_bound_cart_bin(double, int ) const{
1380  /*
1381  // Les alignemenents pour le signe des CL.
1382  double orientation = mp.get_rot_phi() ;
1383  assert ((orientation == 0) || (orientation == M_PI)) ;
1384  int aligne = (orientation == 0) ? 1 : -1 ;
1385 
1386  Vector angular (mp, CON, mp.get_bvect_cart()) ;
1387  Scalar yya (mp) ;
1388  yya = mp.ya ;
1389  Scalar xxa (mp) ;
1390  xxa = mp.xa ;
1391 
1392  angular.set(1) = aligne * om * yya ;
1393  angular.set(2) = - aligne * om * xxa ;
1394  angular.set(3).annule_hard() ;
1395 
1396  angular.set(1).set_spectral_va()
1397  .set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1398  angular.set(2).set_spectral_va()
1399  .set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1400  angular.set(3).set_spectral_va()
1401  .set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1402 
1403  angular.change_triad(mp.get_bvect_spher()) ;
1404 
1405  // HH_tilde
1406  Vector s_tilde = met_gamt.radial_vect() ;
1407 
1408  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1409  hh_tilde.dec_dzpuis(2) ;
1410 
1411  Scalar btilde = b_tilde() - contract(angular, 0, s_tilde.up_down(met_gamt), 0) ;
1412 
1413  // Tangential part of the shift
1414  Vector tmp_vect = btilde * s_tilde ;
1415 
1416  // Value of kss
1417  // --------------
1418 
1419  Scalar kss (mp) ;
1420  kss = - 0.26 ;
1421  kss.std_spectral_base() ;
1422  kss.inc_dzpuis(2) ;
1423 
1424  // Apparent horizon boundary condition
1425  // -----------------------------------
1426 
1427 
1428  // kss frome a fich
1429 
1430  // Construction of the binary
1431  // --------------------------
1432 
1433  char* name_fich = "/home/francois/resu/bin_hor/Test/b11_9x9x8/bin.dat" ;
1434 
1435  FILE* fich_s = fopen(name_fich, "r") ;
1436  Mg3d grid_s (fich_s) ;
1437  Map_af map_un_s (grid_s, fich_s) ;
1438  Map_af map_deux_s (grid_s, fich_s) ;
1439  Bin_hor bin (map_un_s, map_deux_s, fich_s) ;
1440  fclose(fich_s) ;
1441 
1442  // Inititialisation of fields :
1443  // ----------------------------
1444 
1445  bin.set(1).n_comp_import (bin(2)) ;
1446  bin.set(1).psi_comp_import (bin(2)) ;
1447  bin.set(2).n_comp_import (bin(1)) ;
1448  bin.set(2).psi_comp_import (bin(1)) ;
1449  bin.decouple() ;
1450  bin.extrinsic_curvature() ;
1451 
1452  kss = contract(bin(jj).get_k_dd(), 0, 1, bin(jj).get_gam().radial_vect()*
1453  bin(jj).get_gam().radial_vect(), 0, 1) ;
1454 
1455 
1456  cout << "kappa_hor = " << kappa_hor() << endl ;
1457 
1458  // beta^r component
1459  //-----------------
1460  double rho = 5. ; // rho=0 "pure Dirichlet" version
1461  Scalar beta_r_corr = rho * btilde * hh_tilde;
1462  beta_r_corr.inc_dzpuis(2) ;
1463  Scalar nn_KK = nn() * trk() ;
1464  nn_KK.inc_dzpuis(2) ;
1465 
1466  beta_r_corr.set_dzpuis(2) ;
1467  nn_KK.set_dzpuis(2) ;
1468 
1469  Scalar b_tilde_new (mp) ;
1470  b_tilde_new = 2 * contract(s_tilde, 0, btilde.derive_cov(ff), 0)
1471  + beta_r_corr - 3 * nn() * kss + nn_KK ;
1472 
1473  b_tilde_new = b_tilde_new / (hh_tilde * (rho+1)) ;
1474 
1475  b_tilde_new.dec_dzpuis(2) ;
1476 
1477  tmp_vect.set(1) = b_tilde_new * s_tilde(1) + angular(1) ;
1478  tmp_vect.set(2) = b_tilde_new * s_tilde(2) + angular(2) ;
1479  tmp_vect.set(3) = b_tilde_new * s_tilde(3) + angular(3) ;
1480 
1481  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1482 
1483  return tmp_vect ;
1484  */
1485  Vector pipo(mp, CON, mp.get_bvect_cart()) ;
1486  return pipo ;
1487 }
1488 
1489 
1490 // Component x of boundary value of V^i
1491 //-------------------------------------
1492 
1493 const Valeur Isol_hor:: boundary_vv_x(double om)const {
1494 
1495  int nnp = mp.get_mg()->get_np(1) ;
1496  int nnt = mp.get_mg()->get_nt(1) ;
1497 
1498  //Isol_hor boundary conditions
1499 
1500  Valeur lim_x (mp.get_mg()->get_angu()) ;
1501 
1502  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
1503 
1504  Scalar vv_x = vv_bound_cart(om)(1) ;
1505 
1506  for (int k=0 ; k<nnp ; k++)
1507  for (int j=0 ; j<nnt ; j++)
1508  lim_x.set(0, k, j, 0) = vv_x.val_grid_point(1, k, j, 0) ;
1509 
1510  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1511 
1512  return lim_x ;
1513 
1514 
1515 }
1516 
1517 
1518 // Component y of boundary value of V^i
1519 //--------------------------------------
1520 
1521 const Valeur Isol_hor:: boundary_vv_y(double om)const {
1522 
1523  int nnp = mp.get_mg()->get_np(1) ;
1524  int nnt = mp.get_mg()->get_nt(1) ;
1525 
1526  // Isol_hor boundary conditions
1527 
1528  Valeur lim_y (mp.get_mg()->get_angu()) ;
1529 
1530  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
1531 
1532  Scalar vv_y = vv_bound_cart(om)(2) ;
1533 
1534  for (int k=0 ; k<nnp ; k++)
1535  for (int j=0 ; j<nnt ; j++)
1536  lim_y.set(0, k, j, 0) = vv_y.val_grid_point(1, k, j, 0) ;
1537 
1538  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1539 
1540  return lim_y ;
1541 }
1542 
1543 
1544 // Component z of boundary value of V^i
1545 //-------------------------------------
1546 
1547 const Valeur Isol_hor:: boundary_vv_z(double om)const {
1548 
1549  int nnp = mp.get_mg()->get_np(1) ;
1550  int nnt = mp.get_mg()->get_nt(1) ;
1551 
1552  // Isol_hor boundary conditions
1553 
1554  Valeur lim_z (mp.get_mg()->get_angu()) ;
1555 
1556  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
1557 
1558  Scalar vv_z = vv_bound_cart(om)(3) ;
1559 
1560  for (int k=0 ; k<nnp ; k++)
1561  for (int j=0 ; j<nnt ; j++)
1562  lim_z.set(0, k, j, 0) = vv_z.val_grid_point(1, k, j, 0) ;
1563 
1564  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1565 
1566  return lim_z ;
1567 }
1568 
1569 // Component x of boundary value of V^i
1570 //-------------------------------------
1571 
1572 const Valeur Isol_hor:: boundary_vv_x_bin(double om, int jj)const {
1573 
1574  int nnp = mp.get_mg()->get_np(1) ;
1575  int nnt = mp.get_mg()->get_nt(1) ;
1576 
1577  //Isol_hor boundary conditions
1578 
1579  Valeur lim_x (mp.get_mg()->get_angu()) ;
1580 
1581  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
1582 
1583  Scalar vv_x = vv_bound_cart_bin(om, jj)(1) ;
1584 
1585  for (int k=0 ; k<nnp ; k++)
1586  for (int j=0 ; j<nnt ; j++)
1587  lim_x.set(0, k, j, 0) = vv_x.val_grid_point(1, k, j, 0) ;
1588 
1589  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1590 
1591  return lim_x ;
1592 
1593 
1594 }
1595 
1596 // Component y of boundary value of V^i
1597 //--------------------------------------
1598 
1599 const Valeur Isol_hor:: boundary_vv_y_bin(double om, int jj)const {
1600 
1601  int nnp = mp.get_mg()->get_np(1) ;
1602  int nnt = mp.get_mg()->get_nt(1) ;
1603 
1604  // Isol_hor boundary conditions
1605 
1606  Valeur lim_y (mp.get_mg()->get_angu()) ;
1607 
1608  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
1609 
1610  Scalar vv_y = vv_bound_cart_bin(om, jj)(2) ;
1611 
1612  for (int k=0 ; k<nnp ; k++)
1613  for (int j=0 ; j<nnt ; j++)
1614  lim_y.set(0, k, j, 0) = vv_y.val_grid_point(1, k, j, 0) ;
1615 
1616  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1617 
1618  return lim_y ;
1619 }
1620 
1621 
1622 // Component z of boundary value of V^i
1623 //-------------------------------------
1624 
1625 const Valeur Isol_hor:: boundary_vv_z_bin(double om, int jj)const {
1626 
1627  int nnp = mp.get_mg()->get_np(1) ;
1628  int nnt = mp.get_mg()->get_nt(1) ;
1629 
1630  // Isol_hor boundary conditions
1631 
1632  Valeur lim_z (mp.get_mg()->get_angu()) ;
1633 
1634  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
1635 
1636  Scalar vv_z = vv_bound_cart_bin(om, jj)(3) ;
1637 
1638  for (int k=0 ; k<nnp ; k++)
1639  for (int j=0 ; j<nnt ; j++)
1640  lim_z.set(0, k, j, 0) = vv_z.val_grid_point(1, k, j, 0) ;
1641 
1642  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1643 
1644  return lim_z ;
1645 }
1646 }
Bases of the spectral expansions.
Definition: base_val.h:322
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:422
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition: base_val.h:400
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:411
const Valeur boundary_vv_z_bin(double om, int hole=0) const
Component z of boundary value of .
Definition: bound_hor.C:1625
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
Definition: bound_hor.C:1521
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution)
Definition: bound_hor.C:174
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial)
Definition: bound_hor.C:231
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif)
Definition: bound_hor.C:622
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
Definition: bound_hor.C:1493
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
Definition: bound_hor.C:1071
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial)
Definition: bound_hor.C:297
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
Definition: bound_hor.C:1021
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook's boundary condition.
Definition: bound_hor.C:489
const Valeur boundary_b_tilde_Neu() const
Neumann boundary condition for b_tilde.
Definition: bound_hor.C:1169
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition: phys_param.C:136
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
Definition: bound_hor.C:1144
const Valeur boundary_beta_r() const
Component r of boundary value of .
Definition: bound_hor.C:897
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
Definition: bound_hor.C:1547
const Valeur boundary_vv_y_bin(double om, int hole=0) const
Component y of boundary value of .
Definition: bound_hor.C:1599
const Valeur boundary_nn_Dir_lapl(int mer=1) const
Dirichlet boundary condition for N fixing the divergence of the connection form .
Definition: bound_hor.C:693
double boost_z
Boost velocity in z-direction.
Definition: isol_hor.h:275
const Valeur boundary_beta_theta() const
Component theta of boundary value of .
Definition: bound_hor.C:938
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
Definition: bound_hor.C:1155
const Valeur boundary_vv_x_bin(double om, int hole=0) const
Component x of boundary value of .
Definition: bound_hor.C:1572
const Vector vv_bound_cart(double om) const
Vector for boundary conditions in cartesian
Definition: bound_hor.C:1251
double kappa_hor() const
Surface gravity
Definition: phys_param.C:218
const Valeur boundary_beta_z() const
Component z of boundary value of .
Definition: bound_hor.C:1121
double radius
Radius of the horizon in LORENE's units.
Definition: isol_hor.h:266
virtual const Sym_tensor & aa_auto() const
Conformal representation of the traceless part of the extrinsic curvature: Returns the value at the ...
Definition: isol_hor.C:503
const Valeur boundary_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:371
const Vector tradial_vect_hor() const
Vector radial normal tilde.
Definition: phys_param.C:116
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial)
Definition: bound_hor.C:267
const Valeur boundary_beta_phi(double om) const
Component phi of boundary value of .
Definition: bound_hor.C:978
const Vector radial_vect_hor() const
Vector radial normal.
Definition: phys_param.C:95
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif)
Definition: bound_hor.C:586
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
Definition: bound_hor.C:647
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:420
double boost_x
Boost velocity in x-direction.
Definition: isol_hor.h:272
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial)
Definition: bound_hor.C:324
const Valeur boundary_b_tilde_Dir() const
Dirichlet boundary condition for b_tilde.
Definition: bound_hor.C:1212
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution)
Definition: bound_hor.C:203
const Vector vv_bound_cart_bin(double om, int hole=0) const
Vector for boundary conditions in cartesian for binary systems.
Definition: bound_hor.C:1379
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Coord cosp
Definition: map.h:724
Coord sint
Definition: map.h:721
Coord ya
Absolute y coordinate.
Definition: map.h:731
Coord r
r coordinate centered on the grid
Definition: map.h:718
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Coord sinp
Definition: map.h:723
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
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 xa
Absolute x coordinate.
Definition: map.h:730
Coord cost
Definition: map.h:722
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric_flat.C:134
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition: metric.C:362
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
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
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Multi-domain array.
Definition: mtbl.h:118
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition: scalar_deriv.C:461
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:200
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...
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
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
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:507
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Metric & gam() const
Induced metric at the current time step (jtime )
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition: valeur.C:836
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:824
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363
void set_base_t(int base_t)
Sets the expansion basis for functions in all domains.
Definition: valeur.C:849
void set_base_p(int base_p)
Sets the expansion basis for functions in all domains.
Definition: valeur.C:862
Tensor field of valence 1.
Definition: vector.h:188
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition: app_hor.h:64