map_af_poisson2d.C

/*
 * Method of the class Map_af for the resolution of the 2-D Poisson
 *  equation.
 *
 * (see file map.h for documentation).
 */
 
/*
 *   Copyright (c) 2000-2001 Eric Gourgoulhon
 *
 *   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
 *
 */
 
 
char map_af_poisson2d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson2d.C,v 1.6 2014/10/13 08:53:02 j_novak Exp $" ;
 
/*
 * $Id: map_af_poisson2d.C,v 1.6 2014/10/13 08:53:02 j_novak Exp $
 * $Log: map_af_poisson2d.C,v $
 * Revision 1.6  2014/10/13 08:53:02  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.5  2012/09/04 14:53:28  j_novak
 * Replacement of the FORTRAN version of huntm by a C one.
 *
 * Revision 1.4  2012/08/12 17:48:36  p_cerda
 * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
 *
 * Revision 1.3  2002/09/09 13:54:20  e_gourgoulhon
 *
 * Change of name of the Fortran subroutines
 *  poisson2d -> poiss2d
 *  poisson2di -> poiss2di
 * to avoid any clash with Map::poisson2d and Map::poisson2di
 *
 * Revision 1.2  2002/09/09 13:00:39  e_gourgoulhon
 * Modification of declaration of Fortran 77 prototypes for
 * a better portability (in particular on IBM AIX systems):
 * All Fortran subroutine names are now written F77_* and are
 * defined in the new file C++/Include/proto_f77.h.
 *
 * Revision 1.1.1.1  2001/11/20 15:19:27  e_gourgoulhon
 * LORENE
 *
 * Revision 2.1  2000/10/11  15:15:26  eric
 * 1ere version operationnelle.
 *
 * Revision 2.0  2000/10/09  13:47:10  eric
 * *** empty log message ***
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson2d.C,v 1.6 2014/10/13 08:53:02 j_novak Exp $
 *
 */
 
// Header Lorene:
#include "map.h"
#include "cmp.h"
#include "param.h"
#include "proto_f77.h"
 
//*****************************************************************************
 
