scalar.h

/*
 *  Definition of Lorene class Scalar
 *
 */

/*
 *   Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
 *   
 *   Copyright (c) 1999-2000 Jean-Alain Marck (for previous class Cmp)
 *   Copyright (c) 1999-2002 Eric Gourgoulhon (for previous class Cmp)
 *   Copyright (c) 1999-2001 Philippe Grandclement (for previous class Cmp)
 *   Copyright (c) 2000-2002 Jerome Novak (for previous class Cmp)
 *   Copyright (c) 2000-2001 Keisuke Taniguchi (for previous class Cmp)
 *
 *   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 __SCALAR_H_ 
#define __SCALAR_H_ 


/*
 * $Id: scalar.h,v 1.96 2023/08/31 08:27:26 g_servignat Exp $
 * $Log: scalar.h,v $
 * Revision 1.96  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.95  2023/08/28 09:53:33  g_servignat
 * Added ylm filter for Tensor and Scalar in theta + phi directions
 *
 * Revision 1.94  2015/12/18 15:52:51  j_novak
 * New method is_nan() for class Scalar.
 *
 * Revision 1.93  2015/09/10 13:28:05  j_novak
 * New methods for the class Hot_Eos
 *
 * Revision 1.92  2014/10/13 08:52:36  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.91  2013/06/05 15:06:10  j_novak
 * Legendre bases are treated as standard bases, when the multi-grid
 * (Mg3d) is built with BASE_LEG.
 *
 * Revision 1.90  2013/01/11 15:44:54  j_novak
 * Addition of Legendre bases (part 2).
 *
 * Revision 1.89  2012/12/19 13:59:56  j_penner
 * added a few lines to the documentation for scalar_match_tau function
 *
 * Revision 1.88  2012/01/17 15:05:46  j_penner
 * *** empty log message ***
 *
 * Revision 1.87  2012/01/17 10:16:27  j_penner
 * functions added: sarra_filter_r, sarra_filter_r_all_domains, Heaviside
 *
 * Revision 1.86  2011/04/08 13:13:09  e_gourgoulhon
 * Changed the comment of function val_point to indicate specifically the
 * division by r^dzpuis in the compactified external domain.
 *
 * Revision 1.85  2008/09/29 13:23:51  j_novak
 * Implementation of the angular mapping associated with an affine
 * mapping. Things must be improved to take into account the domain index.
 *
 * Revision 1.84  2008/09/22 19:08:01  j_novak
 * New methods to deal with boundary conditions
 *
 * Revision 1.83  2008/05/24 15:05:22  j_novak
 * New method Scalar::match_tau to match the output of an explicit time-marching scheme with the tau method.
 *
 * Revision 1.82  2007/12/21 16:06:16  j_novak
 * Methods to filter Tensor, Vector and Sym_tensor objects.
 *
 * Revision 1.81  2007/10/31 10:33:11  j_novak
 * Added exponential filters to smooth Gibbs-type phenomena.
 *
 * Revision 1.80  2007/06/21 19:56:36  k_taniguchi
 * Introduction of another filtre_r.
 *
 * Revision 1.79  2007/05/06 10:48:08  p_grandclement
 * Modification of a few operators for the vorton project
 *
 * Revision 1.78  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.77  2006/05/26 09:00:09  j_novak
 * New members for multiplication or division by cos(theta).
 *
 * Revision 1.76  2005/11/30 13:48:06  e_gourgoulhon
 * Replaced M_PI/2 by 1.57... in argument list of sol_elliptic_sin_zec
 * (in order not to require the definition of M_PI).
 *
 * Revision 1.75  2005/11/30 11:09:03  p_grandclement
 * Changes for the Bin_ns_bh project
 *
 * Revision 1.74  2005/11/17 15:29:46  e_gourgoulhon
 * Added arithmetics with Mtbl.
 *
 * Revision 1.73  2005/10/25 08:56:34  p_grandclement
 * addition of std_spectral_base in the case of odd functions near the origin
 *
 * Revision 1.72  2005/09/07 13:10:47  j_novak
 * Added a filter setting to zero all mulitpoles in a given range.
 *
 * Revision 1.71  2005/08/26 14:02:38  p_grandclement
 * Modification of the elliptic solver that matches with an oscillatory exterior solution
 * small correction in Poisson tau also...
 *
 * Revision 1.70  2005/08/25 12:14:07  p_grandclement
 * Addition of a new method to solve the scalar Poisson equation, based on a multi-domain Tau-method
 *
 * Revision 1.69  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.68  2005/06/08 12:35:18  j_novak
 * New method for solving divergence-like ODEs.
 *
 * Revision 1.67  2005/05/20 14:42:30  j_novak
 * Added the method Scalar::get_spectral_base().
 *
 * Revision 1.66  2005/04/04 21:28:57  e_gourgoulhon
 * Added argument lambda to method poisson_angu
 * to deal with the generalized angular Poisson equation:
 *    Lap_ang u + lambda u = source.
 *
 * Revision 1.65  2004/12/14 09:09:39  f_limousin
 * Modif. comments.
 *
 * Revision 1.64  2004/11/23 12:41:53  f_limousin
 * Intoduce function poisson_dir_neu(...) to solve a scalar poisson
 * equation with a mixed boundary condition (Dirichlet + Neumann).
 *
 * Revision 1.63  2004/10/11 15:09:00  j_novak
 * The radial manipulation functions take Scalar as arguments, instead of Cmp.
 * Added a conversion operator from Scalar to Cmp.
 * The Cmp radial manipulation function make conversion to Scalar, call to the
 * Map_radial version with a Scalar argument and back.
 *
 * Revision 1.62  2004/08/24 09:14:40  p_grandclement
 * Addition of some new operators, like Poisson in 2d... It now requieres the
 * GSL library to work.
 *
 * Also, the way a variable change is stored by a Param_elliptic is changed and
 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
 * will requiere some modification. (It should concern only the ones about monopoles)
 *
 * Revision 1.61  2004/07/27 08:24:26  j_novak
 * Modif. comments
 *
 * Revision 1.60  2004/07/26 16:02:21  j_novak
 * Added a flag to specify whether the primitive should be zero either at r=0
 * or at r going to infinity.
 *
 * Revision 1.59  2004/07/06 13:36:27  j_novak
 * Added methods for desaliased product (operator |) only in r direction.
 *
 * Revision 1.58  2004/06/22 08:49:57  p_grandclement
 * Addition of everything needed for using the logarithmic mapping
 *
 * Revision 1.57  2004/06/14 15:24:23  e_gourgoulhon
 * Added method primr (radial primitive).
 *
 * Revision 1.56  2004/06/11 14:29:56  j_novak
 * Scalar::multipole_spectrum() is now a const method.
 *
 * Revision 1.55  2004/05/24 14:07:31  e_gourgoulhon
 * Method set_domain now includes a call to del_deriv() for safety.
 *
 * Revision 1.54  2004/05/07 11:26:10  f_limousin
 * New method filtre_r(int*)
 *
 * Revision 1.53  2004/03/22 13:12:43  j_novak
 * Modification of comments to use doxygen instead of doc++
 *
 * Revision 1.52  2004/03/17 15:58:47  p_grandclement
 * Slight modification of sol_elliptic_no_zec
 *
 * Revision 1.51  2004/03/11 12:07:30  e_gourgoulhon
 * Added method visu_section_anim.
 *
 * Revision 1.50  2004/03/08 15:45:38  j_novak
 * Modif. comment
 *
 * Revision 1.49  2004/03/05 15:09:40  e_gourgoulhon
 * Added method smooth_decay.
 *
 * Revision 1.48  2004/03/01 09:57:02  j_novak
 * the wave equation is solved with Scalars. It now accepts a grid with a
 * compactified external domain, which the solver ignores and where it copies
 * the values of the field from one time-step to the next.
 *
 * Revision 1.47  2004/02/27 09:43:58  f_limousin
 * New methods filtre_phi(int) and filtre_theta(int).
 *
 * Revision 1.46  2004/02/26 22:46:26  e_gourgoulhon
 * Added methods derive_cov, derive_con and derive_lie.
 *
 * Revision 1.45  2004/02/21 17:03:49  e_gourgoulhon
 * -- Method "point" renamed "val_grid_point".
 * -- Method "set_point" renamed "set_grid_point".
 *
 * Revision 1.44  2004/02/19 22:07:35  e_gourgoulhon
 * Added argument "comment" in method spectral_display.
 *
 * Revision 1.43  2004/02/11 09:47:44  p_grandclement
 * Addition of a new elliptic solver, matching with the homogeneous solution
 * at the outer shell and not solving in the external domain (more details
 * coming soon ; check your local Lorene dealer...)
 *
 * Revision 1.42  2004/01/28 16:46:22  p_grandclement
 * Addition of the sol_elliptic_fixe_der_zero stuff
 *
 * Revision 1.41  2004/01/28 13:25:40  j_novak
 * The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
 * In the div/mult _r_dzpuis, there is no more default value.
 *
 * Revision 1.40  2004/01/28 10:39:17  j_novak
 * Comments modified.
 *
 * Revision 1.39  2004/01/27 15:10:01  j_novak
 * New methods Scalar::div_r_dzpuis(int) and Scalar_mult_r_dzpuis(int)
 * which replace div_r_inc*. Tried to clean the dzpuis handling.
 * WARNING: no testing at this point!!
 *
 * Revision 1.38  2004/01/23 13:25:44  e_gourgoulhon
 * Added methods set_inner_boundary and set_outer_boundary.
 * Methods set_val_inf and set_val_hor, which are particular cases of
 * the above, have been suppressed.
 *
 * Revision 1.37  2004/01/22 16:10:09  e_gourgoulhon
 * Added (provisory) method div_r_inc().
 *
 * Revision 1.36  2003/12/16 06:32:20  e_gourgoulhon
 * Added method visu_box.
 *
 * Revision 1.35  2003/12/14 21:46:35  e_gourgoulhon
 * Added argument start_dx in visu_section.
 *
 * Revision 1.34  2003/12/11 16:19:38  e_gourgoulhon
 * Added method visu_section for visualization with OpenDX.
 *
 * Revision 1.33  2003/12/11 14:48:47  p_grandclement
 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
 *
 * Revision 1.32  2003/11/13 13:43:53  p_grandclement
 * Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
 *
 * Revision 1.31  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.30  2003/11/04 22:55:50  e_gourgoulhon
 * Added new methods mult_cost(), mult_sint() and div_sint().
 *
 * Revision 1.29  2003/10/29 13:09:11  e_gourgoulhon
 * -- Added integer argument to derivative functions dsdr, etc...
 *    so that one can choose the dzpuis of the result (default=2).
 * -- Change of method name: laplacien --> laplacian.
 *
 * Revision 1.28  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.27  2003/10/20 14:26:02  j_novak
 * New assignement operators.
 *
 * Revision 1.26  2003/10/19 19:46:33  e_gourgoulhon
 * -- Method spectral_display now virtual (from Tensor), list of argument
 *    changed.
 *
 * Revision 1.25  2003/10/17 13:46:14  j_novak
 * The argument is now between 1 and 3 (instead of 0->2)
 *
 * Revision 1.24  2003/10/16 15:23:41  e_gourgoulhon
 * Name of method div_r_ced() changed to div_r_inc2().
 * Name of method div_rsint_ced() changed to div_rsint_inc2().
 *
 * Revision 1.23  2003/10/15 21:10:11  e_gourgoulhon
 * Added method poisson_angu().
 *
 * Revision 1.22  2003/10/15 16:03:35  j_novak
 * Added the angular Laplace operator for Scalar.
 *
 * Revision 1.21  2003/10/15 10:29:05  e_gourgoulhon
 * Added new members p_dsdt and p_stdsdp.
 * Added new methods dsdt(), stdsdp() and div_tant().
 *
 * Revision 1.20  2003/10/13 13:52:39  j_novak
 * Better managment of derived quantities.
 *
 * Revision 1.19  2003/10/10 15:57:27  j_novak
 * Added the state one (ETATUN) to the class Scalar
 *
 * Revision 1.18  2003/10/08 14:24:08  j_novak
 * replaced mult_r_zec with mult_r_ced
 *
 * Revision 1.17  2003/10/06 16:16:02  j_novak
 * New constructor from a Tensor.
 *
 * Revision 1.16  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.15  2003/10/05 21:06:31  e_gourgoulhon
 * - Added new methods div_r_ced() and div_rsint_ced().
 * - Added new virtual method std_spectral_base()
 * - Removed method std_spectral_base_scal()
 *
 * Revision 1.14  2003/10/01 13:02:58  e_gourgoulhon
 * Suppressed the constructor from Map* .
 *
 * Revision 1.13  2003/09/29 12:52:56  j_novak
 * Methods for changing the triad are implemented.
 *
 * Revision 1.12  2003/09/25 09:33:36  j_novak
 * Added methods for integral calculation and various manipulations
 *
 * Revision 1.11  2003/09/25 09:11:21  e_gourgoulhon
 * Added functions for radial operations (divr, etc...)
 *
 * Revision 1.10  2003/09/25 08:55:23  e_gourgoulhon
 * Added members raccord*.
 *
 * Revision 1.9  2003/09/25 08:50:11  j_novak
 * Added the members import
 *
 * Revision 1.8  2003/09/25 08:13:51  j_novak
 * Added method for calculating derivatives
 *
 * Revision 1.7  2003/09/25 07:59:26  e_gourgoulhon
 * Added prototypes for PDE resolutions.
 *
 * Revision 1.6  2003/09/25 07:17:58  j_novak
 * Method asymptot implemented.
 *
 * Revision 1.5  2003/09/24 20:53:38  e_gourgoulhon
 * Added  -- constructor by conversion from a Cmp
 *        -- assignment from Cmp
 *
 * 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 12:01:44  j_novak
 * Added friend functions for math.
 *
 * Revision 1.2  2003/09/24 10:22:01  e_gourgoulhon
 * still in progress...
 *
 * Revision 1.1  2003/09/22 12:50:47  e_gourgoulhon
 * First version: not ready yet!
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/scalar.h,v 1.96 2023/08/31 08:27:26 g_servignat Exp $
 *
 */

