/*
 *  Methods of class Tenseur
 *
 *   (see file tenseur.h for documentation)
 *
 */

/*
 *   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_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $" ;

/*
 * $Id: tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $
 * $Log: tenseur.C,v $
 * Revision 1.15  2014/10/13 08:53:41  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.14  2014/10/06 15:13:18  j_novak
 * Modified #include directives to use c++ syntax.
 *
 * Revision 1.13  2003/03/03 19:37:31  f_limousin
 * Suppression of many assert(verif()).
 *
 * Revision 1.12  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.11  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.10  2002/08/30 13:21:38  j_novak
 * Corrected error in constructor
 *
 * Revision 1.9  2002/08/14 13:46:15  j_novak
 * Derived quantities of a Tenseur can now depend on several Metrique's
 *
 * Revision 1.8  2002/08/13 08:02:45  j_novak
 * Handling of spherical vector/tensor components added in the classes
 * Mg3d and Tenseur. Minor corrections for the class Metconf.
 *
 * Revision 1.7  2002/08/08 15:10:45  j_novak
 * The flag "plat" has been added to the class Metrique to show flat metrics.
 *
 * Revision 1.6  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.5  2002/08/02 15:07:41  j_novak
 * Member function determinant has been added to the class Metrique.
 * A better handling of spectral bases is now implemented for the class Tenseur.
 *
 * Revision 1.4  2002/05/07 07:36:03  e_gourgoulhon
 * Compatibilty with xlC compiler on IBM SP2:
 *    suppressed the parentheses around argument of instruction new:
 * 	e.g.   t = new (Tbl *[nzone])  -->   t = new Tbl*[nzone]
 *
 * Revision 1.3  2002/05/02 15:16:22  j_novak
 * Added functions for more general bi-fluid EOS
 *
 * Revision 1.2  2001/12/04 21:27:54  e_gourgoulhon
 *
 * All writing/reading to a binary file are now performed according to
 * the big endian convention, whatever the system is big endian or
 * small endian, thanks to the functions fwrite_be and fread_be
 *
 * Revision 1.1.1.1  2001/11/20 15:19:30  e_gourgoulhon
 * LORENE
 *
 * Revision 2.21  2001/10/10  13:54:40  eric
 * Modif Joachim: pow(3, *) --> pow(3., *)
 *
 * Revision 2.20  2000/12/20  09:50:08  eric
 * Correction erreur dans operator<< : la sortie doit etre flux et non cout !
 *
 * Revision 2.19  2000/10/12  13:11:23  eric
 * Methode set_std_base(): traitement du cas etat = ETATZERO (return).
 *
 * Revision 2.18  2000/09/13  12:11:40  eric
 * Ajout de la fonction allocate_all().
 *
 * Revision 2.17  2000/05/22  14:40:09  phil
 * ajout de inc_dzpuis et dec_dzpuis
 *
 * Revision 2.16  2000/03/22  09:18:57  eric
 * Traitement du cas ETATZERO dans dec2_dzpuis, inc2_dzpuis et mult_r_zec.
 *
 * Revision 2.15  2000/02/12  11:35:58  eric
 * Modif de la fonction set_std_base : appel de Valeur::set_base plutot
 * que l'affectation directe du membre Valeur::base.
 *
 * Revision 2.14  2000/02/10  18:30:47  eric
 * La fonction set_triad ne fait plus que l'affectation du membre triad.
 *
 * Revision 2.13  2000/02/10  16:11:07  eric
 * Ajout de la fonction change_triad.
 *
 * Revision 2.12  2000/02/09  19:32:39  eric
 * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
 * argument des constructeurs.
 *
 * Revision 2.11  2000/01/24  13:02:39  eric
 * Traitement du cas triad = 0x0 dans la sauvegarde/lecture fichier
 * Constructeur par lecture de fichier: met_depend est desormais initialise
 *  a 0x0.
 *
 * Revision 2.10  2000/01/20  16:02:57  eric
 * Ajout des operator=(double ) et operator=(int ).
 *
 * Revision 2.9  2000/01/20  15:34:39  phil
 * changement traid dans fait_gradient()
 *
 * Revision 2.8  2000/01/14  14:07:26  eric
 * Ajout de la fonction annule.
 *
 * Revision 2.7  2000/01/13  14:10:53  eric
 * Ajout du constructeur par copie d'un Cmp (pour un scalaire)
 * ainsi que l'affectation a un Cmp.
 *
 * Revision 2.6  2000/01/13  13:46:38  eric
 * Ajout du membre p_gradient_spher et des fonctions fait_gradient_spher(),
 *  gradient_spher() pour le calcul du gradient d'un scalaire en
 *  coordonnees spheriques sur la triade spherique associee.
 *
 * Revision 2.5  2000/01/12  13:19:04  eric
 * Les operator::(...) renvoient desormais une reference const sur le c[...]
 * correspondant et non plus un Cmp copie de c[...].
 * (ceci grace a la nouvelle fonction Map::cmp_zero()).
 *
 * Revision 2.4  2000/01/11  11:13:59  eric
 * Changement de nom pour la base vectorielle : base --> triad
 *
 * Revision 2.3  2000/01/10  17:23:07  eric
 * Modif affichage.
 * Methode fait_derive_cov : ajout de
 *   assert( metre.gamma().get_base() == base )
 * Methode set_std_base : ajout de
 *    assert( *base == mp->get_bvect_cart() )
 *
 * Revision 2.2  2000/01/10  15:15:26  eric
 * Ajout du membre base (base vectorielle sur laquelle sont definies
 *   les composantes).
 *
 * Revision 2.1  1999/12/09  12:39:23  phil
 * changement prototypage des derivees
 *
 * Revision 2.0  1999/12/02  17:18:31  phil
 * *** empty log message ***
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur.C,v 1.15 2014/10/13 08:53:41 j_novak Exp $
 *
 */

// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>

// Headers Lorene
#include "tenseur.h"
#include "metrique.h"
#include "utilitaires.h"

			//--------------//
			// Constructors //
			//--------------//
// Consistency check for tensor densities
//---------------------------------------
namespace Lorene {
bool Tenseur::verif() const {
  return ( (poids == 0.) || (metric != 0x0) ) ;
}

void Tenseur::new_der_met() {
    met_depend = new const Metrique*[N_MET_MAX] ;
    p_derive_cov = new Tenseur*[N_MET_MAX] ;
    p_derive_con = new Tenseur*[N_MET_MAX] ;
    p_carre_scal = new Tenseur*[N_MET_MAX] ;
    for (int i=0; i<N_MET_MAX; i++) {
      met_depend[i] = 0x0 ;
    }
    set_der_0x0() ;
}

// Constructor for a scalar field
// ------------------------------
Tenseur::Tenseur (const Map& map, const Metrique* met, double weight) : 
		mp(&map), valence(0), triad(0x0),
		type_indice(0), n_comp(1), etat(ETATNONDEF), poids(weight),
		metric(met) {
    
  //    assert(verif()) ;
    c = new Cmp*[n_comp] ;
    c[0] = 0x0 ;
    new_der_met() ;
}



// Constructor for a scalar field and from a {\tt Cmp} 
// ---------------------------------------------------
Tenseur::Tenseur (const Cmp& ci, const Metrique* met, double weight) : 
		mp(ci.get_mp()), valence(0), triad(0x0),
		type_indice(0), n_comp(1), etat(ci.get_etat()), poids(weight),
		metric(met){
    
    assert(ci.get_etat() != ETATNONDEF) ; 
    assert(verif()) ;
    
    c = new Cmp*[n_comp] ;

    if ( ci.get_etat() != ETATZERO ) {
      assert( ci.get_etat() == ETATQCQ ) ; 
      c[0] = new Cmp(ci) ;
    }
    else {
      c[0] = 0x0 ;
    }
    new_der_met() ;
}

// Standard constructor 
// --------------------
Tenseur::Tenseur(const Map& map, int val, const Itbl& tipe, 
		 const Base_vect& triad_i, const Metrique* met, double weight) 
		: mp(&map), valence(val), triad(&triad_i), type_indice(tipe), 
		   n_comp(int(pow(3., val))), etat(ETATNONDEF), poids(weight),
		metric(met){
		
    // Des verifs :
    assert (valence >= 0) ;
    assert (tipe.get_ndim() == 1) ;
    assert (valence == tipe.get_dim(0)) ;
    for (int i=0 ; i<valence ; i++)
	assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
    assert(verif()) ;
    
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++)
	c[i] = 0x0 ;
    new_der_met() ;
}

// Standard constructor with the triad passed as a pointer
// -------------------------------------------------------
Tenseur::Tenseur(const Map& map, int val, const Itbl& tipe, 
		 const Base_vect* triad_i, const Metrique* met, double weight) 
		: mp(&map), valence(val), triad(triad_i), type_indice(tipe), 
		   n_comp(int(pow(3., val))), etat(ETATNONDEF), poids(weight),
		metric(met){
		
    // Des verifs :
    assert (valence >= 0) ;
    assert (tipe.get_ndim() == 1) ;
    assert (valence == tipe.get_dim(0)) ;
    for (int i=0 ; i<valence ; i++)
	assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
    //   assert(verif()) ;
    
    if (valence == 0) {	    // particular case of a scalar 
	triad = 0x0 ; 
    }   
    
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++)
	c[i] = 0x0 ;
    new_der_met() ;
}




// Standard constructor when all the indices are of the same type
// --------------------------------------------------------------
Tenseur::Tenseur(const Map& map, int val, int tipe, const Base_vect& triad_i,
		 const Metrique* met, double weight) 
		: mp(&map), valence(val), triad(&triad_i), type_indice(val), 
                  n_comp(int(pow(3., val))), etat (ETATNONDEF), poids(weight), 
		  metric(met){
    
    // Des verifs :
    assert (valence >= 0) ;
    assert ((tipe == COV) || (tipe == CON)) ;
    assert(verif()) ;
    type_indice.set_etat_qcq() ;
    type_indice = tipe ;
    
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++)
	c[i] = 0x0 ;
    new_der_met() ;
}	
	
// Copy constructor
// ----------------
Tenseur::Tenseur (const Tenseur& source) : 
    mp(source.mp), valence(source.valence), triad(source.triad), 
    type_indice(source.type_indice), etat (source.etat), poids(source.poids),
    metric(source.metric) {
  
  //   assert(verif()) ;
    
    n_comp = int(pow(3., valence)) ;
        
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++) {
	int place_source = source.donne_place(donne_indices(i)) ;
	if (source.c[place_source] == 0x0)
	    c[i] = 0x0 ;
	else
	    c[i] = new Cmp(*source.c[place_source]) ;
    }

    assert(source.met_depend != 0x0) ;
    assert(source.p_derive_cov != 0x0) ;
    assert(source.p_derive_con != 0x0) ;
    assert(source.p_carre_scal != 0x0) ;
    new_der_met() ;
    
    if (source.p_gradient != 0x0)
	    p_gradient = new Tenseur (*source.p_gradient) ;
    
    if (source.p_gradient_spher != 0x0)
	    p_gradient_spher = new Tenseur (*source.p_gradient_spher) ;

    for (int i=0; i<N_MET_MAX; i++) {
      met_depend[i] = source.met_depend[i] ;
      if (met_depend[i] != 0x0) {
	
	set_dependance (*met_depend[i]) ;
	
	if (source.p_derive_cov[i] != 0x0)
	  p_derive_cov[i] = new Tenseur (*source.p_derive_cov[i]) ;
	if (source.p_derive_con[i] != 0x0)
	  p_derive_con[i] = new Tenseur (*source.p_derive_con[i]) ;
	if (source.p_carre_scal[i] != 0x0)
	    p_carre_scal[i] = new Tenseur (*source.p_carre_scal[i]) ;
      }
    }
}   

