mg3d.C

/*
 * Methods of class Mg3d
 *
 */
 
/*
 *   Copyright (c) 1999-2000 Jean-Alain Marck
 *   Copyright (c) 1999-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
 *
 */
 
 
 
 
/*
 * $Id: mg3d.C,v 1.22 2023/05/24 09:48:50 g_servignat Exp $
 * $Log: mg3d.C,v $
 * Revision 1.22  2023/05/24 09:48:50  g_servignat
 * Added plus_half grid in angular direction for dealiasing
 *
 * Revision 1.21  2020/01/27 11:00:19  j_novak
 * New include <stdexcept> to be compatible with older versions of g++
 *
 * Revision 1.20  2018/12/05 15:03:20  j_novak
 * New Mg3d constructor from a formatted file.
 *
 * Revision 1.19  2016/12/05 16:17:59  j_novak
 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
 *
 * Revision 1.18  2014/10/13 08:53:07  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.17  2014/10/06 15:13:14  j_novak
 * Modified #include directives to use c++ syntax.
 *
 * Revision 1.16  2013/06/05 15:00:26  j_novak
 * Suppression of all classes derived from Grille3d. Now Grille3d is no
 * longer an abstract class. r-samplings are only one of RARE, FIN or
 * UNSURR (FINJAC has been removed). Instead, Mg3d possesses a new member
 * colloc_r[nzone] defining the type of collocation points (spectral
 * bases) in each domain.
 *
 * Revision 1.15  2012/01/17 10:37:42  j_penner
 * added a constructor that treats all domains as type FIN
 *
 * Revision 1.14  2008/08/19 06:42:00  j_novak
 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
 * cast-type operations, and constant strings that must be defined as const char*
 *
 * Revision 1.13  2007/12/11 15:28:15  jl_cornou
 * Jacobi(0,2) polynomials partially implemented
 *
 * Revision 1.12  2006/05/17 13:17:03  j_novak
 * New member g_angu_1dom, the one-domain angular grid associated with the
 * current grid.
 *
 * Revision 1.11  2005/10/07 08:47:21  j_novak
 * Addition of the pointer g_non_axi on a grid, with at least 5 points in the
 * theta direction and 4 in the phi one (for tensor rotations).
 *
 * Revision 1.10  2004/07/06 13:36:29  j_novak
 * Added methods for desaliased product (operator |) only in r direction.
 *
 * Revision 1.9  2003/10/27 16:21:54  e_gourgoulhon
 * Treated the special case nz=1 in the simplified constructor.
 *
 * Revision 1.8  2003/06/20 14:50:15  f_limousin
 * Add the operator==
 *
 * Revision 1.7  2003/06/18 08:45:27  j_novak
 * In class Mg3d: added the member get_radial, returning only a radial grid
 * For dAlembert solver: the way the coefficients of the operator are defined has been changed.
 *
 * Revision 1.6  2002/10/16 14:36:42  j_novak
 * Reorganization of #include instructions of standard C++, in order to
 * use experimental version 3 of gcc.
 *
 * Revision 1.5  2002/05/07 07:36:03  e_gourgoulhon
 * Compatibilty with xlC compiler on IBM SP2:
 *    suppressed the parentheses around argument of instruction new:
 *  e.g.   t = new (Tbl *[nzone])  -->   t = new Tbl*[nzone]
 *
 * Revision 1.4  2001/12/12 09:23:46  e_gourgoulhon
 * Parameter compact added to the simplified constructor of class Mg3d
 *
 * Revision 1.3  2001/12/11 06:48:30  e_gourgoulhon
 * Addition of the simplified constructor
 *
 * Revision 1.2  2001/12/04 21:27:54  e_gourgoulhon
 *
 * All writing/reading to a binary file are now performed according to
 * the big endian convention, whatever the system is big endian or
 * small endian, thanks to the functions fwrite_be and fread_be
 *
 * Revision 1.1.1.1  2001/11/20 15:19:27  e_gourgoulhon
 * LORENE
 *
 * Revision 2.10  2001/05/26  14:50:46  eric
 * *** empty log message ***
 *
 * Revision 2.9  2001/05/26  13:25:59  eric
 * Ajout du membre g_twice (grille double pour le desaliasing)
 * Modif de la declaration de g_angu (pointeur mutable)
 *   g_twice et g_angu ne sont calcules que si necessaire (cad si
 *   on appelle la fonction get_twice() ou get_angu()).
 *
 * Revision 2.8  2000/03/22  13:38:51  eric
 * Remplacement des iendl par endl dans <<
 *
 * Revision 2.7  1999/10/12  15:04:29  eric
 * *** empty log message ***
 *
 * Revision 2.6  1999/10/12  15:03:30  eric
 * *** empty log message ***
 *
 * Revision 2.5  1999/09/30  14:58:16  eric
 * Operator!= declare const
 *
 * Revision 2.4  1999/09/30  14:12:04  eric
 * sauve declaree const.
 *
 * Revision 2.3  1999/09/30  12:52:52  eric
 * Depoussierage.
 * Documentation.
 *
 * Revision 2.2  1999/03/01  14:35:21  eric
 * Modif affichage (operator<<)
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Mg3d/mg3d.C,v 1.22 2023/05/24 09:48:50 g_servignat Exp $
 *
 */
