/*!
  \file mesh.cpp
  \author Olga Abramkina
  \brief Definition of class CMesh.
*/

#include "mesh.hpp"

namespace xios {

/// ////////////////////// Définitions ////////////////////// ///

  CMesh::CMesh(void) :  nbNodesGlo{0}, nbEdgesGlo{0}
            ,  node_start{0}, node_count{0}
            ,  edge_start{0}, edge_count{0}
            ,  nbFaces_{0}, nbNodes_{0}, nbEdges_{0}
            ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false}
            ,  node_lon(), node_lat()
            ,  edge_lon(), edge_lat(), edge_nodes()
            ,  face_lon(), face_lat()
            ,  face_nodes()
            ,  pNodeGlobalIndex{NULL}, pEdgeGlobalIndex{NULL}
  {
  }


  CMesh::~CMesh(void)
  {
    if (pNodeGlobalIndex != NULL) delete pNodeGlobalIndex;
    if (pEdgeGlobalIndex != NULL) delete pEdgeGlobalIndex;
  }

  std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>();
  std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >();

///---------------------------------------------------------------
/*!
 * \fn bool CMesh::getMesh (StdString meshName)
 * Returns a pointer to a mesh. If a mesh has not been created, creates it and adds its name to the list of meshes meshList.
 * \param [in] meshName  The name of a mesh ("name" attribute of a domain).
 * \param [in] nvertex Number of verteces (1 for nodes, 2 for edges, 3 and up for faces).
 */
  CMesh* CMesh::getMesh (StdString meshName, int nvertex)
  {
    CMesh::domainList[meshName].push_back(nvertex);

    if ( CMesh::meshList.begin() != CMesh::meshList.end() )
    {
      for (std::map<StdString, CMesh>::iterator it=CMesh::meshList.begin(); it!=CMesh::meshList.end(); ++it)
      {
        if (it->first == meshName)
          return &meshList[meshName];
        else
        {
          CMesh newMesh;
          CMesh::meshList.insert( make_pair(meshName, newMesh) );
          return &meshList[meshName];
        }
      }
    }
    else
    {
      CMesh newMesh;
      CMesh::meshList.insert( make_pair(meshName, newMesh) );
      return &meshList[meshName];
    }
  }

///----------------------------------------------------------------
  size_t hashPair(size_t first, size_t second)
  {
    HashXIOS<size_t> sizetHash;
    size_t seed = sizetHash(first) + 0x9e3779b9 ;
    seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
    return seed ;
  }

///----------------------------------------------------------------
  size_t hashPairOrdered(size_t first, size_t second)
  {
    size_t seed;
    HashXIOS<size_t> sizetHash;
    if (first < second)
    {
      seed = sizetHash(first) + 0x9e3779b9 ;
      seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
    }
    else
    {
      seed = sizetHash(second) + 0x9e3779b9 ;
      seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
    }
    return seed ;
  }

///----------------------------------------------------------------
/*!
 * \fn size_t generateEdgeIndex(size_t first, size_t second, int rank)
 * Generates an edge index.
 * If the same edge is generated by two processes, each process will have its own edge index.
 * \param [in] first Edge node.
 * \param [in] second Edge node.
 * \param [in] rank MPI process rank.
 */
  size_t generateEdgeIndex(size_t first, size_t second, int rank)
  {
    size_t seed = rank ;
    if (first < second)
    {
      seed = hashPair(seed, first);
      seed = hashPair(seed, second);
    }
    else
    {
      seed = hashPair(seed, second);
      seed = hashPair(seed, first);
    }

    return seed ;
  }

///----------------------------------------------------------------
/*!
 * \fn size_t generateNodeIndex(vector<size_t>& valList, int rank)
 * Generates a node index.
 * If the same node is generated by two processes, each process will have its own node index.
 * \param [in] valList Vector storing four node hashes.
 * \param [in] rank MPI process rank.
 */
  size_t generateNodeIndex(vector<size_t>& valList, int rank)
  {
    // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
    vector<size_t> vec = valList;
    sort (vec.begin(), vec.end());
    size_t seed = rank ;
    int it = 0;
    for(; it != vec.size(); ++it)
    {
       seed = hashPair(seed, vec[it]);
    }
    return seed ;
  }

///----------------------------------------------------------------
/*!
 * \fn size_t CMesh::nodeIndex (double lon, double lat)
 * Returns its index if a node exists; otherwise adds the node and returns -1.
 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
 * \param [in] lon Node longitude in degrees.
 * \param [in] lat Node latitude in degrees ranged from 0 to 360.
 * \return node index if a node exists; -1 otherwise
 */
  size_t CMesh::nodeIndex (double lon, double lat)
  {
    double minBoundLon = 0. ;
    double maxBoundLon = 360. ;
    double minBoundLat = -90 ;
    double maxBoundLat = 90 ;
    double prec=1e-11 ;
    double precLon=prec ;
    double precLat=prec ;

    size_t maxsize_t=numeric_limits<size_t>::max() ;
    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;

    size_t iMinLon=0 ;
    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
    size_t iMinLat=0 ;
    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;

    size_t hash0,hash1,hash2,hash3 ;
    size_t lon0,lon1,lat0,lat1 ;

    lon0=(lon-minBoundLon)/precLon ;
    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
    {
      if (lon0==iMinLon) lon1=iMaxLon ;
      else lon1=lon0-1 ;
    }
    else
    {
      if (lon0==iMaxLon) lon1=iMinLon ;
      else lon1=lon0+1 ;
    }

    lat0=(lat-minBoundLat)/precLat ;
    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
    {
      if (lat0==iMinLat) lat1=lat0 ;
      else lat1=lat0-1 ;
    }
    else
    {
      if (lat0==iMaxLat) lat1=lat0 ;
      else lat1=lat0+1 ;
    }

    hash0=hashPair(lon0,lat0) ;
    hash1=hashPair(lon0,lat1) ;
    hash2=hashPair(lon1,lat0) ;
    hash3=hashPair(lon1,lat1) ;

    boost::unordered_map<size_t, size_t>::iterator end = hashed_map_nodes.end() ;
    size_t mapSize = hashed_map_nodes.size();
    if (hashed_map_nodes.find(hash0)==end && hashed_map_nodes.find(hash1)==end && hashed_map_nodes.find(hash2)==end && hashed_map_nodes.find(hash3)==end)
    {
      hashed_map_nodes[hash0] = mapSize ;
      hashed_map_nodes[hash1] = mapSize + 1;
      hashed_map_nodes[hash2] = mapSize + 2;
      hashed_map_nodes[hash3] = mapSize + 3;
      return -1;
    }
    else
      return ( (hashed_map_nodes[hash0]+1) / 4 );

  } // nodeIndex()

