vector.h

 /*
 *  Definition of Lorene class Vector
 *
 */

/*
 *   Copyright (c) 2003  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 __VECTOR_H_ 
#define __VECTOR_H_ 

/*
 * $Id: vector.h,v 1.42 2014/10/13 08:52:37 j_novak Exp $
 * $Log: vector.h,v $
 * Revision 1.42  2014/10/13 08:52:37  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.41  2008/12/03 10:18:56  j_novak
 * Method 6 is now the default for calls to vector Poisson solver.
 *
 * Revision 1.40  2008/10/29 08:19:08  jl_cornou
 * Typo in the doxygen documentation + spectral bases for pseudo vectors added
 * and curl
 *
 * Revision 1.39  2008/08/27 08:45:21  jl_cornou
 * Implemented routines to solve dirac systems for divergence-free vectors
 *
 * Revision 1.38  2007/12/21 16:06:16  j_novak
 * Methods to filter Tensor, Vector and Sym_tensor objects.
 *
 * Revision 1.37  2007/05/11 11:38:16  n_vasset
 * Added poisson_boundary2.C routine
 *
 * Revision 1.36  2007/01/16 15:05:59  n_vasset
 * New constructor (taking a Scalar in mono-domain angular grid for
 * boundary) for function sol_elliptic_boundary
 *
 * Revision 1.35  2005/06/09 07:56:25  f_limousin
 * Implement a new function sol_elliptic_boundary() and
 * Vector::poisson_boundary(...) which solve the vectorial poisson
 * equation (method 6) with an inner boundary condition.
 *
 * Revision 1.34  2005/02/16 15:00:05  m_forot
 * Add visu_streamlime function
 *
 * Revision 1.33  2005/02/14 13:01:48  j_novak
 * p_eta and p_mu are members of the class Vector. Most of associated functions
 * have been moved from the class Vector_divfree to the class Vector.
 *
 * Revision 1.32  2004/05/25 14:54:01  f_limousin
 * Change arguments of method poisson with parameters.
 *
 * Revision 1.31  2004/05/09 20:54:22  e_gourgoulhon
 * Added method flux (to compute the flux accross a sphere).
 *
 * Revision 1.30  2004/03/29 11:57:13  e_gourgoulhon
 * Added methods ope_killing and ope_killing_conf.
 *
 * Revision 1.29  2004/03/28 21:25:14  e_gourgoulhon
 * Minor modif comments (formula for V^\theta in Vector_divfree).
 *
 * Revision 1.28  2004/03/27 20:59:55  e_gourgoulhon
 * Slight modif comment (doxygen \ingroup).
 *
 * Revision 1.27  2004/03/22 13:12:44  j_novak
 * Modification of comments to use doxygen instead of doc++
 *
 * Revision 1.26  2004/03/10 12:52:19  f_limousin
 * Add a new argument "method" in poisson method.
 *
 * Revision 1.25  2004/03/03 09:07:02  j_novak
 * In Vector::poisson(double, int), the flat metric is taken from the mapping.
 *
 * Revision 1.24  2004/02/26 22:45:44  e_gourgoulhon
 * Added method derive_lie.
 *
 * Revision 1.23  2004/02/22 15:47:45  j_novak
 * Added 2 more methods to solve the vector Poisson equation. Method 1 is not
 * tested yet.
 *
 * Revision 1.22  2004/02/21 16:27:53  j_novak
 * Modif. comments
 *
 * Revision 1.21  2004/02/16 17:40:14  j_novak
 * Added a version of poisson with the flat metric as argument (avoids
 * unnecessary calculations by decompose_div)
 *
 * Revision 1.20  2003/12/14 21:47:24  e_gourgoulhon
 * Added method visu_arrows for visualization through OpenDX.
 *
 * Revision 1.19  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.18  2003/11/03 22:31:02  e_gourgoulhon
 * Class Vector_divfree: parameters of methods set_vr_eta_mu and set_vr_mu
 * are now all references.
 *
 * Revision 1.17  2003/11/03 14:02:17  e_gourgoulhon
 * Class Vector_divfree: the members p_eta and p_mu are no longer declared
 * "const".
 *
 * Revision 1.16  2003/10/20 15:12:17  j_novak
 * New method Vector::poisson
 *
 * Revision 1.15  2003/10/20 14:44:49  e_gourgoulhon
 * Class Vector_divfree: added method poisson().
 *
 * Revision 1.14  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.13  2003/10/17 16:36:53  e_gourgoulhon
 * Class Vector_divfree: Added new methods set_vr_eta_mu and set_vr_mu.
 *
 * Revision 1.12  2003/10/16 21:36:01  e_gourgoulhon
 * Added method Vector_divfree::update_vtvp().
 *
 * Revision 1.11  2003/10/16 15:25:00  e_gourgoulhon
 * Changes in documentation.
 *
 * Revision 1.10  2003/10/16 14:21:33  j_novak
 * The calculation of the divergence of a Tensor is now possible.
 *
 * Revision 1.9  2003/10/13 20:49:26  e_gourgoulhon
 * Corrected typo in the comments.
 *
 * Revision 1.8  2003/10/13 13:52:39  j_novak
 * Better managment of derived quantities.
 *
 * Revision 1.7  2003/10/12 20:33:15  e_gourgoulhon
 * Added new derived class Vector_divfree (preliminary).
 *
 * Revision 1.6  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.5  2003/10/05 21:07:27  e_gourgoulhon
 * Method std_spectral_base() is now virtual.
 *
 * Revision 1.4  2003/10/03 14:08:03  e_gourgoulhon
 * Added constructor from Tensor.
 *
 * Revision 1.3  2003/09/29 13:48:17  j_novak
 * New class Delta.
 *
 * Revision 1.2  2003/09/29 12:52:56  j_novak
 * Methods for changing the triad are implemented.
 *
 * Revision 1.1  2003/09/26 08:07:32  j_novak
 * New class vector
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/vector.h,v 1.42 2014/10/13 08:52:37 j_novak Exp $
 *
 */