// C++ Headers
#include <stdexcept>
 
// C Headers
#include <cstdlib>
#include <cassert>
 
// Lorene headers
#include "grilles.h"
#include "type_parite.h"
#include "utilitaires.h"
 
        //--------------//
        // Multi-grille //
        //--------------//
 
 
//=============================================================================
//    General constructor
//=============================================================================
 
 
namespace Lorene {
Mg3d::Mg3d(int nz, int nbr[], int typr[], int nbt[], int typt, int nbp[], 
       int typp, int* base_r)
    : nzone(nz), type_t(typt), type_p(typp)
{
 
    // Type d'echantillonnage dans chaque zone
    type_r = new int[nz];
    colloc_r = new int[nz] ;
    bool cheb = (base_r == 0x0) ;
    for (int i=0 ; i<nz ; i++) {
    type_r[i] = typr[i];
    colloc_r[i] = cheb ? BASE_CHEB : base_r[i] ;
    }
 
    // Nombre de points
    nr = new int[nz];
    nt = new int[nz];
    np = new int[nz];
    for (int i=0 ; i<nz ; i++) {
    nr[i] = nbr[i] ;
    nt[i] = nbt[i] ;
    np[i] = nbp[i] ;
    }
 
    // Les grilles
    // -----------
    g = new Grille3d*[nz] ;
 
    for (int i=0; i<nz; i++) {
      
      g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p, 
              colloc_r[i]) ;
    }   // fin de la boucle sur les zones
    
    // Pointers on derived grids initiated to 0x0:
    // -------------------------------------------
    
    set_deriv_0x0() ;
    
 
}
 
//=============================================================================
//    Simplified constructor
//=============================================================================
 
Mg3d::Mg3d(int nz, int nbr, int nbt, int nbp, int typt, int typp, 
       bool compact, bool leg)
            : nzone(nz),
              type_t(typt),
              type_p(typp)   {
 
    // Type of r sampling in each domain: 
    type_r = new int[nz];
    colloc_r = new int[nz];
    type_r[0] = RARE ;
    colloc_r[0] = leg ? BASE_LEG : BASE_CHEB ;
    for (int l=1; l<nz-1; l++) {
    type_r[l] = FIN ;
    colloc_r[l] = leg ? BASE_LEG : BASE_CHEB ;
    }
    if (nz > 1) {
      if (compact) {
    type_r[nz-1] = UNSURR ;
    colloc_r[nz-1] = BASE_CHEB ;
      }
      else {
    type_r[nz-1] = FIN ;
    colloc_r[nz-1] = leg ? BASE_LEG : BASE_CHEB ;
      }
    }
 
    // Same number of points in all domains:
    nr = new int[nz];
    nt = new int[nz];
    np = new int[nz];
    for (int l=0 ; l<nz ; l++) {
        nr[l] = nbr ;
        nt[l] = nbt ;
        np[l] = nbp ;
    }
 
    // Les grilles
    // -----------
    g = new Grille3d*[nz] ;
 
    for (int i=0; i<nz; i++) {
      
      g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p, 
              colloc_r[i]) ;
    }   // fin de la boucle sur les zones
    
    // Pointers on derived grids initiated to 0x0:
    // -------------------------------------------
    
    set_deriv_0x0() ;
    
}
 
//=============================================================================
//    Simplified shell constructor
//    Note: This does not handle the nucleus or the CED!
//=============================================================================
 