// Constructor from a symmetric tensor
// -----------------------------------
Tenseur::Tenseur (const Tenseur_sym& source) :
	mp(source.mp), valence(source.valence), triad(source.triad), 
	type_indice(source.type_indice), etat(source.etat), 
	poids(source.poids), metric(source.metric) {
    
    assert(verif()) ;
    n_comp = int(pow(3., valence)) ;
        
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++) {
	int place_source = source.donne_place(donne_indices(i)) ;
	if (source.c[place_source] == 0x0)
	    c[i] = 0x0 ;
	else
	    c[i] = new Cmp(*source.c[place_source]) ;
    }

    assert(source.met_depend != 0x0) ;
    assert(source.p_derive_cov != 0x0) ;
    assert(source.p_derive_con != 0x0) ;
    assert(source.p_carre_scal != 0x0) ;
    new_der_met() ;
    
    if (source.p_gradient != 0x0)
	p_gradient = new Tenseur_sym (*source.p_gradient) ;

    for (int i=0; i<N_MET_MAX; i++) {
      met_depend[i] = source.met_depend[i] ;
      if (met_depend[i] != 0x0) {
	
	set_dependance (*met_depend[i]) ;
	
	if (source.p_derive_cov[i] != 0x0)
	  p_derive_cov[i] = new Tenseur (*source.p_derive_cov[i]) ;
	if (source.p_derive_con[i] != 0x0)
	  p_derive_con[i] = new Tenseur (*source.p_derive_con[i]) ;
	if (source.p_carre_scal[i] != 0x0)
	    p_carre_scal[i] = new Tenseur (*source.p_carre_scal[i]) ;
      }
    }
    
}

// Constructor from a file
// -----------------------
Tenseur::Tenseur(const Map& mapping, const Base_vect& triad_i, FILE* fd, 
		 const Metrique* met)
		 : mp(&mapping), triad(&triad_i), type_indice(fd), 
                   metric(met) {
   
    fread_be(&valence, sizeof(int), 1, fd) ;

    if (valence != 0) {
	Base_vect* triad_fich = Base_vect::bvect_from_file(fd) ; 
	assert( *triad_fich == *triad) ; 
	delete triad_fich ; 
    }
    else{
	triad = 0x0 ; 
    }
    
    fread_be(&n_comp, sizeof(int), 1, fd) ;
    fread_be(&etat, sizeof(int), 1, fd) ;
    
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++)
	c[i] = 0x0 ;
    if (etat == ETATQCQ)
	for (int i=0 ; i<n_comp ; i++)
	    c[i] = new Cmp(*mp, *mp->get_mg(), fd) ;

    new_der_met() ;
    
    if (met == 0x0) 
      poids = 0. ;
    else
      fread_be(&poids, sizeof(double), 1, fd) ;
}


// Constructor from a file for a scalar field
// ------------------------------------------
Tenseur::Tenseur (const Map& mapping, FILE* fd, const Metrique* met) 
		 : mp(&mapping), type_indice(fd), metric(met){
   
    fread_be(&valence, sizeof(int), 1, fd) ;

    assert(valence == 0) ; 
    
    triad = 0x0 ; 
    
    fread_be(&n_comp, sizeof(int), 1, fd) ;
    
    assert(n_comp == 1) ; 

    fread_be(&etat, sizeof(int), 1, fd) ;
    
    c = new Cmp*[n_comp] ;

    if (etat == ETATQCQ) {
	c[0] = new Cmp(*mp, *mp->get_mg(), fd) ;
    }
    else{
	c[0] = 0x0 ; 
    }
    
    new_der_met() ;

    if (met == 0x0) 
      poids = 0. ;
    else
      fread_be(&poids, sizeof(double), 1, fd) ;
}




// Constructor used by the derived classes
// ---------------------------------------
Tenseur::Tenseur (const Map& map, int val, const Itbl& tipe, int compo, 
		const Base_vect& triad_i, const Metrique* met, double weight) :
     mp(&map), valence(val), triad(&triad_i), type_indice(tipe), n_comp(compo), 
	    etat (ETATNONDEF), poids(weight), metric(met) {
     
    // Des verifs :
    assert (valence >= 0) ;
    assert (tipe.get_ndim() == 1) ;
    assert (n_comp > 0) ;   
    assert (valence == tipe.get_dim(0)) ;
    for (int i=0 ; i<valence ; i++)
	assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
    //assert(verif()) ;
    
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++)
	c[i] = 0x0 ;
    
    new_der_met() ;
}

// Constructor used by the derived classes when all the indices are of 
// the same type.
// -------------------------------------------------------------------
Tenseur::Tenseur (const Map& map, int val, int tipe, int compo, 
		const Base_vect& triad_i, const Metrique* met, double weight) :
     mp(&map), valence(val), triad(&triad_i), type_indice(val), n_comp(compo), 
	    etat (ETATNONDEF), poids(weight), metric(met) {
    // Des verifs :
    assert (valence >= 0) ;
    assert (n_comp >= 0) ;
    assert ((tipe == COV) || (tipe == CON)) ;
    //assert(verif()) ;
    type_indice.set_etat_qcq() ;
    type_indice = tipe ;
    
    c = new Cmp*[n_comp] ;
    for (int i=0 ; i<n_comp ; i++)
	c[i] = 0x0 ;

    new_der_met() ;
}

			//--------------//
			//  Destructor  //
			//--------------//


Tenseur::~Tenseur () {
    
    del_t() ;
    delete [] met_depend ;
    delete [] p_derive_cov ;
    delete [] p_derive_con ;
    delete [] p_carre_scal ;
    delete [] c ;
}



