tenseur.h

/*
 *  Definition of Lorene classes Tenseur
 *               Tenseur_sym
 *
 */

/*
 *   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
 *
 */


#ifndef __TENSEUR_H_
#define __TENSEUR_H_


/*
 * $Id: tenseur.h,v 1.17 2010/02/02 13:34:12 e_gourgoulhon Exp $
 * $Log: tenseur.h,v $
 * Revision 1.17  2010/02/02 13:34:12  e_gourgoulhon
 * Marked DEPRECATED (in the documentation).
 *
 * Revision 1.16  2005/08/30 08:35:10  p_grandclement
 * Addition of the Tau version of the vectorial Poisson equation for the Tensors
 *
 * Revision 1.15  2004/12/29 16:24:03  k_taniguchi
 * Addition of the method for vectorial Poisson equations with a multipole
 * falloff condition at the outer boundary.
 *
 * Revision 1.14  2004/12/22 18:25:12  k_taniguchi
 * Change an argument of poisson_vect_falloff
 *
 * Revision 1.13  2004/11/30 20:43:32  k_taniguchi
 * Addition of the method for vectorial Poisson equations with falloff
 * condition at the outer boundary.
 *
 * Revision 1.12  2004/03/22 13:12:43  j_novak
 * Modification of comments to use doxygen instead of doc++
 *
 * Revision 1.11  2003/11/06 12:17:31  r_prix
 * fixed mini-bug in documentation: without explicit argument in function-prototype,
 * doc++ seemed to merge docu of Tensor::operator=(Cmp&) and operator=(Tenseur&)
 *
 * Revision 1.10  2003/06/20 14:23:38  f_limousin
 * Add the functions compare().
 *
 * Revision 1.9  2002/10/16 14:36:29  j_novak
 * Reorganization of #include instructions of standard C++, in order to
 * use experimental version 3 of gcc.
 *
 * Revision 1.8  2002/09/19 14:12:37  e_gourgoulhon
 * Modif documentation for LaTeX compliance.
 *
 * Revision 1.7  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.6  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.5  2002/08/14 13:46:14  j_novak
 * Derived quantities of a Tenseur can now depend on several Metrique's
 *
 * Revision 1.4  2002/08/08 15:10:44  j_novak
 * The flag "plat" has been added to the class Metrique to show flat metrics.
 *
 * Revision 1.3  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.2  2002/06/17 14:05:17  j_novak
 * friend functions are now also declared outside the class definition
 *
 * Revision 1.1.1.1  2001/11/20 15:19:27  e_gourgoulhon
 * LORENE
 *
 * Revision 2.46  2001/08/27  10:03:56  eric
 * Ajout de l'operator% (produit tensoriel avec desaliasing)
 *
 * Revision 2.45  2001/06/19  15:38:00  eric
 * Modif commentaires: mise en conformite Doc++ 3.4.8.
 *
 * Revision 2.44  2001/06/18  13:56:07  novak
 * Ajout de la fonction abs()
 *
 * Revision 2.43  2001/05/29 16:10:43  eric
 * Modif commentaires (mise en conformite Doc++ 3.4.7).
 *
 * Revision 2.42  2001/05/26  15:48:17  eric
 * *** empty log message ***
 *
 * Revision 2.41  2001/05/26  15:42:54  eric
 * Ajout de la fonction flat_scalar_prod_desal (desaliasage)
 *
 * Revision 2.40  2000/10/19  10:37:51  phil
 * *** empty log message ***
 *
 * Revision 2.39  2000/10/19  09:47:15  phil
 * *** empty log message ***
 *
 * Revision 2.38  2000/10/19  09:41:48  phil
 * ajout de inverse_poisson_vect
 *
 * Revision 2.37  2000/10/06  12:37:09  keisuke
 * Add a vectorial Poisson equation Tenseur::poisson_vect_regu.
 *
 * Revision 2.36  2000/09/27  08:52:34  eric
 * Modif commentaires.
 *
 * Revision 2.35  2000/09/13  12:21:36  eric
 * Modif commentaires.
 *
 * Revision 2.34  2000/09/13  12:11:26  eric
 * Ajout de la fonction allocate_all().
 *
 * Revision 2.33  2000/05/22  14:39:11  phil
 * ajout de inc_dzpuis et dec_dzpuis
 *
 * Revision 2.32  2000/04/03  15:18:54  phil
 * suppression de poisson_vect_dirichlet
 *
 * Revision 2.31  2000/03/31  13:29:43  phil
 * poisson_vect_neumann devient poisson_vect_dirichlet
 *
 * Revision 2.30  2000/03/30  16:10:46  phil
 * *** empty log message ***
 *
 * Revision 2.29  2000/02/18  10:45:05  eric
 * Modif commentaires.
 *
 * Revision 2.28  2000/02/11  19:11:55  phil
 * commentaires
 *
 * Revision 2.27  2000/02/10  16:26:48  eric
 * Modif commentaires.
 *
 * Revision 2.26  2000/02/10  16:10:49  eric
 * Ajout de la fonction change_triad.
 *
 * Revision 2.25  2000/02/09  19:29:25  eric
 * MODIF IMPORTANTE: la triade de decomposition est desormais passee en
 * argument des constructeurs.
 *
 * Revision 2.24  2000/02/09  09:53:56  phil
 * ajout de poisson_vect_oohara
 * ,
 *
 * Revision 2.23  2000/02/08  19:04:25  eric
 * Les fonctions arithmetiques ne sont plus amies.
 * Ajout de nouvelles fonctions arithmetiques.
 *
 * Revision 2.22  2000/02/01  15:40:19  eric
 * Ajout de la fonction sqrt
 *
 * Revision 2.21  2000/02/01  14:13:56  eric
 * Modif commentaires.
 * Ajout de la fonction amie flat_scalar_prod.
 *
 * Revision 2.20  2000/01/21  12:48:55  phil
 * changement prototypage de Tenseur::poisson_vect
 *
 * Revision 2.19  2000/01/20  16:02:20  eric
 * Ajout des operator=(double ) et operator=(int ).
 *
 * Revision 2.18  2000/01/20  14:52:08  phil
 * *** empty log message ***
 *
 * Revision 2.17  2000/01/20  14:50:56  phil
 * *** empty log message ***
 *
 * Revision 2.16  2000/01/20  14:39:31  phil
 * *** empty log message ***
 *
 * Revision 2.15  2000/01/20  13:10:56  phil
 * Ajout de Tenseur::poisson_vect (double)
 *
 * Revision 2.14  2000/01/20  11:20:17  phil
 * changement prototypage
 *
 * Revision 2.13  2000/01/20  10:31:30  phil
 * ajout de xksk
 *
 * Revision 2.12  2000/01/14  14:17:47  eric
 * Modif commentaires.
 *
 * Revision 2.11  2000/01/14  14:04:07  eric
 * Ajout de la fonction annule.
 * classe Tenseur: constructeurs pour les classes derivees et fonctions
 *  de gestion memoire (del_t(), ...) declarees protected et non plus
 *  private.
 *
 * Revision 2.10  2000/01/13  14:15:30  eric
 * Modif commentaires.
 *
 * Revision 2.9  2000/01/13  14:10:18  eric
 * Ajout du constructeur par copie d'un Cmp (pour un scalaire)
 * ainsi que l'affectation a un Cmp.
 *
 * Revision 2.8  2000/01/13  13:45:56  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.7  2000/01/12  13:17:49  eric
 * Les operator::(...) renvoie desormais une reference const sur le c[...]
 * correspondant et non plus un Cmp copie de c[...].
 * (ceci grace a Map::cmp_zero()).
 *
 * Revision 2.6  2000/01/11  11:13:16  eric
 * Changement de nom pour la base vectorielle : base --> triad
 *
 * Revision 2.5  2000/01/10  17:22:26  eric
 *  Modif des #include
 *
 * Revision 2.4  2000/01/10  15:14:58  eric
 * Ajout du membre base (base vectorielle sur laquelle sont definies
 *   les composantes).
 *
 * Revision 2.3  1999/12/09  12:39:39  phil
 * changement prototypage des derivees
 *
 * Revision 2.2  1999/12/07  15:24:17  phil
 * ajout include
 *
 * Revision 2.1  1999/12/03  09:37:11  phil
 * *** empty log message ***
 *
 * Revision 2.0  1999/12/02  17:15:32  phil
 * *** empty log message ***
 *
 * Revision 1.1  1999/12/02  17:13:29  phil
 * Initial revision
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/tenseur.h,v 1.17 2010/02/02 13:34:12 e_gourgoulhon Exp $
 *
 */