Mg3d::Mg3d(int nz, int nbr, int nbt, int nbp, int typt, int typp)
            : nzone(nz),
              type_t(typt),
              type_p(typp)   {
 
    // Type of r sampling in each domain: 
    type_r = new int[nz];
    colloc_r = new int[nz] ;
    for (int l=0; l<nz; l++) {
    type_r[l] = FIN ;
    colloc_r[l] = BASE_CHEB ;
    }
 
    // Same number of points in all domains:
    nr = new int[nz];
    nt = new int[nz];
    np = new int[nz];
    for (int l=0 ; l<nz ; l++) {
        nr[l] = nbr ;
        nt[l] = nbt ;
        np[l] = nbp ;
    }
 
    // Les grilles
    // -----------
    g = new Grille3d*[nz] ;
 
    for (int i=0; i<nz; i++) {
 
      g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p, 
              colloc_r[i]) ;
      
    }   // fin de la boucle sur les zones
 
    // Pointers on derived grids initiated to 0x0:
    // -------------------------------------------
 
    set_deriv_0x0() ;
    
}
 
//=============================================================================
//    Constructors from a file
//=============================================================================
 
 
  // From a formatted file...
  //=========================
  Mg3d::Mg3d(const string& filename) {
 
    // Opening of the file
    //-----------------------------------------------------------------------
    ifstream parfile(filename.c_str()) ;
    if (!parfile) {
      string message = "Unable to open file " ;
      message += filename ;
      throw ios_base::failure(message);
    };
 
    string tmp_str = "Definition of the Mg3d" ;
    if (!search_file(parfile, tmp_str)) {
      string mesg = "No data found to contruct an object Mg3d in " ;
      mesg += filename ;
      throw invalid_argument(mesg) ;
    }
    //    parfile >> tmp_str ;
    //cout << tmp_str ;
     parfile.ignore(1000, '\n') ;
 
     // Number of domains
     //---------------------
     parfile >> nzone ; parfile.ignore(1000, '\n') ; // Number of domains
     parfile >> tmp_str ;
     
     // Theta & phi symmetries + number of grid points
     //-----------------------------------------------
     if ( tmp_str.compare("SYM") == 0) {
       type_t = SYM ;
     }
     else if ( tmp_str.compare("NONSYM") == 0)
       {
     type_t = NONSYM ;
       }
     else
       throw
     invalid_argument("Mg3d::Mg3d(string): incorrect value for theta symmetry") ;
     parfile.ignore(1000, '\n') ;
     
     parfile >> tmp_str ;
     if ( tmp_str.compare("SYM") == 0)
       {
     type_p = SYM ;
       }
     else if ( tmp_str.compare("NONSYM") == 0)
       {
     type_p = NONSYM ;
       }
     else
       throw
     invalid_argument("Mg3d::Mg3d(string): incorrect value for phi symmetry") ;
     parfile.ignore(1000, '\n') ;
 
     int nbp ; parfile >> nbp ; parfile.ignore(1000, '\n') ;
     int nbt ; parfile >> nbt ; parfile.ignore(1000, '\n') ;
 
     // Radial number of points and sampling type
     //---------------------------------------------
     nr = new int[nzone] ;
     nt = new int[nzone] ;
     np = new int[nzone] ;
     type_r = new int[nzone] ;
  
     for (int i=0; i<nzone; i++) {
       parfile >> nr[i] >> tmp_str ;
       assert (nr[i] > 0) ;
       if ( tmp_str.compare("nucleus") == 0)
       {
     type_r[i] = RARE ;
       }
       else if ( tmp_str.compare("shell") == 0)
       {
     type_r[i] = FIN ;
       }
       else if ( tmp_str.compare("ced") == 0)
       {
     type_r[i] = UNSURR ;
       }
       else
     throw
       invalid_argument("Mg3d::Mg3d(string): incorrect value for the type of domain") ;
       parfile.ignore(1000, '\n') ;
       nt[i] = nbt ;
       np[i] = nbp ;
     } // end of loop on domains
 
     colloc_r = new int[nzone] ;
     bool legendre ;
     parfile >> legendre ;
     for (int i=0; i<nzone; i++)
       colloc_r[i] = (legendre ? BASE_LEG : BASE_CHEB) ;
    
 
    // Grids
    // -----------
    g = new Grille3d*[nzone] ;
 
    for (int i=0; i<nzone; i++) {
 
      g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p, 
              colloc_r[i]) ;
      
    }   // end of loop on domains
 
    // Pointers on derived grids initiated to 0x0:
    // -------------------------------------------
 
    set_deriv_0x0() ;
         
  }
 
  
