tensor.h

/*
 *  Definition of Lorene classes Tensor and Sym_tensor
 *
 */

/*
 *   Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
 *
 *   Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
 *   Copyright (c) 2000-2001 Eric Gourgoulhon      (for preceding class Tenseur)
 *   Copyright (c) 2002 Jerome Novak               (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
 *
 */


#ifndef __TENSOR_H_ 
#define __TENSOR_H_ 


/*
 * $Id: tensor.h,v 1.63 2023/08/31 08:27:26 g_servignat Exp $
 * $Log: tensor.h,v $
 * Revision 1.63  2023/08/31 08:27:26  g_servignat
 * Added the possibility to filter in the r direction within the ylm filter. An order filtering of 0 means no filtering (for all 3 directions).
 *
 * Revision 1.62  2023/08/28 09:53:33  g_servignat
 * Added ylm filter for Tensor and Scalar in theta + phi directions
 *
 * Revision 1.61  2014/10/13 08:52:37  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.60  2013/06/05 15:43:49  j_novak
 * Suppression of leg_spectral_base()
 *
 * Revision 1.59  2013/01/11 15:44:54  j_novak
 * Addition of Legendre bases (part 2).
 *
 * Revision 1.58  2008/12/05 08:44:02  j_novak
 * New flag to control the "verbosity" of maxabs.
 *
 * Revision 1.57  2007/12/21 16:06:16  j_novak
 * Methods to filter Tensor, Vector and Sym_tensor objects.
 *
 * Revision 1.56  2006/06/07 14:08:58  j_novak
 * New methods set_index_type( / int).
 *
 * Revision 1.55  2005/10/25 08:56:34  p_grandclement
 * addition of std_spectral_base in the case of odd functions near the origin
 *
 * Revision 1.54  2004/07/08 12:21:51  j_novak
 * Replaced tensor::annule_extern_c2 with tensor::annule_extern_cn for a
 * more general transition.
 *
 * Revision 1.53  2004/06/17 06:54:23  e_gourgoulhon
 * Added method annule_extern_c2.
 *
 * Revision 1.52  2004/05/13 21:29:27  e_gourgoulhon
 * Added (external) functions central_value, max_all_domains,
 * min_all_domains and maxabs_all_domains.
 *
 * Revision 1.51  2004/03/24 14:53:39  j_novak
 * Double declarations suppressed
 *
 * Revision 1.50  2004/03/22 13:12:43  j_novak
 * Modification of comments to use doxygen instead of doc++
 *
 * Revision 1.49  2004/02/27 21:12:44  e_gourgoulhon
 * Suppressed function contract_desal (since contract has now the
 * boolean argument "desaliasing").
 *
 * Revision 1.48  2004/02/26 22:44:37  e_gourgoulhon
 * -- constructor of Tensor from Map is now declared explicit.
 * -- class Tensor: added methods compute_derive_lie and derive_lie
 * -- class Tensor_sym: added methods derive_cov, derive_con and derive_lie.
 *
 * Revision 1.47  2004/02/19 22:08:51  e_gourgoulhon
 * Added argument "comment" in method spectral_display,
 * as well as in external functions min, max, maxabs, etc...
 *
 * Revision 1.46  2004/02/18 18:42:41  e_gourgoulhon
 * -- Added methods trace.
 * -- Method scontract suppressed ( since it is the same as trace(int, int) ).
 *
 * Revision 1.45  2004/02/18 15:52:39  e_gourgoulhon
 * -- Added optional argument desaliasing in function contract.
 * -- Added new function contract for double contraction.
 *
 * Revision 1.44  2004/02/16 10:48:06  e_gourgoulhon
 * Added "class Tensor_sym;" at the beginning.
 *
 * Revision 1.43  2004/02/15 21:53:48  e_gourgoulhon
 * Modif. comments: suppressed the mention *** under development ***.
 *
 * Revision 1.42  2004/01/30 12:44:17  e_gourgoulhon
 * Added Tensor_sym operator*(const Tensor_sym&, const Tensor_sym& ).
 *
 * Revision 1.41  2004/01/27 13:05:10  j_novak
 * Removed the method Tensor::mult_r_ced()
 *
 * Revision 1.40  2004/01/19 16:31:40  e_gourgoulhon
 * Added operator()(int, int, int, int) and set(int, int, int, int)
 * for direct access to components of valence 4 tensors.
 *
 * Revision 1.39  2004/01/15 11:09:27  f_limousin
 * Modif in method contract_desal
 *
 * Revision 1.38  2004/01/15 11:00:44  f_limousin
 * Added method contract_desal for the contraction of two tensors with desaliasing
 *
 * Revision 1.37  2004/01/14 11:39:00  f_limousin
 * Added method contract for one tensor
 *
 * Revision 1.36  2004/01/08 09:21:39  e_gourgoulhon
 * Added arithmetics of Tensor_sym.
 * Added arithmetics with Scalar (to solve some ambiguities with respect
 * to the Scalar arithmetics).
 * Added Tensor_sym tensorial product.
 *
 * Revision 1.35  2004/01/04 20:47:37  e_gourgoulhon
 * -- Introduction of new derived class Tensor_sym to store tensor with
 *    two symmetric indices
 * -- Suppression of class Tensor_delta (now a special case of Tensor_sym).
 *
 * Revision 1.34  2003/12/27 14:58:01  e_gourgoulhon
 * Improved documentation. In particular, better description of methods
 * derive_cov(), derive_con() and divergence(), taking into account the
 * new index convention for covariant derivatives.
 *
 * Revision 1.33  2003/12/05 16:41:05  f_limousin
 * Added method operator*
 *
 * Revision 1.32  2003/11/06 14:43:37  e_gourgoulhon
 * Gave a name to const arguments in certain method prototypes (e.g.
 * constructors) to correct a bug of DOC++.
 *
 * Revision 1.31  2003/11/05 15:25:57  e_gourgoulhon
 * Added declaration of external functions:
 * max, min, maxabs, diffrel and diffrelmax.
 *
 * Revision 1.30  2003/11/03 10:58:00  j_novak
 * Suppressed the constructor from a Sym_tensor.
 *
 * Revision 1.29  2003/10/29 11:00:42  e_gourgoulhon
 * Virtual functions dec_dzpuis and inc_dzpuis have now an integer argument to
 *  specify by which amount dzpuis is to be increased.
 * Accordingly virtual methods dec2_dzpuis and inc2_dzpuis have been suppressed.
 *
 * Revision 1.28  2003/10/28 21:21:50  e_gourgoulhon
 * Member function Tensor::contract(int, int) renamed
 *  Tensor::scontract(int, int) in order not to mask
 * the non-member function contract.
 *
 * Revision 1.27  2003/10/27 10:44:00  e_gourgoulhon
 * Declaration of class Sym_tensor is now in file sym_tensor.h.
 *
 * Revision 1.26  2003/10/24 15:00:19  j_novak
 * Forgotten Class declaration... thanks IBM aix!
 *
 * Revision 1.25  2003/10/20 14:26:02  j_novak
 * New assignement operators.
 *
 * Revision 1.24  2003/10/20 09:32:10  j_novak
 * Members p_potential and p_div_free of the Helmholtz decomposition
 * + the method decompose_div(Metric).
 *
 * Revision 1.23  2003/10/19 19:47:31  e_gourgoulhon
 * Introduced new virtual method spectral_display.
 *
 * Revision 1.22  2003/10/16 15:24:30  e_gourgoulhon
 * Name of method annule(int ) changed to annule_domain(int ).
 *
 * Revision 1.21  2003/10/16 14:21:33  j_novak
 * The calculation of the divergence of a Tensor is now possible.
 *
 * Revision 1.20  2003/10/13 13:52:39  j_novak
 * Better managment of derived quantities.
 *
 * Revision 1.19  2003/10/08 14:24:08  j_novak
 * replaced mult_r_zec with mult_r_ced
 *
 * Revision 1.18  2003/10/06 20:48:23  e_gourgoulhon
 * Added methods down and up_down.
 *
 * Revision 1.17  2003/10/06 16:17:29  j_novak
 * Calculation of contravariant derivative and Ricci scalar.
 *
 * Revision 1.16  2003/10/06 15:12:56  e_gourgoulhon
 * Added tensor contraction and raising of index.
 *
 * Revision 1.15  2003/10/06 13:58:45  j_novak
 * The memory management has been improved.
 * Implementation of the covariant derivative with respect to the exact Tensor
 * type.
 *
 * Revision 1.14  2003/10/05 21:07:27  e_gourgoulhon
 * Method std_spectral_base() is now virtual.
 *
 * Revision 1.13  2003/10/03 11:21:45  j_novak
 * More methods for the class Metric
 *
 * Revision 1.12  2003/10/02 15:45:48  j_novak
 * New class Metric
 *
 * Revision 1.11  2003/10/01 15:41:14  e_gourgoulhon
 * class name Delta changed to Tensor_delta.
 *
 * Revision 1.10  2003/10/01 13:03:52  e_gourgoulhon
 * The method get_mp() returns now a reference (and not a pointer)
 * onto a mapping.
 *
 * Revision 1.9  2003/09/29 13:48:17  j_novak
 * New class Delta.
 *
 * Revision 1.8  2003/09/26 14:33:51  j_novak
 * Arithmetic functions for the class Tensor
 *
 * Revision 1.7  2003/09/26 08:05:29  j_novak
 * New class Vector.
 *
 * Revision 1.6  2003/09/25 21:01:50  e_gourgoulhon
 * Improved comments.
 *
 * Revision 1.5  2003/09/25 13:37:38  j_novak
 * Symmetric tensors of valence 2 are now implemented (not tested yet).
 *
 * Revision 1.4  2003/09/24 15:10:54  j_novak
 * Suppression of the etat flag in class Tensor (still present in Scalar)
 *
 * Revision 1.3  2003/09/24 08:46:31  j_novak
 * Added tensor.h and scalar.h to the documentation
 *
 * Revision 1.2  2003/09/23 08:53:11  e_gourgoulhon
 * not ready yet
 *
 * Revision 1.1  2003/09/22 12:50:47  e_gourgoulhon
 * First version: not ready yet!
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/tensor.h,v 1.63 2023/08/31 08:27:26 g_servignat Exp $
 *
 */