#define COV -1
#define CON +1

#define N_MET_MAX 5

// Headers Lorene 
#include "cmp.h"
#include "itbl.h"
#include "base_vect.h"

class Metrique ;
class Tenseur_sym ;

            //---------------------------------//
            //  class Tenseur          //
            //---------------------------------//
            

class Tenseur { 

    // Data : 
    // -----
    protected:
    const Map* const mp ;   
    int valence ;       
    
    const Base_vect* triad ; 

    Itbl type_indice ;  
    
    int n_comp ;    
    int etat ;  
    Cmp** c ;   
    double poids ; 

    const Metrique* metric ;      
    
    // Derived data : 
    // ------------
    protected:
    const Metrique** met_depend ;
    
    mutable Tenseur* p_gradient ;
    
    mutable Tenseur* p_gradient_spher ;
    
    Tenseur** p_derive_cov ;

    Tenseur** p_derive_con ;

    Tenseur** p_carre_scal ;
   
    // Constructors - Destructor :
    // -------------------------
    protected:
    bool verif() const ; 
        
    void new_der_met() ;
    
    public:
    explicit Tenseur (const Map& map, const Metrique* met = 0x0, 
             double weight = 0) ; 

    explicit Tenseur (const Cmp& cmp, const Metrique* met = 0x0, 
             double weight = 0) ; 