///----------------------------------------------------------------
/*!
 * \fn CArray<size_t,1>& CMesh::createHashes (const double longitude, const double latitude)
 * Creates two hash values for each dimension, longitude and latitude.
 * \param [in] longitude Node longitude in degrees.
 * \param [in] latitude Node latitude in degrees ranged from 0 to 360.
 */

  vector<size_t> CMesh::createHashes (const double longitude, const double latitude)
  {
    double minBoundLon = 0. ;
    double maxBoundLon = 360. ;
    double minBoundLat = -90. ;
    double maxBoundLat = 90. ;
    double prec=1e-11 ;
    double precLon=prec ;
    double precLat=prec ;
    double lon = longitude;
    double lat = latitude;

    if (lon > (360.- prec)) lon = 0.;

    size_t maxsize_t=numeric_limits<size_t>::max() ;
    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;

    size_t iMinLon=0 ;
    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
    size_t iMinLat=0 ;
    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;

    vector<size_t> hash(4);
    size_t lon0,lon1,lat0,lat1 ;

    lon0=(lon-minBoundLon)/precLon ;
    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
    {
      if (lon0==iMinLon) lon1=iMaxLon ;
      else lon1=lon0-1 ;
    }
    else
    {
      if (lon0==iMaxLon) lon1=iMinLon ;
      else lon1=lon0+1 ;
    }

    lat0=(lat-minBoundLat)/precLat ;
    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
    {
      if (lat0==iMinLat) lat1=lat0 ;
      else lat1=lat0-1 ;
    }
    else
    {
      if (lat0==iMaxLat) lat1=lat0 ;
      else lat1=lat0+1 ;
    }

    hash[0] = hashPair(lon0,lat0) ;
    hash[1] = hashPair(lon0,lat1) ;
    hash[2] = hashPair(lon1,lat0) ;
    hash[3] = hashPair(lon1,lat1) ;

    return hash;

  } // createHashes

///----------------------------------------------------------------
  std::pair<int,int> make_ordered_pair(int a, int b)
  {
    if ( a < b )
      return std::pair<int,int>(a,b);
    else
      return std::pair<int,int>(b,a);
  }

///----------------------------------------------------------------
/*!
 * \fn void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
 * \param [in] lonvalue  Array of longitudes.
 * \param [in] latvalue  Array of latitudes.
 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
 */
  void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
  {
    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();

    if (nvertex == 1)
    {
      nbNodes_ = lonvalue.numElements();
      node_lon.resizeAndPreserve(nbNodes_);
      node_lat.resizeAndPreserve(nbNodes_);
      for (int nn = 0; nn < nbNodes_; ++nn)
      {
        if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end())
        {
          map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ;
          node_lon(nn) = lonvalue(nn);
          node_lat(nn) = latvalue(nn);
        }
      }
    }
    else if (nvertex == 2)
    {
      nbEdges_ = bounds_lon.shape()[1];

      // Create nodes and edge_node connectivity
      node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes
      node_lat.resizeAndPreserve(nbEdges_*nvertex);
      edge_nodes.resizeAndPreserve(nvertex, nbEdges_);

      for (int ne = 0; ne < nbEdges_; ++ne)
      {
        for (int nv = 0; nv < nvertex; ++nv)
        {
          if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end())
          {
            map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ;
            edge_nodes(nv,ne) = nbNodes_ ;
            node_lon(nbNodes_) = bounds_lon(nv, ne);
            node_lat(nbNodes_) = bounds_lat(nv, ne);
            ++nbNodes_ ;
          }
          else
            edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))];
        }
      }
      node_lon.resizeAndPreserve(nbNodes_);
      node_lat.resizeAndPreserve(nbNodes_);

      // Create edges
      edge_lon.resizeAndPreserve(nbEdges_);
      edge_lat.resizeAndPreserve(nbEdges_);

      for (int ne = 0; ne < nbEdges_; ++ne)
      {
        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
        {
          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
          edge_lon(ne) = lonvalue(ne);
          edge_lat(ne) = latvalue(ne);
        }

      }
      edgesAreWritten = true;
    }
    else
    {
      nbFaces_ = bounds_lon.shape()[1];
  
      // Create nodes and face_node connectivity
      node_lon.resizeAndPreserve(nbFaces_*nvertex);  // Max possible number of nodes
      node_lat.resizeAndPreserve(nbFaces_*nvertex);
      face_nodes.resize(nvertex, nbFaces_);
  
      for (int nf = 0; nf < nbFaces_; ++nf)
      {
        for (int nv = 0; nv < nvertex; ++nv)
        {
          if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end())
          {
            map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ;
            face_nodes(nv,nf) = nbNodes_ ;
            node_lon(nbNodes_) = bounds_lon(nv, nf);
            node_lat(nbNodes_) = bounds_lat(nv ,nf);
            ++nbNodes_ ;
          }
          else
          {
            face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
          }
        }
      }
      node_lon.resizeAndPreserve(nbNodes_);
      node_lat.resizeAndPreserve(nbNodes_);
  
      // Create edges and edge_nodes connectivity
      edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges
      edge_lat.resizeAndPreserve(nbFaces_*nvertex);
      edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex);
      edge_faces.resize(2, nbFaces_*nvertex);
      face_edges.resize(nvertex, nbFaces_);
      face_faces.resize(nvertex, nbFaces_);

      vector<int> countEdges(nbFaces_*nvertex);   // needed in case if edges have been already generated
      vector<int> countFaces(nbFaces_);
      countEdges.assign(nbFaces_*nvertex, 0);
      countFaces.assign(nbFaces_, 0);
      int edge;
      for (int nf = 0; nf < nbFaces_; ++nf)
      {
        for (int nv1 = 0; nv1 < nvertex; ++nv1)
        {
          int nv = 0;
          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
          {
            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ;
            face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
            edge_faces(0,nbEdges_) = nf;
            edge_faces(1,nbEdges_) = -999;
            face_faces(nv1,nf) = 999999;
            edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf);
            edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);
            edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
            ++nbEdges_;
          }
          else
          {
            edge = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
            face_edges(nv1,nf) = edge;
            if (edgesAreWritten)
            {
              edge_faces(countEdges[edge], edge) = nf;
              if (countEdges[edge]==0)
              {
                face_faces(nv1,nf) = 999999;
              }
              else
              {
                int face1 = nf; // = edge_faces(1,edge)
                int face2 = edge_faces(0,edge);
                face_faces(countFaces[face1], face1) =  face2;
                face_faces(countFaces[face2], face2) =  face1;
                ++(countFaces[face1]);
                ++(countFaces[face2]);
              }
            }
            else
            {
              edge_faces(1,edge) = nf;
              int face1 = nf; // = edge_faces(1,edge)
              int face2 = edge_faces(0,edge);
              face_faces(countFaces[face1], face1) =  face2;
              face_faces(countFaces[face2], face2) =  face1;
              ++(countFaces[face1]);
              ++(countFaces[face2]);
            }
            ++(countEdges[edge]);
          }
        }
      }
      edge_nodes.resizeAndPreserve(2, nbEdges_);
      edge_faces.resizeAndPreserve(2, nbEdges_);
      edge_lon.resizeAndPreserve(nbEdges_);
      edge_lat.resizeAndPreserve(nbEdges_);

      // Create faces
      face_lon.resize(nbFaces_);
      face_lat.resize(nbFaces_);
      face_lon = lonvalue;
      face_lat = latvalue;
      facesAreWritten = true;

    } // nvertex > 2
    
  } // createMesh()

