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 char et_rot_extr_curv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_extr_curv.C,v 1.2 2014/10/13 08:52:57 j_novak Exp $" ;
28 
29 
30 /*
31  * $Id: et_rot_extr_curv.C,v 1.2 2014/10/13 08:52:57 j_novak Exp $
32  * $Log: et_rot_extr_curv.C,v $
33  * Revision 1.2 2014/10/13 08:52:57 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
37  * LORENE
38  *
39  * Revision 2.2 2000/11/18 17:14:35 eric
40  * Traitement du cas np=1 (axisymetrie).
41  *
42  * Revision 2.1 2000/10/06 15:07:10 eric
43  * Traitement des cas ETATZERO.
44  *
45  * Revision 2.0 2000/09/18 16:15:45 eric
46  * *** empty log message ***
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_extr_curv.C,v 1.2 2014/10/13 08:52:57 j_novak Exp $
50  *
51  */
52 
53 // Headers Lorene
54 #include "etoile.h"
55 
56 namespace Lorene {
58 
59 
60  // ---------------------------------------
61  // Special treatment for axisymmetric case
62  // ---------------------------------------
63 
64  if ( (mp.get_mg())->get_np(0) == 1) {
65 
66  tkij.set_etat_zero() ; // initialisation
67 
68 
69  // Computation of K_xy
70  // -------------------
71 
72  Cmp dnpdr = nphi().dsdr() ; // d/dr (N^phi)
73  Cmp dnpdt = nphi().srdsdt() ; // 1/r d/dtheta (N^phi)
74 
75  // What follows is valid only for a mapping of class Map_radial :
76  assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
77 
78  if (dnpdr.get_etat() == ETATQCQ) {
79  dnpdr.va = (dnpdr.va).mult_st() ; // multiplication by sin(theta)
80  }
81 
82  if (dnpdt.get_etat() == ETATQCQ) {
83  dnpdt.va = (dnpdt.va).mult_ct() ; // multiplication by cos(theta)
84  }
85 
86  Cmp tmp = dnpdr + dnpdt ;
87 
88  tmp.mult_rsint() ; // multiplication by r sin(theta)
89 
90  if (tmp.get_etat() != ETATZERO) {
91  tkij.set_etat_qcq() ;
92  tkij.set(0,1) = - 0.5 * tmp / nnn() ; // component (x,y)
93  }
94 
95  // Computation of K_yz
96  // -------------------
97 
98  dnpdr = nphi().dsdr() ; // d/dr (N^phi)
99  dnpdt = nphi().srdsdt() ; // 1/r d/dtheta (N^phi)
100 
101  if (dnpdr.get_etat() == ETATQCQ) {
102  dnpdr.va = (dnpdr.va).mult_ct() ; // multiplication by cos(theta)
103  }
104 
105  if (dnpdt.get_etat() == ETATQCQ) {
106  dnpdt.va = (dnpdt.va).mult_st() ; // multiplication by sin(theta)
107  }
108 
109  tmp = dnpdr - dnpdt ;
110 
111  tmp.mult_rsint() ; // multiplication by r sin(theta)
112 
113  if (tmp.get_etat() != ETATZERO) {
114  if (tkij.get_etat() != ETATQCQ) {
115  tkij.set_etat_qcq() ;
116  }
117  tkij.set(1,2) = - 0.5 * tmp / nnn() ; // component (y,z)
118  }
119 
120  // The other components are set to zero
121  // ------------------------------------
122  if (tkij.get_etat() == ETATQCQ) {
123  tkij.set(0,0) = 0 ; // component (x,x)
124  tkij.set(0,2) = 0 ; // component (x,z)
125  tkij.set(1,1) = 0 ; // component (y,y)
126  tkij.set(2,2) = 0 ; // component (z,z)
127  }
128 
129 
130 
131  }
132  else {
133 
134  // ------------
135  // General case
136  // ------------
137 
138  // Gradient (Cartesian components) of the shift
139  // D_j N^i
140 
141  Tenseur dn = shift.gradient() ;
142 
143  // Trace of D_j N^i = divergence of N^i :
144  Tenseur divn = contract(dn, 0, 1) ;
145 
146  if (divn.get_etat() == ETATQCQ) {
147 
148  // Computation of B^{-2} K_{ij}
149  // ----------------------------
150  tkij.set_etat_qcq() ;
151  for (int i=0; i<3; i++) {
152  for (int j=i; j<3; j++) {
153  tkij.set(i, j) = dn(i, j) + dn(j, i) ;
154  }
155  tkij.set(i, i) -= double(2) /double(3) * divn() ;
156  }
157 
158  tkij = - 0.5 * tkij / nnn ;
159 
160  }
161  else{
162  assert( divn.get_etat() == ETATZERO ) ;
163  tkij.set_etat_zero() ;
164  }
165  }
166 
167  // Computation of A^2 K_{ij} K^{ij}
168  // --------------------------------
169 
170  if (tkij.get_etat() == ETATZERO) {
171  ak_car = 0 ;
172  }
173  else {
174  ak_car.set_etat_qcq() ;
175 
176  ak_car.set() = 0 ;
177 
178  for (int i=0; i<3; i++) {
179  for (int j=0; j<3; j++) {
180 
181  ak_car.set() += tkij(i, j) * tkij(i, j) ;
182 
183  }
184  }
185 
186  ak_car = b_car * ak_car ;
187  }
188 
189 }
190 
191 }
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
int get_etat() const
Returns the logical state.
Definition: cmp.h:896
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1507
Tenseur shift
Total shift vector.
Definition: etoile.h:512
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:824
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:116
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
Tenseur ak_car
Scalar .
Definition: etoile.h:1586
int get_etat() const
Returns the logical state.
Definition: tenseur.h:704
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:636
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:461
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition: etoile.h:1567
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542