    Tenseur (const Map& map, int val, const Itbl& tipe, 
         const Base_vect& triad_i, const Metrique* met = 0x0, 
             double weight = 0) ;

    Tenseur (const Map& map, int val, const Itbl& tipe, 
         const Base_vect* triad_i, const Metrique* met = 0x0, 
             double weight = 0) ;

    Tenseur (const Map& map, int val, int tipe, const 
         Base_vect& triad_i, const Metrique* met = 0x0, 
             double weight = 0) ;

    Tenseur (const Tenseur&) ;  

    explicit Tenseur (const Tenseur_sym&) ;
    
    Tenseur (const Map& map, const Base_vect& triad_i, FILE* fich, 
         const Metrique* met = 0x0) ;

    Tenseur (const Map& map, FILE* fich, const Metrique* met = 0x0) ;

    
    protected:
    Tenseur (const Map& map, int val, const Itbl& tipe, int n_comp,
         const Base_vect& triad_i, const Metrique* met = 0x0, 
             double weight = 0) ;

    Tenseur (const Map&, int val, int tipe, int n_comp, 
         const Base_vect& triad_i, const Metrique* met = 0x0, 
             double weight = 0) ;


    public: 

    virtual ~Tenseur() ;    
    
    // Memory management
    // -----------------
    protected:
    void del_t() ;  

    void del_derive_met(int i) const ;

    void del_derive() const ;
    
    void set_der_met_0x0(int i) const ;

    void set_der_0x0() const ;

    // Mutators / assignment
    // ---------------------
    public:
    void set_etat_nondef() ;
    
    void set_etat_zero() ;

    void set_etat_qcq() ;
    
    void allocate_all() ; 

    void change_triad(const Base_vect& new_triad) ; 
    