#define COV -1
#define CON +1

#define N_MET_MAX 5

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

namespace Lorene {
class Scalar ;
class Vector ; 
class Tensor_sym ;
class Sym_tensor ;
class Metric ;

            //-------------------------//
            //       class Tensor      //
            //-------------------------//
            

class Tensor { 

    // Data : 
    // -----
    protected:
    
    const Map* const mp ;   

    int valence ;   
            
    const Base_vect* triad ; 

    Itbl type_indice ;  
    
    int n_comp ;    

    Scalar** cmp ;   


    // Derived data : 
    // ------------
     protected:
    mutable const Metric* met_depend[N_MET_MAX] ; 

    mutable Tensor* p_derive_cov[N_MET_MAX];
    
    mutable Tensor* p_derive_con[N_MET_MAX];

    mutable Tensor* p_divergence[N_MET_MAX];

   
    // Constructors - Destructor :
    // -------------------------
    
    public: 

    Tensor(const Map& map, int val, const Itbl& tipe, 
            const Base_vect& triad_i) ;

    Tensor(const Map& map, int val, const Itbl& tipe, 
            const Base_vect* triad_i) ;

    Tensor(const Map& map, int val, int tipe, 
            const Base_vect& triad_i) ;

    Tensor(const Tensor&) ;  