// Headers Lorene 
#include "valeur.h"
#include "tensor.h"

namespace Lorene {
class Param ; 
class Cmp ;
class Param_elliptic ;

class Scalar : public Tensor {
  
  // Data : 
  // -----
 protected:
  
  int etat ; 
  
  int dzpuis ;  
  
  Valeur va ;       
  
  // Derived data : 
  // ------------
 protected:
  mutable Scalar* p_dsdr ;  

  mutable Scalar* p_srdsdt ;    

  mutable Scalar* p_srstdsdp ;

  mutable Scalar* p_dsdt ;  

  mutable Scalar* p_stdsdp ;    
  
  mutable Scalar* p_dsdx ;  
  
  mutable Scalar* p_dsdy ;  

  mutable Scalar* p_dsdz ;  
  
  mutable Scalar* p_lap ;   
  
  mutable Scalar* p_lapang ;    
   
  mutable Scalar* p_dsdradial ; 

  mutable Scalar* p_dsdrho ;    

  mutable int ind_lap ; 

  mutable Tbl* p_integ ; 
  
  // Constructors - Destructor
  // -------------------------
  
 public:
  
  explicit Scalar(const Map& mpi) ; 

  Scalar(const Tensor& a) ;     

  Scalar(const Scalar& a) ;     
  
