LORENE
et_rot_extr_curv.C
1 /*
2  * Member function of class Etoile_rot to compute the extrinsic curvature
3  */
4 
5 /*
6  * Copyright (c) 2000-2001 Eric Gourgoulhon
7  *
8  * This file is part of LORENE.
9  *
10  * LORENE is free software; you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation; either version 2 of the License, or
13  * (at your option) any later version.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 
27 
28 
29 
30 /*
31  * $Id: et_rot_extr_curv.C,v 1.3 2016/12/05 16:17:54 j_novak Exp $
32  * $Log: et_rot_extr_curv.C,v $
33  * Revision 1.3 2016/12/05 16:17:54 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.2 2014/10/13 08:52:57 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
40  * LORENE
41  *
42  * Revision 2.2 2000/11/18 17:14:35 eric
43  * Traitement du cas np=1 (axisymetrie).
44  *
45  * Revision 2.1 2000/10/06 15:07:10 eric
46  * Traitement des cas ETATZERO.
47  *
48  * Revision 2.0 2000/09/18 16:15:45 eric
49  * *** empty log message ***
50  *
51  *
52  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_extr_curv.C,v 1.3 2016/12/05 16:17:54 j_novak Exp $
53  *
54  */
55 
56 // Headers Lorene
57 #include "etoile.h"
58 
59 namespace Lorene {
61 
62 
63  // ---------------------------------------
64  // Special treatment for axisymmetric case
65  // ---------------------------------------
66 
67  if ( (mp.get_mg())->get_np(0) == 1) {
68 
69  tkij.set_etat_zero() ; // initialisation
70 
71 
72  // Computation of K_xy
73  // -------------------
74 
75  Cmp dnpdr = nphi().dsdr() ; // d/dr (N^phi)
76  Cmp dnpdt = nphi().srdsdt() ; // 1/r d/dtheta (N^phi)
77 
78  // What follows is valid only for a mapping of class Map_radial :
79  assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
80 
81  if (dnpdr.get_etat() == ETATQCQ) {
82  dnpdr.va = (dnpdr.va).mult_st() ; // multiplication by sin(theta)
83  }
84 
85  if (dnpdt.get_etat() == ETATQCQ) {
86  dnpdt.va = (dnpdt.va).mult_ct() ; // multiplication by cos(theta)
87  }
88 
89  Cmp tmp = dnpdr + dnpdt ;
90 
91  tmp.mult_rsint() ; // multiplication by r sin(theta)
92 
93  if (tmp.get_etat() != ETATZERO) {
94  tkij.set_etat_qcq() ;
95  tkij.set(0,1) = - 0.5 * tmp / nnn() ; // component (x,y)
96  }
97 
98  // Computation of K_yz
99  // -------------------
100 
101  dnpdr = nphi().dsdr() ; // d/dr (N^phi)
102  dnpdt = nphi().srdsdt() ; // 1/r d/dtheta (N^phi)
103 
104  if (dnpdr.get_etat() == ETATQCQ) {
105  dnpdr.va = (dnpdr.va).mult_ct() ; // multiplication by cos(theta)
106  }
107 
108  if (dnpdt.get_etat() == ETATQCQ) {
109  dnpdt.va = (dnpdt.va).mult_st() ; // multiplication by sin(theta)
110  }
111 
112  tmp = dnpdr - dnpdt ;
113 
114  tmp.mult_rsint() ; // multiplication by r sin(theta)
115 
116  if (tmp.get_etat() != ETATZERO) {
117  if (tkij.get_etat() != ETATQCQ) {
118  tkij.set_etat_qcq() ;
119  }
120  tkij.set(1,2) = - 0.5 * tmp / nnn() ; // component (y,z)
121  }
122 
123  // The other components are set to zero
124  // ------------------------------------
125  if (tkij.get_etat() == ETATQCQ) {
126  tkij.set(0,0) = 0 ; // component (x,x)
127  tkij.set(0,2) = 0 ; // component (x,z)
128  tkij.set(1,1) = 0 ; // component (y,y)
129  tkij.set(2,2) = 0 ; // component (z,z)
130  }
131 
132 
133 
134  }
135  else {
136 
137  // ------------
138  // General case
139  // ------------
140 
141  // Gradient (Cartesian components) of the shift
142  // D_j N^i
143 
144  Tenseur dn = shift.gradient() ;
145 
146  // Trace of D_j N^i = divergence of N^i :
147  Tenseur divn = contract(dn, 0, 1) ;
148 
149  if (divn.get_etat() == ETATQCQ) {
150 
151  // Computation of B^{-2} K_{ij}
152  // ----------------------------
153  tkij.set_etat_qcq() ;
154  for (int i=0; i<3; i++) {
155  for (int j=i; j<3; j++) {
156  tkij.set(i, j) = dn(i, j) + dn(j, i) ;
157  }
158  tkij.set(i, i) -= double(2) /double(3) * divn() ;
159  }
160 
161  tkij = - 0.5 * tkij / nnn ;
162 
163  }
164  else{
165  assert( divn.get_etat() == ETATZERO ) ;
166  tkij.set_etat_zero() ;
167  }
168  }
169 
170  // Computation of A^2 K_{ij} K^{ij}
171  // --------------------------------
172 
173  if (tkij.get_etat() == ETATZERO) {
174  ak_car = 0 ;
175  }
176  else {
177  ak_car.set_etat_qcq() ;
178 
179  ak_car.set() = 0 ;
180 
181  for (int i=0; i<3; i++) {
182  for (int j=0; j<3; j++) {
183 
184  ak_car.set() += tkij(i, j) * tkij(i, j) ;
185 
186  }
187  }
188 
189  ak_car = b_car * ak_car ;
190  }
191 
192 }
193 
194 }
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Tenseur nnn
Total lapse function.
Definition: etoile.h:512
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1513
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1510
Tenseur shift
Total shift vector.
Definition: etoile.h:515
void extrinsic_curvature()
Computes tkij and ak_car from shift , nnn and b_car .
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:119
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
Tenseur ak_car
Scalar .
Definition: etoile.h:1589
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition: etoile.h:1570
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558