    Tensor(const Map& map, const Base_vect& triad_i, FILE* fich) ;

    protected:
     explicit Tensor(const Map& map) ;

    Tensor(const Map& map, int val, const Itbl& tipe, int n_comp_i,
            const Base_vect& triad_i) ;

    Tensor(const Map& map, int val, int tipe, int n_comp_i, 
         const Base_vect& triad_i) ;


    public: 

    virtual ~Tensor() ; 
    
    // Memory management
    // -----------------
    protected:
    virtual void del_deriv() const ;    

    void set_der_0x0() const ; 

    virtual void del_derive_met(int) const ;

    void set_der_met_0x0(int) const ;

    void set_dependance (const Metric&) const ;

    int get_place_met(const Metric&) const ;
    
    // Mutators / assignment
    // ---------------------
    public:
    virtual void set_etat_nondef() ;
    
    virtual void set_etat_zero() ;

    virtual void set_etat_qcq() ;
    
    virtual void allocate_all() ; 

    virtual void change_triad(const Base_vect& new_triad) ; 
    
    void set_triad(const Base_vect& new_triad) ; 

    
    virtual void operator=(const Tensor&) ;
    
    Scalar& set(const Itbl& ind) ; 
    
    Scalar& set(int i1, int i2) ; 
    
    
    Scalar& set(int i1, int i2, int i3) ; 
    
    Scalar& set(int i1, int i2, int i3, int i4) ; 
    
    void annule_domain(int l) ; 

    virtual void annule(int l_min, int l_max) ; 

         void annule_extern_cn(int l_0, int deg) ;

    virtual void std_spectral_base() ; 
    
    virtual void std_spectral_base_odd() ; 
    
    virtual void dec_dzpuis(int dec = 1) ; 

    virtual void inc_dzpuis(int inc = 1) ; 
    
    virtual void  exponential_filter_r(int lzmin, int lzmax, int p, 
                double alpha= -16.) ;

