map_log_elliptic.C

char map_log_elliptic_C[] = "$Header $" ;

/*
 * $Id: map_log_elliptic.C,v 1.4 2007/01/16 15:08:07 n_vasset Exp $
 * $Log: map_log_elliptic.C,v $
 * Revision 1.4  2007/01/16 15:08:07  n_vasset
 * New constructor, usn Scalar on mono-domain angular grid for boundary,
 * for function sol_elliptic_boundary()
 *
 * Revision 1.3  2005/06/09 07:57:32  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.2  2004/08/24 09:14:42  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.1  2004/06/22 08:49:58  p_grandclement
 * Addition of everything needed for using the logarithmic mapping
 *
 * Revision 1.4  2004/03/17 15:58:48  p_grandclement
 *
 * $Header: /cvsroot/Lorene/C++/Source/Map/map_log_elliptic.C,v 1.4 2007/01/16 15:08:07 n_vasset Exp $
 *
 */

// Header C : 
#include <stdlib.h>
#include <math.h>

// Headers Lorene :
#include "tbl.h"
#include "mtbl_cf.h"
#include "map.h"
#include "param_elliptic.h"
         
            //----------------------------------------------
       //       General elliptic solver
      //----------------------------------------------

void Map_log::sol_elliptic(Param_elliptic& ope_var, const Scalar& source, 
              Scalar& pot) const {
    
  assert(source.get_etat() != ETATNONDEF) ; 
  assert(source.get_mp().get_mg() == mg) ; 
  assert(pot.get_mp().get_mg() == mg) ; 
  
  assert(source.check_dzpuis(2) || source.check_dzpuis(3) || 
     source.check_dzpuis(4)) ; 
  // Spherical harmonic expansion of the source
  // ------------------------------------------
  
  const Valeur& sourva = source.get_spectral_va() ; 
  
  if (sourva.get_etat() == ETATZERO) {
    pot.set_etat_zero() ;
    return ;  
    }
  
  // Spectral coefficients of the source
  assert(sourva.get_etat() == ETATQCQ) ; 
  
  Valeur rho(sourva.get_mg()) ; 
  sourva.coef() ; 
  rho = *(sourva.c_cf) ;    // copy of the coefficients of the source
  
  rho.ylm() ;           // spherical harmonic transforms 
    
  // On met les bonnes bases dans le changement de variable 
  // de ope_var :
  ope_var.var_F.set_spectral_va().coef() ;
  ope_var.var_F.set_spectral_va().ylm() ;
  ope_var.var_G.set_spectral_va().coef() ;
  ope_var.var_G.set_spectral_va().ylm() ;

  // Call to the Mtbl_cf version
  // ---------------------------
  Mtbl_cf resu = elliptic_solver (ope_var, *(rho.c_cf)) ;
  // Final result returned as a Scalar
  // ---------------------------------
  
  pot.set_etat_zero() ;  // to call Scalar::del_t().
  
  pot.set_etat_qcq() ; 
  
  pot.set_spectral_va() = resu ;
  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.       
  
  pot.set_dzpuis(0) ; 
}


         //--------------------------------------------------
         //     General elliptic solver with boundary as Mtbl_cf
         //--------------------------------------------------


void Map_log::sol_elliptic_boundary(Param_elliptic& ope_var, const Scalar& source, 
               Scalar& pot, const Mtbl_cf& bound, 
               double fact_dir, double fact_neu) const {
    
  assert(source.get_etat() != ETATNONDEF) ; 
  assert(source.get_mp().get_mg() == mg) ; 
  assert(pot.get_mp().get_mg() == mg) ; 
  
  assert(source.check_dzpuis(2) || source.check_dzpuis(3) || 
     source.check_dzpuis(4)) ; 
  // Spherical harmonic expansion of the source
  // ------------------------------------------
  
  const Valeur& sourva = source.get_spectral_va() ; 
  
  if (sourva.get_etat() == ETATZERO) {
    pot.set_etat_zero() ;
    return ;  
    }
  
  // Spectral coefficients of the source
  assert(sourva.get_etat() == ETATQCQ) ; 
  
  Valeur rho(sourva.get_mg()) ; 
  sourva.coef() ; 
  rho = *(sourva.c_cf) ;    // copy of the coefficients of the source
  
  rho.ylm() ;           // spherical harmonic transforms 
    
  // On met les bonnes bases dans le changement de variable 
  // de ope_var :
  ope_var.var_F.set_spectral_va().coef() ;
  ope_var.var_F.set_spectral_va().ylm() ;
  ope_var.var_G.set_spectral_va().coef() ;
  ope_var.var_G.set_spectral_va().ylm() ;

  // Call to the Mtbl_cf version
  // ---------------------------
  Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound, 
                       fact_dir, fact_neu) ;
  // Final result returned as a Scalar
  // ---------------------------------
  
  pot.set_etat_zero() ;  // to call Scalar::del_t().
  
  pot.set_etat_qcq() ; 
  
  pot.set_spectral_va() = resu ;
  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.       
  
  pot.set_dzpuis(0) ; 
}

   
         //--------------------------------------------------
         //     General elliptic solver with boundary as Scalar
         //--------------------------------------------------


