vector.C

/*
 *  Methods of class Vector
 *
 *   (see file vector.h for documentation)
 *
 */
 
/*
 *   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 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: vector.C,v 1.31 2016/12/05 16:18:18 j_novak Exp $
 * $Log: vector.C,v $
 * Revision 1.31  2016/12/05 16:18:18  j_novak
 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
 *
 * Revision 1.30  2014/10/13 08:53:44  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.29  2014/10/06 15:13:20  j_novak
 * Modified #include directives to use c++ syntax.
 *
 * Revision 1.28  2008/10/29 14:09:14  jl_cornou
 * Spectral bases for pseudo vectors and curl added
 *
 * Revision 1.27  2008/08/27 08:52:23  jl_cornou
 * Added fonctions for angular potential A
 *
 * Revision 1.26  2007/12/21 16:07:08  j_novak
 * Methods to filter Tensor, Vector and Sym_tensor objects.
 *
 * Revision 1.25  2005/02/14 13:01:50  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.24  2005/01/25 15:37:35  j_novak
 * Solved some dzpuis problem...
 *
 * Revision 1.23  2005/01/12 16:48:23  j_novak
 * Better treatment of the case where all vector components are null in
 * decompose_div .
 *
 * Revision 1.22  2004/10/12 09:58:25  j_novak
 * Better memory management.
 *
 * Revision 1.21  2004/10/11 09:46:31  j_novak
 * Speed improvements.
 *
 * Revision 1.20  2004/05/09 20:55:05  e_gourgoulhon
 * Added method flux.
 *
 * Revision 1.19  2004/03/29 11:57:45  e_gourgoulhon
 * Added methods ope_killing and ope_killing_conf.
 *
 * Revision 1.18  2004/02/26 22:48:50  e_gourgoulhon
 * -- Method divergence: call to Tensor::divergence and cast of the
 *    result.
 * -- Added method derive_lie.
 *
 * Revision 1.17  2004/02/24 09:46:20  j_novak
 * Correction to cope with SGI compiler's warnings.
 *
 * Revision 1.16  2004/02/20 10:53:04  j_novak
 * Added the matching of the potential adapted to the case where the
 * vector is the source of a Poisson equation (dzpuis = 4).
 *
 * Revision 1.15  2004/01/30 10:30:17  j_novak
 * Changed dzpuis handling in Vector::decompose_div (this may be temporary).
 *
 * Revision 1.14  2003/12/30 23:09:47  e_gourgoulhon
 * Change in methods derive_cov() and divergence() to take into account
 *  the change of name: Metric::get_connect() --> Metric::connect().
 *
 * Revision 1.13  2003/12/19 15:18:16  j_novak
 * Shadow variables hunt
 *
 * Revision 1.12  2003/10/29 11:04:34  e_gourgoulhon
 * dec2_dpzuis() replaced by dec_dzpuis(2).
 * inc2_dpzuis() replaced by inc_dzpuis(2).
 *
 * Revision 1.11  2003/10/22 14:24:19  j_novak
 * *** empty log message ***
 *
 * Revision 1.9  2003/10/20 13:00:38  j_novak
 * Memory error corrected
 *
 * Revision 1.8  2003/10/20 09:32:12  j_novak
 * Members p_potential and p_div_free of the Helmholtz decomposition
 * + the method decompose_div(Metric).
 *
 * Revision 1.7  2003/10/16 14:21:37  j_novak
 * The calculation of the divergence of a Tensor is now possible.
 *
 * Revision 1.6  2003/10/13 13:52:40  j_novak
 * Better managment of derived quantities.
 *
 * Revision 1.5  2003/10/06 13:58:48  j_novak
 * The memory management has been improved.
 * Implementation of the covariant derivative with respect to the exact Tensor
 * type.
 *
 * Revision 1.4  2003/10/05 21:14:20  e_gourgoulhon
 * Added method std_spectral_base().
 *
 * Revision 1.3  2003/10/03 14:10:32  e_gourgoulhon
 * Added constructor from Tensor.
 *
 * Revision 1.2  2003/10/03 14:08:46  j_novak
 * Removed old change_trid...
 *
 * Revision 1.1  2003/09/26 08:05:31  j_novak
 * New class Vector.
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector.C,v 1.31 2016/12/05 16:18:18 j_novak Exp $
 *
 */
 
