scalar_manip.C

/*
 *   Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
 *   
 *   Copyright (c) 2000-2001 Philippe Grandclement (Cmp version)
 *
 *   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
 *
 */


 

/*
 * $Id: scalar_manip.C,v 1.23 2025/03/31 08:29:50 j_novak Exp $
 * $Log: scalar_manip.C,v $
 * Revision 1.23  2025/03/31 08:29:50  j_novak
 * Methods to copyt Tensors and Scalars to smaller grids...
 *
 * Revision 1.22  2025/03/28 09:12:50  j_novak
 * News methods to copy Tensors and Scalars to larger grids. To be used for de-aliasing.
 *
 * Revision 1.21  2018/11/16 14:34:37  j_novak
 * Changed minor points to avoid some compilation warnings.
 *
 * Revision 1.20  2016/12/05 16:18:18  j_novak
 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
 *
 * Revision 1.19  2014/10/13 08:53:46  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.18  2014/10/06 15:16:15  j_novak
 * Modified #include directives to use c++ syntax.
 *
 * Revision 1.17  2008/10/03 09:03:52  j_novak
 * Correction of yet another mistake (the array values in physical space was not
 * destroyed).
 *
 * Revision 1.16  2008/09/30 08:35:18  j_novak
 * Correction of forgotten call to coef()
 *
 * Revision 1.15  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.14  2008/09/22 19:08:01  j_novak
 * New methods to deal with boundary conditions
 *
 * Revision 1.13  2007/06/21 20:00:00  k_taniguchi
 * Addition of the method filtre_r (int n, int nz)
 *
 * Revision 1.12  2006/06/28 07:46:41  j_novak
 * Better treatment in the case of a domain set to zero.
 *
 * Revision 1.11  2005/09/07 13:39:10  j_novak
 * *** empty log message ***
 *
 * Revision 1.10  2005/09/07 13:10:48  j_novak
 * Added a filter setting to zero all mulitpoles in a given range.
 *
 * Revision 1.9  2004/11/23 12:47:44  f_limousin
 * Add function filtre_tp(int nn, int nz1, int nz2).
 *
 * Revision 1.8  2004/05/07 11:24:54  f_limousin
 * Implement new method filtre_r (int* nn)
 *
 * Revision 1.7  2004/02/27 09:47:26  f_limousin
 * New methods filtre_phi(int) and filtre_theta(int).
 *
 * Revision 1.6  2004/01/23 13:26:28  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.5  2003/11/13 13:43:55  p_grandclement
 * Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
 *
 * Revision 1.4  2003/10/11 14:47:17  e_gourgoulhon
 * Lines 112 and 145 : replaced "0" by "double(0)" in the comparison
 * statement.
 *
 * Revision 1.3  2003/10/10 15:57:29  j_novak
 * Added the state one (ETATUN) to the class Scalar
 *
 * Revision 1.2  2003/10/08 14:24:09  j_novak
 * replaced mult_r_zec with mult_r_ced
 *
 * Revision 1.1  2003/09/25 09:33:36  j_novak
 * Added methods for integral calculation and various manipulations
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_manip.C,v 1.23 2025/03/31 08:29:50 j_novak Exp $
 *
 */

//standard
#include <cstdlib>
#include <cmath>

// Lorene
#include "tensor.h"
#include "proto.h"
#include "utilitaires.h"

/*
 * Annule tous les l entre l_min et l_max (compris)
 */