void Map_log::sol_elliptic_boundary(Param_elliptic& ope_var, const Scalar& source, 
               Scalar& pot, const Scalar& bound, 
               double fact_dir, double fact_neu) const {
    
  assert(source.get_etat() != ETATNONDEF) ; 
  assert(source.get_mp().get_mg() == mg) ; 
  assert(pot.get_mp().get_mg() == mg) ; 
  
  assert(source.check_dzpuis(2) || source.check_dzpuis(3) || 
     source.check_dzpuis(4)) ; 
  // Spherical harmonic expansion of the source
  // ------------------------------------------
  
  const Valeur& sourva = source.get_spectral_va() ; 
  
  if (sourva.get_etat() == ETATZERO) {
    pot.set_etat_zero() ;
    return ;  
    }
  
  // Spectral coefficients of the source
  assert(sourva.get_etat() == ETATQCQ) ; 
  
  Valeur rho(sourva.get_mg()) ; 
  sourva.coef() ; 
  rho = *(sourva.c_cf) ;    // copy of the coefficients of the source
  
  rho.ylm() ;           // spherical harmonic transforms 
    
  // On met les bonnes bases dans le changement de variable 
  // de ope_var :
  ope_var.var_F.set_spectral_va().coef() ;
  ope_var.var_F.set_spectral_va().ylm() ;
  ope_var.var_G.set_spectral_va().coef() ;
  ope_var.var_G.set_spectral_va().ylm() ;

  // Call to the Mtbl_cf version
  // ---------------------------
    Scalar bbound = bound;
  bbound.set_spectral_va().ylm() ;
  const Map& mapp = bbound.get_mp();
 
  const Mg3d& gri2d = *mapp.get_mg();

  assert(&gri2d == source.get_mp().get_mg()->get_angu_1dom()) ;
  
  Mtbl_cf bound2 (gri2d , bbound.get_spectral_base()) ;
  bound2.annule_hard() ;  

  if (bbound.get_etat() != ETATZERO){ 
 
      int nr = gri2d.get_nr(0) ;
      int nt = gri2d.get_nt(0) ; 
      int np = gri2d.get_np(0) ; 
       
      for(int k=0; k<np+2; k++)
          for (int j=0; j<=nt-1; j++)
          for(int xi=0; xi<= nr-1; xi++)
          {
  bound2.set(0, k , j , xi) = (*bbound.get_spectral_va().c_cf)(0, k, j, xi) ;   
  }
  }  
 
  
  Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound2, 
                       fact_dir, fact_neu) ;
  // Final result returned as a Scalar
  // ---------------------------------
  
  pot.set_etat_zero() ;  // to call Scalar::del_t().
  
  pot.set_etat_qcq() ; 
  
  pot.set_spectral_va() = resu ;
  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.       
  
  pot.set_dzpuis(0) ; 
}

  
            //------------------------------------------------
       //       General elliptic solver with no zec
      //-------------------------------------------------

void Map_log::sol_elliptic_no_zec (Param_elliptic& ope_var, const Scalar& source, 
              Scalar& pot, double val) const {
    
  assert(source.get_etat() != ETATNONDEF) ; 
  assert(source.get_mp().get_mg() == mg) ; 
  assert(pot.get_mp().get_mg() == mg) ; 
  
  assert(source.check_dzpuis(2) || source.check_dzpuis(3) || 
     source.check_dzpuis(4)) ; 
  // Spherical harmonic expansion of the source
  // ------------------------------------------
  
  const Valeur& sourva = source.get_spectral_va() ; 
  
  if (sourva.get_etat() == ETATZERO) {
    pot.set_etat_zero() ;
    return ;  
    }
  
  // Spectral coefficients of the source
  assert(sourva.get_etat() == ETATQCQ) ; 
  
  Valeur rho(sourva.get_mg()) ; 
  sourva.coef() ; 
  rho = *(sourva.c_cf) ;    // copy of the coefficients of the source
  
  rho.ylm() ;           // spherical harmonic transforms 
    
  // On met les bonnes bases dans le changement de variable 
  // de ope_var :
  ope_var.var_F.set_spectral_va().coef() ;
  ope_var.var_F.set_spectral_va().ylm() ;
  ope_var.var_G.set_spectral_va().coef() ;
  ope_var.var_G.set_spectral_va().ylm() ;

  // Call to the Mtbl_cf version
  // ---------------------------
  Mtbl_cf resu = elliptic_solver_no_zec (ope_var, *(rho.c_cf), val) ;
  // Final result returned as a Scalar
  // ---------------------------------
  
  pot.set_etat_zero() ;  // to call Scalar::del_t().
  
  pot.set_etat_qcq() ; 
  
  pot.set_spectral_va() = resu ;
  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.       
  
  pot.set_dzpuis(0) ; 
}

   

---