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
27char 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
56namespace 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) {
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) {
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 // ----------------------------
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 ) ;
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 {
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:446
void mult_rsint()
Multiplication by .
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
void extrinsic_curvature()
Computes tkij and ak_car from shift , nnn and b_car .
Tenseur ak_car
Scalar .
Definition etoile.h:1586
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition etoile.h:1567
Tenseur b_car
Square of the metric factor B.
Definition etoile.h:1507
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Tenseur shift
Total shift vector.
Definition etoile.h:512
Base class for pure radial mappings.
Definition map.h:1536
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64