tensor_calculus.C
/*
* Methods of class Tensor for tensor calculus
*
*
*/
/*
* Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: tensor_calculus.C,v 1.10 2014/10/13 08:53:44 j_novak Exp $
* $Log: tensor_calculus.C,v $
* Revision 1.10 2014/10/13 08:53:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.9 2014/10/06 15:13:20 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.8 2004/02/26 22:49:45 e_gourgoulhon
* Added methods compute_derive_lie and derive_lie.
*
* Revision 1.7 2004/02/18 18:50:07 e_gourgoulhon
* -- Added new methods trace(...).
* -- Tensorial product moved to file tensor_calculus_ext.C, since it is not
* a method of class Tensor.
*
* Revision 1.6 2004/02/18 15:54:23 e_gourgoulhon
* Efficiency improved in method scontract: better handling of work (it is
* now considered as a reference on the relevant component of the result).
*
* Revision 1.5 2003/12/05 16:38:50 f_limousin
* Added method operator*
*
* Revision 1.4 2003/10/28 21:25:34 e_gourgoulhon
* Method contract renamed scontract.
*
* Revision 1.3 2003/10/11 16:47:10 e_gourgoulhon
* Suppressed the call to Ibtl::set_etat_qcq() after the construction
* of the Itbl's, thanks to the new property of the Itbl class.
*
* Revision 1.2 2003/10/06 20:52:22 e_gourgoulhon
* Added methods up, down and up_down.
*
* Revision 1.1 2003/10/06 15:13:38 e_gourgoulhon
* Tensor contraction.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_calculus.C,v 1.10 2014/10/13 08:53:44 j_novak Exp $
*
*/
// Headers C++
#include "headcpp.h"
// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>
// Headers Lorene
#include "tensor.h"
#include "metric.h"
//------------------//
// Trace //
//------------------//
namespace Lorene {
Tensor Tensor::trace(int ind_1, int ind_2) const {
// Les verifications :
assert( (ind_1 >= 0) && (ind_1 < valence) ) ;
assert( (ind_2 >= 0) && (ind_2 < valence) ) ;
assert( ind_1 != ind_2 ) ;
assert( type_indice(ind_1) != type_indice(ind_2) ) ;
// On veut ind_1 < ind_2 :
if (ind_1 > ind_2) {
int auxi = ind_2 ;
ind_2 = ind_1 ;
ind_1 = auxi ;
}
// On construit le resultat :
int val_res = valence - 2 ;
Itbl tipe(val_res) ;
for (int i=0 ; i<ind_1 ; i++)
tipe.set(i) = type_indice(i) ;
for (int i=ind_1 ; i<ind_2-1 ; i++)
tipe.set(i) = type_indice(i+1) ;
for (int i = ind_2-1 ; i<val_res ; i++)
tipe.set(i) = type_indice(i+2) ;
Tensor res(*mp, val_res, tipe, triad) ;
// Boucle sur les composantes de res :
Itbl jeux_indice_source(valence) ;
for (int i=0 ; i<res.get_n_comp() ; i++) {
Itbl jeux_indice_res(res.indices(i)) ;
for (int j=0 ; j<ind_1 ; j++)
jeux_indice_source.set(j) = jeux_indice_res(j) ;
for (int j=ind_1+1 ; j<ind_2 ; j++)
jeux_indice_source.set(j) = jeux_indice_res(j-1) ;
for (int j=ind_2+1 ; j<valence ; j++)
jeux_indice_source.set(j) = jeux_indice_res(j-2) ;
Scalar& work = res.set(jeux_indice_res) ;
work.set_etat_zero() ;
for (int j=1 ; j<=3 ; j++) {
jeux_indice_source.set(ind_1) = j ;
jeux_indice_source.set(ind_2) = j ;
work += (*cmp[position(jeux_indice_source)]) ;
}
}
return res ;
}
Tensor Tensor::trace(int ind1, int ind2, const Metric& gam) const {
// Les verifications :
assert( (ind1 >= 0) && (ind1 < valence) ) ;
assert( (ind2 >= 0) && (ind2 < valence) ) ;
assert( ind1 != ind2 ) ;
if ( type_indice(ind1) != type_indice(ind2) ) {
cout << "Tensor::trace(int, int, const Metric&) : Warning : \n"
<< " the two indices for the trace have opposite types,\n"
<< " hence the metric is useless !\n" ;
return trace(ind1, ind2) ;
}
else {
if ( type_indice(ind1) == COV ) {
return contract(gam.con(), 0, 1, *this, ind1, ind2) ;
}
else{
return contract(gam.cov(), 0, 1, *this, ind1, ind2) ;
}
}
}
Scalar Tensor::trace() const {
// Les verifications :
assert( valence == 2 ) ;
assert( type_indice(0) != type_indice(1) ) ;
Scalar res(*mp) ;
res.set_etat_zero() ;
for (int i=1; i<=3; i++) {
res += operator()(i,i) ;
}
return res ;
}
Scalar Tensor::trace(const Metric& gam) const {
assert( valence == 2 ) ;
if ( type_indice(0) != type_indice(1) ) {
cout << "Tensor::trace(const Metric&) : Warning : \n"
<< " the two indices have opposite types,\n"
<< " hence the metric is useless to get the trace !\n" ;
return trace() ;
}
else {
if ( type_indice(0) == COV ) {
return contract(gam.con(), 0, 1, *this, 0, 1) ;
}
else{
return contract(gam.cov(), 0, 1, *this, 0, 1) ;
}
}
}
//----------------------//
// Index manipulation //
//----------------------//
Tensor Tensor::up(int place, const Metric& met) const {
assert (valence != 0) ; // Aucun interet pour un scalaire...
assert ((place >=0) && (place < valence)) ;
Tensor auxi = Lorene::contract(met.con(), 1, *this, place) ;
// On doit remettre les indices a la bonne place ...
Itbl tipe(valence) ;
for (int i=0 ; i<valence ; i++)
tipe.set(i) = type_indice(i) ;
tipe.set(place) = CON ;
Tensor res(*mp, valence, tipe, triad) ;
Itbl place_auxi(valence) ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl place_res(res.indices(i)) ;
place_auxi.set(0) = place_res(place) ;
for (int j=1 ; j<place+1 ; j++)
place_auxi.set(j) = place_res(j-1) ;
place_res.set(place) = place_auxi(0) ;
for (int j=place+1 ; j<valence ; j++)
place_auxi.set(j) = place_res(j);
res.set(place_res) = auxi(place_auxi) ;
}
return res ;
}
Tensor Tensor::down(int place, const Metric& met) const {
assert (valence != 0) ; // Aucun interet pour un scalaire...
assert ((place >=0) && (place < valence)) ;
Tensor auxi = Lorene::contract(met.cov(), 1, *this, place) ;
// On doit remettre les indices a la bonne place ...
Itbl tipe(valence) ;
for (int i=0 ; i<valence ; i++)
tipe.set(i) = type_indice(i) ;
tipe.set(place) = COV ;
Tensor res(*mp, valence, tipe, triad) ;
Itbl place_auxi(valence) ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl place_res(res.indices(i)) ;
place_auxi.set(0) = place_res(place) ;
for (int j=1 ; j<place+1 ; j++)
place_auxi.set(j) = place_res(j-1) ;
place_res.set(place) = place_auxi(0) ;
for (int j=place+1 ; j<valence ; j++)
place_auxi.set(j) = place_res(j);
res.set(place_res) = auxi(place_auxi) ;
}
return res ;
}
Tensor Tensor::up_down(const Metric& met) const {
Tensor* auxi ;
Tensor* auxi_old = new Tensor(*this) ;
for (int i=0 ; i<valence ; i++) {
if (type_indice(i) == COV) {
auxi = new Tensor( auxi_old->up(i, met) ) ;
}
else{
auxi = new Tensor( auxi_old->down(i, met) ) ;
}
delete auxi_old ;
auxi_old = new Tensor(*auxi) ;
delete auxi ;
}
Tensor result(*auxi_old) ;
delete auxi_old ;
return result ;
}
//-----------------------//
// Lie derivative //
//-----------------------//
// Protected method
//-----------------
void Tensor::compute_derive_lie(const Vector& vv, Tensor& resu) const {
// Protections
// -----------
if (valence > 0) {
assert(vv.get_triad() == triad) ;
assert(resu.get_triad() == triad) ;
}
// Flat metric
// -----------
const Metric_flat* fmet ;
if (valence == 0) {
fmet = &(mp->flat_met_spher()) ;
}
else {
assert( triad != 0x0 ) ;
const Base_vect_spher* bvs =
dynamic_cast<const Base_vect_spher*>(triad) ;
if (bvs != 0x0) {
fmet = &(mp->flat_met_spher()) ;
}
else {
const Base_vect_cart* bvc =
dynamic_cast<const Base_vect_cart*>(triad) ;
if (bvc != 0x0) {
fmet = &(mp->flat_met_cart()) ;
}
else {
cerr << "Tensor::compute_derive_lie : unknown triad type !\n" ;
abort() ;
}
}
}
// Determination of the dzpuis parameter of the input --> dz_in
// ---------------------------------------------------
int dz_in = 0 ;
for (int ic=0; ic<n_comp; ic++) {
int dzp = cmp[ic]->get_dzpuis() ;
assert(dzp >= 0) ;
if (dzp > dz_in) dz_in = dzp ;
}
#ifndef NDEBUG
// Check : do all components have the same dzpuis ?
for (int ic=0; ic<n_comp; ic++) {
if ( !(cmp[ic]->check_dzpuis(dz_in)) ) {
cout << "######## WARNING #######\n" ;
cout << " Tensor::compute_derive_lie: the tensor components \n"
<< " do not have all the same dzpuis ! : \n"
<< " ic, dzpuis(ic), dz_in : " << ic << " "
<< cmp[ic]->get_dzpuis() << " " << dz_in << endl ;
}
}
#endif
// Initialisation to nabla_V T
// ---------------------------
resu = contract(vv, 0, derive_cov(*fmet), valence) ;
// Addition of the terms with derivatives of V (only if valence > 0)
// -------------------------------------------
if (valence > 0) {
const Tensor& dv = vv.derive_cov(*fmet) ; // gradient of V
Itbl ind1(valence) ; // working Itbl to store the indices of resu
Itbl ind0(valence) ; // working Itbl to store the indices of this
Scalar tmp(*mp) ; // working scalar
// loop on all the components of the output tensor:
int ncomp_resu = resu.get_n_comp() ;
for (int ic=0; ic<ncomp_resu; ic++) {
// indices corresponding to the component no. ic in the output tensor
ind1 = resu.indices(ic) ;
tmp = 0 ;
// Loop on the number of indices of this
for (int id=0; id<valence; id++) {
ind0 = ind1 ;
switch( type_indice(id) ) {
case CON : {
for (int k=1; k<=3; k++) {
ind0.set(id) = k ;
tmp -= operator()(ind0) * dv(ind1(id), k) ;
}
break ;
}
case COV : {
for (int k=1; k<=3; k++) {
ind0.set(id) = k ;
tmp += operator()(ind0) * dv(k, ind1(id)) ;
}
break ;
}
default : {
cerr <<
"Tensor::compute_derive_lie: unexpected type of index !\n" ;
abort() ;
break ;
}
} // end of switch on index type
} // end of loop on the number of indices of uu
if (dz_in > 0) tmp.dec_dzpuis() ; // to get the same dzpuis as
// nabla_V T
resu.set(ind1) += tmp ; // Addition to the nabla_V T term
} // end of loop on all the components of the output tensor
} // end of case valence > 0
}
// Public interface
//-----------------
Tensor Tensor::derive_lie(const Vector& vv) const {
Tensor resu(*mp, valence, type_indice, triad) ;
compute_derive_lie(vv, resu) ;
return resu ;
}
}
Lorene::Base_vect_cart
Cartesian vectorial bases (triads).
Definition base_vect.h:201
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
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::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
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::Tensor::trace
Scalar trace() const
Trace on two different type indices for a valence 2 tensor.
Definition tensor_calculus.C:180
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::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::Tensor
Tensor(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i)
Standard constructor.
Definition tensor.C:202
Lorene::Tensor::valence
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...).
Definition tensor.h:298
Lorene::Tensor::get_n_comp
int get_n_comp() const
Returns the number of stored components.
Definition tensor.h:872
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::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
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::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::Tensor::n_comp
int n_comp
Number of stored components, depending on the symmetry.
Definition tensor.h:312
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::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::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::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::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64