sym_tensor.h

/*
 *  Definition of Lorene class Sym_tensor, 
 *  as well as derived classes Sym_tensor_trans and Sym_tensor_tt
 *
 */

/*
 *   Copyright (c) 2003-2004 Eric Gourgoulhon & 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 version 2
 *   as published by the Free Software Foundation.
 *
 *   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 __SYM_TENSOR_H_ 
#define __SYM_TENSOR_H_ 

/*
 * $Id: sym_tensor.h,v 1.50 2019/08/16 08:47:36 j_novak Exp $
 * $Log: sym_tensor.h,v $
 * Revision 1.50  2019/08/16 08:47:36  j_novak
 * *** empty log message ***
 *
 * Revision 1.49  2014/10/13 08:52:36  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.48  2010/10/11 10:23:03  j_novak
 * Removed methods Sym_tensor_trans::solve_hrr() and Sym_tensor_trans::set_WX_det_one(), as they are no longer relevant.
 *
 * Revision 1.47  2008/12/05 08:46:19  j_novak
 * New method Sym_tensor_trans_aux::set_tt_part_det_one.
 *
 * Revision 1.46  2008/12/03 10:18:56  j_novak
 * Method 6 is now the default for calls to vector Poisson solver.
 *
 * Revision 1.45  2008/08/20 14:39:53  n_vasset
 * New Dirac solvers handling degenerate elliptic operators on excised spacetimes.
 *
 * Revision 1.44  2007/12/21 16:06:16  j_novak
 * Methods to filter Tensor, Vector and Sym_tensor objects.
 *
 * Revision 1.43  2007/11/27 15:48:52  n_vasset
 * New member p_tilde_c for class Sym_tensor
 *
 * Revision 1.42  2007/05/04 16:43:50  n_vasset
 * adding of functions sol_Dirac_BC2 and sol_Dirac_A2
 *
 * Revision 1.41  2006/10/24 13:03:17  j_novak
 * New methods for the solution of the tensor wave equation. Perhaps, first
 * operational version...
 *
 * Revision 1.40  2006/08/31 12:13:21  j_novak
 * Added an argument of type Param to Sym_tensor_trans::sol_  rac_A().
 *
 * Revision 1.39  2006/06/20 12:07:13  j_novak
 * Improved execution speed for sol_Dirac_tildeB...
 *
 * Revision 1.38  2006/06/14 10:04:19  j_novak
 * New methods sol_Dirac_l01, set_AtB_det_one and set_AtB_trace_zero.
 *
 * Revision 1.37  2006/06/13 13:30:12  j_novak
 * New members sol_Dirac_A and sol_Dirac_tildeB (see documentation).
 *
 * Revision 1.36  2006/06/12 13:37:23  j_novak
 * Added bounds in l (multipolar momentum) for Sym_tensor_trans::solve_hrr.
 *
 * Revision 1.35  2006/06/12 07:42:28  j_novak
 * Fields A and tilde{B} are defined only for l>1.
 *
 * Revision 1.34  2006/06/12 07:27:18  j_novak
 * New members concerning A and tilde{B}, dealing with the transverse part of the
 * Sym_tensor.
 *
 * Revision 1.33  2005/11/28 14:45:14  j_novak
 * Improved solution of the Poisson tensor equation in the case of a transverse
 * tensor.
 *
 * Revision 1.32  2005/09/16 13:58:10  j_novak
 * New Poisson solver for a Sym_tensor_trans.
 *
 * Revision 1.31  2005/09/07 16:47:42  j_novak
 * Removed method Sym_tensor_trans::T_from_det_one
 * Modified Sym_tensor::set_auxiliary, so that it takes eta/r and mu/r as
 * arguments.
 * Modified Sym_tensor_trans::set_hrr_mu.
 * Added new protected method Sym_tensor_trans::solve_hrr
 *
 * Revision 1.30  2005/04/08 08:22:04  j_novak
 * New methods set_hrr_mu_det_one() and set_WX_det_one(). Not tested yet...
 *
 * Revision 1.29  2005/04/06 15:43:58  j_novak
 * New method Sym_tensor_trans::T_from_det_one(...).
 *
 * Revision 1.28  2005/04/04 15:25:22  j_novak
 * Added new members www, xxx, ttt and the associated methods.
 *
 * Revision 1.27  2005/04/01 14:28:31  j_novak
 * Members p_eta and p_mu are now defined in class Sym_tensor.
 *
 * Revision 1.26  2005/01/03 08:34:58  f_limousin
 * Come back to the previous version.
 *
 * Revision 1.25  2005/01/03 08:15:39  f_limousin
 * The first argument of the function trace_from_det_one(...) is now
 * a Sym_tensor_trans instead of a Sym_tensor_tt (because of a
 * compilation error with some compilators).
 *
 * Revision 1.24  2004/12/28 14:21:46  j_novak
 * Added the method Sym_tensor_trans::trace_from_det_one
 *
 * Revision 1.23  2004/12/28 10:37:22  j_novak
 * Better way of enforcing zero divergence.
 *
 * Revision 1.22  2004/06/14 20:44:44  e_gourgoulhon
 * Added argument method_poisson to Sym_tensor::longit_pot and
 * Sym_tensor::transverse.
 *
 * Revision 1.21  2004/05/25 14:57:20  f_limousin
 * Add parameters in argument of functions transverse, longit_pot,
 * set_tt_trace, tt_part and set_khi_mu for the case of a Map_et.
 *
 * Revision 1.20  2004/05/24 13:44:54  e_gourgoulhon
 * Added parameter dzp to method Sym_tensor_tt::update.
 *
 * Revision 1.19  2004/04/08 16:37:54  e_gourgoulhon
 * Sym_tensor_tt::set_khi_mu: added argument dzp (dzpuis of resulting h^{ij}).
 *
 * Revision 1.18  2004/03/30 14:01:19  j_novak
 * Copy constructors and operator= now copy the "derived" members.
 *
 * Revision 1.17  2004/03/29 16:13:06  j_novak
 * New methods set_longit_trans and set_tt_trace .
 *
 * Revision 1.16  2004/03/22 13:12:43  j_novak
 * Modification of comments to use doxygen instead of doc++
 *
 * Revision 1.15  2004/03/03 13:54:16  j_novak
 * Error in comments corrected.
 *
 * Revision 1.14  2004/03/03 13:16:20  j_novak
 * New potential khi (p_khi) and the functions manipulating it.
 *
 * Revision 1.13  2004/02/26 22:45:13  e_gourgoulhon
 * Added method derive_lie.
 *
 * Revision 1.12  2004/02/18 18:43:22  e_gourgoulhon
 * Method trace() renamed the_trace() in order to avoid
 * any confusion with new method Tensor::trace().
 *
 * Revision 1.11  2004/01/04 20:49:06  e_gourgoulhon
 * Sym_tensor is now a derived class of Tensor_sym.
 * Suppressed methods Sym_tensor::indices and Sym_tensor::position:
 *  they are now implemented at the Tensor_sym level.
 *
 * Revision 1.10  2003/11/27 16:05:11  e_gourgoulhon
 * Changed return value of methods transverse( ) and longit_pot( ).
 *
 * Revision 1.9  2003/11/26 21:56:21  e_gourgoulhon
 * Class Sym_tensor: added the members p_transverse and p_longit_pot,
 * and the associated methods transverse( ), longit_pot( ),
 * del_deriv_met( ) and set_der_met_0x0( ).
 *
 * Revision 1.8  2003/11/07 16:54:23  e_gourgoulhon
 * Added method Sym_tensor_tt::poisson().
 *
 * Revision 1.7  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.6  2003/11/05 15:26:31  e_gourgoulhon
 * Modif documentation.
 *
 * Revision 1.5  2003/11/04 22:57:26  e_gourgoulhon
 * Class Sym_tensor_tt: method set_eta_mu renamed set_rr_eta_mu
 *    method update_tp() renamed update()
 *    added method set_rr_mu.
 *
 * Revision 1.4  2003/11/03 22:29:54  e_gourgoulhon
 * Class Sym_tensor_tt: added functions set_eta_mu and update_tp.
 *
 * Revision 1.3  2003/11/03 17:09:30  e_gourgoulhon
 * Class Sym_tensor_tt: added the methods eta() and mu().
 *
 * Revision 1.2  2003/10/28 21:22:51  e_gourgoulhon
 * Class Sym_tensor_trans: added methods trace() and tt_part().
 *
 * Revision 1.1  2003/10/27 10:45:19  e_gourgoulhon
 * New derived classes Sym_tensor_trans and Sym_tensor_tt.
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/sym_tensor.h,v 1.50 2019/08/16 08:47:36 j_novak Exp $
 *
 */