void Tenseur::del_t() {
    del_derive() ;
    for (int i=0 ; i<n_comp ; i++)
	if (c[i] != 0x0) {
	    delete c[i] ;
	    c[i] = 0x0 ;
	    }
}

void Tenseur::del_derive_met(int j) const {

  assert( (j>=0) && (j<N_MET_MAX) ) ;
  // On gere la metrique ...
  if (met_depend[j] != 0x0) {
    for (int i=0 ; i<N_DEPEND ; i++)
      if (met_depend[j]->dependances[i] == this)
	met_depend[j]->dependances[i] = 0x0 ;
    if (p_derive_cov[j] != 0x0)
      delete p_derive_cov[j] ;
    if (p_derive_con[j] != 0x0)
      delete p_derive_con[j] ;
    if (p_carre_scal[j] != 0x0)
      delete p_carre_scal[j] ;
    set_der_met_0x0(j) ;
  }
}


void Tenseur::del_derive () const {
  for (int i=0; i<N_MET_MAX; i++) 
    del_derive_met(i) ;
  if (p_gradient != 0x0)
    delete p_gradient ;
  if (p_gradient_spher != 0x0)
    delete p_gradient_spher ;
  set_der_0x0() ;
}

void Tenseur::set_der_met_0x0(int i) const {
  met_depend[i] = 0x0 ;
    p_derive_cov[i] = 0x0 ;
    p_derive_con[i] = 0x0 ;
    p_carre_scal[i] = 0x0 ;
}


void Tenseur::set_der_0x0() const {
  for (int i=0; i<N_MET_MAX; i++) 
    set_der_met_0x0(i) ;
  p_gradient = 0x0 ;   
  p_gradient_spher = 0x0 ;   
}

int Tenseur::get_place_met(const Metrique& metre) const {
  int resu = -1 ;
  for (int i=0; i<N_MET_MAX; i++) 
    if (met_depend[i] == &metre) {
      assert(resu == -1) ;
      resu = i ;
    }
  return resu ;
}

void Tenseur::set_dependance (const Metrique& met) const {
    
  int nmet = 0 ;
  bool deja = false ;
  for (int i=0; i<N_MET_MAX; i++) {
    if (met_depend[i] == &met) deja = true ;
    if ((!deja) && (met_depend[i] != 0x0)) nmet++ ;
  }
  if (nmet == N_MET_MAX) {
    cout << "Too many metrics in Tenseur::set_dependances" << endl ;
    abort() ;
  }
  if (!deja) { 
    int conte = 0 ;
    while ((conte < N_DEPEND) && (met.dependances[conte] != 0x0))
      conte ++ ;
    
    if (conte == N_DEPEND) {
      cout << "Too many dependancies in Tenseur::set_dependances " << endl ;
      abort() ;
    }
    else {
      met.dependances[conte] = this ;
      met_depend[nmet] = &met ;
    }
  }
}

void Tenseur::set_etat_qcq() { 
    
    del_derive() ;
    for (int i=0 ; i<n_comp ; i++)
	if (c[i] == 0x0)
	    c[i] = new Cmp(mp) ;
    etat = ETATQCQ ;
}

void Tenseur::set_etat_zero() { 
    del_t() ;
    etat = ETATZERO ;
}

void Tenseur::set_etat_nondef() { 
    del_t() ;
    etat = ETATNONDEF ;
}

// Allocates everything
// --------------------
void Tenseur::allocate_all() {
    
	set_etat_qcq() ; 
	for (int i=0 ; i<n_comp ; i++) {
	    c[i]->allocate_all() ; 
	}
	
} 



void Tenseur::change_triad(const Base_vect& bi) {
    
    bi.change_basis(*this) ; 
    
}

void Tenseur::set_triad(const Base_vect& bi) {
    
    triad = &bi ; 
    
}

void Tenseur::set_poids(double weight) {

    poids = weight ;
}

void Tenseur::set_metric(const Metrique& met) {

    metric = &met ;
}

int Tenseur::donne_place (const Itbl& idx) const {
    
    assert (idx.get_ndim() == 1) ;
    assert (idx.get_dim(0) == valence) ;
    
    for (int i=0 ; i<valence ; i++)
	assert ((idx(i)>=0) && (idx(i)<3)) ;
    int res = 0 ;
    for (int i=0 ; i<valence ; i++)
        res = 3*res+idx(i) ;
    
    return res;
}

Itbl Tenseur::donne_indices (int place) const {
    
    assert ((place >= 0) && (place < n_comp)) ;

    Itbl res(valence) ;
    res.set_etat_qcq() ;
    	    
    for (int i=valence-1 ; i>=0 ; i--) {
	res.set(i) = div(place, 3).rem ;
	place = int((place-res(i))/3) ;
	}
    return res ;
}

void Tenseur::operator=(const Tenseur& t) {
    
    assert (valence == t.valence) ;

    triad = t.triad ; 
    poids = t.poids ;
    metric = t.metric ;

    for (int i=0 ; i<valence ; i++)
	assert (t.type_indice(i) == type_indice(i)) ;
	
    switch (t.etat) {
	case ETATNONDEF: {
	    set_etat_nondef() ;
	    break ;
	}
	
	case ETATZERO: {
	    set_etat_zero() ;
	    break ;
	}
	
	case ETATQCQ: {
	    set_etat_qcq() ;
	    for (int i=0 ; i<n_comp ; i++) {
		int place_t = t.donne_place(donne_indices(i)) ;
		if (t.c[place_t] == 0x0)
		    c[i] = 0x0 ;
		else 
		    *c[i] = *t.c[place_t] ;
	    }
	    break ;
	}
	
	default: {
	    cout << "Unknown state in Tenseur::operator= " << endl ;
	    abort() ;
	    break ;
	    }
    }
}


