LORENE
sym_tensor_tt_etamu.C
1 /*
2  * Methods of class Sym_tensor_tt related to eta and mu
3  *
4  * (see file sym_tensor.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
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 
30 char sym_tensor_tt_etamu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $" ;
31 
32 /*
33  * $Id: sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
34  * $Log: sym_tensor_tt_etamu.C,v $
35  * Revision 1.18 2014/10/13 08:53:44 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.17 2014/10/06 15:13:19 j_novak
39  * Modified #include directives to use c++ syntax.
40  *
41  * Revision 1.16 2006/10/24 13:03:19 j_novak
42  * New methods for the solution of the tensor wave equation. Perhaps, first
43  * operational version...
44  *
45  * Revision 1.15 2005/04/01 14:28:32 j_novak
46  * Members p_eta and p_mu are now defined in class Sym_tensor.
47  *
48  * Revision 1.14 2004/06/04 09:25:58 e_gourgoulhon
49  * Method eta(): eta is no longer computed from h^rr but from khi (in the
50  * case where khi is known).
51  *
52  * Revision 1.13 2004/05/25 15:08:44 f_limousin
53  * Add parameters in argument of the functions update, eta and mu for
54  * the case of a Map_et.
55  *
56  * Revision 1.12 2004/05/24 13:45:29 e_gourgoulhon
57  * Added parameter dzp to method Sym_tensor_tt::update.
58  *
59  * Revision 1.11 2004/05/05 14:24:54 e_gourgoulhon
60  * Corrected a bug in method set_khi_mu: the division of khi by r^2
61  * was ommitted in the case dzp=2 !!!
62  *
63  * Revision 1.10 2004/04/08 16:38:43 e_gourgoulhon
64  * Sym_tensor_tt::set_khi_mu: added argument dzp (dzpuis of resulting h^{ij}).
65  *
66  * Revision 1.9 2004/03/04 09:53:04 e_gourgoulhon
67  * Methods eta(), mu() and upate(): use of Scalar::mult_r_dzpuis and
68  * change of dzpuis behavior of eta and mu.
69  *
70  * Revision 1.8 2004/03/03 13:16:21 j_novak
71  * New potential khi (p_khi) and the functions manipulating it.
72  *
73  * Revision 1.7 2004/02/05 13:44:50 e_gourgoulhon
74  * Major modif. of methods eta(), mu() and update() to treat
75  * any value of dzpuis, thanks to the new definitions of
76  * Scalar::mult_r(), Scalar::dsdr(), etc...
77  *
78  * Revision 1.6 2004/01/28 13:25:41 j_novak
79  * The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
80  * In the div/mult _r_dzpuis, there is no more default value.
81  *
82  * Revision 1.5 2003/11/05 15:28:31 e_gourgoulhon
83  * Corrected error in update.
84  *
85  * Revision 1.4 2003/11/04 23:03:34 e_gourgoulhon
86  * First full version of method update().
87  * Add method set_rr_mu.
88  * Method set_eta_mu ---> set_rr_eta_mu.
89  *
90  * Revision 1.3 2003/11/04 09:35:27 e_gourgoulhon
91  * First operational version of update_tp().
92  *
93  * Revision 1.2 2003/11/03 22:33:36 e_gourgoulhon
94  * Added methods update_tp and set_eta_mu.
95  *
96  * Revision 1.1 2003/11/03 17:08:37 e_gourgoulhon
97  * First version
98  *
99  *
100  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
101  *
102  */
103 
104 // Headers C
105 #include <cstdlib>
106 #include <cassert>
107 
108 // Headers Lorene
109 #include "tensor.h"
110 
111  //--------------//
112  // khi //
113  //--------------//
114 
115 namespace Lorene {
116 const Scalar& Sym_tensor_tt::khi() const {
117 
118  if (p_khi == 0x0) { // a new computation is necessary
119 
120  // All this has a meaning only for spherical components:
121  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
122 
123  // khi is computed from $h^{rr}$ component
124 
125  p_khi = new Scalar(operator()(1,1)) ;
126  p_khi->mult_r() ;
127  p_khi->mult_r() ;
128  }
129 
130  return *p_khi ;
131 
132 }
133 
134 
135  //--------------//
136  // eta //
137  //--------------//
138 
139 
140 const Scalar& Sym_tensor_tt::eta(Param* par) const {
141 
142 
143  if (p_eta == 0x0) { // a new computation is necessary
144 
145 
146  // All this has a meaning only for spherical components:
147  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
148 
149  // eta is computed from the divergence-free condition:
150 
151  int dzp = operator()(1,1).get_dzpuis() ;
152  int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
153 
154  Scalar source_eta(*mp) ;
155 
156  if (p_khi == 0x0) { // eta is computed from h^rr
157  // -------------------------
158 
159  source_eta = - operator()(1,1).dsdr() ;
160 
161  // dhrr contains - dh^{rr}/dr in all domains but the CED,
162  // in the CED: - r^2 dh^{rr}/dr if dzp = 0 (1)
163  // - r^(dzp+1) dh^{rr}/dr if dzp > 0 (2)
164 
165 
166  // Multiplication by r of (-d h^{rr}/dr) (with the same dzpuis as h^{rr})
167  source_eta.mult_r_dzpuis( dzp ) ;
168 
169  // Substraction of the h^rr part and multiplication by r :
170  source_eta -= 3. * operator()(1,1) ;
171 
172  source_eta.mult_r_dzpuis(dzp_resu) ;
173  }
174  else { // eta is computed from khi
175  // ------------------------
176 
177  source_eta = - p_khi->dsdr() ;
178  int diff_dzp = source_eta.get_dzpuis() - dzp_resu ;
179  assert( diff_dzp >= 0 ) ;
180  source_eta.dec_dzpuis(diff_dzp) ;
181 
182  Scalar tmp(*p_khi) ;
183  tmp.div_r_dzpuis(dzp_resu) ;
184 
185  source_eta -= tmp ;
186 
187  }
188 
189 
190  // Resolution of the angular Poisson equation for eta
191  // --------------------------------------------------
192  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
193  p_eta = new Scalar( source_eta.poisson_angu() ) ;
194  }
195  else {
196  Scalar resu (*mp) ;
197  resu = 0. ;
198  mp->poisson_angu(source_eta, *par, resu) ;
199  p_eta = new Scalar( resu ) ;
200  }
201 
202  }
203 
204  return *p_eta ;
205 
206 }
207 
208 
209 
210  //-------------------//
211  // set_rr_eta_mu //
212  //-------------------//
213 
214 
215 void Sym_tensor_tt::set_rr_eta_mu(const Scalar& hrr, const Scalar& eta_i,
216  const Scalar& mu_i) {
217 
218  // All this has a meaning only for spherical components:
219  assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
220 
221  set(1,1) = hrr ; // h^{rr}
222  // calls del_deriv() and therefore delete previous
223  // p_eta and p_mu
224 
225  p_eta = new Scalar( eta_i ) ; // eta
226 
227  p_mu = new Scalar( mu_i ) ; // mu
228 
229  update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
230 
231 }
232 
233  //---------------//
234  // set_rr_mu //
235  //---------------//
236 
237 
238 void Sym_tensor_tt::set_rr_mu(const Scalar& hrr, const Scalar& mu_i) {
239 
240  // All this has a meaning only for spherical components:
241  assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
242 
243  set(1,1) = hrr ; // h^{rr}
244  // calls del_deriv() and therefore delete previous
245  // p_eta and p_mu
246 
247  p_mu = new Scalar( mu_i ) ; // mu
248 
249  eta() ; // computes eta form the divergence-free condition
250 
251  update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
252 
253 }
254 
255  //-------------------//
256  // set_khi_eta_mu //
257  //-------------------//
258 
259 
260 void Sym_tensor_tt::set_khi_eta_mu(const Scalar& khi_i, const Scalar& eta_i,
261  const Scalar& mu_i) {
262 
263  // All this has a meaning only for spherical components:
264  assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
265 
266  set(1,1) = khi_i ;
267  set(1,1).div_r() ;
268  set(1,1).div_r() ; // h^{rr}
269 
270  // calls del_deriv() and therefore delete previous
271  // p_khi, p_eta and p_mu
272 
273  p_khi = new Scalar( khi_i ) ; // khi
274 
275  p_eta = new Scalar( eta_i ) ; // eta
276 
277  p_mu = new Scalar( mu_i ) ; // mu
278 
279  update( khi_i.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
280 
281 }
282 
283  //---------------//
284  // set_khi_mu //
285  //---------------//
286 
287 
288 void Sym_tensor_tt::set_khi_mu(const Scalar& khi_i, const Scalar& mu_i,
289  int dzp, Param* par1, Param* par2, Param* par3) {
290 
291  // All this has a meaning only for spherical components:
292  assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
293 
294  set(1,1) = khi_i ;
295  // calls del_deriv() and therefore delete previous
296  // p_eta and p_mu
297 
298  assert( khi_i.check_dzpuis(0) ) ;
299  if (dzp == 0) {
300  set(1,1).div_r() ;
301  set(1,1).div_r() ; // h^{rr}
302  }
303  else {
304  assert(dzp == 2) ; //## temporary: the other cases are not treated yet
305  set(1,1).div_r_dzpuis(1) ;
306  set(1,1).div_r_dzpuis(2) ;
307  }
308 
309  p_khi = new Scalar ( khi_i ) ; // khi
310 
311  p_mu = new Scalar( mu_i ) ; // mu
312 
313  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
314  eta() ; // computes eta form the divergence-free condition
315  // dzp = 0 ==> eta.dzpuis = 0
316  // dzp = 2 ==> eta.dzpuis = 1
317 
318  update(dzp) ; // all h^{ij}, except for h^{rr}
319  }
320  else {
321  eta(par1) ; // computes eta form the divergence-free condition
322  // dzp = 0 ==> eta.dzpuis = 0
323  // dzp = 2 ==> eta.dzpuis = 1
324 
325  update(dzp, par2, par3) ; // all h^{ij}, except for h^{rr}
326  }
327 
328 }
329 
330 
331 
332  //-------------//
333  // update //
334  //-------------//
335 
336 
337 
338 void Sym_tensor_tt::update(int dzp, Param* par1, Param* par2) {
339 
340  // All this has a meaning only for spherical components:
341  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
342 
343  assert( (p_eta != 0x0) && (p_mu != 0x0) ) ;
344 
345  Itbl idx(2) ;
346  idx.set(0) = 1 ; // r index
347 
348  // h^{r theta} :
349  // ------------
350  idx.set(1) = 2 ; // theta index
351  *cmp[position(idx)] = p_eta->srdsdt() - p_mu->srstdsdp() ;
352 
353  if (dzp == 0) {
354  assert( cmp[position(idx)]->check_dzpuis(2) ) ;
355  cmp[position(idx)]->dec_dzpuis(2) ;
356  }
357 
358  assert( cmp[position(idx)]->check_dzpuis(dzp) ) ;
359 
360  // h^{r phi} :
361  // ------------
362  idx.set(1) = 3 ; // phi index
363  *cmp[position(idx)] = p_eta->srstdsdp() + p_mu->srdsdt() ;
364 
365  if (dzp == 0) {
366  assert( cmp[position(idx)]->check_dzpuis(2) ) ;
367  cmp[position(idx)]->dec_dzpuis(2) ;
368  }
369 
370  assert( cmp[position(idx)]->check_dzpuis(dzp) ) ;
371 
372  // h^{theta phi} and h^{phi phi}
373  // -----------------------------
374 
375  //-------------- Computation of T^theta --> taut :
376 
377  Scalar tautst = operator()(1,2).dsdr() ;
378 
379  // dhrr contains dh^{rt}/dr in all domains but the CED,
380  // in the CED: r^2 dh^{rt}/dr if dzp = 0 (1)
381  // r^(dzp+1) dh^{rt}/dr if dzp > 0 (2)
382 
383  // Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
384  tautst.mult_r_dzpuis( operator()(1,2).get_dzpuis() ) ;
385 
386 
387  // Addition of the remaining parts :
388  tautst += 3 * operator()(1,2) - operator()(1,1).dsdt() ;
389  tautst.mult_sint() ;
390 
391  Scalar tmp = operator()(1,1) ;
392  tmp.mult_cost() ; // h^{rr} cos(th)
393 
394  tautst -= tmp ; // T^th / sin(th)
395 
396  Scalar taut = tautst ;
397  taut.mult_sint() ; // T^th
398 
399 
400  //----------- Computation of T^phi --> taup :
401 
402  Scalar taupst = - operator()(1,3).dsdr() ;
403 
404  // dhrr contains - dh^{rp}/dr in all domains but the CED,
405  // in the CED: - r^2 dh^{rp}/dr if dzp = 0 (3)
406  // - r^(dzp+1) dh^{rp}/dr if dzp > 0 (4)
407 
408  // Multiplication by r of -dh^{rp}/dr (with the same dzpuis than h^{rp})
409  taupst.mult_r_dzpuis( operator()(1,3).get_dzpuis() ) ;
410 
411 
412  // Addition of the remaining part :
413 
414  taupst -= 3 * operator()(1,3) ;
415  taupst.mult_sint() ; // T^ph / sin(th)
416 
417  Scalar taup = taupst ;
418  taup.mult_sint() ; // T^ph
419 
420 
421  //------------------- Computation of F and h^[ph ph}
422 
423  tmp = tautst ;
424  tmp.mult_cost() ; // T^th / tan(th)
425 
426  // dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
427  tmp = taut.dsdt() + tmp + taup.stdsdp() ;
428 
429  Scalar tmp2 (*mp) ;
430  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
431  tmp2 = tmp.poisson_angu() ; // F
432  }
433  else {
434  tmp2 = 0. ;
435  mp->poisson_angu(tmp, *par1, tmp2) ; // F
436  }
437 
438 
439  tmp2.div_sint() ;
440  tmp2.div_sint() ; // h^{ph ph}
441 
442  idx.set(0) = 3 ; // phi index
443  idx.set(1) = 3 ; // phi index
444  *cmp[position(idx)] = tmp2 ; // h^{ph ph} is updated
445 
446 
447  //------------------- Computation of G and h^[th ph}
448 
449  tmp = taupst ;
450  tmp.mult_cost() ; // T^ph / tan(th)
451 
452  // - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
453  tmp = - taut.stdsdp() + taup.dsdt() + tmp ;
454 
455  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
456  tmp2 = tmp.poisson_angu() ; // G
457  }
458  else {
459  tmp2 = 0. ;
460  mp->poisson_angu(tmp, *par2, tmp2) ; // G
461  }
462 
463  tmp2.div_sint() ;
464  tmp2.div_sint() ; // h^{th ph}
465 
466  idx.set(0) = 2 ; // theta index
467  idx.set(1) = 3 ; // phi index
468  *cmp[position(idx)] = tmp2 ; // h^{th ph} is updated
469 
470  // h^{th th} (from the trace-free condition)
471  // ---------
472  idx.set(1) = 2 ; // theta index
473  *cmp[position(idx)] = - operator()(1,1) - operator()(3,3) ;
474 
475 
476  Sym_tensor_trans::del_deriv() ; //## in order not to delete p_eta and p_mu
477 
478 
479 
480 }
481 
482 
483  //-----------------//
484  // set_A_tilde_B //
485  //-----------------//
486 
487 void Sym_tensor_tt::set_A_tildeB(const Scalar& a_in, const Scalar& tb_in,
488  Param* par_bc, Param* par_mat) {
489 
490  Scalar zero(*mp) ;
491  zero.set_etat_zero() ;
492  set_AtB_trace(a_in, tb_in, zero, par_bc, par_mat) ;
493  return ;
494 }
495 
496 
497 
498 
499 
500 
501 }
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
Basic integer array class.
Definition: itbl.h:122
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Affine radial mapping.
Definition: map.h:2027
virtual void poisson_angu(const Scalar &source, Param &par, Scalar &uu, double lambda=0) const =0
Computes the solution of the generalized angular Poisson equation.
Parameter storage.
Definition: param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
const Scalar & srdsdt() const
Returns of *this .
Definition: scalar_deriv.C:145
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
void div_sint()
Division by .
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
void mult_cost()
Multiplication by .
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:200
const Scalar & srstdsdp() const
Returns of *this .
Definition: scalar_deriv.C:177
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:873
void div_r()
Division by r everywhere; dzpuis is not changed.
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
void mult_sint()
Multiplication by .
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void set_AtB_trace(const Scalar &a_in, const Scalar &tb_in, const Scalar &trace, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A , and the trace.
virtual void del_deriv() const
Deletes the derived quantities.
void set_khi_mu(const Scalar &khi_i, const Scalar &mu_i, int dzp=0, Param *par1=0x0, Param *par2=0x0, Param *par3=0x0)
Sets the component , as well as the angular potential (see member p_khi and p_mu ).
void set_rr_mu(const Scalar &hrr, const Scalar &mu_i)
Sets the component , as well as the angular potential (see member p_mu ).
void set_A_tildeB(const Scalar &a_in, const Scalar &tb_in, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A and .
Scalar * p_khi
Field such that the component .
Definition: sym_tensor.h:946
void set_khi_eta_mu(const Scalar &khi_i, const Scalar &eta_i, const Scalar &mu_i)
Sets the component , as well as the angular potentials and (see members p_khi , p_eta and p_mu ).
void update(int dzp, Param *par1=0x0, Param *par2=0x0)
Computes the components , , , and , from and the potentials and .
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
void set_rr_eta_mu(const Scalar &hrr, const Scalar &eta_i, const Scalar &mu_i)
Sets the component , as well as the angular potentials and (see members p_eta and p_mu ).
const Scalar & khi() const
Gives the field such that the component .
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition: sym_tensor.h:274
Scalar * p_eta
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition: sym_tensor.h:260
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition: tensor.C:798
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition: tensor_sym.C:245
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Lorene prototypes.
Definition: app_hor.h:64