namespace Lorene {
class Sym_tensor_trans ;
class Sym_tensor_tt ;




            //---------------------------------//
            //        class Sym_tensor         //
            //---------------------------------//
            
class Sym_tensor : public Tensor_sym {

    // Derived data : 
    // ------------
    protected:
    mutable Sym_tensor_trans* p_transverse[N_MET_MAX] ;

    mutable Vector* p_longit_pot[N_MET_MAX] ;

    mutable Scalar* p_eta ;
    
    mutable Scalar* p_mu ;

    mutable Scalar* p_www ;

    mutable Scalar* p_xxx ;

    mutable Scalar* p_ttt ;

    mutable Scalar* p_aaa ;

    mutable Scalar* p_tilde_b ;

    mutable Scalar* p_tilde_c ;


     



    // Constructors - Destructor :
    // -------------------------
    
    public:
    Sym_tensor(const Map& map, const Itbl& tipe, const Base_vect& triad_i) ;

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

    Sym_tensor(const Sym_tensor& a) ; 

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

    virtual ~Sym_tensor() ;    

      

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

    void set_der_0x0() const ; 

    virtual void del_derive_met(int i) const ;

    void set_der_met_0x0(int i) const ;


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

    virtual void operator=(const Tensor_sym& a) ;

    virtual void operator=(const Tensor& a) ;