// Headers C
#include <cstdlib>
#include <cassert>
#include <cmath>
 
// Headers Lorene
#include "metric.h"
#include "proto.h"
#include "matrice.h"
#include "nbr_spx.h"
 
 
            //--------------//
            // Constructors //
            //--------------//
 
// Standard constructor 
// --------------------
namespace Lorene {
Vector::Vector(const Map& map, int tipe, const Base_vect& triad_i) 
        : Tensor(map, 1, tipe, triad_i) {
        
  set_der_0x0() ;
 
}
 
// Standard constructor with the triad passed as a pointer
// -------------------------------------------------------
Vector::Vector(const Map& map, int tipe, const Base_vect* triad_i) 
        : Tensor(map, 1, tipe, *triad_i) {
        
  set_der_0x0() ;
}
    
// Copy constructor
// ----------------
Vector::Vector (const Vector& source) : 
    Tensor(source) {
  
  assert(valence == 1) ;
  set_der_0x0() ;
 
}   
 
 
// Constructor from a {\tt Tensor}.
//--------------------------------
Vector::Vector(const Tensor& uu) : Tensor(uu) {
 
  assert(valence == 1) ;
  set_der_0x0() ;
 
}
 
 
// Constructor from a file
// -----------------------
Vector::Vector(const Map& mapping, const Base_vect& triad_i, FILE* fd) : 
  Tensor(mapping, triad_i, fd) {
   
  assert ( (valence == 1) && (n_comp == 3) ) ;
  set_der_0x0() ;
 
}
 
 
            //--------------//
            //  Destructor  //
            //--------------//
 
 
Vector::~Vector () {
 
  Vector::del_deriv() ;
 
}
 
 
            //-------------------//
            // Memory managment  //
            //-------------------//
 
void Vector::del_deriv() const {
 
  for (int i=0; i<N_MET_MAX; i++) 
    del_derive_met(i) ;
  
  if (p_A != 0x0) delete p_A ;
  if (p_eta != 0x0) delete p_eta ; 
  if (p_mu != 0x0) delete p_mu ; 
  set_der_0x0() ;
  Tensor::del_deriv() ;
 
}
 
void Vector::set_der_0x0() const {
 
  for (int i=0; i<N_MET_MAX; i++) 
    set_der_met_0x0(i) ;
  p_A = 0x0 ;
  p_eta = 0x0 ; 
  p_mu = 0x0 ; 
 
}
 
void Vector::del_derive_met(int j) const {
 
  assert( (j>=0) && (j<N_MET_MAX) ) ;
 
  if (met_depend[j] != 0x0) {
    if (p_potential[j] != 0x0)
      delete p_potential[j] ;
    if (p_div_free[j] != 0x0)
      delete p_div_free[j] ;
    
    set_der_met_0x0(j) ;
    
    Tensor::del_derive_met(j) ;
  }
}
 
void Vector::set_der_met_0x0(int i) const {
 
  assert( (i>=0) && (i<N_MET_MAX) ) ;
 
  p_potential[i] = 0x0 ;
  p_div_free[i] = 0x0 ;
 
}
 
void Vector::operator=(const Vector& t) {
    
    triad = t.triad ; 
 
    assert(t.type_indice(0) == type_indice(0)) ;
 
    for (int i=0 ; i<3 ; i++) {
      *cmp[i] = *t.cmp[i] ;
    }
    del_deriv() ;
}
 
void Vector::operator=(const Tensor& t) {
    
    assert (t.valence == 1) ;
 
    triad = t.triad ; 
 
    assert(t.type_indice(0) == type_indice(0)) ;
 
    for (int i=0 ; i<3 ; i++) {
      *cmp[i] = *t.cmp[i] ;
    }
    del_deriv() ;
}
 
 
 
// Affectation d'une composante :
Scalar& Vector::set(int index) {
    
  assert ( (index>=1) && (index<=3) ) ;
 
  del_deriv() ;
  
  return *cmp[index - 1] ;
}
 
const Scalar& Vector::operator()(int index) const {
    
  assert ( (index>=1) && (index<=3) ) ;
 
  return *cmp[index - 1] ;
 
}
 
 
// Sets the standard spectal bases of decomposition for each component
 
void Vector::std_spectral_base() {
 
    Base_val** bases = 0x0 ;
 
    if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
 
        // Cartesian case
        bases = mp->get_mg()->std_base_vect_cart() ;
 
    }
    else {
        // Spherical case
        assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
        bases = mp->get_mg()->std_base_vect_spher() ;
    }
        
