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