///----------------------------------------------------------------
/*!
 * \fn void CMesh::createMeshEpsilon(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
 * \param [in] comm
 * \param [in] lonvalue  Array of longitudes.
 * \param [in] latvalue  Array of latitudes.
 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
 */
  void CMesh::createMeshEpsilon(const ep_lib::MPI_Comm& comm,
                                const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
  {

    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
    int mpiRank, mpiSize;
    MPI_Comm_rank(comm, &mpiRank);
    MPI_Comm_size(comm, &mpiSize);
    double prec = 1e-11;  // used in calculations of edge_lon/lat

    if (nvertex == 1)
    {
      nbNodes_ = lonvalue.numElements();
      node_lon.resize(nbNodes_);
      node_lat.resize(nbNodes_);
      node_lon = lonvalue;
      node_lat = latvalue;

      // Global node indexes
      vector<size_t> hashValues(4);
      CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
      for (size_t nn = 0; nn < nbNodes_; ++nn)
      {
        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn));
        for (size_t nh = 0; nh < 4; ++nh)
        {
          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes_ + nn);
        }
      }
      pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
      nodesAreWritten = true;
    }

    else if (nvertex == 2)
    {
      nbEdges_ = bounds_lon.shape()[1];
      edge_lon.resize(nbEdges_);
      edge_lat.resize(nbEdges_);
      edge_lon = lonvalue;
      edge_lat = latvalue;
      edge_nodes.resize(nvertex, nbEdges_);

      // For determining the global edge index
      size_t nbEdgesOnProc = nbEdges_;
      size_t nbEdgesAccum;
      MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
      nbEdgesAccum -= nbEdges_;

      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo;
      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;

      // Case (1): node indexes have been generated by domain "nodes"
      if (nodesAreWritten)
      {
        vector<size_t> hashValues(4);
        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4);
        for (int ne = 0; ne < nbEdges_; ++ne)      // size_t doesn't work with CArray<double, 2>
        {
          for (int nv = 0; nv < nvertex; ++nv)
          {
            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
            for (int nh = 0; nh < 4; ++nh)
            {
              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
            }
          }
        }

        // Recuperating the node global indexing and writing edge_nodes
        // Creating map edgeHash2IdxGlo <hash, idxGlo>
        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
        size_t nodeIdxGlo1, nodeIdxGlo2;
        for (int ne = 0; ne < nbEdges_; ++ne)
        {
          for (int nv = 0; nv < nvertex; ++nv)
          {
            int nh = 0;
            it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
            // The loop below is needed in case if a hash generated by domain "edges" differs
            // from that generated by domain "nodes" because of possible precision issues
            while (it == nodeHash2IdxGlo.end())
            {
              ++nh;
              it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
            }
            edge_nodes(nv,ne) = it->second[0];
            if (nv ==0)
              nodeIdxGlo1 = it->second[0];
            else
              nodeIdxGlo2 = it->second[0];
          }
          size_t edgeIdxGlo = nbEdgesAccum + ne;
          edgeHash2IdxGlo[ hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2) ].push_back(edgeIdxGlo);
        }
      } // nodesAreWritten


      // Case (2): node indexes have not been generated previously
      else
      {
        // (2.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
        vector<size_t> hashValues(4);
        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4);
        for (int ne = 0; ne < nbEdges_; ++ne)
        {
          for (int nv = 0; nv < nvertex; ++nv)
          {
            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
            for (int nh = 0; nh < 4; ++nh)
            {
              if (nodeHash2Idx[hashValues[nh]].size() == 0)
              {
                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues, mpiRank));
                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
              }
              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
            }
          }
        }

        // (2.2) Generating global node indexes
        // The ownership criterion: priority of the process holding the smaller index
        // Maps generated in this step are:
        // nodeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
        // nodeIdx2IdxMin = <idx, idxMin>
        // nodeIdx2IdxGlo = <idxMin, idxMin>

        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
        CArray<size_t,1> nodeIdxMinList(nbEdges_*nvertex);

        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
        size_t iIdxMin = 0;

        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
        {
          size_t idx = (it->second)[0];
          size_t idxMin = (it->second)[0];
          for (int i = 2; i < (it->second).size();)
          {
            if (mpiRank == (it->second)[i+1])
            {
              idx = (it->second)[i];
            }
            if ((it->second)[i] < idxMin)
            {
                idxMin = (it->second)[i];
                (it->second)[i] = (it->second)[i-2];
            }
            i += 2;
          }
          (it->second)[0] = idxMin;
          if (nodeIdx2IdxMin.count(idx) == 0)
          {
            nodeIdx2IdxMin[idx].push_back(idxMin);
            if (idx == idxMin)
              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
            nodeIdxMinList(iIdxMin) = idxMin;
            ++iIdxMin;
          }
        }
        nodeIdxMinList.resizeAndPreserve(iIdxMin);
        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
        // nodeIdx2IdxGlo holds global indexes only for nodes owned by a process
        // nodeIdxMin2IdxGlo holds global indexes for all nodes generated by a process

        // (2.3) Saving variables: node_lon, node_lat, edge_nodes
        // Creating map nodeHash2IdxGlo <hash, idxGlo>
        // Creating map edgeHash2IdxGlo <hash, idxGlo>
        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
        node_count = dhtNodeIdxGlo.getIndexCount();
        node_start = dhtNodeIdxGlo.getIndexStart();
        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
        node_lon.resize(node_count);
        node_lat.resize(node_count);
        vector <size_t> edgeNodes;
        size_t idxGlo = 0;

        for (int ne = 0; ne < nbEdges_; ++ne)
        {
          for (int nv = 0; nv < nvertex; ++nv)
          {
            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
            size_t myIdx = generateNodeIndex(hashValues, mpiRank);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxMin.find(myIdx);
            size_t ownerIdx = (itIdx->second)[0];

            if (myIdx == ownerIdx)
            {
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdx2IdxGlo.find(myIdx);
              idxGlo = (itIdxGlo->second)[0];
//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
              node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
              node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
            }
            else
            {
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdxMin2IdxGlo.find(ownerIdx);
              idxGlo = (itIdxGlo->second)[0];
            }

            edge_nodes(nv,ne) = idxGlo;
            for (int nh = 0; nh < 4; ++nh)
              nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo);
            edgeNodes.push_back(idxGlo);
          }
          if (edgeNodes[0] != edgeNodes[1])
          {
            size_t edgeIdxGlo = nbEdgesAccum + ne;
            edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(edgeIdxGlo);
          }
          edgeNodes.clear();
        }
        pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
      } // !nodesAreWritten

      pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm);
      edgesAreWritten = true;
    } //nvertex = 2

    else
    {
      nbFaces_ = bounds_lon.shape()[1];
      face_lon.resize(nbFaces_);
      face_lat.resize(nbFaces_);
      face_lon = lonvalue;
      face_lat = latvalue;
      face_nodes.resize(nvertex, nbFaces_);
      face_edges.resize(nvertex, nbFaces_);

      // For determining the global face index
      size_t nbFacesOnProc = nbFaces_;
      size_t nbFacesAccum;
      MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
      nbFacesAccum -= nbFaces_;

      // Case (1): edges have been previously generated
      if (edgesAreWritten)
      {
        // (1.1) Recuperating node global indexing and saving face_nodes
        vector<size_t> hashValues(4);
        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv = 0; nv < nvertex; ++nv)
            {
              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
              for (int nh = 0; nh < 4; ++nh)
                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
            }
        }
        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
        size_t nEdge = 0;
        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            int nh1 = 0;
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            while (it1 == nodeHash2IdxGlo.end())
            {
              ++nh1;
              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            }
            int nh2 = 0;
            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
            while (it2 == nodeHash2IdxGlo.end())
            {
              ++nh2;
              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
            }
            face_nodes(nv1,nf) = it1->second[0];
            if (it1->second[0] != it2->second[0])
            {
              edgeHashList(nEdge) = hashPairOrdered(it1->second[0], it2->second[0]);
              ++nEdge;
            }
          }
        }
        edgeHashList.resizeAndPreserve(nEdge);

        // (1.2) Recuperating edge global indexing and saving face_edges
        pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash;
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
        size_t iIdx = 0;

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            int nh1 = 0;
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            while (it1 == nodeHash2IdxGlo.end())
            {
              ++nh1;
              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            }
            int nh2 = 0;
            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
            while (it2 == nodeHash2IdxGlo.end())
            {
              ++nh2;
              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
            }
            if (it1->second[0] != it2->second[0])
            {
              size_t faceIdxGlo = nbFacesAccum + nf;
              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
              itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
              size_t edgeIdxGlo = itEdgeHash->second[0];
              face_edges(nv1,nf) = edgeIdxGlo;
              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
              {
                edgeIdxGloList(iIdx) = edgeIdxGlo;
                ++iIdx;
              }
              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
              edgeHash2Rank[edgeHash].push_back(itEdgeHash->first);
              edgeHash2Rank[edgeHash].push_back(mpiRank);
            }
            else
            {
              face_edges(nv1,nf) = 999999;
            }
          }
        }
        edgeIdxGloList.resizeAndPreserve(iIdx);

        // (1.3) Saving remaining variables edge_faces and face_faces

        // Establishing edge ownership
        // The ownership criterion: priority of the process with smaller rank
        CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm);
        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo;

        // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>>
        // edgeHashMin2IdxGlo = <edgeHashMin, idxGlo>
        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
        {
          vector <size_t> edgeInfo = it->second;
          if (edgeInfo.size() == 4)                       // two processes generate the same edge
            if (edgeInfo[1] > edgeInfo[3])
              edgeInfo[1] = edgeInfo[3];
          if (edgeInfo[1] == mpiRank)
            edgeHashMin2IdxGlo[it->first].push_back(edgeInfo[0]);
        }

        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeHashMin2IdxGlo, comm);
        edge_count = dhtEdgeIdxGlo.getIndexCount();
        edge_start = dhtEdgeIdxGlo.getIndexStart();

        edge_faces.resize(2, edge_count);
        face_faces.resize(nvertex, nbFaces_);

        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2;

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            int nh1 = 0;
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            while (it1 == nodeHash2IdxGlo.end())
            {
              ++nh1;
              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            }
            int nh2 = 0;
            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
            while (it2 == nodeHash2IdxGlo.end())
            {
              ++nh2;
              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
            }

            if (it1->second[0] != it2->second[0])
            {
              size_t faceIdxGlo = nbFacesAccum + nf;
              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
              itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
              if (itEdgeHash != edgeHashMin2IdxGlo.end())
              {
                int edgeIdxGlo = itEdgeHash->second[0];
                itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                int face1 = itFace1->second[0];
                if (itFace1->second.size() == 1)
                {
                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
                  edge_faces(1, edgeIdxGlo - edge_start) = -999;
                  face_faces(nv1, nf) = 999999;
                }
                else
                {
                  int face2 = itFace1->second[1];
                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
                  edge_faces(1, edgeIdxGlo - edge_start) = face2;
                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
                }
              } // edge owner
              else
              {
                itEdgeHash = edgeHashMin2IdxGlo.find(edgeHash);
                size_t edgeIdxGlo = itEdgeHash->second[0];
                itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                int face1 = itFace1->second[0];
                int face2 = itFace1->second[1];
                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
              } // not an edge owner
            } // node1 != node2
            else
            {
              face_faces(nv1, nf) = 999999;
            }
          }
        }
      } // edgesAreWritten

      // Case (2): nodes have been previously generated
      else if (nodesAreWritten)
      {
        // (2.1) Generating nodeHashList
        vector<size_t> hashValues(4);
        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv = 0; nv < nvertex; ++nv)
            {
              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
              for (int nh = 0; nh < 4; ++nh)
                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
            }
        }

        // (2.2) Recuperating node global indexing and saving face_nodes
        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            int nh1 = 0;
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            while (it1 == nodeHash2IdxGlo.end())
            {
              ++nh1;
              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            }
            int nh2 = 0;
            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
            while (it2 == nodeHash2IdxGlo.end())
            {
              ++nh2;
              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
            }
            face_nodes(nv1,nf) = it1->second[0];
            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(it1->second[0], it2->second[0], mpiRank));
            edgeHash2Idx[edgeHash].push_back(mpiRank);
            edgeHashList(nf*nvertex + nv1) = edgeHash;
          }
        }

        // (2.2) Generating global edge indexes
        // The ownership criterion: priority of the process holding the smaller index
        // Maps generated in this step are:
        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
        // edgeIdx2IdxMin = = <idx, idxMin>
        // edgeIdx2IdxGlo = <idxMin, idxGlo>

        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();

        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
        CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex);
        size_t iIdxMin = 0;

        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
        {
          size_t idxMin = (it->second)[0];
          size_t idx = (it->second)[0];
          for (int i = 2; i < (it->second).size();)
          {
            if (mpiRank == (it->second)[i+1])
            {
              idx = (it->second)[i];
            }
            if ((it->second)[i] < idxMin)
            {
                idxMin = (it->second)[i];
                (it->second)[i] = (it->second)[i-2];
            }
            i += 2;
          }
          (it->second)[0] = idxMin;
          if (edgeIdx2IdxMin.count(idx) == 0)
          {
            edgeIdx2IdxMin[idx].push_back(idxMin);
            if (idx == idxMin)
              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
            edgeIdxMinList(iIdxMin) = idxMin;
            ++iIdxMin;
          }
        }
        edgeIdxMinList.resizeAndPreserve(iIdxMin);
        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
        // edgeIdx2IdxGlo holds global indexes only for edges owned by a process
        // edgeIdxMin2IdxGlo holds global indexes for all edges generated by a process

        // (2.3) Saving variables: edge_lon, edge_lat, face_edges
        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
        edge_count = dhtEdgeIdxGlo.getIndexCount();
        edge_start = dhtEdgeIdxGlo.getIndexStart();
        edge_lon.resize(edge_count);
        edge_lat.resize(edge_count);
        edge_nodes.resize(2, edge_count);
        face_edges.resize(nvertex, nbFaces_);

        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
        size_t iIdx = 0;

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            // Getting global indexes of edge's nodes
            int nh1 = 0;
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            while (it1 == nodeHash2IdxGlo.end())
            {
              ++nh1;
              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            }
            int nh2 = 0;
            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
            while (it2 == nodeHash2IdxGlo.end())
            {
              ++nh2;
              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
            }
            // Getting edge global index
            size_t nodeIdxGlo1 = it1->second[0];
            size_t nodeIdxGlo2 = it2->second[0];
            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
              size_t ownerIdx = (itIdx->second)[0];
              int edgeIdxGlo = 0;
              size_t faceIdxGlo = nbFacesAccum + nf;

              if (myIdx == ownerIdx)
              {
                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
                edgeIdxGlo = (itIdxGlo->second)[0];
                double edgeLon;
                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
                if (diffLon < (180.- prec))
                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
                else if (diffLon > (180.+ prec))
                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
                else
                  edgeLon = 0.;
                edge_lon(edgeIdxGlo - edge_start) = edgeLon;
                edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
              }
              else
              {
                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
                edgeIdxGlo = (itIdxGlo->second)[0];
              }
              face_edges(nv1,nf) = edgeIdxGlo;
              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
              {
                edgeIdxGloList(iIdx) = edgeIdxGlo;
                ++iIdx;
              }
              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
            } // nodeIdxGlo1 != nodeIdxGlo2
            else
            {
              face_edges(nv1,nf) = 999999;
            }
          }
        }
        edgeIdxGloList.resizeAndPreserve(iIdx);

        // (2.4) Saving remaining variables edge_faces and face_faces
        edge_faces.resize(2, edge_count);
        face_faces.resize(nvertex, nbFaces_);

        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            // Getting global indexes of edge's nodes
            int nh1 = 0;
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            while (it1 == nodeHash2IdxGlo.end())
            {
              ++nh1;
              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
            }
            int nh2 = 0;
            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
            while (it2 == nodeHash2IdxGlo.end())
            {
              ++nh2;
              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
            }
            size_t nodeIdxGlo1 = it1->second[0];
            size_t nodeIdxGlo2 = it2->second[0];

            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
            size_t ownerIdx = (itIdx->second)[0];
            size_t faceIdxGlo = nbFacesAccum + nf;

            if (myIdx == ownerIdx)
            {
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
              int edgeIdxGlo = (itIdxGlo->second)[0];
              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
              int face1 = it1->second[0];
              if (it1->second.size() == 1)
              {
                edge_faces(0, edgeIdxGlo - edge_start) = face1;
                edge_faces(1, edgeIdxGlo - edge_start) = -999;
                face_faces(nv1, nf) = 999999;
              }
              else
              {
                int face2 = it1->second[1];
                edge_faces(0, edgeIdxGlo - edge_start) = face1;
                edge_faces(1, edgeIdxGlo - edge_start) = face2;
                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
              }
            }
            else
            {
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
              size_t edgeIdxGlo = (itIdxGlo->second)[0];
              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
              int face1 = it1->second[0];
              int face2 = it1->second[1];
              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
            }
          }
        }
      } // nodesAreWritten

      // Case (3): Neither nodes nor edges have been previously generated
      else
      {
        // (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
        vector<size_t> hashValues(4);
        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
        size_t iHash = 0;
        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv = 0; nv < nvertex; ++nv)
          {
            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
            for (int nh = 0; nh < 4; ++nh)
            {
              if (nodeHash2Idx.count(hashValues[nh])==0)
              {
                nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
                nodeHashList(iHash) = hashValues[nh];
                ++iHash;
              }
            }
          }
        }
        nodeHashList.resizeAndPreserve(iHash);

        // (3.2) Generating global node indexes
        // The ownership criterion: priority of the process holding the smaller index
        // Maps generated in this step are:
        // nodeHash2Info = <hash, [[idxMin, rankMin], [idx1, rank1], [idx3, rank3]..]>
        // nodeIdx2IdxMin = = <idx, idxMin>
        // nodeIdx2IdxGlo = <idxMin, idxGlo>

        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();

        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
        CArray<size_t,1> nodeIdxMinList(nbFaces_*nvertex*4);
        size_t iIdxMin = 0;

        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
        {
          size_t idxMin = (it->second)[0];
          size_t idx = (it->second)[0];
          for (int i = 2; i < (it->second).size();)
          {
            if (mpiRank == (it->second)[i+1])
            {
              idx = (it->second)[i];
            }
            if ((it->second)[i] < idxMin)
            {
                idxMin = (it->second)[i];
                (it->second)[i] = (it->second)[i-2];
            }
            i += 2;
          }
          (it->second)[0] = idxMin;
          if (nodeIdx2IdxMin.count(idx) == 0)
          {
            nodeIdx2IdxMin[idx].push_back(idxMin);
            if (idx == idxMin)
              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
            nodeIdxMinList(iIdxMin) = idxMin;
            ++iIdxMin;
          }
        }

        nodeIdxMinList.resizeAndPreserve(iIdxMin);
        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();

        // (3.3) Saving node data: node_lon, node_lat, and face_nodes
        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
        node_count = dhtNodeIdxGlo.getIndexCount();
        node_start = dhtNodeIdxGlo.getIndexStart();
        node_lon.resize(node_count);
        node_lat.resize(node_count);
        size_t nodeIdxGlo1 = 0;
        size_t nodeIdxGlo2 = 0;
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
        size_t nEdgeHash = 0;

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxMin.find(myNodeIdx1);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxMin.find(myNodeIdx2);
            size_t ownerNodeIdx = (itNodeIdx1->second)[0];

            if (myNodeIdx1 == ownerNodeIdx)
            {
              itNodeIdxGlo1 = nodeIdx2IdxGlo.find(myNodeIdx1);
              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
            }
            else
            {
              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find(ownerNodeIdx);
              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
            }
            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]);
            nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
              edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank));
              edgeHash2Idx[edgeHash].push_back(mpiRank);
              edgeHashList(nEdgeHash) = edgeHash;
              ++nEdgeHash;
            }
            face_nodes(nv1,nf) = nodeIdxGlo1;
          }
        }
        edgeHashList.resizeAndPreserve(nEdgeHash);

        // (3.4) Generating global edge indexes
        // Maps generated in this step are:
        // edgeIdx2IdxMin = = <idx, idxMin>
        // edgeIdx2IdxGlo = <idxMin, idxGlo>

        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>

        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
        CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex);
        iIdxMin = 0;

        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
        {
          size_t idxMin = (it->second)[0];
          size_t idx = (it->second)[0];

          for (int i = 2; i < (it->second).size();)
          {
            if (mpiRank == (it->second)[i+1])
            {
              idx = (it->second)[i];
            }
            if ((it->second)[i] < idxMin)
            {
                idxMin = (it->second)[i];
                (it->second)[i] = (it->second)[i-2];
            }
            i += 2;
          }
          (it->second)[0] = idxMin;
          if (edgeIdx2IdxMin.count(idx) == 0)
          {
            edgeIdx2IdxMin[idx].push_back(idxMin);
            if (idx == idxMin)
              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
            edgeIdxMinList(iIdxMin) = idxMin;
            ++iIdxMin;
          }
        }
        edgeIdxMinList.resizeAndPreserve(iIdxMin);
        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();

        // (3.5) Saving variables: edge_lon, edge_lat, face_edges
        // Creating map edgeIdxGlo2Face <idxGlo, face>

        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
        edge_count = dhtEdgeIdxGlo.getIndexCount();
        edge_start = dhtEdgeIdxGlo.getIndexStart();
        edge_lon.resize(edge_count);
        edge_lat.resize(edge_count);
        edge_nodes.resize(2, edge_count);
        face_edges.resize(nvertex, nbFaces_);

        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
        size_t nEdge = 0;

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            // Getting global indexes of edge's nodes
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));

            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
            it1 = nodeIdx2IdxMin.find(myNodeIdx1);
            it2 = nodeIdx2IdxMin.find(myNodeIdx2);
            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];

            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
              size_t ownerIdx = (itIdx->second)[0];
              int edgeIdxGlo;
              size_t faceIdxGlo = nbFacesAccum + nf;

              if (myIdx == ownerIdx)
              {
                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
                edgeIdxGlo = (itIdxGlo->second)[0];
                double edgeLon;
                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
                if (diffLon < (180.- prec))
                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
                else if (diffLon > (180.+ prec))
                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
                else
                  edgeLon = 0.;
                edge_lon(edgeIdxGlo - edge_start) = edgeLon;
                edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
              }
              else
              {
                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
                edgeIdxGlo = (itIdxGlo->second)[0];
              }
              face_edges(nv1,nf) = edgeIdxGlo;
              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
              {
                edgeIdxGloList(nEdge) = edgeIdxGlo;
                ++nEdge;
              }
              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
            } // nodeIdxGlo1 != nodeIdxGlo2
            else
            {
              face_edges(nv1,nf) = 999999;
            }
          }
        }
        edgeIdxGloList.resizeAndPreserve(nEdge);

        // (3.6) Saving remaining variables edge_faces and face_faces
        edge_faces.resize(2, edge_count);
        face_faces.resize(nvertex, nbFaces_);

        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();

        for (int nf = 0; nf < nbFaces_; ++nf)
        {
          for (int nv1 = 0; nv1 < nvertex; ++nv1)
          {
            // Getting global indexes of edge's nodes
            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));

            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
            if (myNodeIdx1 != myNodeIdx2)
            {
              it1 = nodeIdx2IdxMin.find(myNodeIdx1);
              it2 = nodeIdx2IdxMin.find(myNodeIdx2);
              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
              itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
              size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
              size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];

              size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
              size_t ownerIdx = (itIdx->second)[0];
              size_t faceIdxGlo = nbFacesAccum + nf;

              if (myIdx == ownerIdx)
              {
                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
                int edgeIdxGlo = (itIdxGlo->second)[0];
                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                int face1 = it1->second[0];
                if (it1->second.size() == 1)
                {
                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
                  edge_faces(1, edgeIdxGlo - edge_start) = -999;
                  face_faces(nv1, nf) = 999999;
                }
                else
                {
                  size_t face2 = it1->second[1];
                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
                  edge_faces(1, edgeIdxGlo - edge_start) = face2;
                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
                }
              } // myIdx == ownerIdx
              else
              {
                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
                size_t edgeIdxGlo = (itIdxGlo->second)[0];
                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                int face1 = it1->second[0];
                int face2 = it1->second[1];
                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
              } // myIdx != ownerIdx
            } // myNodeIdx1 != myNodeIdx2
            else
              face_faces(nv1, nf) = 999999;
          }
        }
      }
      facesAreWritten = true;
    } // nvertex >= 3

  } // createMeshEpsilon

  ///----------------------------------------------------------------
  /*!
   * \fn  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
                                              const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                              CArray<int, 2>& nghbFaces)
   * Finds neighboring cells of a local domain for node-type of neighbors.
   * \param [in] comm
   * \param [in] face_idx Array with global indexes.
   * \param [in] bounds_lon Array of boundary longitudes.
   * \param [in] bounds_lat Array of boundary latitudes.
   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
   */

  void CMesh::getGloNghbFacesNodeType(const ep_lib::MPI_Comm& comm, const CArray<int, 1>& face_idx,
                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                               CArray<int, 2>& nghbFaces)
  {
    int nvertex = bounds_lon.rows();
    int nbFaces = bounds_lon.shape()[1];
    nghbFaces.resize(2, nbFaces*10);    // some estimate on max number of neighbouring cells

    int mpiRank, mpiSize;
    MPI_Comm_rank(comm, &mpiRank);
    MPI_Comm_size(comm, &mpiSize);

    // (1) Generating unique node indexes
    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
    vector<size_t> hashValues(4);
    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
    size_t iIdx = 0;
    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv = 0; nv < nvertex; ++nv)
      {
        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
        for (int nh = 0; nh < 4; ++nh)
        {
          if (nodeHash2Idx.count(hashValues[nh])==0)
         {
            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
            nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
           nodeHashList(iIdx) = hashValues[nh];
           ++iIdx;
         }
        }
      }
    }
    nodeHashList.resizeAndPreserve(iIdx);

    // (1.2) Generating node indexes
    // The ownership criterion: priority of the process holding the smaller index
    // Maps generated in this step are:
    // nodeHash2Info = <hash, idx1, idx2, idx3....>
    // nodeIdx2IdxMin = <idx, idxMin>
    // idxMin is a unique node identifier

    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
    dhtNodeHash.computeIndexInfoMapping(nodeHashList);
    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();

    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;

    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
    {
      size_t idxMin = (it->second)[0];
      size_t idx = (it->second)[0];
      for (int i = 2; i < (it->second).size();)
      {
        if (mpiRank == (it->second)[i+1])
        {
          idx = (it->second)[i];
        }
        if ((it->second)[i] < idxMin)
        {
          idxMin = (it->second)[i];
          (it->second)[i] = (it->second)[i-2];
          (it->second)[i+1] = (it->second)[i-1];
        }
        i += 2;
      }
      (it->second)[0] = idxMin;
      if (nodeIdx2IdxMin.count(idx) == 0)
      {
        nodeIdx2IdxMin[idx].push_back(idxMin);
      }
    }

    // (2) Creating maps nodeIdxMin2Face = <nodeIdxMin, [face1, face2, ...]>
    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdxMin2Face;
    CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4);

    size_t nNode = 0;

    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv = 0; nv < nvertex; ++nv)
      {
        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
        it = nodeIdx2IdxMin.find(myNodeIdx);
        size_t nodeIdxMin = (it->second)[0];
        size_t faceIdx = face_idx(nf);
        if (nodeIdxMin2Face.count(nodeIdxMin) == 0)
        {
          nodeIdxMinList(nNode) = nodeIdxMin;
          ++nNode;
        }
        nodeIdxMin2Face[nodeIdxMin].push_back(faceIdx);
        nodeIdxMin2Face[nodeIdxMin].push_back(mpiRank);
      }
    }
    nodeIdxMinList.resizeAndPreserve(nNode);

    // (3) Face_face connectivity

    // nodeIdxMin2Info = <nodeIdxMin, [face1, face2,...]>
    CClientClientDHTSizet dhtNode2Face (nodeIdxMin2Face, comm);
    dhtNode2Face.computeIndexInfoMapping(nodeIdxMinList);
    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2Info = dhtNode2Face.getInfoIndexMap();
    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map

    int nbNghb = 0;
    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode;

    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv = 0; nv < nvertex; ++nv)
      {
        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
        itNode = nodeIdx2IdxMin.find(myNodeIdx);
        size_t nodeIdxMin = (itNode->second)[0];

        itNode = nodeIdxMin2Info.find(nodeIdxMin);
        for (int i = 0; i < itNode->second.size();)
        {
          size_t face = itNode->second[i];
          size_t rank = itNode->second[i+1];
          if (rank != mpiRank)
            if (mapFaces.count(face) == 0)
            {
              nghbFaces(0, nbNghb) = face;
              nghbFaces(1, nbNghb) = rank;
              ++nbNghb;
              mapFaces[face].push_back(face);
            }
          i += 2;
        }
      }
    }
    nghbFaces.resizeAndPreserve(2, nbNghb);
  } // getGloNghbFacesNodeType

  ///----------------------------------------------------------------
  /*!
   * \fn  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                               CArray<int, 2>& nghbFaces)
   * Finds neighboring cells of a local domain for edge-type of neighbors.
   * \param [in] comm
   * \param [in] face_idx Array with global indexes.
   * \param [in] bounds_lon Array of boundary longitudes.
   * \param [in] bounds_lat Array of boundary latitudes.
   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
   */

  void CMesh::getGloNghbFacesEdgeType(const ep_lib::MPI_Comm& comm, const CArray<int, 1>& face_idx,
                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                               CArray<int, 2>& nghbFaces)
  {
    int nvertex = bounds_lon.rows();
    int nbFaces = bounds_lon.shape()[1];
    nghbFaces.resize(2, nbFaces*10);    // estimate of max number of neighbouring cells

    int mpiRank, mpiSize;
    MPI_Comm_rank(comm, &mpiRank);
    MPI_Comm_size(comm, &mpiSize);

    // (1) Generating unique node indexes
    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
    vector<size_t> hashValues(4);
    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
    size_t iIdx = 0;
    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv = 0; nv < nvertex; ++nv)
      {
        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
        for (int nh = 0; nh < 4; ++nh)
        {
          if (nodeHash2Idx.count(hashValues[nh])==0)
          {
            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
            nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
            nodeHashList(iIdx) = hashValues[nh];
            ++iIdx;
          }
        }
      }
    }
    nodeHashList.resizeAndPreserve(iIdx);

    // (1.2) Generating node indexes
    // The ownership criterion: priority of the process holding the smaller index
    // Maps generated in this step are:
    // nodeHash2Info = <hash, idx1, idx2, idx3....>
    // nodeIdx2IdxMin = <idx, idxMin>
    // idxMin is a unique node identifier

    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
    dhtNodeHash.computeIndexInfoMapping(nodeHashList);
    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();

    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;

    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
    {
      size_t idxMin = (it->second)[0];
      size_t idx = (it->second)[0];
      for (int i = 2; i < (it->second).size();)
      {
        if (mpiRank == (it->second)[i+1])
        {
          idx = (it->second)[i];
        }
        if ((it->second)[i] < idxMin)
        {
          idxMin = (it->second)[i];
          (it->second)[i] = (it->second)[i-2];
          (it->second)[i+1] = (it->second)[i-1];
        }
        i += 2;
      }
      (it->second)[0] = idxMin;
      if (nodeIdx2IdxMin.count(idx) == 0)
      {
        nodeIdx2IdxMin[idx].push_back(idxMin);
      }
    }

    // (2) Creating map edgeHash2Face = <edgeHash, [[face1, rank1], [face2, rank2]]>, where rank1 = rank2 = ...

    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2, it;
    CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Face;
    CArray<size_t,1> edgeHashList(nbFaces*nvertex);

    size_t nEdge = 0;

    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv1 = 0; nv1 < nvertex; ++nv1)
      {
        // Getting indexes of edge's nodes
        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
        it1 = nodeIdx2IdxMin.find(myNodeIdx1);
        it2 = nodeIdx2IdxMin.find(myNodeIdx2);
        size_t nodeIdxMin1 = (it1->second)[0];
        size_t nodeIdxMin2 = (it2->second)[0];
        size_t faceIdx = face_idx(nf);

        if (nodeIdxMin1 != nodeIdxMin2)
        {
          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
          if (edgeHash2Face.count(edgeHash) == 0)
          {
            edgeHashList(nEdge) = edgeHash;
            ++nEdge;
          }
          edgeHash2Face[edgeHash].push_back(faceIdx);
          edgeHash2Face[edgeHash].push_back(mpiRank);
        } // nodeIdxMin1 != nodeIdxMin2
      }
    }
    edgeHashList.resizeAndPreserve(nEdge);

    // (3) Face_face connectivity

    int nbNghb = 0;
    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode1, itNode2;

    // edgeHash2Info = <edgeHash, [[face1, rank1], [face2, rank2]]>
    CClientClientDHTSizet dhtEdge2Face (edgeHash2Face, comm);
    dhtEdge2Face.computeIndexInfoMapping(edgeHashList);
    CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdge2Face.getInfoIndexMap();
    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map

    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv1 = 0; nv1 < nvertex; ++nv1)
      {
        // Getting indexes of edge's nodes
        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));

        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
        itNode1 = nodeIdx2IdxMin.find(myNodeIdx1);
        itNode2 = nodeIdx2IdxMin.find(myNodeIdx2);
        size_t nodeIdxMin1 = (itNode1->second)[0];
        size_t nodeIdxMin2 = (itNode2->second)[0];

        if (nodeIdxMin1 != nodeIdxMin2)
        {
          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
          it1 = edgeHash2Info.find(edgeHash);

          for (int i = 0; i < it1->second.size();)
          {
            size_t face = it1->second[i];
            size_t rank = it1->second[i+1];
            if (rank != mpiRank)
              if (mapFaces.count(face) == 0)
              {
                nghbFaces(0, nbNghb) = face;
                nghbFaces(1, nbNghb) = rank;
                ++nbNghb;
                mapFaces[face].push_back(face);
              }
            i += 2;
          }
        } // nodeIdxMin1 != nodeIdxMin2
      }
    }
    nghbFaces.resizeAndPreserve(2, nbNghb);
  } // getGloNghbFacesEdgeType

  ///----------------------------------------------------------------
  /*!
   * \fn void getGlobalNghbFaces (const int nghbType, const MPI_Comm& comm, const CArray<int, 1>& face_idx,
                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                  CArray<size_t, 1>& nghbFaces)
   * Finds neighboring faces owned by other procs.
   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
   * \param [in] comm
   * \param [in] face_idx Array with global indexes.
   * \param [in] bounds_lon Array of boundary longitudes.
   * \param [in] bounds_lat Array of boundary latitudes.
   * \param [out] nghbFaces 2D array containing neighboring faces and owner ranks.
   */

  void CMesh::getGlobalNghbFaces(const int nghbType, const ep_lib::MPI_Comm& comm,
                                 const CArray<int, 1>& face_idx,
                                 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                 CArray<int, 2>& nghbFaces)
  {
    if (nghbType == 0)
      getGloNghbFacesNodeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
    else
      getGloNghbFacesEdgeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
  } // getGlobalNghbFaces

  ///----------------------------------------------------------------
  /*!
   * \fn void getLocalNghbFaces (const int nghbType,
                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                  CArray<size_t, 1>& nghbFaces)
   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
   * \param [in] bounds_lon Array of boundary longitudes.
   * \param [in] bounds_lat Array of boundary latitudes.
   * \param [out] nghbFaces 1D array containing neighboring faces.
   */

  void CMesh::getLocalNghbFaces(const int nghbType, const CArray<int, 1>& face_idx,
                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
  {
    if (nghbType == 0)
      getLocNghbFacesNodeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
    else
      getLocNghbFacesEdgeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
  } // getLocalNghbFaces

  ///----------------------------------------------------------------
  /*!
   * \fn void getLocNghbFacesNodeType (const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                       CArray<int, 2>& nghbFaces)
   * \param [in] face_idx Array with local face indexing.
   * \param [in] bounds_lon Array of boundary longitudes.
   * \param [in] bounds_lat Array of boundary latitudes.
   * \param [out] nghbFaces 2D array containing neighboring faces.
   * \param [out] nbNghbFaces Array containing number of neighboring faces.
   */

  void CMesh::getLocNghbFacesNodeType (const CArray<int, 1>& face_idx,
                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
  {
    int nvertex = bounds_lon.rows();
    int nbFaces = bounds_lon.shape()[1];
    int nbNodes = 0;
    nbNghbFaces.resize(nbFaces);
    nbNghbFaces = 0;

    // nodeToFaces connectivity
    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces;
    for (int nf = 0; nf < nbFaces; ++nf)
      for (int nv = 0; nv < nvertex; ++nv)
      {
        size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0];
        nodeToFaces[nodeHash].push_back(face_idx(nf));
      }

    // faceToFaces connectivity
    boost::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
    int maxNb = 20;                            // some assumption on the max possible number of neighboring cells
    faceToFaces.resize(maxNb, nbFaces);
    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
    for (it = nodeToFaces.begin(); it != nodeToFaces.end(); ++it)
    {
      int size = it->second.size();
      for (int i = 0; i < (size-1); ++i)
      {
        int face1 = it->second[i];
        for (int j = i+1; j < size; ++j)
        {
          int face2 = it->second[j];
          if (face2 != face1)
          {
            int hashFace = hashPairOrdered(face1, face2);
            if (mapFaces.count(hashFace) == 0)
            {
              faceToFaces(nbNghbFaces(face1), face1) = face2;
              faceToFaces(nbNghbFaces(face2), face2) = face1;
              ++nbNghbFaces(face1);
              ++nbNghbFaces(face2);
              mapFaces[hashFace] = hashFace;
            }
          }
        }
      }
    }
  } //getLocNghbFacesNodeType


  ///----------------------------------------------------------------
  /*!
   * \fn void getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
   *                                   const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
   *                                   CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
   * \param [in] face_idx Array with local face indexing.
   * \param [in] bounds_lon Array of boundary longitudes.
   * \param [in] bounds_lat Array of boundary latitudes.
   * \param [out] nghbFaces 2D array containing neighboring faces.
   * \param [out] nbNghbFaces Array containing number of neighboring faces.
   */

  void CMesh::getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
  {
    int nvertex = bounds_lon.rows();
    int nbFaces = bounds_lon.shape()[1];
    int nbNodes = 0;
    int nbEdges = 0;
    nbNghbFaces.resize(nbFaces);
    nbNghbFaces = 0;

    // faceToNodes connectivity
    CArray<double, 2> faceToNodes (nvertex, nbFaces);

    boost::unordered_map <pair<double,double>, int> mapNodes;

    for (int nf = 0; nf < nbFaces; ++nf)
      for (int nv = 0; nv < nvertex; ++nv)
      {
        if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end())
        {
          mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes;
          faceToNodes(nv,nf) = nbNodes ;
          ++nbNodes ;
        }
        else
          faceToNodes(nv,nf) = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
      }

    // faceToFaces connectivity
    boost::unordered_map <pair<int,int>, int> mapEdges;
    faceToFaces.resize(nvertex, nbFaces);
    CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible

    for (int nf = 0; nf < nbFaces; ++nf)
    {
      for (int nv1 = 0; nv1 < nvertex; ++nv1)
      {
        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
        int face = face_idx(nf);
        int node1 = faceToNodes(nv1,face);
        int node2 = faceToNodes(nv2,face);
        if (node1 != node2)
        {
          if (mapEdges.find(make_ordered_pair (node1, node2)) == mapEdges.end())
          {
            mapEdges[make_ordered_pair (node1, node2)] = nbEdges;
            edgeToFaces(0,nbEdges) = face;
            ++nbEdges;
          }
          else
          {
            int edge = mapEdges[make_ordered_pair (node1, node2)];
            edgeToFaces(1, edge) = face;
            int face1 = face;
            int face2 = edgeToFaces(0,edge);
            faceToFaces(nbNghbFaces(face1), face1) =  face2;
            faceToFaces(nbNghbFaces(face2), face2) =  face1;
            ++nbNghbFaces(face1);
            ++nbNghbFaces(face2);
          }
        } // node1 != node2
      } // nv
    } // nf

  } //getLocNghbFacesEdgeType


} // namespace xios