tenseur_operateur.C
/*
* Copyright (c) 1999-2001 Philippe Grandclement
* Copyright (c) 2000-2001 Eric Gourgoulhon
* Copyright (c) 2002 Jerome Novak
*
* 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 tenseur_operateur_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_operateur.C,v 1.10 2014/10/13 08:53:42 j_novak Exp $" ;
/*
* $Id: tenseur_operateur.C,v 1.10 2014/10/13 08:53:42 j_novak Exp $
* $Log: tenseur_operateur.C,v $
* Revision 1.10 2014/10/13 08:53:42 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.9 2014/10/06 15:13:19 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.8 2004/05/27 07:17:19 p_grandclement
* Correction of some shadowed variables
*
* Revision 1.7 2003/06/20 14:53:38 f_limousin
* Add the function contract_desal()
*
* Revision 1.6 2003/03/03 19:38:41 f_limousin
* Suppression of an assert on a metric associated with a tensor.
*
* Revision 1.5 2002/10/16 14:37:14 j_novak
* Reorganization of #include instructions of standard C++, in order to
* use experimental version 3 of gcc.
*
* Revision 1.4 2002/09/10 13:44:17 j_novak
* The method "manipule" of one indice has been removed for Tenseur_sym objects
* (the result cannot be a Tenseur_sym).
* The method "sans_trace" now computes the traceless part of a Tenseur (or
* Tenseur_sym) of valence 2.
*
* Revision 1.3 2002/09/06 14:49:25 j_novak
* Added method lie_derive for Tenseur and Tenseur_sym.
* Corrected various errors for derive_cov and arithmetic.
*
* Revision 1.2 2002/08/07 16:14:11 j_novak
* class Tenseur can now also handle tensor densities, this should be transparent to older codes
*
* Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
* LORENE
*
* Revision 2.11 2001/08/27 10:04:21 eric
* Ajout de l'operator% (produit tensoriel avec desaliasing)
*
* Revision 2.10 2001/05/26 15:43:17 eric
* Ajout de la fonction flat_scalar_prod_desal (desaliasage)
*
* Revision 2.9 2000/10/06 15:08:40 eric
* Traitement des cas ETATZERO dans contract et flat_scal_prod
*
* Revision 2.8 2000/09/13 09:43:29 eric
* Modif skxk : appel des nouvelles fonctions Valeur::mult_cp() et
* Valeur::mult_sp() pour la multiplication par cos(phi) et sin(phi).
*
* Revision 2.7 2000/02/09 19:32:11 eric
* MODIF IMPORTANTE: la triade de decomposition est desormais passee en
* argument des constructeurs.
*
* Revision 2.6 2000/02/01 14:14:25 eric
* Ajout de la fonction amie flat_scalar_prod.
*
* Revision 2.5 2000/01/21 12:59:18 phil
* ajout de skxk
*
* Revision 2.4 2000/01/14 09:29:43 eric
* *** empty log message ***
*
* Revision 2.3 2000/01/13 17:22:37 phil
* la fonction contraction de deux tenseurs ne passe plus par produit tensoriel
*
* Revision 2.2 2000/01/11 11:14:29 eric
* Changement de nom pour la base vectorielle : base --> triad
*
* Revision 2.1 2000/01/10 17:25:15 eric
* Gestion des bases vectorielles (triades de decomposition).
*
* Revision 2.0 1999/12/02 17:19:06 phil
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_operateur.C,v 1.10 2014/10/13 08:53:42 j_novak Exp $
*
*/
// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>
// Headers Lorene
#include "tenseur.h"
#include "metrique.h"
namespace Lorene {
Tenseur operator*(const Tenseur& t1, const Tenseur& t2) {
assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
assert (t1.mp == t2.mp) ;
int val_res = t1.valence + t2.valence ;
double poids_res = t1.poids + t2.poids ;
poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
const Metrique* met_res = 0x0 ;
if (poids_res != 0.) {
// assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
if (t1.metric != 0x0) met_res = t1.metric ;
else met_res = t2.metric ;
}
// cas scalaire :
if (val_res == 0) {
Tenseur scal(*t1.mp, met_res, poids_res) ;
// cas ou un des deux est nul :
if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
scal.set_etat_zero() ;
else {
scal.set_etat_qcq() ;
scal.set() = t1() * t2() ;
}
return scal ;
}
else {
Itbl tipe (val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<t1.valence ; i++)
tipe.set(i) = t1.type_indice(i) ;
for (int i=0 ; i<t2.valence ; i++)
tipe.set(i+t1.valence) = t2.type_indice(i) ;
if ( (t1.valence != 0) && (t2.valence != 0) ) {
assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
}
const Base_vect* triad_res ;
if (t1.valence != 0) {
triad_res = t1.get_triad() ;
}
else{
triad_res = t2.get_triad() ;
}
Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
res.set_etat_zero() ;
else {
res.set_etat_qcq() ;
Itbl jeux_indice_t1 (t1.valence) ;
jeux_indice_t1.set_etat_qcq() ;
Itbl jeux_indice_t2 (t2.valence) ;
jeux_indice_t2.set_etat_qcq() ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl jeux_indice_res(res.donne_indices(i)) ;
for (int j=0 ; j<t1.valence ; j++)
jeux_indice_t1.set(j) = jeux_indice_res(j) ;
for (int j=0 ; j<t2.valence ; j++)
jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
res.set(jeux_indice_res) = t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
}
}
return res ;
}
}
//------------------------------------//
// Tensorial product wiht desaliasing //
//------------------------------------//
Tenseur operator%(const Tenseur& t1, const Tenseur& t2) {
assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
assert (t1.mp == t2.mp) ;
int val_res = t1.valence + t2.valence ;
double poids_res = t1.poids + t2.poids ;
poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
const Metrique* met_res = 0x0 ;
if (poids_res != 0.) {
// assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
if (t1.metric != 0x0) met_res = t1.metric ;
else met_res = t2.metric ;
}
// cas scalaire :
if (val_res == 0) {
Tenseur scal(*t1.mp, met_res, poids_res) ;
// cas ou un des deux est nul :
if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
scal.set_etat_zero() ;
else {
scal.set_etat_qcq() ;
scal.set() = t1() % t2() ;
}
return scal ;
}
else {
Itbl tipe (val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<t1.valence ; i++)
tipe.set(i) = t1.type_indice(i) ;
for (int i=0 ; i<t2.valence ; i++)
tipe.set(i+t1.valence) = t2.type_indice(i) ;
if ( (t1.valence != 0) && (t2.valence != 0) ) {
assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
}
const Base_vect* triad_res ;
if (t1.valence != 0) {
triad_res = t1.get_triad() ;
}
else{
triad_res = t2.get_triad() ;
}
Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
if ((t1.etat == ETATZERO) || (t2.etat == ETATZERO))
res.set_etat_zero() ;
else {
res.set_etat_qcq() ;
Itbl jeux_indice_t1 (t1.valence) ;
jeux_indice_t1.set_etat_qcq() ;
Itbl jeux_indice_t2 (t2.valence) ;
jeux_indice_t2.set_etat_qcq() ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl jeux_indice_res(res.donne_indices(i)) ;
for (int j=0 ; j<t1.valence ; j++)
jeux_indice_t1.set(j) = jeux_indice_res(j) ;
for (int j=0 ; j<t2.valence ; j++)
jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
res.set(jeux_indice_res) = t1(jeux_indice_t1) %
t2(jeux_indice_t2) ;
}
}
return res ;
}
}
Tenseur contract(const Tenseur& source, int ind_1, int ind_2) {
// Les verifications :
assert (source.etat != ETATNONDEF) ;
assert ((ind_1 >= 0) && (ind_1 < source.valence)) ;
assert ((ind_2 >= 0) && (ind_2 < source.valence)) ;
assert (source.type_indice(ind_1) != source.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 = source.valence - 2 ;
Itbl tipe (val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<ind_1 ; i++)
tipe.set(i) = source.type_indice(i) ;
for (int i=ind_1 ; i<ind_2-1 ; i++)
tipe.set(i) = source.type_indice(i+1) ;
for (int i = ind_2-1 ; i<val_res ; i++)
tipe.set(i) = source.type_indice(i+2) ;
Tenseur res(*source.mp, val_res, tipe, source.triad, source.metric,
source.poids) ;
// Cas particulier d'une source nulle
if (source.etat == ETATZERO) {
res.set_etat_zero() ;
return res ;
}
res.set_etat_qcq() ;
Cmp work(source.mp) ;
// Boucle sur les composantes de res :
Itbl jeux_indice_source(source.valence) ;
jeux_indice_source.set_etat_qcq() ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl jeux_indice_res (res.donne_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<source.valence ; j++)
jeux_indice_source.set(j) = jeux_indice_res(j-2) ;
work.set_etat_zero() ;
for (int j=0 ; j<3 ; j++) {
jeux_indice_source.set(ind_1) = j ;
jeux_indice_source.set(ind_2) = j ;
work = work + source(jeux_indice_source) ;
}
res.set(jeux_indice_res) = work ;
}
return res ;
}
Tenseur contract (const Tenseur& t1, int ind1, const Tenseur& t2, int ind2) {
assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
// Verifs :
assert ((ind1>=0) && (ind1<t1.valence)) ;
assert ((ind2>=0) && (ind2<t2.valence)) ;
assert (*(t1.mp) == *(t2.mp)) ;
// Contraction possible ?
if ( (t1.valence != 0) && (t2.valence != 0) ) {
assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
}
assert (t1.type_indice(ind1) != t2.type_indice(ind2)) ;
int val_res = t1.valence + t2.valence - 2;
double poids_res = t1.poids + t2.poids ;
poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
const Metrique* met_res = 0x0 ;
if (poids_res != 0.) {
// assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
if (t1.metric != 0x0) met_res = t1.metric ;
else met_res = t2.metric ;
}
Itbl tipe(val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<ind1 ; i++)
tipe.set(i) = t1.type_indice(i) ;
for (int i=ind1 ; i<t1.valence-1 ; i++)
tipe.set(i) = t1.type_indice(i+1) ;
for (int i=t1.valence-1 ; i<t1.valence+ind2-1 ; i++)
tipe.set(i) = t2.type_indice(i-t1.valence+1) ;
for (int i = t1.valence+ind2-1 ; i<val_res ; i++)
tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
const Base_vect* triad_res = (val_res == 0) ? 0x0 : t1.get_triad() ;
Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
// Cas particulier ou l'un des deux tenseurs est nul
if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
res.set_etat_zero() ;
return res ;
}
res.set_etat_qcq() ;
Cmp work(t1.mp) ;
// Boucle sur les composantes de res :
Itbl jeux_indice_t1(t1.valence) ;
Itbl jeux_indice_t2(t2.valence) ;
jeux_indice_t1.set_etat_qcq() ;
jeux_indice_t2.set_etat_qcq() ;
for (int comp=0 ; comp<res.n_comp ; comp++) {
Itbl jeux_indice_res (res.donne_indices(comp)) ;
for (int i=0 ; i<ind1 ; i++)
jeux_indice_t1.set(i) = jeux_indice_res(i) ;
for (int i=ind1+1 ; i<t1.valence ; i++)
jeux_indice_t1.set(i) = jeux_indice_res(i-1) ;
for (int i=0 ; i<ind2 ; i++)
jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-1) ;
for (int i=ind2+1 ; i<t2.valence ; i++)
jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
work.set_etat_zero() ;
for (int j=0 ; j<3 ; j++) {
jeux_indice_t1.set(ind1) = j ;
jeux_indice_t2.set(ind2) = j ;
work = work + t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
}
res.set(jeux_indice_res) = work ;
}
return res ;
}
Tenseur contract_desal (const Tenseur& t1, int ind1, const Tenseur& t2, int ind2) {
assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
// Verifs :
assert ((ind1>=0) && (ind1<t1.valence)) ;
assert ((ind2>=0) && (ind2<t2.valence)) ;
assert (t1.mp == t2.mp) ;
// Contraction possible ?
if ( (t1.valence != 0) && (t2.valence != 0) ) {
assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
}
assert (t1.type_indice(ind1) != t2.type_indice(ind2)) ;
int val_res = t1.valence + t2.valence - 2;
double poids_res = t1.poids + t2.poids ;
poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
const Metrique* met_res = 0x0 ;
if (poids_res != 0.) {
// assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
if (t1.metric != 0x0) met_res = t1.metric ;
else met_res = t2.metric ;
}
Itbl tipe(val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<ind1 ; i++)
tipe.set(i) = t1.type_indice(i) ;
for (int i=ind1 ; i<t1.valence-1 ; i++)
tipe.set(i) = t1.type_indice(i+1) ;
for (int i=t1.valence-1 ; i<t1.valence+ind2-1 ; i++)
tipe.set(i) = t2.type_indice(i-t1.valence+1) ;
for (int i = t1.valence+ind2-1 ; i<val_res ; i++)
tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
const Base_vect* triad_res = (val_res == 0) ? 0x0 : t1.get_triad() ;
Tenseur res(*t1.mp, val_res, tipe, triad_res, met_res, poids_res) ;
// Cas particulier ou l'un des deux tenseurs est nul
if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
res.set_etat_zero() ;
return res ;
}
res.set_etat_qcq() ;
Cmp work(t1.mp) ;
// Boucle sur les composantes de res :
Itbl jeux_indice_t1(t1.valence) ;
Itbl jeux_indice_t2(t2.valence) ;
jeux_indice_t1.set_etat_qcq() ;
jeux_indice_t2.set_etat_qcq() ;
for (int comp=0 ; comp<res.n_comp ; comp++) {
Itbl jeux_indice_res (res.donne_indices(comp)) ;
for (int i=0 ; i<ind1 ; i++)
jeux_indice_t1.set(i) = jeux_indice_res(i) ;
for (int i=ind1+1 ; i<t1.valence ; i++)
jeux_indice_t1.set(i) = jeux_indice_res(i-1) ;
for (int i=0 ; i<ind2 ; i++)
jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-1) ;
for (int i=ind2+1 ; i<t2.valence ; i++)
jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
work.set_etat_zero() ;
for (int j=0 ; j<3 ; j++) {
jeux_indice_t1.set(ind1) = j ;
jeux_indice_t2.set(ind2) = j ;
work = work + t1(jeux_indice_t1)%t2(jeux_indice_t2) ;
}
res.set(jeux_indice_res) = work ;
}
return res ;
}
Tenseur manipule(const Tenseur& t1, const Metrique& met, int place) {
assert (t1.etat != ETATNONDEF) ;
assert (met.get_etat() != ETATNONDEF) ;
int valen = t1.valence ;
assert (valen != 0) ; // Aucun interet pour un scalaire...
assert ((place >=0) && (place < valen)) ;
Itbl tipe (valen) ;
tipe.set_etat_qcq() ;
tipe.set(0) = -t1.type_indice(place) ;
for (int i=1 ; i<place+1 ; i++)
tipe.set(i) = t1.type_indice(i-1) ;
for (int i=place+1 ; i<valen ; i++)
tipe.set(i) = t1.type_indice(i) ;
Tenseur auxi(*t1.mp, valen, tipe, t1.triad) ;
if (t1.type_indice(place) == COV)
auxi = contract (met.con(), 1, t1, place) ;
else
auxi = contract (met.cov(), 1, t1, place) ;
// On doit remettre les indices a la bonne place ...
for (int i=0 ; i<valen ; i++)
tipe.set(i) = t1.type_indice(i) ;
tipe.set(place) *= -1 ;
Tenseur res(*t1.mp, valen, tipe, t1.triad, auxi.metric, auxi.poids) ;
res.set_etat_qcq() ;
Itbl place_auxi(valen) ;
place_auxi.set_etat_qcq() ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl place_res (res.donne_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<valen ; j++)
place_auxi.set(j) = place_res(j);
res.set(place_res) = auxi(place_auxi) ;
}
return res ;
}
Tenseur manipule (const Tenseur& t1, const Metrique& met) {
Tenseur* auxi ;
Tenseur* auxi_old = new Tenseur(t1) ;
for (int i=0 ; i<t1.valence ; i++) {
auxi = new Tenseur(manipule(*auxi_old, met, i)) ;
delete auxi_old ;
auxi_old = new Tenseur(*auxi) ;
delete auxi ;
}
Tenseur result(*auxi_old) ;
delete auxi_old ;
return result ;
}
Tenseur skxk(const Tenseur& source) {
// Verification
assert (source.valence > 0) ;
assert (source.etat != ETATNONDEF) ;
assert (*source.triad == source.mp->get_bvect_cart()) ;
// Le resultat :
int val_res = source.valence-1 ;
Itbl tipe (val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<val_res ; i++)
tipe.set(i) = source.type_indice(i) ;
Tenseur res (*source.mp, val_res, tipe, source.triad, source.metric,
source.poids) ;
if (source.etat == ETATZERO)
res.set_etat_zero() ;
else {
res.set_etat_qcq() ;
for (int i=0 ; i<res.n_comp ; i++) {
Itbl indices_res (res.donne_indices(i)) ;
Itbl indices_so (val_res+1) ;
indices_so.set_etat_qcq() ;
for (int j=0 ; j<val_res ; j++)
indices_so.set(j) = indices_res(j) ;
// x S_x
// -----
indices_so.set(val_res) = 0 ;
Cmp resu(source(indices_so)) ;
resu.mult_r() ; // Multipl. by r
// What follows is valid only for a mapping of class Map_radial :
assert( dynamic_cast<const Map_radial*>(source.get_mp()) != 0x0) ;
resu.va = (resu.va).mult_st() ; // Multipl. by sin(theta)
resu.va = (resu.va).mult_cp() ; // Multipl. by cos(phi)
// y S_y
// -----
indices_so.set(val_res) = 1 ;
Cmp auxiliaire (source(indices_so)) ;
auxiliaire.mult_r() ; // Multipl. by r
auxiliaire.va = (auxiliaire.va).mult_st() ; // Multipl. by sin(theta)
auxiliaire.va = (auxiliaire.va).mult_sp() ; // Multipl. by sin(phi)
resu = resu + auxiliaire ;
// z S_z
// -----
indices_so.set(val_res) = 2 ;
auxiliaire = source(indices_so) ;
auxiliaire.mult_r() ; // Multipl. by r
auxiliaire.va = (auxiliaire.va).mult_ct() ; // Multipl. by cos(theta)
resu = resu + auxiliaire ;
res.set(indices_res) = resu ;
// The End
// -------
}
}
return res ;
}
Tenseur flat_scalar_prod(const Tenseur& t1, const Tenseur& t2) {
assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
// Verifs :
assert (t1.mp == t2.mp) ;
// Contraction possible ?
if ( (t1.valence != 0) && (t2.valence != 0) ) {
assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
}
int val_res = t1.valence + t2.valence - 2;
double poids_res = t1.poids + t2.poids ;
poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
const Metrique* met_res = 0x0 ;
if (poids_res != 0.) {
assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
if (t1.metric != 0x0) met_res = t1.metric ;
else met_res = t2.metric ;
}
Itbl tipe(val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<t1.valence - 1 ; i++)
tipe.set(i) = t1.type_indice(i) ;
for (int i = t1.valence-1 ; i<val_res ; i++)
tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
Tenseur res(*t1.mp, val_res, tipe, t1.triad, met_res, poids_res) ;
// Cas particulier ou l'un des deux tenseurs est nul
if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
res.set_etat_zero() ;
return res ;
}
res.set_etat_qcq() ;
Cmp work(t1.mp) ;
// Boucle sur les composantes de res :
Itbl jeux_indice_t1(t1.valence) ;
Itbl jeux_indice_t2(t2.valence) ;
jeux_indice_t1.set_etat_qcq() ;
jeux_indice_t2.set_etat_qcq() ;
for (int ir=0 ; ir<res.n_comp ; ir++) { // Boucle sur les composantes
// du resultat
// Indices du resultat correspondant a la position ir :
Itbl jeux_indice_res = res.donne_indices(ir) ;
// Premiers indices correspondant dans t1 :
for (int i=0 ; i<t1.valence - 1 ; i++) {
jeux_indice_t1.set(i) = jeux_indice_res(i) ;
}
// Derniers indices correspondant dans t2 :
for (int i=1 ; i<t2.valence ; i++) {
jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
}
work.set_etat_zero() ;
// Sommation sur le dernier indice de t1 et le premier de t2 :
for (int j=0 ; j<3 ; j++) {
jeux_indice_t1.set(t1.valence - 1) = j ;
jeux_indice_t2.set(0) = j ;
work = work + t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
}
res.set(jeux_indice_res) = work ;
} // fin de la boucle sur les composantes du resultat
return res ;
}
Tenseur flat_scalar_prod_desal(const Tenseur& t1, const Tenseur& t2) {
assert ((t1.etat != ETATNONDEF) && (t2.etat != ETATNONDEF)) ;
// Verifs :
assert (t1.mp == t2.mp) ;
// Contraction possible ?
if ( (t1.valence != 0) && (t2.valence != 0) ) {
assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
}
int val_res = t1.valence + t2.valence - 2;
double poids_res = t1.poids + t2.poids ;
poids_res = (fabs(poids_res) < 1.e-10 ? 0. : poids_res) ;
const Metrique* met_res = 0x0 ;
if (poids_res != 0.) {
assert((t1.metric != 0x0) || (t2.metric != 0x0)) ;
if (t1.metric != 0x0) met_res = t1.metric ;
else met_res = t2.metric ;
}
Itbl tipe(val_res) ;
tipe.set_etat_qcq() ;
for (int i=0 ; i<t1.valence - 1 ; i++)
tipe.set(i) = t1.type_indice(i) ;
for (int i = t1.valence-1 ; i<val_res ; i++)
tipe.set(i) = t2.type_indice(i-t1.valence+2) ;
Tenseur res(*t1.mp, val_res, tipe, t1.triad, met_res, poids_res) ;
// Cas particulier ou l'un des deux tenseurs est nul
if ( (t1.etat == ETATZERO) || (t2.etat == ETATZERO) ) {
res.set_etat_zero() ;
return res ;
}
res.set_etat_qcq() ;
Cmp work(t1.mp) ;
// Boucle sur les composantes de res :
Itbl jeux_indice_t1(t1.valence) ;
Itbl jeux_indice_t2(t2.valence) ;
jeux_indice_t1.set_etat_qcq() ;
jeux_indice_t2.set_etat_qcq() ;
for (int ir=0 ; ir<res.n_comp ; ir++) { // Boucle sur les composantes
// du resultat
// Indices du resultat correspondant a la position ir :
Itbl jeux_indice_res = res.donne_indices(ir) ;
// Premiers indices correspondant dans t1 :
for (int i=0 ; i<t1.valence - 1 ; i++) {
jeux_indice_t1.set(i) = jeux_indice_res(i) ;
}
// Derniers indices correspondant dans t2 :
for (int i=1 ; i<t2.valence ; i++) {
jeux_indice_t2.set(i) = jeux_indice_res(t1.valence+i-2) ;
}
work.set_etat_zero() ;
// Sommation sur le dernier indice de t1 et le premier de t2 :
for (int j=0 ; j<3 ; j++) {
jeux_indice_t1.set(t1.valence - 1) = j ;
jeux_indice_t2.set(0) = j ;
work = work + t1(jeux_indice_t1) % t2(jeux_indice_t2) ;
}
res.set(jeux_indice_res) = work ;
} // fin de la boucle sur les composantes du resultat
return res ;
}
Tenseur lie_derive (const Tenseur& t, const Tenseur& x, const Metrique* met)
{
assert ( (t.get_etat() != ETATNONDEF) && (x.get_etat() != ETATNONDEF) ) ;
assert(x.get_valence() == 1) ;
assert(x.get_type_indice(0) == CON) ;
assert(x.get_poids() == 0.) ;
assert(t.get_mp() == x.get_mp()) ;
int val = t.get_valence() ;
double poids = t.get_poids() ;
Itbl tipe (val+1) ;
tipe.set_etat_qcq() ;
tipe.set(0) = COV ;
Itbl tipx(2) ;
tipx.set_etat_qcq() ;
tipx.set(0) = COV ;
tipx.set(1) = CON ;
for (int i=0 ; i<val ; i++)
tipe.set(i+1) = t.get_type_indice(i) ;
Tenseur dt(*t.get_mp(), val+1, tipe, t.get_triad(), t.get_metric(), poids) ;
Tenseur dx(*x.get_mp(), 2, tipx, x.get_triad()) ;
if (met == 0x0) {
dx = x.gradient() ;
dt = t.gradient() ;
}
else {
dx = x.derive_cov(*met) ;
dt = t.derive_cov(*met) ;
}
Tenseur resu(contract(x,0,dt,0)) ;
Tenseur* auxi ;
if ((val!=0)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO)) {
assert(t.get_triad()->identify() == x.get_triad()->identify()) ;
for (int i=0 ; i<val ; i++) {
if (t.get_type_indice(i) == COV) {
auxi = new Tenseur(contract(t,i,dx,1)) ;
Itbl indices_aux(val) ;
indices_aux.set_etat_qcq() ;
for (int j=0 ; j<resu.get_n_comp() ; j++) {
Itbl indices (resu.donne_indices(j)) ;
indices_aux.set(val-1) = indices(i) ;
for (int idx=0 ; idx<val-1 ; idx++)
if (idx<i)
indices_aux.set(idx) = indices(idx) ;
else
indices_aux.set(idx) = indices(idx+1) ;
resu.set(indices) += (*auxi)(indices_aux) ;
}
}
else {
auxi = new Tenseur(contract(t,i,dx,0)) ;
Itbl indices_aux(val) ;
indices_aux.set_etat_qcq() ;
//On range comme il faut :
for (int j=0 ; j<resu.get_n_comp() ; j++) {
Itbl indices (resu.donne_indices(j)) ;
indices_aux.set(val-1) = indices(i) ;
for (int idx=0 ; idx<val-1 ; idx++)
if (idx<i)
indices_aux.set(idx) = indices(idx) ;
else
indices_aux.set(idx) = indices(idx+1) ;
resu.set(indices) -= (*auxi)(indices_aux) ;
}
}
delete auxi ;
}
}
if ((poids != 0.)&&(t.get_etat()!=ETATZERO)&&(x.get_etat()!=ETATZERO))
resu = resu + poids*contract(dx,0,1)*t ;
return resu ;
}
Tenseur sans_trace(const Tenseur& t, const Metrique& metre)
{
assert(t.get_etat() != ETATNONDEF) ;
assert(metre.get_etat() != ETATNONDEF) ;
assert(t.get_valence() == 2) ;
Tenseur resu(t) ;
if (resu.get_etat() == ETATZERO) return resu ;
assert(resu.get_etat() == ETATQCQ) ;
int t0 = t.get_type_indice(0) ;
int t1 = t.get_type_indice(1) ;
Itbl mix(2) ;
mix.set_etat_qcq() ;
mix.set(0) = (t0 == t1 ? -t0 : t0) ;
mix.set(1) = t1 ;
Tenseur tmp(*t.get_mp(), 2, mix, *t.get_triad(), t.get_metric(),
t.get_poids()) ;
if (t0 == t1)
tmp = manipule(t, metre, 0) ;
else
tmp = t ;
Tenseur trace(contract(tmp, 0, 1)) ;
if (t0 == t1) {
switch (t0) {
case COV : {
resu = resu - 1./3.*trace * metre.cov() ;
break ;
}
case CON : {
resu = resu - 1./3.*trace * metre.con() ;
break ;
}
default :
cout << "Erreur bizarre dans sans_trace!" << endl ;
abort() ;
break ;
}
}
else {
Tenseur_sym delta(*t.get_mp(), 2, mix, *t.get_triad(),
t.get_metric(), t.get_poids()) ;
delta.set_etat_qcq() ;
for (int i=0; i<3; i++)
for (int j=i; j<3; j++)
delta.set(i,j) = (i==j ? 1 : 0) ;
resu = resu - trace/3. * delta ;
}
resu.set_std_base() ;
return resu ;
}
}
Lorene::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Lorene::Base_vect::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Base_vect the object belongs to.
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::mult_r
void mult_r()
Multiplication by r everywhere.
Definition cmp_r_manip.C:91
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition itbl.C:261
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Map_radial
Base class for pure radial mappings.
Definition map.h:1536
Lorene::Tenseur_sym
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition tenseur.h:1253
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
Lorene::Tenseur::donne_indices
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
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Lorene::Tenseur::get_type_indice
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
Lorene::Tenseur::metric
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition tenseur.h:325
Lorene::Tenseur::mp
const Map *const mp
Reference mapping.
Definition tenseur.h:306
Lorene::Tenseur::derive_cov
const Tenseur & derive_cov(const Metrique &met) const
Returns the covariant derivative of *this , with respect to met .
Definition tenseur.C:1554
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
Lorene::Tenseur::Tenseur
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition tenseur.C:209
Lorene::Tenseur::get_mp
const Map * get_mp() const
Returns pointer on the mapping.
Definition tenseur.h:699
Lorene::Tenseur::get_poids
double get_poids() const
Returns the weight.
Definition tenseur.h:738
Lorene::Tenseur::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tenseur.h:312
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates).
Definition tenseur.C:1542
Lorene::Tenseur::n_comp
int n_comp
Number of components, depending on the symmetry.
Definition tenseur.h:320
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
Lorene::Tenseur::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
Lorene::Tenseur::type_indice
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
Lorene::Tenseur::valence
int valence
Valence.
Definition tenseur.h:307
Lorene::Tenseur::etat
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition tenseur.h:321
Lorene::Tenseur::poids
double poids
For tensor densities: the weight.
Definition tenseur.h:323
Lorene::Tenseur::get_valence
int get_valence() const
Returns the valence.
Definition tenseur.h:710
Lorene::Tenseur::get_metric
const Metrique * get_metric() const
Returns a pointer on the metric defining the conformal factor for tensor densities.
Definition tenseur.h:745
Lorene::Tenseur::get_n_comp
int get_n_comp() const
Returns the number of components.
Definition tenseur.h:713
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Lorene::operator*
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Definition base_val_mult.C:102
Lorene::operator%
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition cmp_arithm.C:364
Lorene::skxk
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Definition tenseur_operateur.C:580
Lorene::flat_scalar_prod
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Definition tenseur_operateur.C:653
Lorene::sans_trace
Tenseur sans_trace(const Tenseur &tens, const Metrique &metre)
Computes the traceless part of a Tenseur of valence 2.
Definition tenseur_operateur.C:897
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene::flat_scalar_prod_desal
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Definition tenseur_operateur.C:735
Lorene::lie_derive
Tenseur lie_derive(const Tenseur &t, const Tenseur &x, const Metrique *=0x0)
Lie Derivative of t with respect to x .
Definition tenseur_operateur.C:816
Lorene::manipule
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
Definition tenseur_operateur.C:509
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:726