  Scalar(const Cmp& a) ;    
  
  Scalar(const Map&, const Mg3d&, FILE* ) ;         
  
  virtual ~Scalar() ;           
  
  
  // Memory management
  // -----------------
 protected:
  void del_t() ;            
  virtual void del_deriv() const;       
  void set_der_0x0() const;     
  
 public:
  
  virtual void set_etat_nondef() ;   
  
  virtual void set_etat_zero() ;        
  
  virtual void set_etat_qcq() ;     
  
  void set_etat_one() ;     
  
  virtual void allocate_all() ; 
  
  void annule_hard() ;
  
  // Extraction of information
  // -------------------------
    public:
  int get_etat() const {return etat;} ; 
  
  int get_dzpuis() const {return dzpuis;} ; 
  
  bool dz_nonzero() const ; 
    
  bool check_dzpuis(int dzi) const ; 

  bool is_nan(bool verbose=false) const ;

  // Assignment
  // -----------
 public: 
  void operator=(const Scalar& a) ; 
  
  virtual void operator=(const Tensor& a) ; 

  void operator=(const Cmp& a) ;    
  void operator=(const Valeur& a) ; 
  void operator=(const Mtbl& a) ;    
  void operator=(double ) ;  
  void operator=(int ) ;         
  
  // Conversion oprators
  //---------------------
  operator Cmp() const ; 