    for (int i=0 ; i<3 ; i++) {
        cmp[i]->set_spectral_base( *bases[i] ) ; 
    }
        
    for (int i=0 ; i<3 ; i++) {
        delete bases[i] ;
    }
    delete [] bases ;
 
 
}
 
// Sets the standard spectral bases of decomposition for each component for a pseudo vector
 
void Vector::pseudo_spectral_base() {
 
    Base_val** bases = 0x0 ;
 
    if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
 
        // Cartesian case
        bases = mp->get_mg()->pseudo_base_vect_cart() ;
 
    }
    else {
        // Spherical case
        assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
        bases = mp->get_mg()->pseudo_base_vect_spher() ;
    }
        
    for (int i=0 ; i<3 ; i++) {
        cmp[i]->set_spectral_base( *bases[i] ) ; 
    }
        
    for (int i=0 ; i<3 ; i++) {
        delete bases[i] ;
    }
    delete [] bases ;
 
 
}
 
 
 
                    //-------------------------------//
                    //    Computational methods      //
                    //-------------------------------//
                    
 
const Scalar& Vector::divergence(const Metric& metre) const {
  
    const Scalar* pscal = 
        dynamic_cast<const Scalar*>( &(Tensor::divergence(metre)) ) ;
 
    assert(pscal != 0x0) ;
 
    return *pscal ;
}
 
 
Vector Vector::derive_lie(const Vector& vv) const {
 
    Vector resu(*mp, type_indice(0), triad) ; 
    
    compute_derive_lie(vv, resu) ;
    
    return resu ; 
    
}
 
const Vector_divfree Vector::curl() const {
 
    if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
    const Metric_flat& metc = mp->flat_met_cart() ;
    Vector_divfree resu(*mp, mp->get_bvect_cart(), metc) ;
    resu.set(1)= cmp[3]->dsdy() - cmp[2]->dsdz();
    resu.set(2)= cmp[1]->dsdz() - cmp[3]->dsdx();
    resu.set(3)= cmp[2]->dsdx() - cmp[1]->dsdy();
    resu.pseudo_spectral_base();
    return resu ;
    }
    else    {
        assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
        const Metric_flat& mets = mp->flat_met_spher() ;
        Vector_divfree resu(*mp, mp->get_bvect_spher(), mets);
        Scalar tmp = *cmp[3] ;
        tmp.div_tant();
        tmp += cmp[3]->dsdt();
        tmp.div_r();
        resu.set(1) = tmp - cmp[2]->srstdsdp() ;
        tmp = *cmp[3] ;
        tmp.mult_r();
        tmp = tmp.dsdr();
        tmp.div_r();
        resu.set(2) = cmp[1]->srstdsdp() - tmp ;
        tmp = *cmp[2];
        tmp.mult_r();   
        resu.set(3) = tmp.dsdr() - cmp[1]->dsdt() ;
        resu.set(3).div_r(); 
        resu.pseudo_spectral_base();
        return resu ;
        }
 
}
 
 
Sym_tensor Vector::ope_killing(const Metric& gam) const {
 
    Sym_tensor resu(*mp, type_indice(0), *triad) ; 
 
    if (type_indice(0) == CON ) {
        for (int i=1; i<=3; i++) {
            for (int j=i; j<=3; j++) {
                resu.set(i,j) = derive_con(gam)(i,j) + derive_con(gam)(j,i)  ;
            }
        }
    }
    else {
        for (int i=1; i<=3; i++) {
            for (int j=i; j<=3; j++) {
                resu.set(i,j) = derive_cov(gam)(i,j) + derive_cov(gam)(j,i)  ;
            }
        }
    }
    
    return resu ; 
 
} 
 
 
Sym_tensor Vector::ope_killing_conf(const Metric& gam) const {
 
    Sym_tensor resu(*mp, type_indice(0), *triad) ; 
 
    if (type_indice(0) == CON ) {
        for (int i=1; i<=3; i++) {
            for (int j=i; j<=3; j++) {
                resu.set(i,j) = derive_con(gam)(i,j) + derive_con(gam)(j,i)  
                                - 0.6666666666666666* divergence(gam) 
                                            * gam.con()(i,j) ;
            }
        }
    }
    else {
        for (int i=1; i<=3; i++) {
            for (int j=i; j<=3; j++) {
                resu.set(i,j) = derive_cov(gam)(i,j) + derive_cov(gam)(j,i)  
                                - 0.6666666666666666* derive_con(gam).trace() 
                                            * gam.cov()(i,j) ;
            }
        }
    }
 
    return resu ; 
 
} 
 
 
 
 
const Scalar& Vector::potential(const Metric& metre) const {
 
  set_dependance(metre) ;
  int j = get_place_met(metre) ;
  assert ((j>=0) && (j<N_MET_MAX)) ;
  if (p_potential[j] == 0x0) {
    decompose_div(metre) ;
  }
 
  return *p_potential[j] ;
}
 