void Tenseur::operator=(const Cmp& ci) {
    
    assert (valence == 0) ;
    poids = 0. ;
    metric = 0x0 ;

    switch (ci.get_etat()) {
	case ETATNONDEF: {
	    set_etat_nondef() ;
	    break ;
	}
	
	case ETATZERO: {
	    set_etat_zero() ;
	    break ;
	}
	
	case ETATQCQ: {
	    set_etat_qcq() ;
	    *(c[0]) = ci ; 
	    break ;
	}
	
	default: {
	    cout << "Unknown state in Tenseur::operator= " << endl ;
	    abort() ;
	    break ;
	}
    }
}

void Tenseur::operator=(double x) {

    poids = 0. ;
    metric = 0x0 ;
    if (x == double(0)) {
	set_etat_zero() ;
    }
    else {
	assert(valence == 0) ; 
	set_etat_qcq() ;
	*(c[0]) = x ; 
    }

}

void Tenseur::operator=(int x) {

    poids = 0. ;
    metric = 0x0 ;
    if (x == 0) {
	set_etat_zero() ;
    }
    else {
	assert(valence == 0) ; 
	set_etat_qcq() ;
	*(c[0]) = x ; 
    }

}


// Affectation d'un scalaire ...
Cmp& Tenseur::set () {
    
    del_derive() ;
    assert(etat == ETATQCQ) ;
    assert (valence == 0) ;  
    return *c[0] ;
}


// Affectation d'un vecteur :
Cmp& Tenseur::set (int ind) {
    
    del_derive() ;
    assert(valence == 1) ;
    assert (etat == ETATQCQ) ;
    assert ((ind >= 0) && (ind < 3)) ;
    
    return *c[ind] ;
}

// Affectation d'un tenseur d'ordre 2 :
Cmp& Tenseur::set (int ind1, int ind2) {
    
    del_derive() ;
    assert (valence == 2) ;
    assert (etat == ETATQCQ) ;
    assert ((ind1 >= 0) && (ind1 < 3)) ;
    assert ((ind2 >= 0) && (ind2 < 3)) ;
    
    Itbl ind (valence) ;
    ind.set_etat_qcq() ;
    ind.set(0) = ind1 ;
    ind.set(1) = ind2 ;
    
    int place = donne_place(ind) ;
    
    return *c[place] ;
}

// Affectation d'un tenseur d'ordre 3 :
Cmp& Tenseur::set (int ind1, int ind2, int ind3) {
    
    del_derive() ;
    assert (valence == 3) ;
    assert (etat == ETATQCQ) ;
    assert ((ind1 >= 0) && (ind1 < 3)) ;
    assert ((ind2 >= 0) && (ind2 < 3)) ;
    assert ((ind3 >= 0) && (ind3 < 3)) ;
    
    Itbl indices(valence) ;
    indices.set_etat_qcq() ;
    indices.set(0) = ind1 ;
    indices.set(1) = ind2 ;
    indices.set(2) = ind3 ;
    int place = donne_place(indices) ;
 
    return *c[place] ;
}

// Affectation cas general
Cmp& Tenseur::set(const Itbl& indices) {
    
    assert (indices.get_ndim() == 1) ;
    assert (indices.get_dim(0) == valence) ;
    
    del_derive() ;
    assert (etat == ETATQCQ) ;
    for (int i=0 ; i<valence ; i++)
	assert ((indices(i)>=0) && (indices(i)<3)) ;
	
    int place = donne_place(indices) ;
    
    return *c[place] ;
}

// Annulation dans des domaines
void Tenseur::annule(int l) {
    
    annule(l, l) ;     
}

void Tenseur::annule(int l_min, int l_max) {
    
    // Cas particulier: annulation globale : 
    if ( (l_min == 0) && (l_max == mp->get_mg()->get_nzone()-1) ) {
	set_etat_zero() ;
	return ; 
    }
    
    assert( etat != ETATNONDEF ) ; 
    
    if ( etat == ETATZERO ) {
	return ;		// rien n'a faire si c'est deja zero
    }
    else {
	assert( etat == ETATQCQ ) ;	// sinon...
	
	// Annulation des composantes:
	for (int i=0 ; i<n_comp ; i++) {
	    c[i]->annule(l_min, l_max) ; 
	}
	
	// Annulation des membres derives
	if (p_gradient != 0x0) p_gradient->annule(l_min, l_max) ;
	if (p_gradient_spher != 0x0) p_gradient_spher->annule(l_min, l_max) ;
	for (int j=0; j<N_MET_MAX; j++) {
	  if (p_derive_cov[j] != 0x0) p_derive_cov[j]->annule(l_min, l_max) ;
	  if (p_derive_con[j] != 0x0) p_derive_con[j]->annule(l_min, l_max) ;
	  if (p_carre_scal[j] != 0x0) p_carre_scal[j]->annule(l_min, l_max) ;
	}

    }
    
}