namespace Lorene {
void Scalar::annule_l (int l_min, int l_max, bool ylm_output) {

    assert (etat != ETATNONDEF) ;
    assert (l_min <= l_max) ;
    assert (l_min >= 0) ;
    if (etat == ETATZERO )
    return ;

    if (etat == ETATUN) {
    if (l_min == 0) set_etat_zero() ;
    else return ;
    }

    va.ylm() ;
    Mtbl_cf& m_coef = *va.c_cf ;
    const Base_val& base = va.base ; 
    int l_q, m_q, base_r ;
    for (int lz=0; lz<mp->get_mg()->get_nzone(); lz++)
    if (m_coef(lz).get_etat() != ETATZERO)
        for (int k=0; k<mp->get_mg()->get_np(lz)+1; k++)
        for (int j=0; j<mp->get_mg()->get_nt(lz); j++)
            for (int i=0; i<mp->get_mg()->get_nr(lz); i++) {
            base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
            if ((l_min <= l_q) && (l_q<= l_max)) 
                m_coef.set(lz, k, j, i) = 0 ;
            }
    if (va.c != 0x0) {
    delete va.c ;
    va.c = 0x0 ;
    }
    if (!ylm_output)
    va.ylm_i() ;
    
    return ;
}

/*
 * Annule les n derniers coefficients en r dans la derniere zone
 */
 
void Scalar::filtre (int n) {
    
    assert (etat != ETATNONDEF) ;
    if ( (etat == ETATZERO) || (etat == ETATUN) )
    return ;
    
    int nz = mp->get_mg()->get_nzone() ;
    int np = mp->get_mg()->get_np(nz-1) ;
    int nt = mp->get_mg()->get_nt(nz-1) ;
    int nr = mp->get_mg()->get_nr(nz-1) ;
    
    del_deriv() ;
    
    va.coef() ;
    va.set_etat_cf_qcq() ;

    
    for (int k=0 ; k<np+1 ; k++)
    if (k!=1)
        for (int j=0 ; j<nt ; j++)
        for (int i=nr-1 ; i>nr-1-n ; i--)
            va.c_cf->set(nz-1, k, j, i) = 0 ;
}


/*
 * Annule les n derniers coefficients en r dans toutes les zones
 */
 
void Scalar::filtre_r (int* nn) {
    assert (etat != ETATNONDEF) ;
    if ( (etat == ETATZERO) || (etat == ETATUN) )
    return ;
    
    del_deriv() ;
    
    va.coef() ;
    va.set_etat_cf_qcq() ;
    int nz = mp->get_mg()->get_nzone() ;
    int* nr = new int[nz];
    int* nt = new int[nz];
    int* np = new int[nz];
    for (int l=0; l<=nz-1; l++) {
    nr[l] = mp->get_mg()->get_nr(l) ; 
    nt[l] = mp->get_mg()->get_nt(l) ; 
    np[l] = mp->get_mg()->get_np(l) ; 
    }
 
    for (int l=0; l<=nz-1; l++)
      for (int k=0 ; k<np[l]+1 ; k++)
    if (k!=1)
      for (int j=0 ; j<nt[l] ; j++)
        for (int i=nr[l]-1; i>nr[l]-1-nn[l] ; i--)
          va.c_cf->set(l, k, j, i) = 0 ;
    
    if (va.c != 0x0) {
      delete va.c ;
      va.c = 0x0 ;
    }
    
}


/*
 * Annule les n derniers coefficients en r dans zone nz
 */

void Scalar::filtre_r (int n, int nz) {
    assert (etat != ETATNONDEF) ;
    if ( (etat == ETATZERO) || (etat == ETATUN) )
    return ;
    
    del_deriv() ;
    
    va.coef() ;
    va.set_etat_cf_qcq() ;
    int nr = mp->get_mg()->get_nr(nz) ; 
    int nt = mp->get_mg()->get_nt(nz) ; 
    int np = mp->get_mg()->get_np(nz) ; 

    for (int k=0 ; k<np+1 ; k++)
        if (k!=1)
        for (int j=0 ; j<nt ; j++)
            for (int i=nr-1; i>nr-1-n ; i--)
            va.c_cf->set(nz, k, j, i) = 0 ;

    if (va.c != 0x0) {
        delete va.c ;
    va.c = 0x0 ;
    }

}


/*
 * Annule les n derniers coefficients en phi dans zone nz
 */
 
void Scalar::filtre_phi (int n, int nz) {
    assert (etat != ETATNONDEF) ;
    if ( (etat == ETATZERO) || (etat == ETATUN) )
    return ;
    
    del_deriv() ;
    
    va.coef() ;
    va.set_etat_cf_qcq() ;
    int np = mp->get_mg()->get_np(nz) ;
    int nt = mp->get_mg()->get_nt(nz) ;
    int nr = mp->get_mg()->get_nr(nz) ;
    
    for (int k=np+1-n ; k<np+1 ; k++)
    for (int j=0 ; j<nt ; j++)
        for (int i=0 ; i<nr ; i++)
        va.c_cf->set(nz, k, j, i) = 0 ;

}


void Scalar::filtre_tp(int nn, int nz1, int nz2) {

    va.filtre_tp(nn, nz1, nz2) ;

}


  /* Sets the value of the {\tt Scalar} at the inner boundary of a given 
   * domain. 
   * @param l [input] domain index
   * @param x [input] (constant) value at the inner boundary of domain no. {\tt l}
   */

void Scalar::set_inner_boundary(int l0, double x0) {
    
    assert (etat != ETATNONDEF) ;
    if (etat == ETATZERO) {
        if (x0 == double(0)) return ;
        else annule_hard() ;
    }
    
    if (etat == ETATUN) {
        if (x0 == double(1)) return ;
        else etat = ETATQCQ ;
    }
    
    del_deriv() ;
    
    int nt = mp->get_mg()->get_nt(l0) ;
    int np = mp->get_mg()->get_np(l0) ;
    
    va.coef_i() ;
    va.set_etat_c_qcq() ;
    
    for (int k=0 ; k<np ; k++)
    for (int j=0 ; j<nt ; j++)
        va.set(l0, k, j, 0) = x0 ;
}