const Vector_divfree& Vector::div_free(const Metric& metre) const {
 
  set_dependance(metre) ;
  int j = get_place_met(metre) ;
  assert ((j>=0) && (j<N_MET_MAX)) ;
  if (p_div_free[j] == 0x0) {
    decompose_div(metre) ;
  }
 
  return *p_div_free[j] ;
}
 
void Vector::decompose_div(const Metric& metre) const {
 
  assert( type_indice(0) == CON ) ; //Only for contravariant vectors...
 
  set_dependance(metre) ;
  int ind =  get_place_met(metre) ;
  assert ((ind>=0) && (ind<N_MET_MAX)) ;
 
  if ( (p_potential[ind] != 0x0) && (p_div_free[ind] != 0x0) ) 
    return ; // Nothing to do ...
 
  else {
    if (p_div_free[ind] != 0x0)
      delete p_div_free[ind] ;  
    if (p_potential[ind] != 0x0) 
      delete p_potential[ind] ;
 
    const Mg3d* mmg = mp->get_mg() ;
    int nz = mmg->get_nzone() ;
 
    int dzp = cmp[0]->get_dzpuis() ;
#ifndef NDEBUG
    bool dz_zero = cmp[0]->check_dzpuis(0) ;
    bool dz_four = cmp[0]->check_dzpuis(4) ;
#endif
    assert( dz_zero || dz_four) ;
    assert (cmp[1]->check_dzpuis(dzp)) ;
    assert (cmp[2]->check_dzpuis(dzp)) ;
 
    Scalar dive = divergence(metre) ;
 
    if (dive.get_etat() == ETATZERO) {
    p_potential[ind] = new Scalar(*mp) ;
    p_potential[ind]->set_etat_zero() ;
    p_potential[ind]->set_dzpuis(dzp) ;
    }
    else {
        //No matching of the solution at this point
    p_potential[ind] = new Scalar(dive.poisson()) ; 
 
    if (dzp == 4) {
        assert (mmg->get_type_r(nz-1) == UNSURR) ;
        // Let's now do the matching ... in the case dzpuis = 4
        const Map_af* mapping = dynamic_cast<const Map_af*>(mp) ;
        assert (mapping != 0x0) ;
        Valeur& val = p_potential[ind]->set_spectral_va() ;
        val.ylm() ;
        Mtbl_cf& mtc = *val.c_cf ;
        if (nz > 1) {
        int np = mmg->get_np(0) ;
        int nt = mmg->get_nt(0) ;
#ifndef NDEBUG
        for (int lz=0; lz<nz; lz++) {
            assert (mmg->get_np(lz) == np) ;
            assert (mmg->get_nt(lz) == nt) ;
        }
#endif
        int lmax = 0 ;
        for (int k=0; k<np+1; k++) 
            for (int j=0; j<nt; j++) 
            if ( nullite_plm(j, nt, k, np, val.base)) {
                int m_quant, l_quant, base_r ;
                donne_lm (nz, 0, j, k, val.base, m_quant, 
                      l_quant, base_r) ; 
                lmax = (l_quant > lmax ? l_quant : lmax) ;
            }
        Scalar** ri = new Scalar*[lmax+1] ;
        Scalar** rmi = new Scalar*[lmax+1] ;
        Scalar erre(*mp) ;
        erre = mp->r ;
        for (int l=0; l<=lmax; l++) {
            ri[l] = new Scalar(*mp) ;
            rmi[l] = new Scalar(*mp) ;
            if (l == 0) *(ri[l]) = 1. ;
            else *(ri[l]) = pow(erre, l) ;
            ri[l]->annule_domain(nz - 1) ;
            ri[l]->std_spectral_base() ; //exact base in r will be set later
            if (l==2) *(rmi[l]) = 1 ;
            else *(rmi[l]) = pow(erre, 2-l) ;
            rmi[l]->annule(0,nz-2) ;
            Scalar tmp = pow(erre, -l-1) ;
            tmp.annule_domain(nz-1) ;
            tmp.annule_domain(0) ;
            *(rmi[l]) += tmp ;
            rmi[l]->set_dzpuis(3) ;
            rmi[l]->std_spectral_base() ;//exact base in r will be set later
        }
 
        for (int k=0; k<np+1; k++) {
            for (int j=0; j<nt; j++) {
            if ( nullite_plm(j, nt, k, np, val.base)) {
                int m_quant, l_quant, base_r, l, c ;
                donne_lm (nz, 0, j, k, val.base, m_quant, l_quant, base_r) ; 
                bool lzero = (l_quant == 0) || (l_quant == 1) ;
                int nb_hom_ced  = (l_quant < 2 ? 0 : 1) ; 
          
                int taille = 2*(nz-1) - 1 + nb_hom_ced ;
                Tbl deuz(taille) ;
                deuz.set_etat_qcq() ;
                Matrice systeme(taille,taille) ;
                systeme.set_etat_qcq() ;
                for (l=0; l<taille; l++) 
                for (c=0; c<taille; c++) systeme.set(l,c) = 0 ;
                for (l=0; l<taille; l++) deuz.set(l) = 0 ;
          
                //---------
                // Nucleus
                //---------
                assert(mmg->get_type_r(0) == RARE) ;
                assert ((base_r == R_CHEBP)||(base_r == R_CHEBI)) ;
                ri[l_quant]->set_spectral_va().set_base_r(0, base_r) ;
      
                double alpha = mapping->get_alpha()[0] ;
                int nr = mmg->get_nr(0) ;
                Tbl partn(nr) ;
                partn.set_etat_qcq() ;
                for (int i=0; i<nr; i++)
                partn.set(i) = mtc(0, k, j, i) ;
                l=0 ; c=0 ;
 
                systeme.set(l,c) += pow(alpha, double(l_quant)) ;
 
                deuz.set(l) -= val1_dern_1d(0, partn, base_r) ;
                l++ ;
 
                if (!lzero) {
                systeme.set(l,c) += l_quant * pow(alpha, double(l_quant - 1)) ;
   
                deuz.set(l) -= val1_dern_1d(1, partn, base_r) / alpha ;
                }
 
                //----------
                //  Shells
                //----------
 
                for (int lz=1; lz<nz-1; lz++) {
                alpha = mapping->get_alpha()[lz] ;
                double beta = mapping->get_beta()[lz] ;
                double xm1 = beta - alpha ;
                double xp1 = alpha + beta ;
                nr = mmg->get_nr(lz) ;
                Tbl parts(nr) ;
                parts.set_etat_qcq() ;
                for (int i=0; i<nr; i++)
                    parts.set(i) = mtc(lz, k, j, i) ;
  
                //Function at x = -1
                l-- ; 
                c = 2*lz - 1 ;
                systeme.set(l,c) -= pow(xm1, double(l_quant)) ;
                c++;
                systeme.set(l,c) -= pow(xm1, double(-l_quant-1)) ;
          
                deuz.set(l) += valm1_dern_1d(0, parts, R_CHEB) ;
        
                if ((lz>1) || (!lzero)) {
                    //First derivative at x=-1
                    l++ ;
                    c-- ;
                    systeme.set(l,c) -=  l_quant*pow(xm1, double(l_quant - 1)) ;
                    c++;
                    systeme.set(l,c) -= (-l_quant - 1)*
                    pow(xm1, double(-l_quant - 2)) ;
 
                    deuz.set(l) += valm1_dern_1d(1, parts, R_CHEB) / alpha ;
                }
 
                //Function at x = 1
                l++ ; 
                c-- ;
                systeme.set(l,c) += pow(xp1, double(l_quant)) ;
                c++;
                systeme.set(l,c) += pow(xp1, double(-l_quant-1)) ;
 
                deuz.set(l) -= val1_dern_1d(0, parts, R_CHEB) ;
        
                //First derivative at x = 1
                l++ ; 
                c-- ;
                systeme.set(l,c) +=  l_quant*pow(xp1, double(l_quant - 1)) ;
                c++;
                systeme.set(l,c) += (-l_quant - 1)*
                    pow(xp1, double(-l_quant - 2)) ;
        
                deuz.set(l) -= val1_dern_1d(1, parts, R_CHEB) / alpha ;
        
                } //End of the loop on different shells
 
                //-------------------------------
                //  Compactified external domain
                //-------------------------------
                assert(mmg->get_type_r(nz-1) == UNSURR) ;
                nr = mmg->get_nr(nz-1) ;
                Tbl partc(nr) ;
                partc.set_etat_qcq() ;
                for (int i=0; i<nr; i++)
                partc.set(i) = mtc(nz-1, k, j, i) ;
 
                alpha = mapping->get_alpha()[nz-1] ;
                double beta =  mapping->get_beta()[nz-1] ;
                double xm1 = beta - alpha ; // 1 / r_left
                double part0, part1 ;
                part0 = valm1_dern_1d(0, partc, R_CHEB) ;
                part1 = xm1*valm1_dern_1d(1, partc, R_CHEB) / alpha ;
                assert (p_potential[ind]->get_dzpuis() == 3) ;
        
                //Function at x = -1
                l--;
                if (!lzero) {
                c++;
                systeme.set(l,c) -= pow(xm1, double(l_quant+1)) ;
                }
                deuz.set(l) += part0*xm1*xm1*xm1 ;
 
                // First derivative at x = -1
                l++ ;
                if (!lzero) {
                systeme.set(l,c) -= (-l_quant - 1)*
                    pow(xm1, double(l_quant + 2)) ;
                }
                if ( (nz > 2) || (!lzero))
                deuz.set(l) += -xm1*xm1*xm1*xm1*(3*part0 + part1) ;
 
                //--------------------------------------
                //   Solution of the linear system
                //--------------------------------------
      
                int inf = 1 + nb_hom_ced;
                int sup = 3 - nb_hom_ced;
                systeme.set_band(sup, inf) ;
                systeme.set_lu() ;
                Tbl facteur(systeme.inverse(deuz)) ;
                ri[l_quant]->set_spectral_va().coef() ;
                rmi[l_quant]->set_spectral_va().coef() ;
 
                //Linear combination in the nucleus
                nr = mmg->get_nr(0) ;
                for (int i=0; i<nr; i++) 
                mtc.set(0, k, j, i) += 
         facteur(0)*(*(ri[l_quant]->get_spectral_va().c_cf))(0, 0, 0, i) ;
        
                //Linear combination in the shells
                for (int lz=1; lz<nz-1; lz++) {
                nr = mmg->get_nr(lz) ;
                for (int i=0; i<nr; i++) 
                    mtc.set(lz, k, j, i) += facteur(2*lz - 1)*
           (*(ri[l_quant]->get_spectral_va().c_cf))(lz, 0, 0, i) ;
                for (int i=0; i<nr; i++) 
                    mtc.set(lz, k, j, i) += facteur(2*lz)*
           (*(rmi[l_quant]->get_spectral_va().c_cf))(lz, 0, 0, i) ;
                }
 
                //Linear combination in the CED
                nr = mmg->get_nr(nz-1) ;
                if (!lzero) {
                for (int i=0; i<nr; i++) 
                    mtc.set(nz-1, k, j, i) += 
                    facteur(taille - 1)*
          (*(rmi[l_quant]->get_spectral_va().c_cf))(nz-1, 0, 0, i) ;
                }
            } //End of nullite_plm ...
 
            } //End of j/theta loop   
        } //End of k/phi loop 
 
        for (int l=0; l<=lmax; l++) {
            delete ri[l] ;
            delete rmi[l] ;
        }
        delete [] ri ;
        delete [] rmi ;
 
        } //End of the case of more than one domain
 
        val.ylm_i() ;
 
    } //End of the case dzp = 4
    }
    p_div_free[ind] = new Vector_divfree(*mp, *triad, metre) ;
 
    Vector gradient = p_potential[ind]->derive_con(metre) ;
    if (dzp != 4)  gradient.dec_dzpuis(2) ;
 
    *p_div_free[ind] = ( *this - gradient ) ;
 
  }
  
}
 
 
 
