LORENE
tenseur_sym_operateur.C
1/*
2 * Copyright (c) 1999-2001 Philippe Grandclement
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 * Copyright (c) 2002 Jerome Novak
5 *
6 * This file is part of LORENE.
7 *
8 * LORENE is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * LORENE is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with LORENE; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 *
22 */
23
24
25char tenseur_sym_operateur_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym_operateur.C,v 1.8 2014/10/13 08:53:43 j_novak Exp $" ;
26
27/*
28 * $Id: tenseur_sym_operateur.C,v 1.8 2014/10/13 08:53:43 j_novak Exp $
29 * $Log: tenseur_sym_operateur.C,v $
30 * Revision 1.8 2014/10/13 08:53:43 j_novak
31 * Lorene classes and functions now belong to the namespace Lorene.
32 *
33 * Revision 1.7 2014/10/06 15:13:19 j_novak
34 * Modified #include directives to use c++ syntax.
35 *
36 * Revision 1.6 2003/03/03 19:41:34 f_limousin
37 * Suppression of an assert on a metric associated with a tensor.
38 *
39 * Revision 1.5 2002/10/16 14:37:15 j_novak
40 * Reorganization of #include instructions of standard C++, in order to
41 * use experimental version 3 of gcc.
42 *
43 * Revision 1.4 2002/09/10 13:44:17 j_novak
44 * The method "manipule" of one indice has been removed for Tenseur_sym objects
45 * (the result cannot be a Tenseur_sym).
46 * The method "sans_trace" now computes the traceless part of a Tenseur (or
47 * Tenseur_sym) of valence 2.
48 *
49 * Revision 1.3 2002/09/06 14:49:26 j_novak
50 * Added method lie_derive for Tenseur and Tenseur_sym.
51 * Corrected various errors for derive_cov and arithmetic.
52 *
53 * Revision 1.2 2002/08/07 16:14:12 j_novak
54 * class Tenseur can now also handle tensor densities, this should be transparent to older codes
55 *
56 * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
57 * LORENE
58 *
59 * Revision 2.2 2000/02/09 19:32:22 eric
60 * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
61 * argument des constructeurs.
62 *
63 * Revision 2.1 2000/01/11 11:15:08 eric
64 * Gestion de la base vectorielle (triad).
65 *
66 * Revision 2.0 1999/12/02 17:19:02 phil
67 * *** empty log message ***
68 *
69 *
70 * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym_operateur.C,v 1.8 2014/10/13 08:53:43 j_novak Exp $
71 *
72 */
73
74// Headers C
75#include <cstdlib>
76#include <cassert>
77#include <cmath>
78
79// Headers Lorene
80#include "tenseur.h"
81#include "metrique.h"
82
83namespace Lorene {
85
86 assert ((t1.get_etat() != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
87 assert (t1.get_mp() == t2.mp) ;
88
89 int val_res = t1.get_valence() + t2.valence ;
90 double poids_res = t1.get_poids() + t2.poids ;
91 poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
92 const Metrique* met_res = 0x0 ;
93 if (poids_res != 0.) {
94 // assert((t1.get_metric() != 0x0) || (t2.metric != 0x0)) ;
95 if (t1.get_metric() != 0x0) met_res = t1.get_metric() ;
96 else met_res = t2.metric ;
97 }
98
99 Itbl tipe (val_res) ;
100 tipe.set_etat_qcq() ;
101 for (int i=0 ; i<t1.get_valence() ; i++)
102 tipe.set(i) = t1.get_type_indice(i) ;
103 for (int i=0 ; i<t2.valence ; i++)
104 tipe.set(i+t1.get_valence()) = t2.type_indice(i) ;
105
106
107 if ( t1.get_valence() != 0 ) {
108 assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
109 }
110
111 Tenseur_sym res(*t1.get_mp(), val_res, tipe, *(t2.get_triad()),
113
114
115 if ((t1.get_etat() == ETATZERO) || (t2.etat == ETATZERO))
116 res.set_etat_zero() ;
117 else {
118 res.set_etat_qcq() ;
119 Itbl jeux_indice_t1 (t1.get_valence()) ;
120 jeux_indice_t1.set_etat_qcq() ;
121 Itbl jeux_indice_t2 (t2.valence) ;
122 jeux_indice_t2.set_etat_qcq() ;
123
124 for (int i=0 ; i<res.n_comp ; i++) {
125 Itbl jeux_indice_res(res.donne_indices(i)) ;
126 for (int j=0 ; j<t1.get_valence() ; j++)
128 for (int j=0 ; j<t2.valence ; j++)
129 jeux_indice_t2.set(j) = jeux_indice_res(j+t1.get_valence()) ;
130
132 }
133 }
134 return res ;
135}
136
137Tenseur_sym manipule (const Tenseur_sym& t1, const Metrique& met) {
138
139 Tenseur* auxi ;
140 Tenseur* auxi_old = new Tenseur(t1) ;
141
142 for (int i=0 ; i<t1.valence ; i++) {
143 auxi = new Tenseur(manipule(*auxi_old, met, i)) ;
144 delete auxi_old ;
145 auxi_old = new Tenseur(*auxi) ;
146 delete auxi ;
147 }
148
150 delete auxi_old ;
151 return result ;
152}
153
155 const Metrique* met)
156{
157 assert ( (t.get_etat() != ETATNONDEF) && (x.get_etat() != ETATNONDEF) ) ;
158 assert(x.get_valence() == 1) ;
159 assert(x.get_type_indice(0) == CON) ;
160 assert(x.get_poids() == 0.) ;
161 assert(t.get_mp() == x.get_mp()) ;
162
163 int val = t.get_valence() ;
164 double poids = t.get_poids() ;
165 Itbl tipe (val+1) ;
166 tipe.set_etat_qcq() ;
167 tipe.set(0) = COV ;
168 Itbl tipx(2) ;
169 tipx.set_etat_qcq() ;
170 tipx.set(0) = COV ;
171 tipx.set(1) = CON ;
172 for (int i=0 ; i<val ; i++)
173 tipe.set(i+1) = t.get_type_indice(i) ;
174 Tenseur_sym dt(*t.get_mp(), val+1, tipe, *t.get_triad(), t.get_metric(),
175 poids) ;
176 Tenseur dx(*x.get_mp(), 2, tipx, x.get_triad()) ;
177 if (met == 0x0) {
178 dx = x.gradient() ;
179 dt = t.gradient() ;
180 }
181 else {
182 dx = x.derive_cov(*met) ;
183 dt = t.derive_cov(*met) ;
184 }
185 Tenseur_sym resu(contract(x,0,dt,0)) ;
186 Tenseur* auxi ;
187 if ((val!=0)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO)) {
188 assert(t.get_triad()->identify() == x.get_triad()->identify()) ;
189
190 for (int i=0 ; i<val ; i++) {
191 if (t.get_type_indice(i) == COV) {
192 auxi = new Tenseur(contract(t,i,dx,1)) ;
193
194 Itbl indices_aux(val) ;
195 indices_aux.set_etat_qcq() ;
196 for (int j=0 ; j<resu.get_n_comp() ; j++) {
197
198 Itbl indices (resu.donne_indices(j)) ;
199 indices_aux.set(val-1) = indices(i) ;
200 for (int idx=0 ; idx<val-1 ; idx++)
201 if (idx<i)
202 indices_aux.set(idx) = indices(idx) ;
203 else
204 indices_aux.set(idx) = indices(idx+1) ;
205
206 resu.set(indices) += (*auxi)(indices_aux) ;
207 }
208 }
209 else {
210 auxi = new Tenseur(contract(t,i,dx,0)) ;
211
212 Itbl indices_aux(val) ;
213 indices_aux.set_etat_qcq() ;
214
215 //On range comme il faut :
216 for (int j=0 ; j<resu.get_n_comp() ; j++) {
217
218 Itbl indices (resu.donne_indices(j)) ;
219 indices_aux.set(val-1) = indices(i) ;
220 for (int idx=0 ; idx<val-1 ; idx++)
221 if (idx<i)
222 indices_aux.set(idx) = indices(idx) ;
223 else
224 indices_aux.set(idx) = indices(idx+1) ;
225 resu.set(indices) -= (*auxi)(indices_aux) ;
226 }
227 }
228 delete auxi ;
229 }
230 }
231 if ((poids != 0.)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO))
232 resu = resu + poids*contract(dx,0,1)*t ;
233
234 return resu ;
235}
236
237Tenseur_sym sans_trace(const Tenseur_sym& t, const Metrique& metre)
238{
239 assert(t.get_etat() != ETATNONDEF) ;
240 assert(metre.get_etat() != ETATNONDEF) ;
241 assert(t.get_valence() == 2) ;
242
243 Tenseur_sym resu(t) ;
244 if (resu.get_etat() == ETATZERO) return resu ;
245 assert(resu.get_etat() == ETATQCQ) ;
246
247 int t0 = t.get_type_indice(0) ;
248 int t1 = t.get_type_indice(1) ;
249 Itbl mix(2) ;
250 mix.set_etat_qcq() ;
251 mix.set(0) = (t0 == t1 ? -t0 : t0) ;
252 mix.set(1) = t1 ;
253
254 Tenseur tmp(*t.get_mp(), 2, mix, *t.get_triad(), t.get_metric(),
255 t.get_poids()) ;
256 if (t0 == t1)
257 tmp = manipule(t, metre, 0) ;
258 else
259 tmp = t ;
260
261 Tenseur trace(contract(tmp, 0, 1)) ;
262
263 if (t0 == t1) {
264 switch (t0) {
265 case COV : {
266 resu = resu - 1./3.*trace * metre.cov() ;
267 break ;
268 }
269 case CON : {
270 resu = resu - 1./3.*trace * metre.con() ;
271 break ;
272 }
273 default :
274 cout << "Erreur bizarre dans sans_trace!" << endl ;
275 abort() ;
276 break ;
277 }
278 }
279 else {
280 Tenseur_sym delta(*t.get_mp(), 2, mix, *t.get_triad(),
281 t.get_metric(), t.get_poids()) ;
282 delta.set_etat_qcq() ;
283 for (int i=0; i<3; i++)
284 for (int j=i; j<3; j++)
285 delta.set(i,j) = (i==j ? 1 : 0) ;
286 resu = resu - trace/3. * delta ;
287 }
288 resu.set_std_base() ;
289 return resu ;
290}
291}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition tenseur.h:1253
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
const Tenseur & derive_cov(const Metrique &met) const
Returns the covariant derivative of *this , with respect to met .
Definition tenseur.C:1554
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
const Map * get_mp() const
Returns pointer on the mapping.
Definition tenseur.h:699
double get_poids() const
Returns the weight.
Definition tenseur.h:738
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
int get_valence() const
Returns the valence.
Definition tenseur.h:710
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition tenseur.h:745
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Tenseur lie_derive(const Tenseur &t, const Tenseur &x, const Metrique *=0x0)
Lie Derivative of t with respect to x .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur sans_trace(const Tenseur &tens, const Metrique &metre)
Computes the traceless part of a Tenseur of valence 2.
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
Lorene prototypes.
Definition app_hor.h:64