namespace Lorene {
class Vector_divfree ;

            //-------------------------//
            //       class Vector      //
            //-------------------------//
            

class Vector: public Tensor {

    // Derived data : 
    // ------------
    protected:
        mutable Scalar* p_potential[N_MET_MAX] ;

    mutable Vector_divfree* p_div_free[N_MET_MAX] ;

    mutable Scalar* p_eta ;
    
    mutable Scalar* p_mu ;

    mutable Scalar* p_A ;
    
   // Constructors - Destructor
    // -------------------------
    public:
    Vector(const Map& map, int tipe, const Base_vect& triad_i) ;

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

    Vector(const Vector& a) ;       

    Vector(const Tensor& a) ;   

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

    virtual ~Vector() ;         

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


    // Mutators / assignment
    // ---------------------
    public:
    virtual void change_triad(const Base_vect& ) ; 

    virtual void operator=(const Vector& a) ;   
    
    virtual void operator=(const Tensor& a) ;   

    void set_vr_eta_mu(const Scalar& vr_i, const Scalar& eta_i,
        const Scalar& mu_i) ; 

    void decompose_div(const Metric&) const ;

    const Scalar& potential(const Metric& ) const ;
    
    const Vector_divfree& div_free(const Metric& ) const;

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

    
    // Accessors
    // ---------
    public:
    Scalar& set(int ) ; 

    const Scalar& operator()(int ) const; 

    virtual int position(const Itbl& idx) const {
      assert (idx.get_ndim() == 1) ;
      assert (idx.get_dim(0) == 1) ;
      assert ((idx(0) >= 1) && (idx(0) <= 3)) ;

      return (idx(0) - 1) ;
    
    } ;

    virtual Itbl indices(int place) const {
      assert((place>=0) && (place<3)) ;
      
      Itbl res(1) ;
      res = place + 1;
      return res ;

    };
    
    virtual void std_spectral_base() ; 

    virtual void pseudo_spectral_base() ; 

    virtual const Scalar& eta() const ;
    
    virtual const Scalar& mu() const ;


    virtual const Scalar& A() const ;

    
    void update_vtvp() ;

    
    // Differential operators/ PDE solvers
    // -----------------------------------
    public:
    const Scalar& divergence(const Metric&) const ; 

    const Vector_divfree curl() const ;