  // Access to individual elements
  // -----------------------------
    public:
  
  const Valeur& get_spectral_va() const {return va;} ; 
  
  Valeur& set_spectral_va() {return va;} ; 
  
  Tbl& set_domain(int l) {
    assert(etat == ETATQCQ) ;
    del_deriv() ; 
    return va.set(l) ;
  };
  
  const Tbl& domain(int l) const {
    assert( (etat == ETATQCQ) || (etat == ETATUN) ) ;
    return va(l) ;
  };
  
  
  double val_grid_point(int l, int k, int j, int i) const {
    assert(etat != ETATNONDEF) ;
    if (etat == ETATZERO) {
      double zero = 0. ;
      return zero ; 
    }
    else {
      if (etat == ETATUN) {
      double one = 1. ;
      return one ;
      }
      else{         
    return va(l, k, j, i) ;
      }
    }
  };
  
  double val_point(double r, double theta, double phi) const ; 
  
  
  double& set_grid_point(int l, int k, int j, int i) {
    assert(etat == ETATQCQ) ;
    return va.set(l, k, j, i) ;
  };
  
    
  virtual void annule(int l_min, int l_max) ; 
  
  void set_inner_boundary(int l, double x) ;
    
  void set_outer_boundary(int l, double x) ;

  Tbl multipole_spectrum () const ;