    virtual void exponential_filter_ylm(int lzmin, int lzmax, int p,
                double alpha= -16.) ;

    virtual void exponential_filter_ylm_phi(int lzmin, int lzmax, int p_r, int p_tet, int p_phi,
                double alpha= -16.) ;

    // Computational methods
    // ---------------------
    
    protected: 
        void compute_derive_lie(const Vector& v, Tensor& resu) const ; 

    
    public:
    const Tensor& derive_cov(const Metric& gam) const ; 

    const Tensor& derive_con(const Metric& gam) const ; 

    const Tensor& divergence(const Metric& gam) const ; 


        Tensor derive_lie(const Vector& v) const ; 

    Tensor up(int ind, const Metric& gam) const ; 

    Tensor down(int ind, const Metric& gam) const ; 

    Tensor up_down(const Metric& gam) const ; 

    Tensor trace(int ind1, int ind2) const ; 

    Tensor trace(int ind1, int ind2, const Metric& gam) const ; 

    Scalar trace() const ; 

    Scalar trace(const Metric& gam) const ; 


    // Accessors
    // ---------
        public:
    virtual int position(const Itbl& ind) const ;

    virtual Itbl indices(int pos) const ;
    
    public:
    const Map& get_mp() const {return *mp ;} ; 

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

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

    Itbl get_index_type() const {return type_indice ; } ;

    int& set_index_type(int i) {return type_indice.set(i) ;};

    Itbl& set_index_type() {return type_indice ; } ;

    
    const Scalar& operator()(const Itbl& ind) const ; 

    const Scalar& operator()(int i1, int i2) const ; 

    const Scalar& operator()(int i1, int i2, int i3) const ; 
    
    const Scalar& operator()(int i1, int i2, int i3, int i4) const ; 
    
    // Member arithmetics
    // ------------------
    public:
    void operator+=(const Tensor &) ;           
    void operator-=(const Tensor &) ;           

    // Outputs
    // -------
    public:
    virtual void sauve(FILE *) const ;      

    virtual void spectral_display(const char* comment = 0x0, 
                            double threshold = 1.e-7, int precision = 4, 
                ostream& ostr = cout) const ;

    friend ostream& operator<<(ostream& , const Tensor & ) ;
    

    // Friend classes 
    // ---------------
    friend class Scalar ;
    friend class Vector ;
    friend class Sym_tensor ;
    friend class Tensor_sym ; 
    friend class Metric ;
  
    // Mathematical operators
    // ----------------------
    
    friend Scalar operator+(const Tensor&, const Scalar&) ;
    friend Scalar operator+(const Scalar&, const Tensor&) ;
    friend Scalar operator-(const Tensor&, const Scalar&) ;
    friend Scalar operator-(const Scalar&, const Tensor&) ;
    friend Tensor operator*(const Tensor&, const Tensor&) ; 
    friend Tensor_sym operator*(const Tensor&, const Tensor_sym&) ; 
    friend Tensor_sym operator*(const Tensor_sym&, const Tensor&) ;     
    friend Tensor_sym operator*(const Tensor_sym&, const Tensor_sym&) ;     
   
};



            //-------------------------//
            //    class Tensor_sym     //
            //-------------------------//

class Tensor_sym : public Tensor { 

    // Data : 
    // -----
    protected:
    
        int id_sym1 ;
        
        int id_sym2 ;
        
         
    // Constructors - Destructor :
    // -------------------------
    
    public: 

    Tensor_sym(const Map& map, int val, const Itbl& tipe, 
            const Base_vect& triad_i, int index_sym1, 
                        int index_sym2) ;

    Tensor_sym(const Map& map, int val, int tipe, const Base_vect& triad_i,
                    int index_sym1, int index_sym2) ;

    Tensor_sym(const Map& map, int tipe0, int tipe1, int tipe2, 
                   const Base_vect& triad_i,
                   int index_sym1, int index_sym2) ;

    Tensor_sym(const Tensor_sym& a) ;  

    Tensor_sym(const Map& map, const Base_vect& triad_i, FILE* fich) ;

    public: 

    virtual ~Tensor_sym() ; 
    
    // Mutators / assignment
    // ---------------------
    public:
    
    virtual void operator=(const Tensor_sym& a) ;
    
    virtual void operator=(const Tensor& a) ; 
    

    // Accessors
    // ---------
        public:
        int sym_index1() const {return id_sym1;} ;
        
        int sym_index2() const {return id_sym2;} ;
        