    void set_triad(const Base_vect& new_triad) ; 
    void set_poids(double weight) ; 

    void set_metric(const Metrique& met) ;
    
    virtual void operator=(const Tenseur& tens) ; 
    
    void operator=(const Cmp& field) ; 
    
    void operator=(double ) ;   

    void operator=(int ) ;  

    Cmp& set () ;  
    Cmp& set (int) ; 
    Cmp& set (int, int) ; 
    Cmp& set (int, int, int) ; 
    Cmp& set (const Itbl&) ; 
        
    void annule(int l) ; 

    void annule(int l_min, int l_max) ; 

    void set_std_base() ; 
    
    void dec_dzpuis() ; 
    void inc_dzpuis() ; 
    void dec2_dzpuis() ;    
    void inc2_dzpuis() ;    
    void mult_r_zec() ; 
    
    Tenseur inverse_poisson_vect (double lambda) const ;
    
    // Accessors
    // ---------
    public:
    virtual int donne_place (const Itbl& idx) const ;

    virtual Itbl donne_indices (int place) const ;
    
    const Map* get_mp() const {return mp ;} ; 

    const Base_vect* get_triad() const {return triad;} ; 
    
    int get_etat() const {return etat ;} ;
 
    int get_valence() const {return valence ; } ;

    int get_n_comp() const {return n_comp ;} ; 
    
    int get_type_indice (int i) const {return type_indice(i) ;};

    Itbl get_type_indice () const {return type_indice ; } ;

    double get_poids() const {return poids ; } ; 

    const Metrique* get_metric() const {return metric ; } ; 
    
    const Cmp& operator()() const ; 
    const Cmp& operator()(int) const ; 
    const Cmp& operator()(int, int) const ; 
    const Cmp& operator()(int, int, int) const ; 
    const Cmp& operator()(const Itbl&) const ; 
    
    // Outputs
    // -------
    public:
    void sauve(FILE *) const ;      
    friend ostream& operator<<(ostream& , const Tenseur & ) ;
    
    // Computation of derived members
    // ------------------------------
    protected:
    virtual void fait_gradient () const ;

    void fait_gradient_spher () const ;

    virtual void fait_derive_cov (const Metrique& met, int i) const ;

    virtual void fait_derive_con (const Metrique&, int i) const ;

    void fait_carre_scal (const Metrique&, int i) const ;
    
    void set_dependance (const Metrique& met) const ;

    int get_place_met(const Metrique& metre) const ;
    
    // Differential operators
    // ----------------------
    public:
    const Tenseur& gradient() const ; 
    
    const Tenseur& gradient_spher() const ; 
    
    const Tenseur& derive_cov (const Metrique& met) const ;

    const Tenseur& derive_con (const Metrique&) const ; 

    const Tenseur& carre_scal (const Metrique&) const ;
    
    // Resolution d'EDP :
    void poisson_vect(double lambda, Param& par, Tenseur& shift, Tenseur& vect
            , Tenseur& scal) const ;
    
    /*
    * Same as poisson_vect with a Tau method
    **/
    void poisson_vect_tau(double lambda, Param& par, Tenseur& shift, Tenseur& vect
            , Tenseur& scal) const ;
            
    void poisson_vect_falloff(double lambda, Param& par, Tenseur& shift, 
                  Tenseur& vect, Tenseur& scal, int* k_falloff) const ;

    void poisson_vect_ylm(double lambda, Param& para, Tenseur& shift,
              Tenseur& vecteur, Tenseur& scalaire, int nylm,
              double* intvec) const ;

    Tenseur poisson_vect(double lambda, Tenseur& vect , Tenseur& scal ) const ;
    
    /*
    * Same as poisson_vect with a Tau method
    **/
    Tenseur poisson_vect_tau(double lambda, Tenseur& vect , Tenseur& scal ) const ;
            
    Tenseur poisson_vect_falloff(double lambda, Tenseur& vect , 
                 Tenseur& scal, int* k_falloff ) const ;