  Tbl tbl_out_bound(int l_dom, bool leave_ylm = false) ;
  
  Tbl tbl_in_bound(int n, bool leave_ylm = false) ;
  
  Scalar scalar_out_bound(int n, bool leave_ylm = false) ;
  
  // Differential operators and others
  // ---------------------------------
 public:
  const Scalar& dsdr() const ; 
 
  const Scalar& srdsdt() const ; 
  
  const Scalar& srstdsdp() const ; 
  
  const Scalar& dsdt() const ; 
  
  const Scalar& dsdradial() const ; 
  
  const Scalar& dsdrho() const ;
 
  const Scalar& stdsdp() const ; 
  
  const Scalar& dsdx() const ;  
  
  const Scalar& dsdy() const ;  
  
  const Scalar& dsdz() const ;  
  
  const Scalar& deriv(int i) const ;    
  
    const Vector& derive_cov(const Metric& gam) const ; 


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

    Scalar derive_lie(const Vector& v) const ; 


  const Scalar& laplacian(int ced_mult_r = 4) const ; 
  
  const Scalar& lapang() const ; 
  
  void div_r() ;    
 
  void div_r_dzpuis(int ced_mult_r) ; 
  
  void div_r_ced() ;

  void mult_r() ;  
  
  void mult_r_dzpuis(int ced_mult_r) ; 
  
  void mult_r_ced() ;
  
  void mult_rsint() ;   
  
  void mult_rsint_dzpuis(int ced_mult_r) ; 

  void div_rsint() ;    
  