/*
 * Construction a partir d'un fichier.
 * Cette facon de faire est abominable. Cependant je ne vois pas comment
 * faire autrement... j.a.
 */
Mg3d::Mg3d(FILE* fd, bool read_base)
{
    // Lecture sur le fichier
    fread_be(&nzone, sizeof(int), 1, fd) ;      // nzone
    nr = new int[nzone] ;
    fread_be(nr, sizeof(int), nzone, fd) ;      // nr
    nt = new int[nzone] ;
    fread_be(nt, sizeof(int), nzone, fd) ;      // nt
    np = new int[nzone] ;
    fread_be(np, sizeof(int), nzone, fd) ;      // np
    type_r = new int[nzone] ;
    fread_be(type_r, sizeof(int), nzone, fd) ;  // type_r
    fread_be(&type_t, sizeof(int), 1, fd) ; // type_t
    fread_be(&type_p, sizeof(int), 1, fd) ; // type_p
    colloc_r = new int[nzone] ;
    if (read_base)
      fread_be(colloc_r, sizeof(int), nzone, fd) ; // colloc_r
 
    // Les grilles
    // -----------
 
    g = new Grille3d*[nzone] ;
    for (int i=0; i<nzone; i++) {
      if (!read_base) colloc_r[i] = BASE_CHEB ;
      g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p, 
              colloc_r[i]) ;
      
    }   // fin de la boucle sur les zones
 
    // Pointers on derived grids initiated to 0x0:
    // -------------------------------------------
 
    set_deriv_0x0() ;
 
}
 
// Destructeur
// -----------
Mg3d::~Mg3d() {
 
    del_deriv() ;   // Deletes the derived quantities
 
    delete [] nr ;
    delete [] nt ;
    delete [] np ;
    delete [] type_r ;
    delete [] colloc_r ;
    for (int i=0 ; i<nzone ; i++) {
    delete g[i] ;
    }
    delete [] g ;
 
}
 
//==================================================================
//  Write in a file
//==================================================================
 
void Mg3d::sauve(FILE* fd, bool save_base) const {  
        fwrite_be(&nzone, sizeof(int), 1, fd) ; // nzone
        fwrite_be(nr, sizeof(int), nzone, fd) ; // nr
        fwrite_be(nt, sizeof(int), nzone, fd) ; // nt
        fwrite_be(np, sizeof(int), nzone, fd) ; // np
        fwrite_be(type_r, sizeof(int), nzone, fd) ; // type_r
        fwrite_be(&type_t, sizeof(int), 1, fd) ;    // type_t
        fwrite_be(&type_p, sizeof(int), 1, fd) ;    // type_p
        if (save_base) {
          fwrite_be(colloc_r, sizeof(int), nzone, fd) ; // colloc_r
        }
        else 
          for (int l=0; l<nzone; l++) 
        if (colloc_r[l] != BASE_CHEB) {
          cout << "Mg3d::sauve(FILE*, bool) : " << endl ;
          cout << "The multi-grid is not with Chebyshev basis!!" << endl ;
          cout << "Consider setting the 'save_base' flaf to 'true'!!" 
               << endl ;
          arrete() ;
          
        }
}
 
        //--------------------------//
        // Surcharge des operateurs //
        //--------------------------//
 
// Operateur <<
ostream& operator<<(ostream& o, const Mg3d& g) {
    const char* tr[3] ;
    tr[FIN] = "FIN" ; tr[RARE] = "RARE" ; tr[UNSURR] = "UNSURR" ;
    const char* tang[2] ;
    tang[NONSYM] = "NONSYM" ; tang[SYM] = "SYM" ;
    const char* tbase[3] ;
    tbase[BASE_CHEB] = "Chebyshev" ; tbase[BASE_LEG] = "Legendre" ; 
    tbase[BASE_JAC02] = "Jacobi(0,2)" ;
    o << "Number of domains: " << g.nzone << endl ;
    for (int i=0 ; i< g.nzone ; i++) {
    o << "  Domain #" << i << ": "
      << "nr = " << g.nr[i] << ", " << tr[g.type_r[i]] << "; "
      << "nt = " << g.nt[i] << ", " << tang[g.type_t] << "; "
      << "np = " << g.np[i] << ", " << tang[g.type_p] << "; "
      << "Collocation points type : " << tbase[g.colloc_r[i]] << endl ;
    }
    o << endl ;
    return o ;
}
 