    Tenseur poisson_vect_ylm(double lambda, Tenseur& vecteur, 
                 Tenseur& scalaire, int nylm, double* intvec) const ;

    void poisson_vect_oohara(double lambda, Param& par, Tenseur& shift, 
                    Tenseur& scal) const ;
       
    /*
    * Same as poisson_vect_oohara with a Tau method
    **/
    void poisson_vect_oohara_tau(double lambda, Param& par, Tenseur& shift, 
                    Tenseur& scal) const ;
            
    Tenseur poisson_vect_oohara(double lambda, Tenseur& scal) const ;
   /*
    * Same as poisson_vect_oohara with a Tau method
    **/
    Tenseur poisson_vect_oohara_tau(double lambda, Tenseur& scal) const ;
            
    void poisson_vect_regu(int k_div, int nzet, double unsgam1,
                           double lambda, Param& par, Tenseur& shift,
                           Tenseur& vect, Tenseur& scal) const ;
   
    

    void compare( const Tenseur& tens, const char* name ) ;
    void compare( FILE* fich, const char* name_i ) ;
    


    // Friend classes 
    // ---------------
    friend class Tenseur_sym ;
    friend class Metrique ;
    
    // Mathematical operators
    // ----------------------
    
    friend Tenseur operator* (const Tenseur&, const Tenseur&) ; 
    friend Tenseur operator% (const Tenseur&, const Tenseur&) ; 
    friend Tenseur contract(const Tenseur&, int id1, int id2) ;
    friend Tenseur contract(const Tenseur&, int id1, const Tenseur&, 
                int id2) ;
    friend Tenseur contract_desal(const Tenseur&, int id1, const Tenseur&, 
                int id2) ;
    friend Tenseur flat_scalar_prod(const Tenseur& t1, const Tenseur& t2) ;
    friend Tenseur flat_scalar_prod_desal(const Tenseur& t1, 
                      const Tenseur& t2) ;
    friend Tenseur manipule(const Tenseur&, const Metrique&, int idx) ;
    friend Tenseur manipule(const Tenseur&, const Metrique&) ;
    friend Tenseur skxk (const Tenseur&) ;
    friend Tenseur lie_derive(const Tenseur& , const Tenseur& , 
                  const Metrique* ) ;

};



Tenseur operator*(const Tenseur&, const Tenseur&) ; 

Tenseur operator%(const Tenseur&, const Tenseur&) ; 

Tenseur contract(const Tenseur&, int id1, int id2) ;

Tenseur contract(const Tenseur&, int id1, const Tenseur&, int id2) ;

Tenseur flat_scalar_prod(const Tenseur& t1, const Tenseur& t2) ;
    
Tenseur flat_scalar_prod_desal(const Tenseur& t1, const Tenseur& t2) ;
    
Tenseur manipule(const Tenseur&, const Metrique&, int idx) ;

Tenseur manipule(const Tenseur&, const Metrique&) ;
    
Tenseur skxk (const Tenseur&) ;

Tenseur lie_derive (const Tenseur& t, const Tenseur& x, const Metrique* = 0x0);

Tenseur sans_trace(const Tenseur& tens, const Metrique& metre) ;


Tenseur operator+(const Tenseur& ) ;            
Tenseur operator-(const Tenseur& ) ;            
Tenseur operator+(const Tenseur&, const Tenseur &) ;    

Tenseur operator+(const Tenseur&, double ) ;        

Tenseur operator+(double, const Tenseur& ) ;        

Tenseur operator+(const Tenseur&, int ) ;       

Tenseur operator+(int, const Tenseur& ) ;       

Tenseur operator-(const Tenseur &, const Tenseur &) ;   

Tenseur operator-(const Tenseur&, double ) ;        

Tenseur operator-(double, const Tenseur& ) ;        

Tenseur operator-(const Tenseur&, int ) ;       

Tenseur operator-(int, const Tenseur& ) ;       

