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 
31 char 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 
93 namespace Lorene {
94 Tensor Tensor::trace(int ind_1, int ind_2) const {
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 ) ;
100  assert( type_indice(ind_1) != type_indice(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 
125  Itbl jeux_indice_source(valence) ;
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++)
132  jeux_indice_source.set(j) = jeux_indice_res(j) ;
133  for (int j=ind_1+1 ; j<ind_2 ; j++)
134  jeux_indice_source.set(j) = jeux_indice_res(j-1) ;
135  for (int j=ind_2+1 ; j<valence ; j++)
136  jeux_indice_source.set(j) = jeux_indice_res(j-2) ;
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 ;
144  work += (*cmp[position(jeux_indice_source)]) ;
145  }
146 
147  }
148 
149  return res ;
150 }
151 
152 
153 Tensor 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 
197 Scalar 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 
225 Tensor 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 
241  Tensor res(*mp, valence, tipe, triad) ;
242 
243  Itbl place_auxi(valence) ;
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 
265 Tensor 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 
280  Tensor res(*mp, valence, tipe, triad) ;
281 
282  Itbl place_auxi(valence) ;
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 
305 Tensor Tensor::up_down(const Metric& met) const {
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 
325  Tensor result(*auxi_old) ;
326  delete auxi_old ;
327 
328  return result ;
329 }
330 
331 
332  //-----------------------//
333  // Lie derivative //
334  //-----------------------//
335 
336 // Protected method
337 //-----------------
338 
339 void Tensor::compute_derive_lie(const Vector& vv, Tensor& resu) const {
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 
405  resu = contract(vv, 0, derive_cov(*fmet), valence) ;
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 
481 Tensor Tensor::derive_lie(const Vector& vv) const {
482 
483  Tensor resu(*mp, valence, type_indice, triad) ;
484 
485  compute_derive_lie(vv, resu) ;
486 
487  return resu ;
488 
489 }
490 
491 
492 
493 
494 
495 
496 
497 
498 
499 
500 
501 }
Lorene::Tensor::up
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Definition: tensor_calculus.C:225
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Lorene::Metric
Metric for tensor calculation.
Definition: metric.h:90
Lorene::Tensor
Tensor handling.
Definition: tensor.h:288
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Lorene::Tensor::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition: tensor.C:798
Lorene::Base_vect_cart
Cartesian vectorial bases (triads).
Definition: base_vect.h:201
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map::flat_met_spher
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:321
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
Lorene::Tensor::up_down
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Definition: tensor_calculus.C:305
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition: metric.h:261
Lorene::Itbl
Basic integer array class.
Definition: itbl.h:122
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Scalar::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition: scalar_r_manip.C:418
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Lorene::Tensor::position
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
Lorene::Tensor::down
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
Definition: tensor_calculus.C:265
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
Lorene::Tensor::derive_lie
Tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: tensor_calculus.C:481
Lorene::Tensor::trace
Scalar trace() const
Trace on two different type indices for a valence 2 tensor.
Definition: tensor_calculus.C:180
Lorene::Map::flat_met_cart
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:331
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Lorene::Tensor::get_n_comp
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.h:872
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Tensor::Tensor
Tensor(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i)
Standard constructor.
Definition: tensor.C:202
Lorene::Tensor::type_indice
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
Lorene::Tensor::indices
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition: tensor.C:539
Lorene::Tensor::valence
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition: tensor.h:298
Lorene::Tensor::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
Lorene::Tensor::n_comp
int n_comp
Number of stored components, depending on the symmetry.
Definition: tensor.h:312
Lorene::Tensor::compute_derive_lie
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...
Definition: tensor_calculus.C:339
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279