    virtual int position(const Itbl& ind) const ;

    virtual Itbl indices(int pos) const ;
        
            
    // Outputs
    // -------
    public:
    virtual void sauve(FILE *) const ;      
    

    // Tensor calculus
    // ---------------
    public:
    
    const Tensor_sym& derive_cov(const Metric& gam) const ; 

    const Tensor_sym& derive_con(const Metric& gam) const ; 

        Tensor_sym derive_lie(const Vector& v) const ; 


    // Mathematical operators
    // ----------------------
    
    friend Tensor_sym operator*(const Tensor&, const Tensor_sym&) ; 
    friend Tensor_sym operator*(const Tensor_sym&, const Tensor&) ; 
    
}; 



Tensor operator*(const Tensor& a, const Tensor& b) ; 

Tensor_sym operator*(const Tensor& a, const Tensor_sym& b) ; 

Tensor_sym operator*(const Tensor_sym& a, const Tensor& b) ; 

Tensor_sym operator*(const Tensor_sym& a, const Tensor_sym& b) ; 


Tensor contract(const Tensor& t1, int ind1, const Tensor& t2, int ind2, 
        bool desaliasing = false) ;

Tensor contract(const Tensor& t1, int ind_i1, int ind_j1, 
                const Tensor& t2, int ind_i2, int ind_j2,
                bool desaliasing = false) ;


Tensor contract(const Tensor& t1, int ind1, int ind2) ;


Tbl max(const Tensor& aa, const char* comment = 0x0, ostream& ost = cout) ; 


Tbl min(const Tensor& aa, const char* comment = 0x0, ostream& ost = cout) ; 

Tbl maxabs(const Tensor& aa, const char* comment = 0x0, ostream& ost = cout, 
       bool verb = true) ; 


Tbl diffrel(const Tensor& aa, const Tensor& bb, const char* comment = 0x0,
            ostream& ost = cout) ; 

Tbl diffrelmax(const Tensor& aa, const Tensor& bb, const char* comment = 0x0,
               ostream& ost = cout) ; 

Tbl central_value(const Tensor& aa, const char* comment = 0x0, ostream& ost = cout) ; 

Tbl max_all_domains(const Tensor& aa, int l_excluded = -1, const char* comment = 0x0, 
    ostream& ost = cout) ; 


Tbl min_all_domains(const Tensor& aa, int l_excluded = -1, const char* comment = 0x0, 
    ostream& ost = cout) ; 

Tbl maxabs_all_domains(const Tensor& aa, int l_excluded = -1, const char* comment = 0x0, 
    ostream& ost = cout, bool verb = true) ; 



Tensor operator+(const Tensor& ) ;          
Tensor operator-(const Tensor& ) ;          
Tensor operator+(const Tensor& a, const Tensor& b) ;    

Scalar operator+(const Tensor& a, const Scalar& b) ;    

Scalar operator+(const Scalar& a, const Tensor& b) ;    

Tensor operator-(const Tensor& a, const Tensor& b) ;    

Scalar operator-(const Tensor& a, const Scalar& b) ;    

Scalar operator-(const Scalar& a, const Tensor& b) ;    

Tensor operator*(const Scalar& a , const Tensor& b) ;   
Tensor operator*(const Tensor& a, const Scalar& b) ;    
Tensor operator*(double , const Tensor&) ;              
Tensor operator* (const Tensor&, double) ;              
Tensor operator*(int, const Tensor &) ;                 
Tensor operator*(const Tensor&, int) ;                 
Tensor operator/(const Tensor&, const Scalar&) ;       
Tensor operator/(const Tensor&, double) ;              
Tensor operator/(const Tensor&, int) ;                 


Tensor_sym operator+(const Tensor_sym&) ;  

Tensor_sym operator-(const Tensor_sym&) ;  

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

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

Tensor_sym operator*(const Scalar& a, const Tensor_sym& b) ;   

Tensor_sym operator*(const Tensor_sym& a, const Scalar& b) ;  

Tensor_sym operator*(double, const Tensor_sym&) ;  

Tensor_sym operator*(const Tensor_sym&, double) ;  

Tensor_sym operator*(int, const Tensor_sym&) ;  

Tensor_sym operator*(const Tensor_sym&, int) ;  

Tensor_sym operator/(const Tensor_sym&, const Scalar&) ;  

Tensor_sym operator/(const Tensor_sym&, double) ;  

Tensor_sym operator/(const Tensor_sym&, int) ;  

}
#include "scalar.h"

#include "vector.h"

#include "sym_tensor.h"


#endif