double Vector::flux(double radius, const Metric& met) const {
 
    assert(type_indice(0) == CON) ; 
    
    const Map_af* mp_af = dynamic_cast<const Map_af*>(mp) ; 
    if (mp_af == 0x0) {
        cerr << 
        "Vector::flux : the case of a mapping outside the class Map_af\n"
            << " is not implemented yet !" << endl ; 
        abort() ; 
    } 
    
    const Metric_flat* ff = dynamic_cast<const Metric_flat*>(&met) ;
    if (ff == 0x0) {
        cerr << 
        "Vector::flux : the case of a non flat metric is not implemented yet !"
             << endl ; 
        abort() ; 
    }
    
    const Base_vect_cart* bcart = dynamic_cast<const Base_vect_cart*>(triad) ; 
    Vector* vspher = 0x0 ; 
    if (bcart != 0x0) { // switch to spherical components:
        vspher = new Vector(*this) ; 
        vspher->change_triad(mp->get_bvect_spher()) ;
    }
    else assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ; 
 
    const Vector* vv = (bcart != 0x0) ? vspher : this ; 
 
    const Scalar& vr = vv->operator()(1) ; 
    
    double resu ; 
    if (radius == __infinity) {
        resu = mp_af->integrale_surface_infini(vr) ; 
    }
    else {
        resu = mp_af->integrale_surface(vr, radius) ;
    }
    
    return resu ; 
}
 
void Vector::exponential_filter_r(int lzmin, int lzmax, int p, 
              double alpha) {
    if( triad->identify() == (mp->get_bvect_cart()).identify() ) 
    for (int i=0; i<n_comp; i++)
        cmp[i]->exponential_filter_r(lzmin, lzmax, p, alpha) ;
    else {
    assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
    Scalar vr_tmp = operator()(1) ; 
    vr_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
    Scalar eta_tmp = eta() ;
    eta_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
    Scalar mu_tmp = mu() ;
    mu_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
    set_vr_eta_mu(vr_tmp, eta_tmp, mu_tmp) ;
    }
}
 
void Vector::exponential_filter_ylm(int lzmin, int lzmax, int p, 
              double alpha) {
    if( triad->identify() == (mp->get_bvect_cart()).identify() ) 
    for (int i=0; i<n_comp; i++)
        cmp[i]->exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
    else {
    assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
    Scalar vr_tmp = operator()(1) ; 
    vr_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
    Scalar eta_tmp = eta() ;
    eta_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
    Scalar mu_tmp = mu() ;
    mu_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
    set_vr_eta_mu(vr_tmp, eta_tmp, mu_tmp) ;
    }
}
}