  /* Sets the value of the {\tt Scalar} at the outer boundary of a given 
   * domain. 
   * @param l [input] domain index
   * @param x [input] (constant) value at the outer boundary of domain no. {\tt l}
   */

void Scalar::set_outer_boundary(int l0, double x0) {
    
    assert (etat != ETATNONDEF) ;
    if (etat == ETATZERO) {
        if (x0 == double(0)) return ;
        else annule_hard() ;
    }
    
    if (etat == ETATUN) {
        if (x0 == double(1)) return ;
        else etat = ETATQCQ ;
    }
    
    del_deriv() ;
    
    int nrm1 = mp->get_mg()->get_nr(l0) - 1 ;
    int nt = mp->get_mg()->get_nt(l0) ;
    int np = mp->get_mg()->get_np(l0) ;
    
    va.coef_i() ;
    va.set_etat_c_qcq() ;
    
    for (int k=0 ; k<np ; k++)
    for (int j=0 ; j<nt ; j++)
        va.set(l0, k, j, nrm1) = x0 ;
}

/*
 * Permet de fixer la decroissance du cmp a l infini en viurant les 
 * termes en 1/r^n
 */
void Scalar::fixe_decroissance (int puis) {
    
    if (puis<dzpuis)
    return ;
    else {
    
    int nbre = puis-dzpuis ;
    
    // le confort avant tout ! (c'est bien le confort ...)
    int nz = mp->get_mg()->get_nzone() ;
    int np = mp->get_mg()->get_np(nz-1) ;
    int nt = mp->get_mg()->get_nt(nz-1) ;
    int nr = mp->get_mg()->get_nr(nz-1) ;
    
    const Map_af* map  = dynamic_cast<const Map_af*>(mp) ;
    if (map == 0x0) {
        cout << "Le mapping doit etre affine" << endl ;
        abort() ;
    }
    
    double alpha = map->get_alpha()[nz-1] ;
    
    Scalar courant (*this) ;
    
    va.coef() ;
    va.set_etat_cf_qcq() ;
    
    for (int conte=0 ; conte<nbre ; conte++) {
        
        int base_r = courant.va.base.get_base_r(nz-1) ;
        
        courant.va.coef() ;
        
        // On calcul les coefficients de 1/r^conte
        double* coloc = new double [nr] ;
        int * deg = new int[3] ;
        deg[0] = 1 ; 
        deg[1] = 1 ;
        deg[2] = nr ;
            
        for (int i=0 ; i<nr ; i++)
        coloc[i] =pow(alpha, double(conte))*
            pow(-1-cos(M_PI*i/(nr-1)), double(conte)) ;
            
        cfrcheb(deg, deg, coloc, deg, coloc) ;
        
        for (int k=0 ; k<np+1 ; k++)
        if (k != 1)
        for (int j=0 ; j<nt ; j++) {
            
            // On doit determiner le coefficient du truc courant :
            double* coef = new double [nr] ;
            double* auxi = new double[1] ;
            for (int i=0 ; i<nr ; i++)
            coef[i] = (*courant.va.c_cf)(nz-1, k, j, i) ;
            switch (base_r) {
            case R_CHEBU :
            som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
            break ;
            default :
            som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
            break ;
            }
            
            // On modifie le cmp courant :
            courant.va.coef() ;
            courant.va.set_etat_cf_qcq() ;
            courant.va.c_cf->set(nz-1, k, j, 0) -= *auxi ;  
            
            for (int i=0 ; i<nr ; i++)
                this->va.c_cf->set(nz-1, k, j, i) -= *auxi * coloc[i] ;

              
            delete [] coef ;
            delete [] auxi ;
        }
        delete [] coloc ;
        delete [] deg ;
        
        courant.mult_r_ced() ;
    }
    }
}

Tbl Scalar::tbl_out_bound(int l_zone, bool output_ylm) {
    va.coef() ;
    if (output_ylm) va.ylm() ;

    int np = mp->get_mg()->get_np(l_zone) ;
    int nt = mp->get_mg()->get_nt(l_zone) ;
    Tbl resu(np+2, nt) ;
    if (etat == ETATZERO) resu.set_etat_zero() ;
    else {
    assert(etat == ETATQCQ) ;
    resu.set_etat_qcq() ;
    for (int k=0; k<np+2; k++)
        for (int j=0; j<nt; j++)
        resu.set(k, j) = va.c_cf->val_out_bound_jk(l_zone, j, k) ;
    }
    return resu ;
}

Tbl Scalar::tbl_in_bound(int l_zone, bool output_ylm) {
    assert(mp->get_mg()->get_type_r(l_zone) != RARE) ;
    va.coef() ;
    if (output_ylm) va.ylm() ;

    int np = mp->get_mg()->get_np(l_zone) ;
    int nt = mp->get_mg()->get_nt(l_zone) ;
    Tbl resu(np+2, nt) ;
    if (etat == ETATZERO) resu.set_etat_zero() ;
    else {
    assert(etat == ETATQCQ) ;
    resu.set_etat_qcq() ;
    for (int k=0; k<np+2; k++)
        for (int j=0; j<nt; j++)
        resu.set(k, j) = va.c_cf->val_in_bound_jk(l_zone, j, k) ;
    }
    return resu ;
}

Scalar Scalar::scalar_out_bound(int l_zone, bool output_ylm) {
    va.coef() ;
    if (output_ylm) va.ylm() ;

    Scalar resu(mp->mp_angu(l_zone)) ;
    resu.std_spectral_base() ;
    Base_val base = resu.get_spectral_base() ;
    base.set_base_t(va.base.get_base_t(l_zone)) ;
    resu.set_spectral_base(base) ;
    if (etat == ETATZERO) resu.set_etat_zero() ;
    else {
    assert(etat == ETATQCQ) ;
    resu.annule_hard() ;
    int np = mp->get_mg()->get_np(l_zone) ;
    int nt = mp->get_mg()->get_nt(l_zone) ;
    for (int k=0; k<np+2; k++)
        for (int j=0; j<nt; j++)
        resu.set_spectral_va().c_cf->set(0, k, j, 0) 
            = va.c_cf->val_out_bound_jk(l_zone, j, k) ;
    delete resu.set_spectral_va().c ;
    resu.set_spectral_va().c = 0x0 ;
    }
    return resu ;
}


