LORENE
tensor_calculus.C
1/*
2 * Methods of class Tensor for tensor calculus
3 *
4 *
5 */
6
7/*
8 * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
9 *
10 * Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
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 tensor_calculus_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_calculus.C,v 1.10 2014/10/13 08:53:44 j_novak Exp $" ;
32
33/*
34 * $Id: tensor_calculus.C,v 1.10 2014/10/13 08:53:44 j_novak Exp $
35 * $Log: tensor_calculus.C,v $
36 * Revision 1.10 2014/10/13 08:53:44 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.9 2014/10/06 15:13:20 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.8 2004/02/26 22:49:45 e_gourgoulhon
43 * Added methods compute_derive_lie and derive_lie.
44 *
45 * Revision 1.7 2004/02/18 18:50:07 e_gourgoulhon
46 * -- Added new methods trace(...).
47 * -- Tensorial product moved to file tensor_calculus_ext.C, since it is not
48 * a method of class Tensor.
49 *
50 * Revision 1.6 2004/02/18 15:54:23 e_gourgoulhon
51 * Efficiency improved in method scontract: better handling of work (it is
52 * now considered as a reference on the relevant component of the result).
53 *
54 * Revision 1.5 2003/12/05 16:38:50 f_limousin
55 * Added method operator*
56 *
57 * Revision 1.4 2003/10/28 21:25:34 e_gourgoulhon
58 * Method contract renamed scontract.
59 *
60 * Revision 1.3 2003/10/11 16:47:10 e_gourgoulhon
61 * Suppressed the call to Ibtl::set_etat_qcq() after the construction
62 * of the Itbl's, thanks to the new property of the Itbl class.
63 *
64 * Revision 1.2 2003/10/06 20:52:22 e_gourgoulhon
65 * Added methods up, down and up_down.
66 *
67 * Revision 1.1 2003/10/06 15:13:38 e_gourgoulhon
68 * Tensor contraction.
69 *
70 *
71 * $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_calculus.C,v 1.10 2014/10/13 08:53:44 j_novak Exp $
72 *
73 */
74
75// Headers C++
76#include "headcpp.h"
77
78// Headers C
79#include <cstdlib>
80#include <cassert>
81#include <cmath>
82
83// Headers Lorene
84#include "tensor.h"
85#include "metric.h"
86
87
88 //------------------//
89 // Trace //
90 //------------------//
91
92
93namespace Lorene {
95
96 // Les verifications :
97 assert( (ind_1 >= 0) && (ind_1 < valence) ) ;
98 assert( (ind_2 >= 0) && (ind_2 < valence) ) ;
99 assert( ind_1 != ind_2 ) ;
101
102 // On veut ind_1 < ind_2 :
103 if (ind_1 > ind_2) {
104 int auxi = ind_2 ;
105 ind_2 = ind_1 ;
106 ind_1 = auxi ;
107 }
108
109 // On construit le resultat :
110 int val_res = valence - 2 ;
111
112 Itbl tipe(val_res) ;
113
114 for (int i=0 ; i<ind_1 ; i++)
115 tipe.set(i) = type_indice(i) ;
116 for (int i=ind_1 ; i<ind_2-1 ; i++)
117 tipe.set(i) = type_indice(i+1) ;
118 for (int i = ind_2-1 ; i<val_res ; i++)
119 tipe.set(i) = type_indice(i+2) ;
120
121 Tensor res(*mp, val_res, tipe, triad) ;
122
123 // Boucle sur les composantes de res :
124
126
127 for (int i=0 ; i<res.get_n_comp() ; i++) {
128
129 Itbl jeux_indice_res(res.indices(i)) ;
130
131 for (int j=0 ; j<ind_1 ; j++)
133 for (int j=ind_1+1 ; j<ind_2 ; j++)
135 for (int j=ind_2+1 ; j<valence ; j++)
137
138 Scalar& work = res.set(jeux_indice_res) ;
139 work.set_etat_zero() ;
140
141 for (int j=1 ; j<=3 ; j++) {
142 jeux_indice_source.set(ind_1) = j ;
143 jeux_indice_source.set(ind_2) = j ;
145 }
146
147 }
148
149 return res ;
150}
151
152
153Tensor Tensor::trace(int ind1, int ind2, const Metric& gam) const {
154
155 // Les verifications :
156 assert( (ind1 >= 0) && (ind1 < valence) ) ;
157 assert( (ind2 >= 0) && (ind2 < valence) ) ;
158 assert( ind1 != ind2 ) ;
159
160 if ( type_indice(ind1) != type_indice(ind2) ) {
161 cout << "Tensor::trace(int, int, const Metric&) : Warning : \n"
162 << " the two indices for the trace have opposite types,\n"
163 << " hence the metric is useless !\n" ;
164
165 return trace(ind1, ind2) ;
166 }
167 else {
168 if ( type_indice(ind1) == COV ) {
169 return contract(gam.con(), 0, 1, *this, ind1, ind2) ;
170 }
171 else{
172 return contract(gam.cov(), 0, 1, *this, ind1, ind2) ;
173 }
174 }
175
176}
177
178
179
181
182 // Les verifications :
183 assert( valence == 2 ) ;
184 assert( type_indice(0) != type_indice(1) ) ;
185
186 Scalar res(*mp) ;
187 res.set_etat_zero() ;
188
189 for (int i=1; i<=3; i++) {
190 res += operator()(i,i) ;
191 }
192
193 return res ;
194}
195
196
197Scalar Tensor::trace(const Metric& gam) const {
198
199 assert( valence == 2 ) ;
200
201 if ( type_indice(0) != type_indice(1) ) {
202 cout << "Tensor::trace(const Metric&) : Warning : \n"
203 << " the two indices have opposite types,\n"
204 << " hence the metric is useless to get the trace !\n" ;
205
206 return trace() ;
207 }
208 else {
209 if ( type_indice(0) == COV ) {
210 return contract(gam.con(), 0, 1, *this, 0, 1) ;
211 }
212 else{
213 return contract(gam.cov(), 0, 1, *this, 0, 1) ;
214 }
215 }
216
217}
218
219
220 //----------------------//
221 // Index manipulation //
222 //----------------------//
223
224
225Tensor Tensor::up(int place, const Metric& met) const {
226
227 assert (valence != 0) ; // Aucun interet pour un scalaire...
228 assert ((place >=0) && (place < valence)) ;
229
230
231 Tensor auxi = Lorene::contract(met.con(), 1, *this, place) ;
232
233 // On doit remettre les indices a la bonne place ...
234
235 Itbl tipe(valence) ;
236
237 for (int i=0 ; i<valence ; i++)
238 tipe.set(i) = type_indice(i) ;
239 tipe.set(place) = CON ;
240
242
244
245 for (int i=0 ; i<res.n_comp ; i++) {
246
247 Itbl place_res(res.indices(i)) ;
248
249 place_auxi.set(0) = place_res(place) ;
250 for (int j=1 ; j<place+1 ; j++)
251 place_auxi.set(j) = place_res(j-1) ;
252 place_res.set(place) = place_auxi(0) ;
253
254 for (int j=place+1 ; j<valence ; j++)
255 place_auxi.set(j) = place_res(j);
256
257 res.set(place_res) = auxi(place_auxi) ;
258 }
259
260 return res ;
261
262}
263
264
265Tensor Tensor::down(int place, const Metric& met) const {
266
267 assert (valence != 0) ; // Aucun interet pour un scalaire...
268 assert ((place >=0) && (place < valence)) ;
269
270 Tensor auxi = Lorene::contract(met.cov(), 1, *this, place) ;
271
272 // On doit remettre les indices a la bonne place ...
273
274 Itbl tipe(valence) ;
275
276 for (int i=0 ; i<valence ; i++)
277 tipe.set(i) = type_indice(i) ;
278 tipe.set(place) = COV ;
279
281
283
284 for (int i=0 ; i<res.n_comp ; i++) {
285
286 Itbl place_res(res.indices(i)) ;
287
288 place_auxi.set(0) = place_res(place) ;
289 for (int j=1 ; j<place+1 ; j++)
290 place_auxi.set(j) = place_res(j-1) ;
291 place_res.set(place) = place_auxi(0) ;
292
293 for (int j=place+1 ; j<valence ; j++)
294 place_auxi.set(j) = place_res(j);
295
296 res.set(place_res) = auxi(place_auxi) ;
297 }
298
299 return res ;
300
301}
302
303
304
306
307 Tensor* auxi ;
308 Tensor* auxi_old = new Tensor(*this) ;
309
310 for (int i=0 ; i<valence ; i++) {
311
312 if (type_indice(i) == COV) {
313 auxi = new Tensor( auxi_old->up(i, met) ) ;
314 }
315 else{
316 auxi = new Tensor( auxi_old->down(i, met) ) ;
317 }
318
319 delete auxi_old ;
320 auxi_old = new Tensor(*auxi) ;
321 delete auxi ;
322
323 }
324
326 delete auxi_old ;
327
328 return result ;
329}
330
331
332 //-----------------------//
333 // Lie derivative //
334 //-----------------------//
335
336// Protected method
337//-----------------
338
340
341
342 // Protections
343 // -----------
344 if (valence > 0) {
345 assert(vv.get_triad() == triad) ;
346 assert(resu.get_triad() == triad) ;
347 }
348
349
350 // Flat metric
351 // -----------
352
353 const Metric_flat* fmet ;
354
355 if (valence == 0) {
356 fmet = &(mp->flat_met_spher()) ;
357 }
358 else {
359 assert( triad != 0x0 ) ;
360
361 const Base_vect_spher* bvs =
362 dynamic_cast<const Base_vect_spher*>(triad) ;
363 if (bvs != 0x0) {
364 fmet = &(mp->flat_met_spher()) ;
365 }
366 else {
367 const Base_vect_cart* bvc =
368 dynamic_cast<const Base_vect_cart*>(triad) ;
369 if (bvc != 0x0) {
370 fmet = &(mp->flat_met_cart()) ;
371 }
372 else {
373 cerr << "Tensor::compute_derive_lie : unknown triad type !\n" ;
374 abort() ;
375 }
376 }
377 }
378
379 // Determination of the dzpuis parameter of the input --> dz_in
380 // ---------------------------------------------------
381 int dz_in = 0 ;
382 for (int ic=0; ic<n_comp; ic++) {
383 int dzp = cmp[ic]->get_dzpuis() ;
384 assert(dzp >= 0) ;
385 if (dzp > dz_in) dz_in = dzp ;
386 }
387
388#ifndef NDEBUG
389 // Check : do all components have the same dzpuis ?
390 for (int ic=0; ic<n_comp; ic++) {
391 if ( !(cmp[ic]->check_dzpuis(dz_in)) ) {
392 cout << "######## WARNING #######\n" ;
393 cout << " Tensor::compute_derive_lie: the tensor components \n"
394 << " do not have all the same dzpuis ! : \n"
395 << " ic, dzpuis(ic), dz_in : " << ic << " "
396 << cmp[ic]->get_dzpuis() << " " << dz_in << endl ;
397 }
398 }
399#endif
400
401
402 // Initialisation to nabla_V T
403 // ---------------------------
404
406
407
408 // Addition of the terms with derivatives of V (only if valence > 0)
409 // -------------------------------------------
410
411 if (valence > 0) {
412
413 const Tensor& dv = vv.derive_cov(*fmet) ; // gradient of V
414
415 Itbl ind1(valence) ; // working Itbl to store the indices of resu
416 Itbl ind0(valence) ; // working Itbl to store the indices of this
417 Scalar tmp(*mp) ; // working scalar
418
419 // loop on all the components of the output tensor:
420
421 int ncomp_resu = resu.get_n_comp() ;
422
423 for (int ic=0; ic<ncomp_resu; ic++) {
424
425 // indices corresponding to the component no. ic in the output tensor
426 ind1 = resu.indices(ic) ;
427
428 tmp = 0 ;
429
430 // Loop on the number of indices of this
431 for (int id=0; id<valence; id++) {
432
433 ind0 = ind1 ;
434
435 switch( type_indice(id) ) {
436
437 case CON : {
438 for (int k=1; k<=3; k++) {
439 ind0.set(id) = k ;
440 tmp -= operator()(ind0) * dv(ind1(id), k) ;
441 }
442 break ;
443 }
444
445 case COV : {
446 for (int k=1; k<=3; k++) {
447 ind0.set(id) = k ;
448 tmp += operator()(ind0) * dv(k, ind1(id)) ;
449 }
450 break ;
451 }
452
453 default : {
454 cerr <<
455 "Tensor::compute_derive_lie: unexpected type of index !\n" ;
456 abort() ;
457 break ;
458 }
459
460 } // end of switch on index type
461
462 } // end of loop on the number of indices of uu
463
464
465 if (dz_in > 0) tmp.dec_dzpuis() ; // to get the same dzpuis as
466 // nabla_V T
467
468 resu.set(ind1) += tmp ; // Addition to the nabla_V T term
469
470
471 } // end of loop on all the components of the output tensor
472
473 } // end of case valence > 0
474
475}
476
477
478// Public interface
479//-----------------
480
482
484
486
487 return resu ;
488
489}
490
491
492
493
494
495
496
497
498
499
500
501}
Cartesian vectorial bases (triads).
Definition base_vect.h:201
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Flat metric for tensor calculation.
Definition metric.h:261
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Scalar trace() const
Trace on two different type indices for a valence 2 tensor.
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor.C:525
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition tensor.h:298
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
void compute_derive_lie(const Vector &v, Tensor &resu) const
Computes the Lie derivative of this with respect to some vector field v (protected method; the public...
int n_comp
Number of stored components, depending on the symmetry.
Definition tensor.h:312
Tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Itbl type_indice
1D array of integers (class Itbl ) of size valence containing the type of each index: COV for a cov...
Definition tensor.h:310
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64