// Operateur !=
bool Mg3d::operator!=(const Mg3d & titi) const {
 
    if (nzone != titi.nzone) return true ;   // C'est vrai que c'est faux...
 
    for (int i=0 ; i<nzone ; i++) {
    if (nr[i] != titi.nr[i]) return true ;
    if (nt[i] != titi.nt[i]) return true ;
    if (np[i] != titi.np[i]) return true ;
 
    if (type_r[i] != titi.type_r[i]) return true ;
    if (colloc_r[i] != titi.colloc_r[i]) return true ;
    }
 
    if (type_t != titi.type_t) return true ;
    if (type_p != titi.type_p) return true ;
 
    // C'est faux que c'est vrai...
    return false ;
}
 
 
            //----------------------------------//
            // Management of derived quantities //
            //----------------------------------//
 
void Mg3d::del_deriv() const {
 
    if (g_angu != 0x0) delete g_angu ;
    if (g_angu_1dom != 0x0) delete g_angu_1dom ;
    if (g_radial != 0x0) delete g_radial ;
    if (g_twice != 0x0) delete g_twice ;
    if (g_plus_half != 0x0) delete g_plus_half ;
    if (g_plus_half_angu != 0x0) delete g_plus_half_angu ;
    if (g_non_axi != 0x0) delete g_non_axi ;
 
    set_deriv_0x0() ;
 
}
 
void Mg3d::set_deriv_0x0() const {
 
    g_angu = 0x0 ;
    g_angu_1dom = 0x0 ;
    g_radial = 0x0 ;
    g_twice = 0x0 ;
    g_plus_half = 0x0 ;
    g_plus_half_angu = 0x0 ;
    g_non_axi = 0x0 ;
}
 
 
                //--------------//
                // Angular grid //
                //--------------//
 
const Mg3d* Mg3d::get_angu() const {
 
    if (g_angu == 0x0) {      // The construction is required
 
    int* nbr_angu = new int[nzone] ;
    for (int i=0 ; i<nzone ; i++) {
        nbr_angu[i] = 1 ;
    }
    g_angu = new Mg3d(nzone, nbr_angu, type_r, nt, type_t, np, type_p,
              colloc_r) ;
    delete [] nbr_angu ;
    }
 
    return g_angu ;
 
}
    
                //-----------------------------//
                // Angular grid for one domain //
                //-----------------------------//
 
const Mg3d* Mg3d::get_angu_1dom() const {
 
    if (g_angu_1dom == 0x0) {     // The construction is required
    int* nbr_angu = new int(1) ;
    int* nbt_angu = new int(nt[0]) ;
    int* nbp_angu = new int(np[0]) ;
    int* type_r_angu = new int(FIN) ;
    
    g_angu_1dom = new Mg3d(1, nbr_angu, type_r_angu, nbt_angu, type_t, 
                   nbp_angu, type_p) ;
    delete nbr_angu ;
    delete nbt_angu ;
    delete nbp_angu ;
    delete type_r_angu ;
    }
 
    return g_angu_1dom ;
 
}
    
                //--------------//
                //  Radial grid //
                //--------------//
 
const Mg3d* Mg3d::get_radial() const {
 
    if (g_radial == 0x0) {    // The construction is required
 
    int* nbr_radial = new int[nzone] ;
    for (int i=0 ; i<nzone ; i++) {
        nbr_radial[i] = 1 ;
    }
    g_radial = new Mg3d(nzone, nr, type_r, nbr_radial, SYM, nbr_radial, 
                SYM, colloc_r) ;
    delete [] nbr_radial ;
    }
 
    return g_radial ;
 
}
    
          //--------------------------------------//
          // Grid with twice the number of points //
          //--------------------------------------//
 