  void Scalar::copy_coefs_from_small_grid(const Scalar& scal_small_grid) {

    assert(&scal_small_grid != this) ;

    const Base_val& b_i = scal_small_grid.get_spectral_base() ;

    const Mg3d& mgl = *mp->get_mg() ;
    int nz = mgl.get_nzone( );
    
    const Mg3d& mgs = *scal_small_grid.get_mp().get_mg() ;
    
    // Some grid compatibility checks ...
    //-------------------------------------
    assert(scal_small_grid.etat != ETATNONDEF) ;
    assert(mgs.get_nzone() == nz) ;
#ifndef NDEBUG
    for (int l=0; l<nz; l++) {
      assert(mgs.get_nr(l) <= mgl.get_nr(l)) ;
      assert(mgs.get_nt(l) <= mgl.get_nt(l)) ;
      assert(mgs.get_np(l) <= mgl.get_np(l)) ;
      assert(mgs.get_type_r(l) == mgl.get_type_r(l)) ;
    }
    assert(mgs.get_type_t() == mgl.get_type_t()) ;
    assert(mgs.get_type_p() == mgl.get_type_p()) ;
#endif

    if (scal_small_grid.get_etat() == ETATZERO) {     // Nothing to do...
      set_etat_zero() ;
      return ;
    }
    else {  // Copy must be done ...
      set_etat_qcq() ;
      va.set_etat_cf_qcq() ;
      scal_small_grid.get_spectral_va().coef() ;
      const Mtbl_cf& c_s = *(scal_small_grid.get_spectral_va().c_cf) ;
      assert( c_s.get_etat() == ETATQCQ ) ;
      Mtbl_cf& c_l = *(va.c_cf) ;
      c_l.set_etat_qcq() ;
      c_l.base = b_i ;

      for (int l=0; l<nz; l++) {
    if (c_s.t[l]->get_etat() == ETATZERO) {
      c_l.t[l]->set_etat_zero() ;
    }
    else {
      int nrl = mgl.get_nr(l) ;
      int nrs = mgs.get_nr(l) ;
      int ntl = mgl.get_nt(l) ;
      int nts = mgs.get_nt(l) ;
      int npl = mgl.get_np(l)+2 ;
      int nps = mgs.get_np(l)+2 ;
      c_l.t[l]->set_etat_qcq() ;
      double* p_out = c_l.t[l]->t ;
      double* p_in = c_s.t[l]->t ;
      for (int k=0; k<nps; k++) {
        for (int j=0; j<nts; j++) {
          for (int i=0; i<nrs; i++) {
        *p_out = *p_in ;
        p_out++ ; p_in++ ;
          }
          for (int i=nrs; i<nrl; i++) {
        *p_out = 0. ; p_out++ ;
          } // loop on i (r)
        } // loop on j (theta)
        for (int j=nts; j<ntl; j++) {
          for (int i=0; i<nrl; i++) {
        *p_out = 0. ; p_out++ ;
          }       
        } // loop on j (theta)
      } // loop on k (phi)
      for (int k=nps; k<npl; k++) {
        for (int j=0; j<ntl; j++) {
          for (int i=0; i<nrl; i++) {
        *p_out = 0. ; p_out++ ;
          }
        } // loop on j (theta)
      } // loop on k (phi)
    } // case etat != ETATZERO 
      } // loop on l (domains)
      if (scal_small_grid.get_etat() == ETATUN) etat = ETATUN ;
    } // case when copy must be done
    dzpuis = scal_small_grid.get_dzpuis() ;
    set_spectral_base(b_i ) ;
  }