    Vector derive_lie(const Vector& v) const ; 
        
     Sym_tensor ope_killing(const Metric& gam) const ; 
     
     Sym_tensor ope_killing_conf(const Metric& gam) const ; 

    Vector poisson(double lambda, int method = 6) const ;
     
    Vector poisson(double lambda, const Metric_flat& met_f, int method = 6) const ;
    
    Vector poisson(const double lambda, Param& par,
           int method = 6) const ;
    
    void poisson_boundary(double lambda, const Mtbl_cf& limit_vr, 
              const Mtbl_cf& limit_eta, const Mtbl_cf& limit_mu, 
              int num_front, double fact_dir, double fact_neu,
              Vector& resu) const ;
 

    void poisson_boundary2(double lam, Vector& resu, Scalar boundvr, 
               Scalar boundeta, Scalar boundmu, double dir_vr, 
               double neum_vr, double dir_eta, double neum_eta, 
               double dir_mu, double neum_mu ) const ;

 
    Vector poisson_dirichlet(double lambda, const Valeur& limit_vr, 
              const Valeur& limit_vt, const Valeur& limit_vp, 
              int num_front) const ;

    Vector poisson_neumann(double lambda, const Valeur& limit_vr, 
              const Valeur& limit_vt, const Valeur& limit_vp, 
              int num_front) const ;

    Vector poisson_robin(double lambda, const Valeur& limit_vr, 
             const Valeur& limit_vt, const Valeur& limit_vp, 
             double fact_dir, double fact_neu, 
             int num_front) const ;
 

    double flux(double radius, const Metric& met) const ; 

    void poisson_block(double lambda, Vector& resu) const ;

        // Graphics
        // --------
 public:
    void visu_arrows(double xmin, double xmax, double ymin, double ymax,
    double zmin, double zmax, const char* title0 = 0x0, 
    const char* filename0 = 0x0, bool start_dx = true, int nx = 8, int ny = 8, 
    int nz = 8) const ;    

    void visu_streamline(double xmin, double xmax, double ymin, double ymax,
    double zmin, double zmax, const char* title0 = 0x0, 
    const char* filename0 = 0x0, bool start_dx = true, int nx = 8, int ny = 8, 
    int nz = 8) const ;   
        
     
 
};



            //---------------------------------//
            //       class Vector_divfree      //
            //---------------------------------//
            

class Vector_divfree: public Vector {

    // Data : 
    // -----
    protected:
    const Metric* const met_div ; 
    
    // Constructors - Destructor
    // -------------------------
    public:
    Vector_divfree(const Map& map, const Base_vect& triad_i, 
        const Metric& met) ;

    Vector_divfree(const Vector_divfree& ) ;       

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

    virtual ~Vector_divfree() ;         

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

    void set_der_0x0() const ; 


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

    public:
    void operator=(const Vector_divfree& a) ;   
    
    virtual void operator=(const Vector& a) ;   
    
    virtual void operator=(const Tensor& a) ;   
    
    void set_vr_mu(const Scalar& vr_i, const Scalar& mu_i) ; 
    
    void set_vr_eta_mu(const Scalar& vr_i, const Scalar& eta_i, const Scalar& mu_i) ;

    void set_A_mu(const Scalar& A_i, const Scalar& mu_i, const Param* par_bc) ; 
    
    // Computational methods
    // ---------------------
    public:
    virtual const Scalar& eta() const ;
    
    Vector_divfree poisson() const ; 

    void update_etavr() ;
    
    Vector_divfree poisson(Param& par) const ; 
    protected:
    
    void sol_Dirac_A(const Scalar& aaa, Scalar& eta, Scalar& vr,
              const Param* par_bc = 0x0) const ;    

    void sol_Dirac_A_tau(const Scalar& aaa, Scalar& eta, Scalar& vr,
              const Param* par_bc = 0x0) const ;

    void sol_Dirac_A_poisson(const Scalar& aaa, Scalar& eta, Scalar& vr,
              const Param* par_bc = 0x0) const ;

    void sol_Dirac_A_1z(const Scalar& aaa, Scalar& eta, Scalar& vr,
              const Param* par_bc = 0x0) const ;    
    
};





}
#endif