  void div_rsint_dzpuis(int ced_mult_r) ; 
  
  void mult_cost() ;   

  void div_cost() ;    

  void mult_sint() ;   

  void div_sint() ;    

  void div_tant() ;    
  
  Scalar primr(bool null_infty = true) const  ;         

  double integrale() const ; 
    
  const Tbl& integrale_domains() const ; 
    
  virtual void dec_dzpuis(int dec = 1) ; 

  virtual void inc_dzpuis(int inc = 1) ; 
    
  virtual void change_triad(const Base_vect& new_triad) ; 
    
  void filtre (int n) ;
    
  void filtre_r (int* nn) ;

  void filtre_r (int n, int nzone) ;

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

  void sarra_filter_r(int lzmin, int lzmax, double p, 
                double alpha= -1E-16) ;

  void exp_filter_r_all_domains(Scalar &ss, int p, double alpha=-16.) ;

  void sarra_filter_r_all_domains(double p, double alpha=1E-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.) ;

  void annule_l (int l_min, int l_max, bool ylm_output= false ) ;

  void filtre_phi (int n, int zone) ;
 
  void filtre_tp(int nn, int nz1, int nz2) ;
  
   
  void fixe_decroissance (int puis) ;

    void smooth_decay(int k, int n) ; 

  void raccord(int n) ;
    
  void raccord_c1_zec(int puis, int nbre, int lmax) ;

  void raccord_externe(int puis, int nbre, int lmax) ;

  void match_tau(Param& par_bc, Param* par_mat=0x0) ;

  void match_tau_star(Param& par_bc, Param* par_mat=0x0) ;

  void match_tau_shell(Param& par_bc, Param* par_mat=0x0) ;

  void match_collocation(Param& par_bc, Param* par_mat=0x0) ;

  // 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 Scalar & ) ;  
  
    void visu_section(const char section_type, double aa, double umin, double umax, double vmin,
        double vmax, const char* title = 0x0, const char* filename = 0x0,
        bool start_dx = true, int nu = 200, int nv = 200) const ;   

    void visu_section(const Tbl& plane, double umin, double umax, double vmin,
        double vmax, const char* title = 0x0, const char* filename = 0x0,
        bool start_dx = true, int nu = 200, int nv = 200) const ;   

   void visu_section_anim(const char section_type, double aa, double umin, 
        double umax, double vmin, double vmax, int jtime, double ttime, 
        int jgraph = 1, const char* title = 0x0, const char* filename_root = 0x0, 
        bool start_dx = false, int nu = 200, int nv = 200) const ;   

    void visu_box(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 = 40, int ny = 40, 
        int nz = 40) const ;      
        


  // Member arithmetics
  // ------------------
 public:
  void operator+=(const Scalar &) ;         
  void operator-=(const Scalar &) ;         
  void operator*=(const Scalar &) ;         

  // Manipulation of spectral bases
  // ------------------------------    
  virtual void std_spectral_base() ;     
  
  virtual void std_spectral_base_odd() ;     
  
  void set_spectral_base(const Base_val& ) ;     

  const Base_val& get_spectral_base( ) const {return va.base ;} ;    

  void set_dzpuis(int ) ; 

  Valeur** asymptot(int n, const int flag = 0) const ; 
    

  // PDE resolution 
  // --------------
 public:
  Scalar poisson() const ;

  void poisson(Param& par, Scalar& uu) const ;
      
  Scalar poisson_tau() const ;
   
  void poisson_tau(Param& par, Scalar& uu) const ;

  Scalar poisson_dirichlet (const Valeur& limite, int num) const ;
    
  Scalar poisson_neumann   (const Valeur&, int) const ;


  Scalar poisson_dir_neu  (const Valeur& limite , int num, 
               double fact_dir, double fact_neu) const ;

  Scalar poisson_frontiere_double   (const Valeur&, const Valeur&, int) const ;

  void poisson_regular(int k_div, int nzet, double unsgam1, Param& par,
               Scalar& uu, Scalar& uu_regu, Scalar& uu_div,
               Tensor& duu_div,
               Scalar& source_regu, Scalar& source_div) const ;

