LORENE
tenseur_sym.C
1/*
2 * Methods of class Tenseur_sym
3 *
4 * (see file tenseur.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 1999-2001 Philippe Grandclement
10 * Copyright (c) 2000-2001 Eric Gourgoulhon
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30
31char tenseur_sym_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym.C,v 1.8 2014/10/13 08:53:42 j_novak Exp $" ;
32
33/*
34 * $Id: tenseur_sym.C,v 1.8 2014/10/13 08:53:42 j_novak Exp $
35 * $Log: tenseur_sym.C,v $
36 * Revision 1.8 2014/10/13 08:53:42 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.7 2014/10/06 15:13:19 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.6 2003/03/03 19:39:58 f_limousin
43 * Modification of an assert to have a check on a triad and not only on a pointer.
44 *
45 * Revision 1.5 2002/10/16 14:37:14 j_novak
46 * Reorganization of #include instructions of standard C++, in order to
47 * use experimental version 3 of gcc.
48 *
49 * Revision 1.4 2002/09/06 14:49:25 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.3 2002/08/14 13:46:15 j_novak
54 * Derived quantities of a Tenseur can now depend on several Metrique's
55 *
56 * Revision 1.2 2002/08/07 16:14:11 j_novak
57 * class Tenseur can now also handle tensor densities, this should be transparent to older codes
58 *
59 * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
60 * LORENE
61 *
62 * Revision 2.3 2001/10/10 13:55:23 eric
63 * Modif Joachim: pow(3, *) --> pow(3., *)
64 *
65 * Revision 2.2 2000/02/09 19:33:38 eric
66 * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
67 * argument des constructeurs.
68 *
69 * Revision 2.1 2000/01/11 11:14:41 eric
70 * Gestion de la base vectorielle (triad).
71 *
72 * Revision 2.0 1999/12/02 17:18:37 phil
73 * *** empty log message ***
74 *
75 *
76 * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_sym.C,v 1.8 2014/10/13 08:53:42 j_novak Exp $
77 *
78 */
79
80// Headers C
81#include <cstdlib>
82#include <cassert>
83#include <cmath>
84
85// Headers Lorene
86#include "tenseur.h"
87#include "metrique.h"
88
89 //--------------//
90 // Constructors //
91 //--------------//
92
93// Standard constructor
94// --------------------
95namespace Lorene {
96Tenseur_sym::Tenseur_sym(const Map& map, int val, const Itbl& tipe,
97 const Base_vect& triad_i, const Metrique* met,
98 double weight)
99 : Tenseur(map, val, tipe, int(pow(3., val-2)) * 6, triad_i,
100 met, weight) {
101
102 assert (val >= 2) ;
103}
104
105// Standard constructor when all the indices are of the same type
106// --------------------------------------------------------------
107Tenseur_sym::Tenseur_sym(const Map& map, int val, int tipe,
108 const Base_vect& triad_i, const Metrique* met,
109 double weight)
110 : Tenseur(map, val, tipe, int(pow(3., val-2)) * 6, triad_i,
111 met, weight) {
112
113 assert (val >= 2) ;
114}
115
116// Copy constructor
117// ----------------
118Tenseur_sym::Tenseur_sym (const Tenseur_sym& source) :
119 Tenseur (*source.mp, source.valence, source.type_indice,
120 int(pow(3., source.valence-2)*6), *(source.triad), source.metric,
121 source.poids) {
122
123 assert (valence >= 2) ;
124 for (int i=0 ; i<n_comp ; i++) {
125 int place_source = source.donne_place(donne_indices(i)) ;
126 if (source.c[place_source] == 0x0)
127 c[i] = 0x0 ;
128 else
129 c[i] = new Cmp (*source.c[place_source]) ;
130 }
131 etat = source.etat ;
132 assert(source.met_depend != 0x0) ;
133 assert(source.p_derive_cov != 0x0) ;
134 assert(source.p_derive_con != 0x0) ;
135 assert(source.p_carre_scal != 0x0) ;
136
137 if (source.p_gradient != 0x0)
138 p_gradient = new Tenseur_sym (*source.p_gradient) ;
139
140 for (int i=0; i<N_MET_MAX; i++) {
141 met_depend[i] = source.met_depend[i] ;
142 if (met_depend[i] != 0x0) {
143
145
146 if (source.p_derive_cov[i] != 0x0)
147 p_derive_cov[i] = new Tenseur (*source.p_derive_cov[i]) ;
148 if (source.p_derive_con[i] != 0x0)
149 p_derive_con[i] = new Tenseur (*source.p_derive_con[i]) ;
150 if (source.p_carre_scal[i] != 0x0)
151 p_carre_scal[i] = new Tenseur (*source.p_carre_scal[i]) ;
152 }
153 }
154}
155
156
157// Constructor from a Tenseur
158// --------------------------
159Tenseur_sym::Tenseur_sym (const Tenseur& source) :
160 Tenseur (*source.mp, source.valence, source.type_indice,
161 int(pow(3., source.valence-2)*6), *(source.triad), source.metric,
162 source.poids) {
163
164 assert (valence >= 2) ;
165
166 for (int i=0 ; i<n_comp ; i++) {
167 int place_source = source.donne_place(donne_indices(i)) ;
168 if (source.c[place_source] == 0x0)
169 c[i] = 0x0 ;
170 else
171 c[i] = new Cmp (*source.c[place_source]) ;
172 }
173
174 etat = source.etat ;
175
176 assert(source.met_depend != 0x0) ;
177 assert(source.p_derive_cov != 0x0) ;
178 assert(source.p_derive_con != 0x0) ;
179 assert(source.p_carre_scal != 0x0) ;
180
181 if (source.p_gradient != 0x0)
182 p_gradient = new Tenseur (*source.p_gradient) ;
183
184}
185
186
187// Constructor from a file
188// -----------------------
189Tenseur_sym::Tenseur_sym(const Map& map, const Base_vect& triad_i, FILE* fd,
190 const Metrique* met)
191 : Tenseur(map, triad_i, fd, met) {
192
193 assert (valence >= 2) ;
194 assert (n_comp == int(pow(3., valence-2))*6) ;
195}
196
197 //--------------//
198 // Destructor //
199 //--------------//
200
202
203
204
205
206
208
209 assert (idx.get_ndim() == 1) ;
210 assert (idx.get_dim(0) == valence) ;
211 for (int i=0 ; i<valence ; i++)
212 assert ((idx(i) >= 0) && (idx(i) < 3)) ;
213
214
215 // Gestion des deux derniers indices :
216 int last = idx(valence-1) ;
217 int lastm1 = idx(valence-2) ;
218 if (last < lastm1) {
219 int auxi = last ;
220 last = lastm1 ;
221 lastm1 = auxi ;
222 }
223
224 int place_fin ;
225 switch (lastm1) {
226 case 0 : {
227 place_fin = last ;
228 break ;
229 }
230 case 1 : {
231 place_fin = 2+last ;
232 break ;
233 }
234 case 2 : {
235 place_fin = 5 ;
236 break ;
237 }
238 default : {
239 abort() ;
240 }
241 }
242
243 int res = 0 ;
244 for (int i=0 ; i<valence-2 ; i++)
245 res = 3*res+idx(i) ;
246
247 res = 6*res + place_fin ;
248
249 return res ;
250}
251
253 Itbl res(valence) ;
254 res.set_etat_qcq() ;
255 assert ((place>=0) && (place<n_comp)) ;
256
257 int reste = div(place, 6).rem ;
258 place = int((place-reste)/6) ;
259
260 for (int i=valence-3 ; i>=0 ; i--) {
261 res.set(i) = div(place, 3).rem ;
262 place = int((place-res(i))/3) ;
263 }
264
265 if (reste<3) {
266 res.set(valence-2) = 0 ;
267 res.set(valence-1) = reste ;
268 }
269
270 if ((reste>2) && (reste<5)) {
271 res.set(valence-2) = 1 ;
272 res.set(valence-1) = reste - 2 ;
273 }
274
275 if (reste == 5) {
276 res.set(valence-2) = 2 ;
277 res.set(valence-1) = 2 ;
278 }
279
280 return res ;
281}
282
284
285 assert (valence == t.get_valence()) ;
286
287 triad = t.triad ;
288 poids = t.poids ;
289 metric = t.metric ;
290
291 for (int i=0 ; i<valence ; i++)
292 assert (type_indice(i) == t.type_indice(i)) ;
293
294 switch (t.get_etat()) {
295 case ETATNONDEF: {
297 break ;
298 }
299
300 case ETATZERO: {
301 set_etat_zero() ;
302 break ;
303 }
304
305 case ETATQCQ: {
306 set_etat_qcq() ;
307 for (int i=0 ; i<n_comp ; i++) {
309 if (t.c[place_t] == 0x0)
310 c[i] = 0x0 ;
311 else
312 *c[i] = *t.c[place_t] ;
313 }
314 break ;
315 }
316
317 default: {
318 cout << "Unknown state in Tenseur_sym::operator= " << endl ;
319 abort() ;
320 break ;
321 }
322 }
323}
324
326
327 assert (etat != ETATNONDEF) ;
328
329 if (p_gradient != 0x0)
330 return ;
331 else {
332
333 // Construction du resultat :
334 Itbl tipe (valence+1) ;
335 tipe.set_etat_qcq() ;
336 tipe.set(0) = COV ;
337 for (int i=0 ; i<valence ; i++)
338 tipe.set(i+1) = type_indice(i) ;
339
340 // Vectorial basis
341 // ---------------
342 assert(*triad == mp->get_bvect_cart()) ;
343
345 mp->get_bvect_cart(), metric, poids) ;
346
347
348 if (etat == ETATZERO)
350 else {
352 // Boucle sur les indices :
353
355 indices_source.set_etat_qcq() ;
356
357 for (int j=0 ; j<p_gradient->n_comp ; j++) {
358
360
361 for (int m=0 ; m<valence ; m++)
362 indices_source.set(m) = indices_res(m+1) ;
363
365 (*this)(indices_source).deriv(indices_res(0)) ;
366 }
367 }
368 }
369}
370
371
372void Tenseur_sym::fait_derive_cov (const Metrique& metre, int ind) const {
373
374 assert (etat != ETATNONDEF) ;
375 assert (valence != 0) ;
376
377 if (p_derive_cov[ind] != 0x0)
378 return ;
379 else {
381
382 if ((valence != 0) && (etat != ETATZERO)) {
383
384 assert( *(metre.gamma().get_triad()) == *triad ) ;
385
386 Tenseur* auxi ;
387 for (int i=0 ; i<valence ; i++) {
388
389 if (type_indice(i) == COV) {
390 auxi = new Tenseur(contract(metre.gamma(), 0,(*this), i)) ;
391
393 indices_gamma.set_etat_qcq() ;
394 //On range comme il faut :
395 for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
396
397 Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
398 indices_gamma.set(0) = indices(0) ;
399 indices_gamma.set(1) = indices(i+1) ;
400 for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
401 if (idx<=i+1)
402 indices_gamma.set(idx) = indices(idx-1) ;
403 else
404 indices_gamma.set(idx) = indices(idx) ;
405
406 p_derive_cov[ind]->set(indices) -= (*auxi)(indices_gamma) ;
407 }
408 }
409 else {
410 auxi = new Tenseur(contract(metre.gamma(), 1, (*this), i)) ;
411
413 indices_gamma.set_etat_qcq() ;
414
415 //On range comme il faut :
416 for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
417
418 Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
419 indices_gamma.set(0) = indices(i+1) ;
420 indices_gamma.set(1) = indices(0) ;
421 for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
422 if (idx<=i+1)
423 indices_gamma.set(idx) = indices(idx-1) ;
424 else
425 indices_gamma.set(idx) = indices(idx) ;
426 p_derive_cov[ind]->set(indices) += (*auxi)(indices_gamma) ;
427 }
428 }
429 delete auxi ;
430 }
431 }
432 if ((poids != 0.)&&(etat != ETATZERO))
434 poids*contract(metre.gamma(), 0, 2) * (*this) ;
435
436 }
437}
438
439
440
441void Tenseur_sym::fait_derive_con (const Metrique& metre, int ind) const {
442
443 if (p_derive_con[ind] != 0x0)
444 return ;
445 else {
446 // On calcul la derivee covariante :
447 if (valence != 0)
449 (contract(metre.con(), 1, derive_cov(metre), 0)) ;
450
451 else
453 (contract(metre.con(), 1, gradient(), 0)) ;
454 }
455}
456}
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Base class for coordinate mappings.
Definition map.h:670
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition tenseur.h:1253
virtual void fait_gradient() const
Calculates, if needed, the gradient of *this .
virtual void fait_derive_cov(const Metrique &met, int i) const
Calculates, if needed, the covariant derivative of *this , with respect to met .
virtual void operator=(const Tenseur &)
Assignment from a Tenseur .
virtual int donne_place(const Itbl &idx) const
Returns the position in the Cmp 1-D array c of a component given by its indices.
virtual void fait_derive_con(const Metrique &, int i) const
Calculates, if needed, the contravariant derivative of *this , with respect to met .
virtual ~Tenseur_sym()
Destructor.
virtual Itbl donne_indices(int place) const
Returns the indices of a component given by its position in the Cmp 1-D array c .
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
virtual Itbl donne_indices(int place) const
Returns the indices of a component given by its position in the Cmp 1-D array c .
Definition tenseur.C:704
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Tenseur ** p_derive_con
Array of pointers on the contravariant derivatives of *this with respect to the corresponding metri...
Definition tenseur.h:365
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition tenseur.h:325
Cmp ** c
The components.
Definition tenseur.h:322
const Map *const mp
Reference mapping.
Definition tenseur.h:306
const Metrique ** met_depend
Array of pointers on the Metrique 's used to calculate derivatives members.
Definition tenseur.h:337
virtual int donne_place(const Itbl &idx) const
Returns the position in the Cmp 1-D array c of a component given by its indices.
Definition tenseur.C:690
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
Tenseur ** p_carre_scal
Array of pointers on the scalar squares of *this with respect to the corresponding metric in *met_d...
Definition tenseur.h:372
friend Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tenseur.h:312
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
int n_comp
Number of components, depending on the symmetry.
Definition tenseur.h:320
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
Itbl type_indice
Array of size valence contening the type of each index, COV for a covariant one and CON for a contrav...
Definition tenseur.h:318
int valence
Valence.
Definition tenseur.h:307
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition tenseur.h:321
double poids
For tensor densities: the weight.
Definition tenseur.h:323
Tenseur ** p_derive_cov
Array of pointers on the covariant derivatives of *this with respect to the corresponding metric in...
Definition tenseur.h:358
int get_valence() const
Returns the valence.
Definition tenseur.h:710
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Definition tenseur.C:650
Tenseur * p_gradient
Pointer on the gradient of *this .
Definition tenseur.h:343
void set_dependance(const Metrique &met) const
To be used to describe the fact that the derivatives members have been calculated with met .
Definition tenseur.C:608
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene prototypes.
Definition app_hor.h:64