Tenseur operator*(const Tenseur&, double ) ;        

Tenseur operator*(double, const Tenseur& ) ;        

Tenseur operator*(const Tenseur&, int ) ;       

Tenseur operator*(int, const Tenseur& ) ;       

Tenseur operator/(const Tenseur& a, const Tenseur& b) ; 

Tenseur operator/(const Tenseur&, double ) ;    

Tenseur operator/(double, const Tenseur &) ;        

Tenseur operator/(const Tenseur&, int ) ;       

Tenseur operator/(int, const Tenseur &) ;       

Tenseur exp(const Tenseur& ) ;      
Tenseur log(const Tenseur& ) ;      
Tenseur sqrt(const Tenseur& ) ;     
Tenseur abs(const Tenseur& ) ;      
Tenseur pow(const Tenseur&, int ) ; 
Tenseur pow(const Tenseur&, double ) ;  

            //---------------------------------//
            //  class Tenseur_sym      //
            //---------------------------------//
            
class Tenseur_sym : public Tenseur {

    // Constructors - Destructor :
    // -------------------------
    
    public:
    Tenseur_sym (const Map& map, int val, const Itbl& tipe, 
             const Base_vect& triad_i, const Metrique* met = 0x0,
             double weight = 0) ;

    Tenseur_sym (const Map& map, int val, int tipe, 
             const Base_vect& triad_i, const Metrique* met = 0x0,
             double weight = 0) ;

    Tenseur_sym (const Tenseur_sym&) ; 

    explicit Tenseur_sym (const Tenseur&) ;
    
    Tenseur_sym (const Map& map, const Base_vect& triad_i, FILE* fich,
             const Metrique* met = 0x0) ;

    virtual ~Tenseur_sym() ;    
    
    // Mutators / assignment
    // ---------------------
    public:
    virtual void operator= (const Tenseur&) ;
    

    // Accessors
    // ---------
    public:
    virtual int donne_place (const Itbl& idx) const ;

    virtual Itbl donne_indices (int place) const ;
        
    // Computation of derived members
    // ------------------------------
    protected:
    virtual void fait_gradient () const ;

    virtual void fait_derive_cov (const Metrique& met, int i) const ;

    virtual void fait_derive_con (const Metrique&, int i) const ;
    
    // Mathematical operators
    // ----------------------
    friend Tenseur_sym operator* (const Tenseur&, const Tenseur_sym&) ; 
    friend Tenseur_sym manipule(const Tenseur_sym&, const Metrique&) ;
    friend Tenseur lie_derive (const Tenseur& , const Tenseur& , 
                const Metrique* );
 
} ;

Tenseur_sym operator* (const Tenseur&, const Tenseur_sym&) ; 

Tenseur_sym manipule(const Tenseur_sym&, const Metrique&) ;

Tenseur_sym lie_derive (const Tenseur_sym& t, const Tenseur& x, 
                const Metrique* = 0x0);

Tenseur_sym sans_trace(const Tenseur_sym& tens, const Metrique& metre) ;

Tenseur_sym operator+(const Tenseur_sym& ) ;    
Tenseur_sym operator-(const Tenseur_sym& ) ;    

Tenseur_sym operator+(const Tenseur_sym&, const Tenseur_sym &) ;    

Tenseur_sym operator-(const Tenseur_sym &, const Tenseur_sym &) ;   

Tenseur_sym operator*(const Tenseur_sym&, double ) ;        

Tenseur_sym operator*(double, const Tenseur_sym& ) ;        

Tenseur_sym operator*(const Tenseur_sym&, int ) ;       

Tenseur_sym operator*(int, const Tenseur_sym& ) ;       

Tenseur_sym operator/(const Tenseur_sym& a, const Tenseur& b) ; 

Tenseur_sym operator/(const Tenseur_sym&, double ) ;  

Tenseur_sym operator/(const Tenseur_sym&, int ) ;     

#endif

---