namespace Lorene {
 
void Map_af::poisson2d(const Cmp& source_mat, const Cmp& source_quad, 
               Param& par, Cmp& uu) const {
    
    assert(source_mat.get_etat() != ETATNONDEF) ; 
    assert(source_quad.get_etat() != ETATNONDEF) ; 
    assert(source_mat.get_mp()->get_mg() == mg) ; 
    assert(source_quad.get_mp()->get_mg() == mg) ; 
    assert(uu.get_mp()->get_mg() == mg) ; 
 
    assert( source_quad.check_dzpuis(4) ) ; 
    
    int mpsymm = uu.get_mp()->get_mg()->get_type_t();
 
    double& lambda = par.get_double_mod(0) ; 
 
    // Special case of a vanishing source 
    // ----------------------------------
 
    if (    (source_mat.get_etat() == ETATZERO) 
     && (source_quad.get_etat() == ETATZERO) ) {
    
    uu = 0 ; 
    lambda = 1 ; 
    return ; 
    }
 
    // ---------------------------------------------------------------------
    // Preparation of the parameters for the Fortran subroutines F77_poisson2d
    //  and F77_poisson2di
    // ---------------------------------------------------------------------
    
    int nz = mg->get_nzone() ;
    int np1 = 1 ;       // Axisymmetry enforced
    int nt = mg->get_nt(0) ; 
    int nt2 = 0 ;
 
    switch ( mpsymm ){
    case SYM: {
      nt2 = 2*nt - 1 ;  // Number of points for the theta sampling
      break;                //  in [0,Pi], instead of [0,Pi/2]
    }
    case NONSYM: {
      nt2 = nt;
      break;
    }
    }
      
 
    // Array NDL
    // ---------
    int* ndl = new int[nz+4] ; 
    ndl[0] = nz ; 
    for (int l=0; l<nz; l++) {
    ndl[1+l] = mg->get_nr(l) ; 
    }
    ndl[1+nz] = nt2 ; 
    ndl[2+nz] = np1 ; 
    ndl[3+nz] = nz ; 
    
    // Array INDD
    // ----------
    int* indd = new int[nz] ; 
    for (int l=0; l<nz; l++) {
    switch ( mg->get_type_r(l) ) {
        case RARE : {
        indd[l] = 0 ; 
        break ; 
        }
        case FIN : {
        indd[l] = 1 ; 
        break ; 
        }
        case UNSURR : {
        indd[l] = 2 ; 
        break ; 
        }
        default : {
        cout << "Map_af::poisson2d: unknown type_r !" << endl ; 
        abort() ; 
        break ; 
        }
    }
    }
    
    // Parameters NDR, NDT, NDP and NDZ
    // --------------------------------
    int nrmax = 0 ;
    for (int l=0; l<nz ; l++) {
    nrmax = ( ndl[1+l] > nrmax ) ? ndl[1+l] : nrmax ; 
    }
    int ndr = nrmax + 5 ;   // Le +5 est impose par les routines de resolution
                // de l'equation de Poisson du type gr2p3s_
    int ndt = nt2 + 2 ;
    int ndp = np1 + 2 ;
    int ndz = nz ; 
    
    // Array ERRE
    // ----------
    
    double* erre = new double [ndz*ndr] ; 
    
    for (int l=0; l<nz; l++) {
    for (int i=0; i<ndl[1+l]; i++) {
        double xr = mg->get_grille3d(l)->x[i] ;
        erre[ ndr*l + i ] = alpha[l] * xr + beta[l]  ;
    }
    }
 
    // Arrays containing the data
    // --------------------------
    
    int ndrt = ndr*ndt ; 
    int ndrtp = ndr*ndt*ndp ; 
    int taille = ndrtp*ndz ;
 
    double* tsou_m = new double[ taille ] ;    
    double* tsou_q = new double[ taille ] ;    
    double* tuu = new double[ taille ] ;    
    
    // Initialisation to zero :
    for (int i=0; i<taille; i++) {
    tsou_m[i] = 0 ; 
    tsou_q[i] = 0 ; 
    tuu[i] = 0 ; 
    }
    
    // Copy of source_mat into tsou_m
    // ------------------------------
    const Valeur& va_m = source_mat.va ; 
    assert(va_m.get_etat() == ETATQCQ) ; 
    va_m.coef_i() ; 
    const Mtbl* s_m = va_m.c ; 
    assert(s_m->get_etat() == ETATQCQ) ; 
    
    Base_val base_s = va_m.base ; 
 
    for (int l=0; l<nz; l++) {
    int nr = mg->get_nr(l) ; 
    int nrt = nr*nt ;  
    if (s_m->t[l]->get_etat() == ETATZERO) {
        for (int k=0; k<np1; k++) {
        for (int j=0; j<nt; j++) {
            for (int i=0; i<nr; i++) {
            tsou_m[ndrtp*l + ndrt*k + ndr*j + i] = 0 ;
            // point symetrique par rapport au plan theta = pi/2 :
            if ( mpsymm == SYM ) tsou_m[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = 0 ;         
            }
        }
        }
        
    }
    else {
        assert( s_m->t[l]->get_etat() == ETATQCQ ) ; 
        for (int k=0; k<np1; k++) {
        for (int j=0; j<nt; j++) {
            for (int i=0; i<nr; i++) {
            double xx = s_m->t[l]->t[nrt*k + nr*j + i] ;
            tsou_m[ndrtp*l + ndrt*k + ndr*j + i] = xx ;
            // point symetrique par rapport au plan theta = pi/2 :
            if ( mpsymm == SYM ) tsou_m[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = xx ;            
            }
        }
        }
    }  // End of case etat != ETATZERO
    }  
 
    // Copy of source_quad into tsou_q
    // -------------------------------
 
    if (source_quad.get_etat() != ETATZERO) {
 
    const Valeur& va_q = source_quad.va ; 
    assert(va_q.get_etat() == ETATQCQ) ; 
    va_q.coef_i() ; 
    const Mtbl* s_q = va_q.c ;  
 
    assert( va_q.base == base_s ) ; 
 
    assert(s_q->get_etat() == ETATQCQ) ; 
 
    for (int l=0; l<nz; l++) {
        int nr = mg->get_nr(l) ; 
        int nrt = nr*nt ;  
        if (s_q->t[l]->get_etat() == ETATZERO) {
        for (int k=0; k<np1; k++) {
            for (int j=0; j<nt; j++) {
            for (int i=0; i<nr; i++) {
            tsou_q[ndrtp*l + ndrt*k + ndr*j + i] = 0 ;
            // point symetrique par rapport au plan theta = pi/2 :
            if ( mpsymm == SYM ) tsou_q[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = 0 ;         
            }
            }
        }
        
        }
        else {
        assert( s_q->t[l]->get_etat() == ETATQCQ ) ; 
        for (int k=0; k<np1; k++) {
            for (int j=0; j<nt; j++) {
            for (int i=0; i<nr; i++) {
            double xx = s_q->t[l]->t[nrt*k + nr*j + i] ;
            tsou_q[ndrtp*l + ndrt*k + ndr*j + i] = xx ;
            // point symetrique par rapport au plan theta = pi/2 :
            if ( mpsymm == SYM ) tsou_q[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = xx ;            
            }
            }
        }
        }  // End of case s_q->t[l]->etat != ETATZERO
    }
    
    } // End of case source_quad.etat != ETATZERO  
 
    //-----------------------------------------------------------
    //  Call of the Fortran subroutine poisson2d_ or poisson2di_
    //-----------------------------------------------------------
 
    int base_t = (va_m.base).get_base_t(0) ; 
 
    Base_val base_uu(nz) ;  // Output spectral bases
 
    switch (base_t) {
 
        case T_COS :
    case T_COS_P : {
    
        double lambda0 ; 
        
        F77_poiss2d(ndl, &ndr, &ndt, &ndp, indd, erre, tsou_m, tsou_q, 
               &lambda0, tuu) ;
              
        base_uu = base_s ;      // output bases = input bases
               
        lambda = lambda0 ; 
        break ;         
    }
        case T_SIN :
    case T_SIN_I : {
 
        double* tsou = new double[taille] ; 
        for (int i=0; i<taille; i++) {
        tsou[i] = tsou_m[i] + tsou_q[i] ; 
        }
 
        F77_poiss2di(ndl, &ndr, &ndt, &ndp, indd, erre, tsou, tuu) ;
        
        base_uu = base_s ;      // output bases = input bases
 
        lambda = double(1) ; 
        
        delete [] tsou ; 
        break ;         
    }
    
    default : {
        cout << "Map_af::poisson2d : unkown theta basis !" << endl ;
        cout << "  basis : " << hex << base_t << endl ; 
        abort() ; 
        break ; 
    }
    }
    
    //-------------------------------
    // Copy of tuu into uu
    //-------------------------------
 
    uu.allocate_all() ; 
    (uu.va).set_etat_c_qcq() ;  // To make sure that the coefficients are
                //  deleted
    
    for (int l=0; l<nz; l++) {
    int nr = mg->get_nr(l) ; 
    for (int k=0; k<mg->get_np(l); k++) {
        for (int j=0; j<nt; j++) {
        for (int i=0; i<nr; i++) {
            uu.set(l, k, j, i) = tuu[ndrtp*l + ndr*j + i]  ;
        }
        }
    }
    }
 
    (uu.va).set_base( base_uu ) ;   // Bases for spectral expansions
    
    uu.set_dzpuis(0) ; 
   
    // Cleaning
    // --------
    
    delete [] ndl ; 
    delete [] indd ; 
    delete [] erre ; 
    delete [] tsou_m ; 
    delete [] tsou_q ; 
    delete [] tuu ; 
    
    
}
 
 
}