  Tbl test_poisson(const Scalar& uu, ostream& ostr, 
           bool detail = false) const ;  

    Scalar poisson_angu(double lambda =0) const ;

  Scalar avance_dalembert(Param& par, const Scalar& fJm1, const Scalar& source) 
    const ;

  Scalar sol_elliptic(Param_elliptic& params) const ;
 
  Scalar sol_elliptic_boundary(Param_elliptic& params, const Mtbl_cf& bound,
                   double fact_dir, double fact_neu) const ;
 
  Scalar sol_elliptic_boundary(Param_elliptic& params, const Scalar& bound,
                   double fact_dir, double fact_neu) const ;
 

  Scalar sol_elliptic_2d(Param_elliptic&) const ;
 
  Scalar sol_elliptic_pseudo_1d(Param_elliptic&) const ;

  Scalar sol_elliptic_no_zec(Param_elliptic& params, double val = 0) const ;
  
  Scalar sol_elliptic_only_zec(Param_elliptic& params, double val) const ;

  Scalar sol_elliptic_sin_zec(Param_elliptic& params, double* coefs, double* phases) const ;

   
  Scalar sol_elliptic_fixe_der_zero(double val, 
                    Param_elliptic& params) const ;
  

  Scalar sol_divergence(int n) const ;
    
  // Import from other mapping 
  // -------------------------

  void import(const Scalar& ci) ;    

  void import_symy(const Scalar& ci) ;   

  void import_asymy(const Scalar& ci) ;  

  void import(int nzet, const Scalar& ci) ;  

  void import_symy(int nzet, const Scalar& ci) ;     

  void import_asymy(int nzet, const Scalar& ci) ;    

 protected:
  void import_gal(int nzet, const Scalar& ci) ;  

  void import_align(int nzet, const Scalar& ci) ;    

  void import_anti(int nzet, const Scalar& ci) ;     

  void import_align_symy(int nzet, const Scalar& ci) ;   

  void import_anti_symy(int nzet, const Scalar& ci) ;    

  void import_align_asymy(int nzet, const Scalar& ci) ;  

  void import_anti_asymy(int nzet, const Scalar& ci) ;   
    

  friend Scalar operator-(const Scalar& ) ;         
  friend Scalar operator+(const Scalar&, const Scalar &) ;  
  friend Scalar operator+(const Scalar&, const Mtbl&) ; 
  friend Scalar operator+(const Scalar&, double ) ;     
  friend Scalar operator-(const Scalar &, const Scalar &) ;
  friend Scalar operator-(const Scalar&, const Mtbl&) ; 
  friend Scalar operator-(const Scalar&, double ) ;     
  friend Scalar operator*(const Scalar &, const Scalar &) ;
  friend Scalar operator%(const Scalar &, const Scalar &) ;
  friend Scalar operator|(const Scalar &, const Scalar &) ;
  friend Scalar operator*(const Mtbl&, const Scalar &) ;        
  friend Scalar operator*(double, const Scalar &) ;     
  friend Scalar operator/(const Scalar &, const Scalar &) ;
  friend Scalar operator/(const Scalar &, const Mtbl &) ;
  friend Scalar operator/(const Mtbl &, const Scalar &) ;
  friend Scalar operator/(const Scalar&, double ) ;        
  friend Scalar operator/(double, const Scalar &) ;

  friend Scalar sin(const Scalar& ) ;
  friend Scalar cos(const Scalar& ) ;
  friend Scalar tan(const Scalar& ) ;
  friend Scalar asin(const Scalar& ) ;
  friend Scalar acos(const Scalar& ) ;
  friend Scalar atan(const Scalar& ) ;
  friend Scalar exp(const Scalar& ) ;   
  friend Scalar Heaviside(const Scalar& ) ; 
  friend Scalar log(const Scalar& ) ;   
  friend Scalar log10(const Scalar& ) ; 
  friend Scalar sqrt(const Scalar& ) ;  
  friend Scalar racine_cubique (const Scalar& ) ;
  friend Scalar pow(const Scalar& , int ) ; 
  friend Scalar pow(const Scalar& , double ) ; 
  friend Scalar abs(const Scalar& ) ;   