// Exctraction :
const Cmp& Tenseur::operator()() const {
    
    assert(valence == 0) ;
    
    if (etat == ETATQCQ) return *c[0] ;	    // pour la performance,
					    // ce cas est traite en premier,
					    // en dehors du switch
    switch (etat) {
	
	case ETATNONDEF : {
	    cout << "Undefined Tensor in Tenseur::operator() ..." << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	    
	case ETATZERO : {
	    return mp->cmp_zero() ;
	    }
	    
	    
	default : {
	    cout <<"Unknown state in Tenseur::operator()" << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	}
}


const Cmp& Tenseur::operator() (int indice) const {
    
    assert ((indice>=0) && (indice<3)) ;
    assert(valence == 1) ;
    
    if (etat == ETATQCQ) return *c[indice] ;	 // pour la performance,
					    // ce cas est traite en premier,
					    // en dehors du switch
    switch (etat) {
	
	case ETATNONDEF : {
	    cout << "Undefined Tensor in Tenseur::operator(int) ..." << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	    
	case ETATZERO : {
	    return mp->cmp_zero() ;
	    }
	    	    
	default : {
	    cout <<"Unknown state in Tenseur::operator(int)" << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	}
}

const Cmp& Tenseur::operator() (int indice1, int indice2) const {
    
    assert ((indice1>=0) && (indice1<3)) ;
    assert ((indice2>=0) && (indice2<3)) ;
    assert(valence == 2) ;
    
    if (etat == ETATQCQ) {		// pour la performance,
	    Itbl idx(2) ;		// ce cas est traite en premier,
	    idx.set_etat_qcq() ;	// en dehors du switch
	    idx.set(0) = indice1 ;
	    idx.set(1) = indice2 ;
	    return *c[donne_place(idx)] ;	
    }

    switch (etat) {
	
	case ETATNONDEF : {
	    cout << "Undefined Tensor in Tenseur::operator(int, int) ..." << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	    
	case ETATZERO : {
	    return mp->cmp_zero() ;
	    }
	   	    
	default : {
	    cout <<"Unknown state in Tenseur::operator(int, int)" << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	} 
}

const Cmp& Tenseur::operator() (int indice1, int indice2, int indice3) const {
    
    assert ((indice1>=0) && (indice1<3)) ;
    assert ((indice2>=0) && (indice2<3)) ;
    assert ((indice3>=0) && (indice3<3)) ;
    assert(valence == 3) ;
    
    if (etat == ETATQCQ) {		// pour la performance,
	    Itbl idx(3) ;		// ce cas est traite en premier,
	    idx.set_etat_qcq() ;	// en dehors du switch
	    idx.set(0) = indice1 ;
	    idx.set(1) = indice2 ;
	    idx.set(2) = indice3 ;
	    return *c[donne_place(idx)] ;
    }

    switch (etat) {
	
	case ETATNONDEF : {
	    cout << "Undefined Tensor in Tenseur::operator(int, int, int) ..." << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	    
	case ETATZERO : {
	    return mp->cmp_zero() ;
	    }
	    	    
	default : {
	    cout <<"Unknown state in Tenseur::operator(int, int, int)" << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	}
}


const Cmp& Tenseur::operator() (const Itbl& indices) const {
    
    assert (indices.get_ndim() == 1) ;
    assert (indices.get_dim(0) == valence) ;
    for (int i=0 ; i<valence ; i++)
	assert ((indices(i)>=0) && (indices(i)<3)) ;
    
    if (etat == ETATQCQ) {		    // pour la performance,
	return *c[donne_place(indices)]	;   // ce cas est traite en premier,
    }					    // en dehors du switch

    switch (etat) {
	
	case ETATNONDEF : {
	    cout << "Undefined Tensor in Tenseur::operator(const Itbl&) ..." << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	    
	case ETATZERO : {
	    return mp->cmp_zero() ;
	    }
	    
	default : {
	    cout <<"Unknown state in Tenseur::operator(const Itbl& )" << endl ;
	    abort() ;
	    return *c[0] ;  // bidon pour satisfaire le compilateur
	    }
	}
    
}

// Gestion de la ZEC :
void Tenseur::dec_dzpuis() {
    
    if (etat == ETATZERO) {
	return ; 
    }
    
    assert(etat == ETATQCQ) ;
   
    for (int i=0 ; i<n_comp ; i++)
	if (c[i] != 0x0)
	    c[i]->dec_dzpuis() ;
}

void Tenseur::inc_dzpuis() {
    
    if (etat == ETATZERO) {
	return ; 
    }

    assert(etat == ETATQCQ) ;
    
    for (int i=0 ; i<n_comp ; i++)
	if (c[i] != 0x0)
	    c[i]->inc_dzpuis() ;
}

void Tenseur::dec2_dzpuis() {
    
    if (etat == ETATZERO) {
	return ; 
    }
    
    assert(etat == ETATQCQ) ;
   
    for (int i=0 ; i<n_comp ; i++)
	if (c[i] != 0x0)
	    c[i]->dec2_dzpuis() ;
}

void Tenseur::inc2_dzpuis() {
    
    if (etat == ETATZERO) {
	return ; 
    }

    assert(etat == ETATQCQ) ;
    
    for (int i=0 ; i<n_comp ; i++)
	if (c[i] != 0x0)
	    c[i]->inc2_dzpuis() ;
}

void Tenseur::mult_r_zec() {
    
    if (etat == ETATZERO) {
	return ; 
    }

    assert(etat == ETATQCQ) ;
    
    for (int i=0 ; i<n_comp ; i++) 
    	if (c[i] != 0x0)
	    c[i]->mult_r_zec() ;
}

// Gestion des bases spectrales (valence <= 2)
void Tenseur::set_std_base() {
    
    if (etat == ETATZERO) {
	return ; 
    }
    
    assert(etat == ETATQCQ) ;
    switch (valence) {
	
	case 0 : {
	    c[0]->std_base_scal() ;
	    break ;
	}
	    
	case 1 : {

	  if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {

	    Base_val** bases = mp->get_mg()->std_base_vect_cart() ;

	    for (int i=0 ; i<3 ; i++)
		(c[i]->va).set_base( *bases[i] ) ;
	    for (int i=0 ; i<3 ; i++)
		delete bases[i] ;
	    delete [] bases ;
	  }
	  else {
	    assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
	    Base_val** bases = mp->get_mg()->std_base_vect_spher() ;

	    for (int i=0 ; i<3 ; i++)
		(c[i]->va).set_base( *bases[i] ) ;
	    for (int i=0 ; i<3 ; i++)
		delete bases[i] ;
	    delete [] bases ;
	  }
	  break ;
	    
	}
	    
	case 2 : {

	  if( triad->identify() == (mp->get_bvect_cart()).identify() ) {

	    Base_val** bases = mp->get_mg()->std_base_vect_cart() ;
	    
	    Itbl indices (2) ;
	    indices.set_etat_qcq() ;
	    for (int i=0 ; i<n_comp ; i++) {   
		indices = donne_indices(i) ;
		(c[i]->va).set_base( (*bases[indices(0)]) * 
				     (*bases[indices(1)]) ) ;
	    }
	    for (int i=0 ; i<3 ; i++)
		delete bases[i] ;
	    delete [] bases ;
	  }
	  else {
	    assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
	    Base_val** bases = mp->get_mg()->std_base_vect_spher() ;
	    
	    Itbl indices (2) ;
	    indices.set_etat_qcq() ;
	    for (int i=0 ; i<n_comp ; i++) {   
		indices = donne_indices(i) ;
		(c[i]->va).set_base( (*bases[indices(0)]) * 
				     (*bases[indices(1)]) ) ;
	    }
	    for (int i=0 ; i<3 ; i++)
		delete bases[i] ;
	    delete [] bases ;
	  }
	  break ;
	}
	   
	    default : {
		cout << 
		"Tenseur::set_std_base() : the case valence = " << valence
		 << " is not treated !" << endl ;
		abort() ;
		break ;
	    }
	}
}

 // Le cout :
ostream& operator<<(ostream& flux, const Tenseur &source ) {
    
    flux << "Valence : " << source.valence << endl ;
    if (source.get_poids() != 0.) 
      flux << "Tensor density of weight " << source.poids << endl ;

    if (source.get_triad() != 0x0) {
	flux << "Vectorial basis (triad) on which the components are defined :" 
	     << endl ; 
	flux << *(source.get_triad()) << endl ;
    }
    
    if (source.valence != 0)
	flux << "Type of the indices : " << endl ;
    for (int i=0 ; i<source.valence ; i++) {
	flux << "Index " << i << " : " ;
	if (source.type_indice(i) == CON)
	    flux << " contravariant." << endl ;
	else
	    flux << " covariant." << endl ;
	}
    
    switch (source.etat) {
	
	case ETATZERO : {
	    flux << "Null Tenseur. " << endl ;
	    break ;
	    }
	
	case ETATNONDEF : {
	    flux << "Undefined Tenseur. " << endl ;
	    break ;
	    }
	
	case ETATQCQ : {
	    for (int i=0 ; i<source.n_comp ; i++) {

		Itbl num_indices (source.donne_indices(i)) ;
		flux << "Component " ;
		
		if (source.valence != 0) {
		for (int j=0 ; j<source.valence ; j++)
		    flux << "  " << num_indices(j) ;
		    }
		else
		    flux << "  " << 0 ;
		flux << " : " << endl ;
		flux << "-------------" << endl ; 


		if (source.c[i] != 0x0)
		    flux << *source.c[i] << endl ;
		else
		    flux << "Unknown component ... " << endl ;
		    
	    }
	    break ;
	    }
	default : {
	    cout << "Unknown case in operator<< (ostream&, const Tenseur&)" << endl ;
	    abort() ;
	    break ;
	}
    }
    
    flux << " -----------------------------------------------------" << endl ;
    return flux ;
}

void Tenseur::sauve(FILE* fd) const {
    
    type_indice.sauve(fd) ;	// type des composantes
    fwrite_be(&valence, sizeof(int), 1, fd) ;    // la valence
    
    if (valence != 0) {
	triad->sauve(fd) ;	    // Vectorial basis
    }
    
    fwrite_be(&n_comp, sizeof(int), 1, fd) ; // nbre composantes
    fwrite_be(&etat, sizeof(int), 1, fd) ; // etat
   
    if (etat == ETATQCQ)
	for (int i=0 ; i<n_comp ; i++)
	    c[i]->sauve(fd) ;
    if (poids != 0.)
      fwrite_be(&poids, sizeof(double), 1, fd) ; //poids, si pertinent
}


void Tenseur::fait_gradient () const {
    
    assert (etat != ETATNONDEF) ;
    
    if (p_gradient != 0x0)
	return ;
    else {
 
	// Construction du resultat :
	Itbl tipe (valence+1) ;
	tipe.set_etat_qcq() ;
	tipe.set(0) = COV ;
	for (int i=0 ; i<valence ; i++)
	    tipe.set(i+1) = type_indice(i) ;
	
	// Vectorial basis
	// ---------------
	if (valence != 0) {
	  //	    assert (*triad == mp->get_bvect_cart()) ;
	}

	p_gradient = new Tenseur(*mp, valence+1, tipe, mp->get_bvect_cart(), 
				 metric, poids) ;
    
	if (etat == ETATZERO)
	    p_gradient->set_etat_zero() ;
	else {
	    p_gradient->set_etat_qcq() ;
	    // Boucle sur les indices :
	
	    Itbl indices_source(valence) ;
	    indices_source.set_etat_qcq() ;
	
	    for (int j=0 ; j<p_gradient->n_comp ; j++) {    
		
		Itbl indices_res(p_gradient->donne_indices(j)) ;

		for (int m=0 ; m<valence ; m++)
		    indices_source.set(m) = indices_res(m+1) ;
		 
		p_gradient->set(indices_res) =  
		    (*this)(indices_source).deriv(indices_res(0)) ;
		}
	  }
    }
}

void Tenseur::fait_gradient_spher () const {
    
    assert (etat != ETATNONDEF) ;
    
    if (p_gradient_spher != 0x0)
	return ;
    else {
 
	// Construction du resultat :
	
	if (valence != 0) {
	    cout << 
	    "Tenseur::fait_gradient_spher : the valence must be zero !" 
	    << endl ; 
	    abort() ; 
	}

	p_gradient_spher = new Tenseur(*mp, 1, COV, mp->get_bvect_spher(),
				       metric, poids) ;
    
	if (etat == ETATZERO) {
	    p_gradient_spher->set_etat_zero() ;
	}
	else {
	    assert( etat == ETATQCQ ) ; 
	    p_gradient_spher->set_etat_qcq() ;

	    p_gradient_spher->set(0) = c[0]->dsdr() ;	    // d/dr 
	    p_gradient_spher->set(1) = c[0]->srdsdt() ;	    // 1/r d/dtheta
	    p_gradient_spher->set(2) = c[0]->srstdsdp() ;   // 1/(r sin(theta))
							    //		  d/dphi

	}
    }
}


void Tenseur::fait_derive_cov (const Metrique& metre, int ind) const {
    
  assert (etat != ETATNONDEF) ;
  assert (valence != 0) ;
  
  if (p_derive_cov[ind] != 0x0)
    return ;
  else {
    
    p_derive_cov[ind] = new Tenseur (gradient()) ;
    
    if ((valence != 0) && (etat != ETATZERO)) {

      
      //    assert( *(metre.gamma().get_triad()) == *triad ) ; 
      Tenseur* auxi ;
      for (int i=0 ; i<valence ; i++) {
	
	if (type_indice(i) == COV) {
	  auxi = new Tenseur(contract(metre.gamma(), 0,(*this), i)) ;
	  
	  Itbl indices_gamma(p_derive_cov[ind]->valence) ;
	  indices_gamma.set_etat_qcq() ;
	  //On range comme il faut :
	  for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
	    
	    Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
	    indices_gamma.set(0) = indices(0) ;
	    indices_gamma.set(1) = indices(i+1) ;
	    for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
	      if (idx<=i+1)
		indices_gamma.set(idx) = indices(idx-1) ;
	      else
		indices_gamma.set(idx) = indices(idx) ;
	    
	      p_derive_cov[ind]->set(indices) -= (*auxi)(indices_gamma) ;
	  }
	}   
	else {
	  auxi = new Tenseur(contract(metre.gamma(), 1, (*this), i)) ;
	  
	  Itbl indices_gamma(p_derive_cov[ind]->valence) ;
	  indices_gamma.set_etat_qcq() ;
	  
	  //On range comme il faut :
	  for (int j=0 ; j<p_derive_cov[ind]->n_comp ; j++) {
	    
	    Itbl indices (p_derive_cov[ind]->donne_indices(j)) ;
	    indices_gamma.set(0) = indices(i+1) ;
	    indices_gamma.set(1) = indices(0) ;
	    for (int idx=2 ; idx<p_derive_cov[ind]->valence ; idx++)
	      if (idx<=i+1)
		indices_gamma.set(idx) = indices(idx-1) ;
	      else
		indices_gamma.set(idx) = indices(idx) ;
	    p_derive_cov[ind]->set(indices) += (*auxi)(indices_gamma) ;
	  }
	}
	delete auxi ;
      }
    }
    if ((poids != 0.)&&(etat != ETATZERO)) 
      *p_derive_cov[ind] = *p_derive_cov[ind] - 
	poids*contract(metre.gamma(), 0, 2) * (*this) ;
    
  }
}



void Tenseur::fait_derive_con (const Metrique& metre, int ind) const {
    
  if (p_derive_con[ind] != 0x0)
    return ;
  else {
    // On calcul la derivee covariante :
    if (valence != 0)
      p_derive_con[ind] = new Tenseur
	(contract(metre.con(), 1, derive_cov(metre), 0)) ;
    
    else {
      Tenseur grad (gradient()) ;
      grad.change_triad( *(metre.con().get_triad()) ) ;
      p_derive_con[ind] = new Tenseur (contract(metre.con(), 1, grad, 0)) ;
    }
  }
}

void Tenseur::fait_carre_scal (const Metrique& met, int ind) const {
    
    if (p_carre_scal[ind] != 0x0)
	return ;
    else {
	assert (valence != 0) ;   // A ne pas appeler sur un scalaire ;
       
	// On bouge tous les indices :
	Tenseur op_t(manipule(*this, met)) ;
   
	Tenseur* auxi = new Tenseur(contract(*this, 0, op_t, 0)) ;
	Tenseur* auxi_old ;
    
	// On contracte tous les indices restant :
	for (int indice=1 ; indice<valence ; indice++) {
	    auxi_old = new Tenseur(contract(*auxi, 0, valence-indice)) ;
	    delete auxi ;
	    auxi = new Tenseur(*auxi_old) ;
	    delete auxi_old ;
	}
	p_carre_scal[ind] = new Tenseur (*auxi) ;
	delete auxi ;
    }
}
    
const Tenseur& Tenseur::gradient () const {
    if (p_gradient == 0x0)
	fait_gradient() ;
    return *p_gradient ;
}

const Tenseur& Tenseur::gradient_spher() const {
    if (p_gradient_spher == 0x0)
	fait_gradient_spher() ;
    return *p_gradient_spher ;
}

const Tenseur& Tenseur::derive_cov (const Metrique& metre) const {
    
    if (valence == 0)
	return gradient() ;
    else {
	set_dependance(metre) ;
	int j = get_place_met(metre) ;
	assert(j!=-1) ;
	if (p_derive_cov[j] == 0x0)
	    fait_derive_cov (metre,j) ;
	return *p_derive_cov[j] ;
    }
}

const Tenseur& Tenseur::derive_con (const Metrique& metre) const {
    set_dependance(metre) ;
    int j = get_place_met(metre) ;
    assert(j!=-1) ;
    if (p_derive_con[j] == 0x0)
	fait_derive_con (metre, j) ;
    return *p_derive_con[j] ;
}

const Tenseur& Tenseur::carre_scal (const Metrique& metre) const {
    set_dependance(metre) ;
    int j = get_place_met(metre) ;
    assert(j!=-1) ;
    if (p_carre_scal[j] == 0x0)
	fait_carre_scal (metre, j) ;
    return *p_carre_scal[j] ;
}
}