    void set_longit_trans( const Vector& v, const Sym_tensor_trans& a) ;

    void set_auxiliary( const Scalar& trr, const Scalar& eta_over_r, const
                Scalar& mu_over_r, const Scalar& www, const Scalar&
                xxx, const Scalar& ttt ) ;

    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.) ;

    // Computation of derived members
    // ------------------------------
    public:


    const Vector& divergence(const Metric&) const ; 

        Sym_tensor derive_lie(const Vector& v) const ; 

    const Sym_tensor_trans& transverse(const Metric& gam, Param* par = 0x0,
                int method_poisson = 6) const ; 

    const Vector& longit_pot(const Metric& gam, Param* par = 0x0,
                int method_poisson = 6) const ; 
    
    virtual const Scalar& eta(Param* par = 0x0) const ;

    const Scalar& mu(Param* par = 0x0) const ;

    const Scalar& www() const ;

    const Scalar& xxx() const ;

    const Scalar& ttt() const ;

    const Scalar& compute_A(bool output_ylm = true, Param* par = 0x0) const ;

    const Scalar& compute_tilde_B(bool output_ylm = true, Param* par = 0x0) const ;

    Scalar compute_tilde_B_tt(bool output_ylm = true, Param* par = 0x0) const ;

    const Scalar& compute_tilde_C(bool output_ylm = true, Param* par = 0x0) const ;





 protected:
    Scalar get_tilde_B_from_TT_trace(const Scalar& tilde_B_tt_in, const Scalar&
        trace) const ;
    
    // Mathematical operators
    // ----------------------
 protected:
    Sym_tensor* inverse() const ;

    // Friend classes
    //-----------------
    friend class Metric ;
 
} ;


            //---------------------------------//
            //    class Sym_tensor_trans       //
            //---------------------------------//
            

class Sym_tensor_trans: public Sym_tensor {

    // Data : 
    // -----
    protected:
    const Metric* const met_div ; 
    
    mutable Scalar* p_trace ; 
    
    mutable Sym_tensor_tt* p_tt ;
    
    // Constructors - Destructor
    // -------------------------
    public:
    Sym_tensor_trans(const Map& map, const Base_vect& triad_i, 
        const Metric& met) ;

    Sym_tensor_trans(const Sym_tensor_trans& ) ;       

    Sym_tensor_trans(const Map& map, const Base_vect& triad_i, 
        const Metric& met, FILE* fich) ;

    virtual ~Sym_tensor_trans() ;           

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

    void set_der_0x0() const ; 


    // Accessors
    // ---------
        public:
    const Metric& get_met_div() const {return *met_div ; } ;

    // Mutators / assignment
    // ---------------------

    public:
    virtual void operator=(const Sym_tensor_trans& a) ; 
    
    virtual void operator=(const Sym_tensor& a) ;   
    
    virtual void operator=(const Tensor_sym& a) ;

    virtual void operator=(const Tensor& a) ;   
    
    void set_tt_trace(const Sym_tensor_tt& a, const Scalar& h, 
              Param* par = 0x0) ;

    // Computational methods
    // ---------------------
    const Scalar& the_trace() const ; 
    
    const Sym_tensor_tt& tt_part(Param* par = 0x0) const ; 

 protected:
    void sol_Dirac_A(const Scalar& aaa, Scalar& tilde_mu, Scalar& xxx,
             const Param* par_bc = 0x0) const ;


    void sol_Dirac_tilde_B(const Scalar& tilde_b, const Scalar& hh, Scalar& hrr,
                   Scalar& tilde_eta, Scalar& www, Param* par_bc=0x0,
                   Param* par_mat=0x0) const ;

    void sol_Dirac_l01(const Scalar& hh, Scalar& hrr, Scalar& tilde_eta,
               Param* par_mat) const ;