const Mg3d* Mg3d::get_twice() const {
 
    if (g_twice == 0x0) {     // The construction is required
 
    int* nbr = new int[nzone] ;
    int* nbt = new int[nzone] ;
    int* nbp = new int[nzone] ;
 
    for (int l=0; l<nzone; l++) {
        if (nr[l] == 1) {
        nbr[l] = 1 ;
        }
        else {
        nbr[l] = 2*nr[l] - 1 ;
        }
 
        if (nt[l] == 1) {
        nbt[l] = 1 ;
        }
        else {
        nbt[l] = 2*nt[l] - 1 ;
        }
    
        if (np[l] == 1) {
        nbp[l] = 1 ;
        }
        else {
        nbp[l] = 2*np[l] ;
        }
    }
 
    g_twice = new Mg3d(nzone, nbr, type_r, nbt, type_t, nbp, type_p, colloc_r) ;
 
    delete [] nbr ;
    delete [] nbt ;
    delete [] nbp ;
 
    }
 
    return g_twice ;
 
}
    
 
          //--------------------------------------//
          // Grid with 50% additional points in r //
          //--------------------------------------//
 
const Mg3d* Mg3d::plus_half() const {
 
  if (g_plus_half == 0x0) {   // The construction is required
 
    int* nbr = new int[nzone] ;
    
    for (int l=0; l<nzone; l++) {
      if (nr[l] == 1) 
    nbr[l] = 1 ;
      else 
    nbr[l] = (3*nr[l])/2  ;
    }
    
    g_plus_half = new Mg3d(nzone, nbr, type_r, nt, type_t, np, type_p, colloc_r) ;
 
    delete [] nbr ;
 
 
  }
 
  return g_plus_half ;
 
}
 
      //-------------------------------------------//
          // Grid with 50% additional points in angles //
          //-------------------------------------------//
 
const Mg3d* Mg3d::plus_half_angu() const {
 
  if (g_plus_half_angu == 0x0) {      // The construction is required
 
    int* nbt = new int[nzone] ;
    int* nbp = new int[nzone] ;
    
    for (int l=0; l<nzone; l++) {
      if (nt[l] == 1) 
    nbt[l] = 1 ;
      else {
  int nt_bis = (3*(nt[l]))/2;
  if (nt_bis % 2 == 0) nt_bis += 1 ;
    nbt[l] = nt_bis ;        
      }
      
      if (np[l] == 1) 
    nbp[l] = 1 ;
      else 
    nbp[l] = 3*np[l]/2 ;
    }
    
    g_plus_half_angu = new Mg3d(nzone, nr, type_r, nbt, type_t, nbp, type_p, colloc_r) ;
 
    delete [] nbt ;
    delete [] nbp ;
 
 
  }
 
  return g_plus_half_angu ;
 
}
    
          //----------------------------------------------//
          // Grid for rotations (5/4 points in theta/phi) //
          //----------------------------------------------//
 
const Mg3d* Mg3d::get_non_axi() const {
 
  if (g_non_axi == 0x0) {     // The construction is required
 
    int* nbt = new int[nzone] ;
    int* nbp = new int[nzone] ;
    
    for (int l=0; l<nzone; l++) {
      if (nt[l] < 5)   
      nbt[l] = 5 ;
      else
      nbt[l] = nt[l] ;
      if (np[l] < 4)
      nbp[l] = 4 ;
      else
      nbp[l] = np[l] ;
    }
    
    g_non_axi = new Mg3d(nzone, nr, type_r, nbt, type_t, nbp, type_p, colloc_r) ;
 
    delete [] nbt ;
    delete [] nbp ;
 
 
  }
 
  return g_non_axi ;
 
}
    
 
bool Mg3d::operator==(const Mg3d& mgi) const {
  
  bool resu = true ;
 
  if (mgi.get_nzone() != nzone) {
    resu = false ;
  }
  else {
    for (int i=0; i<nzone; i++) {
      if (mgi.get_nr(i) != nr[i]) resu = false ;
      if (mgi.get_np(i) != np[i]) resu = false ;
      if (mgi.get_nt(i) != nt[i]) resu = false ;
      if (mgi.get_type_r(i) != type_r[i]) resu =false ;
      if (mgi.get_colloc_r(i) != colloc_r[i]) resu = false ;
    }
  }
  
  if (mgi.get_type_t() != type_t) resu = false ;
  if (mgi.get_type_p() != type_p) resu = false ;
 
  return resu ;
 
}
}