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
29char 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
173namespace 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
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
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
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
647const 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
978const 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
1021const 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
1071const 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
1251const 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
1379const 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
1493const 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
1521const 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
1547const 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
1572const 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
1599const 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
1625const 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
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 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
Coord sinp
Definition map.h:723
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
Coord cost
Definition map.h:722
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
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 .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
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...
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
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
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
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition valeur.h:480
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition valeur.C:836
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
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