 public:


        
    void sol_Dirac_Abound(const Scalar& aaa, Scalar& tilde_mu, Scalar& x_new,
              Scalar bound_mu, const Param* par_bc);
 
        
    void sol_Dirac_A2(const Scalar& aaa, Scalar& tilde_mu, Scalar& x_new,
              Scalar bound_mu, const Param* par_bc);       

    void sol_Dirac_BC2(const Scalar& bb, const Scalar& cc, const Scalar& hh, 
                Scalar& hrr, Scalar& tilde_eta, Scalar& ww, Scalar bound_eta,double dir, double neum, double rhor, Param* par_bc, Param* par_mat); 
 

    void sol_Dirac_BC3(const Scalar& bb, const Scalar& hh, 
                Scalar& hrr, Scalar& tilde_eta, Scalar& ww, Scalar bound_hrr, Scalar bound_eta, Param* par_bc, Param* par_mat); 



     // Solving the electric system for l=0 and l=1 only (simpler case), with boundary conditions imposed by the degenerate elliptic system.
   
    void sol_Dirac_l01_bound(const Scalar& hh, Scalar& hrr, Scalar& tilde_eta, Scalar& bound_hrr, Scalar& bound_eta, Param* par_mat) ;

    // Provisory: just for compilation, to be removed
    void sol_Dirac_l01_2(const Scalar& hh, Scalar& hrr, Scalar& tilde_eta, Param* par_mat) ;


    void trace_from_det_one(const Sym_tensor_tt& htt, 
                double precis = 1.e-14, int it_max = 100) ;

    void set_hrr_mu_det_one(const Scalar& hrr, const Scalar& mu_in,
                double precis = 1.e-14, int it_max = 100) ;

    void set_tt_part_det_one(const Sym_tensor_tt& hijtt, const 
                 Scalar* h_prev = 0x0, Param* par_mat = 0x0, 
                 double precis = 1.e-14, int it_max = 100) ;

    void set_AtBtt_det_one(const Scalar& a_in, const Scalar& tbtt_in, 
                   const Scalar* h_prev = 0x0, Param* par_bc = 0x0,
                   Param* par_mat = 0x0, double precis = 1.e-14, 
                   int it_max = 100) ;

    void set_AtB_trace(const Scalar& a_in, const Scalar& tb_in, const 
               Scalar& trace, Param* par_bc = 0x0, Param* par_mat = 0x0) ;

    Sym_tensor_trans poisson(const Scalar* h_guess = 0x0) const ; 
} ; 
    

            //------------------------------//
            //    class Sym_tensor_tt       //
            //------------------------------//
            

class Sym_tensor_tt: public Sym_tensor_trans {

    // Data : 
    // -----

    protected:
    mutable Scalar* p_khi ;
    
    
    // Constructors - Destructor
    // -------------------------
    public:
    Sym_tensor_tt(const Map& map, const Base_vect& triad_i, 
        const Metric& met) ;

    Sym_tensor_tt(const Sym_tensor_tt& ) ;       

    Sym_tensor_tt(const Map& map, const Base_vect& triad_i, 
        const Metric& met, FILE* fich) ;

    virtual ~Sym_tensor_tt() ;          

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

    void set_der_0x0() const ; 


    // Mutators / assignment
    // ---------------------

    public:
    virtual void operator=(const Sym_tensor_tt& a) ;    
    
    virtual void operator=(const Sym_tensor_trans& a) ; 
    
    virtual void operator=(const Sym_tensor& a) ;   
    
    virtual void operator=(const Tensor_sym& a) ;

    virtual void operator=(const Tensor& a) ;   
    
    void set_rr_eta_mu(const Scalar& hrr, const Scalar& eta_i, 
                        const Scalar& mu_i) ; 

    void set_rr_mu(const Scalar& hrr, const Scalar& mu_i) ; 
    
    
    void set_khi_eta_mu(const Scalar& khi_i, const Scalar& eta_i, 
                        const Scalar& mu_i) ; 
        
    void set_khi_mu(const Scalar& khi_i, const Scalar& mu_i, int dzp = 0,
            Param* par1 = 0x0, Param* par2 = 0x0, 
            Param* par3 = 0x0) ; 

    void set_A_tildeB(const Scalar& a_in, const Scalar& tb_in, Param* par_bc = 0x0,
              Param* par_mat = 0x0) ;

    // Computational methods
    // ---------------------
    
    public:
    const Scalar& khi() const ;
    
    virtual const Scalar& eta(Param* par = 0x0) const ;

    protected:
    void update(int dzp, Param* par1 = 0x0, Param* par2 = 0x0) ;

    public:
    Sym_tensor_tt poisson(int dzfin = 2) const ; 
    


} ; 
    



}
#endif