  friend double totalmax(const Scalar& ) ;   
  friend double totalmin(const Scalar& ) ;   
  friend Tbl max(const Scalar& ) ;   
  friend Tbl min(const Scalar& ) ;   
  friend Tbl norme(const Scalar& ) ;   
  friend Tbl diffrel(const Scalar& a, const Scalar& b) ; 
  friend Tbl diffrelmax(const Scalar& a, const Scalar& b) ; 

};

ostream& operator<<(ostream& , const Scalar & ) ;   

// Prototypage de l'arithmetique
Scalar operator+(const Scalar& ) ;          
Scalar operator-(const Scalar& ) ;          
Scalar operator+(const Scalar&, const Scalar &) ;   
Scalar operator+(const Scalar&, const Mtbl&) ;  
Scalar operator+(const Mtbl&, const Scalar&) ;  
Scalar operator+(const Scalar&, double ) ;      
Scalar operator+(double, const Scalar& ) ;      
Scalar operator+(const Scalar&, int ) ;     
Scalar operator+(int, const Scalar& ) ;     
Scalar operator-(const Scalar &, const Scalar &) ;  
Scalar operator-(const Scalar&, const Mtbl&) ;  
Scalar operator-(const Mtbl&, const Scalar&) ;  
Scalar operator-(const Scalar&, double ) ;      
Scalar operator-(double, const Scalar& ) ;      
Scalar operator-(const Scalar&, int ) ;     
Scalar operator-(int, const Scalar& ) ;     
Scalar operator*(const Scalar &, const Scalar &) ;  

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

Scalar operator|(const Scalar &, const Scalar &) ;

Scalar operator*(const Mtbl&, const Scalar&) ; 
Scalar operator*(const Scalar&, const Mtbl&) ; 
    
Scalar operator*(const Scalar&, double ) ;      
Scalar operator*(double, const Scalar &) ;      
Scalar operator*(const Scalar&, int ) ;     
Scalar operator*(int, const Scalar& ) ;     
Scalar operator/(const Scalar &, const Scalar &) ;  
Scalar operator/(const Scalar&, double ) ;      
Scalar operator/(double, const Scalar &) ;      
Scalar operator/(const Scalar&, int ) ;     
Scalar operator/(int, const Scalar &) ;     
Scalar operator/(const Scalar &, const Mtbl&) ; 
Scalar operator/(const Mtbl&, const Scalar &) ; 
    

Scalar sin(const Scalar& ) ;        
Scalar cos(const Scalar& ) ;        
Scalar tan(const Scalar& ) ;        
Scalar asin(const Scalar& ) ;       
Scalar acos(const Scalar& ) ;       
Scalar atan(const Scalar& ) ;       
Scalar exp(const Scalar& ) ;        
Scalar Heaviside(const Scalar& ) ;      
Scalar log(const Scalar& ) ;        
Scalar log10(const Scalar& ) ;  
Scalar sqrt(const Scalar& ) ;       
Scalar racine_cubique (const Scalar& ) ;        
Scalar pow(const Scalar& , int ) ;  
Scalar pow(const Scalar& , double ) ; 
Scalar abs(const Scalar& ) ;        

double totalmax(const Scalar& ) ;   

double totalmin(const Scalar& ) ;   

Tbl max(const Scalar& ) ;   

Tbl min(const Scalar& ) ;   

Tbl norme(const Scalar& ) ;   

Tbl diffrel(const Scalar& a, const Scalar& b) ; 

Tbl diffrelmax(const Scalar& a, const Scalar& b) ; 

void exp_filter_ylm_all_domains(Scalar& ss, int p, double alpha=-16.) ;

void exp_filter_ylm_all_domains(Scalar& ss, int p_tet, int p_phi, double alpha=-16.) ;

}
#endif