  void Scalar::copy_coefs_from_large_grid(const Scalar& scal_large_grid) {

    assert(&scal_large_grid != this) ;

    const Base_val& b_i = scal_large_grid.get_spectral_base() ;

    const Mg3d& mgs = *mp->get_mg() ;
    int nz = mgs.get_nzone( );

    const Mg3d& mgl = *scal_large_grid.get_mp().get_mg() ;
    
    // Some grid compatibility checks ...
    //-------------------------------------
#ifndef NDEBUG
    assert(scal_large_grid.etat != ETATNONDEF) ;
    assert(mgl.get_nzone() == nz) ;    
    for (int l=0; l<nz; l++) {
      assert(mgs.get_nr(l) <= mgl.get_nr(l)) ;
      assert(mgs.get_nt(l) <= mgl.get_nt(l)) ;
      assert(mgs.get_np(l) <= mgl.get_np(l)) ;
      assert(mgs.get_type_r(l) == mgl.get_type_r(l)) ;
    }
    assert(mgs.get_type_t() == mgl.get_type_t()) ;
    assert(mgs.get_type_p() == mgl.get_type_p()) ;
#endif

    if (scal_large_grid.get_etat() == ETATZERO) {     // Nothing to do...
      set_etat_zero() ;
      return ;
    }
    else {  // Copy must be done ...
      set_etat_qcq() ;
      va.set_etat_cf_qcq() ;
      scal_large_grid.get_spectral_va().coef() ;
      const Mtbl_cf& c_l = *(scal_large_grid.get_spectral_va().c_cf) ;
      assert( c_l.get_etat() == ETATQCQ ) ;
      Mtbl_cf& c_s = *(va.c_cf) ;
      c_s.set_etat_qcq() ;
      c_s.base = b_i ;

      for (int l=0; l<nz; l++) {
    if (c_l.t[l]->get_etat() == ETATZERO) {
      c_s.t[l]->set_etat_zero() ;
    }
    else {
      int nrs = mgs.get_nr(l) ;
      int nts = mgs.get_nt(l) ;
      int nps = mgs.get_np(l)+1 ;
      int nrl = mgl.get_nr(l) ;
      int dif_r = nrl - nrs ;
      int dif_t = mgl.get_nt(l) - nts ;
      c_s.t[l]->set_etat_qcq() ;
      double* p_out = c_s.t[l]->t ;
      double* p_in = c_l.t[l]->t ;
      for (int k=0; k<nps; k++) {
        for (int j=0; j<nts; j++) {
          for (int i=0; i<nrs; i++) {
        *p_out = *p_in ;
        p_out++ ; p_in++ ;
          }
          p_in += dif_r ;
        } // loop on j (theta)
        p_in += dif_t*nrl ;
      } // loop on k (phi)
      // The last coef. in phi is set to 0
      for (int j=0; j<nts; j++) {
        for (int i=0; i<nrs; i++) {
          *p_out = 0. ; p_out++ ;
        }
      } // loop on j (theta)
    } // case etat != ETATZERO 
      } // loop on l (domains)
      if (scal_large_grid.get_etat() == ETATUN) etat = ETATUN ;
    } // case when copy must be done
    dzpuis = scal_large_grid.get_dzpuis() ;
    set_spectral_base(b_i ) ;
  }


} // End of namespace Lorene