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

#include "mesh.hpp"
#include <boost/functional/hash.hpp>
//#include <unordered_map>

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 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 generateNodeIndex(vector<size_t>& valList)
   * Generates a node index unique for all processes.
   * \param [in] valList Vector storing four node hashes.
   */
    size_t generateNodeIndex(vector<size_t>& valList)
    {
      // 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 = vec[0] ;
      int it = 1;
      for(; it != vec.size(); ++it)
      {
         seed = hashPair(seed, vec[it]);
      }
      return seed ;
    }

///----------------------------------------------------------------
/*!
 * \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 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
      unsigned long nbEdgesOnProc = nbEdges_;
      unsigned long 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);
        int nbHash = 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));
            for (int nh = 0; nh < 4; ++nh)
            {
              if (nodeHash2Idx[hashValues[nh]].size() == 0)
              {
                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues));
                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
                nodeHashList(nbHash) = hashValues[nh];
                ++nbHash;
              }
            }
          }
        }
        nodeHashList.resizeAndPreserve(nbHash);

        // (2.2) Generating global node indexes
        // The ownership criterion: priority of the process of smaller index
        // Maps generated in this step are:
        // Maps generated in this step are:
         // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
         // nodeIdx2Idx = <idx, <rankOwner, idx>>

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


         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
         CArray<size_t,1> nodeIdxList(nbEdges_*nvertex*4);
         size_t nIdx = 0;

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

         // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
         // Solution: global node indexing by hand.
         // Maps modified in this step:
         // nodeIdx2Idx = <idx, idxGlo>
         unsigned long nodeCount = nodeIdx2Idx.size();
         unsigned long nodeStart, nbNodes;
         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
         int nNodes = nodeStart;
         MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
         nbNodesGlo = nNodes;

         nodeStart -= nodeCount;
         node_start = nodeStart;
         node_count = nodeCount;
         CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
         size_t count = 0;

         for (int ne = 0; ne < nbEdges_; ++ne)
         {
           for (int nv = 0; nv < nvertex; ++nv)
           {
             vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
             size_t nodeIdx = generateNodeIndex(hashValues);
             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
             if (it != nodeIdx2Idx.end())
             {
               if (dummyMap.count(nodeIdx) == 0)
               {
                 dummyMap[nodeIdx].push_back(nodeIdx);
                 (it->second)[1] = node_start + count;
                 ++count;
               }
             }
           }
         }

         CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
         dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();

        // (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);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxGlo.find(myIdx);
            idxGlo = (itIdx->second)[1];

            if (mpiRank == (itIdx->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);
            }
            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
      unsigned long nbFacesOnProc = nbFaces_;
      unsigned long 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> edgeIdxList(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)
              {
                edgeIdxList(iIdx) = edgeIdxGlo;
                ++iIdx;
              }
              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
              edgeHash2Rank[edgeHash].push_back(mpiRank);
              edgeHash2Rank[edgeHash].push_back(itEdgeHash->second[0]);
            }
            else
            {
              face_edges(nv1,nf) = 999999;
            }
          }
        }
        edgeIdxList.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();

        // edgeHash2Info = <edgeHash, < rank1, idxGlo, rank2, idxGlo>>
        int edgeCount = 0;
        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
        {
          vector <size_t> edgeInfo = it->second;
          if (edgeInfo[0] == mpiRank)
          {
            ++edgeCount;
          }
        }

        unsigned long edgeStart, nbEdges;
        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
        int nEdges = edgeStart;
        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
        nbEdgesGlo = nEdges;

        // edges to be splitted equally between procs
        if ( (nbEdgesGlo % mpiSize) == 0)
        {
          edge_count = nbEdgesGlo/mpiSize;
          edge_start = mpiRank*edge_count;
        }
        else
        {
          if (mpiRank == (mpiSize - 1) )
          {
            edge_count = nbEdgesGlo/mpiSize;
            edge_start = mpiRank*(nbEdgesGlo/mpiSize + 1);
          }
          else
          {
            edge_count = nbEdgesGlo/mpiSize + 1;
            edge_start = mpiRank*edge_count;
          }
        }
        CArray<size_t,1> edgeIdxGloList(edge_count);
        for (int i = 0; i < edge_count; ++i)
        {
          edgeIdxGloList(i) = i + edge_start;
        }

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


        edge_faces.resize(2, edge_count);
        for (int i = 0; i < edge_count; ++i)
        {
          CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdxGlo2FaceIdx.find(i + edge_start);
          int indexGlo = it->first;
          vector<size_t> faces = it->second;
          int face1 = faces[0];
          edge_faces(0, indexGlo - edge_start) = face1;
          if (faces.size() == 2)
          {
            int face2 = faces[1];
            edge_faces(1, indexGlo - edge_start) = face2;
          }
          else
          {
            edge_faces(1, indexGlo - edge_start) = -999;
          }
        }

        size_t tmp;
        vector<size_t> tmpVec;
        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2FaceIdx.begin(); it != edgeIdx2FaceIdx.end(); it++)
        {
          tmp = it->first;
          tmpVec = it->second;
          tmp++;
        }

        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2, itIndex;
        face_faces.resize(nvertex, nbFaces_);
        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 = edgeHash2Info.find(edgeHash);
              int edgeIdxGlo = (itEdgeHash->second)[1];

              if ( (itEdgeHash->second)[0] == mpiRank)
              {
                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
                int face1 = itFace1->second[0];
                if (itFace1->second.size() == 1)
                {
                  face_faces(nv1, nf) = 999999;
                }
                else
                {
                  int face2 = itFace1->second[1];
                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
                }
              } // edge owner
              else
              {
                itFace1 = edgeIdx2FaceIdx.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);
        int nEdgeHash = 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];
            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
            if (edgeHash2Idx.count(edgeHash) == 0)
            {
              edgeHash2Idx[edgeHash].push_back(edgeHash);
              edgeHash2Idx[edgeHash].push_back(mpiRank);
              edgeHashList(nEdgeHash) = edgeHash;
              ++nEdgeHash;
            }
          }
        }
        edgeHashList.resizeAndPreserve(nEdgeHash);

        // (2.3) Generating global edge indexes
        // The ownership criterion: priority of the process with smaller rank
        // Maps generated in this step are:
        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>

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

        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;

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

          for (int i = 2; i < (it->second).size();)
          {
            if ((it->second)[i+1] < rankMin)
            {
              rankMin = (it->second)[i+1];
              idx = (it->second)[i];
              (it->second)[i+1] = (it->second)[i-1];
            }
            i += 2;
          }
          if (edgeIdx2Idx.count(idx) == 0)
          {
            if (mpiRank == rankMin)
            {
              edgeIdx2Idx[idx].push_back(rankMin);
              edgeIdx2Idx[idx].push_back(idx);
            }
          }
        }

        unsigned long edgeCount = edgeIdx2Idx.size();
        unsigned long edgeStart, nbEdges;
        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
        int nEdges = edgeStart;
        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
        nbEdgesGlo = nEdges;

        edgeStart -= edgeCount;
        edge_start = edgeStart;
        edge_count = edgeCount;
        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
        int count = 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));
            }
            size_t nodeIdxGlo1 = it1->second[0];
            size_t nodeIdxGlo2 = it2->second[0];

            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
              if (it != edgeIdx2Idx.end())
              {
                if (dummyEdgeMap.count(edgeIdx) == 0)
                {
                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
                  (it->second)[1] = edge_start + count;
                  ++count;
                }
              }
            }
          }
        }

        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();

        // (2.4) Saving variables: edge_lon, edge_lat, face_edges
        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 = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
              int edgeIdxGlo = (itIdx->second)[1];
              size_t faceIdxGlo = nbFacesAccum + nf;

              if (mpiRank == (itIdx->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;
              }
              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.5) 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;
        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx;

        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 = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
            itIdx = edgeIdx2IdxGlo.find(myIdx);
            size_t faceIdxGlo = nbFacesAccum + nf;
            int edgeIdxGlo = (itIdx->second)[1];

            if (mpiRank == (itIdx->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
            {
              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);
            size_t nodeIndex = generateNodeIndex(hashValues);
            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 with smaller rank.
        // With any other criterion it is not possible to have consistent node indexing for different number of procs.
        // Maps generated in this step are:
        // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
        // nodeIdx2Idx = <idx, <rankOwner, idx>>

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

        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
        CArray<size_t,1> nodeIdxList(nbFaces_*nvertex*4);
        size_t nIdx = 0;

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

//        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2Idx, comm);
        // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
        // Solution: global node indexing by hand.
        // Maps modified in this step:
        // nodeIdx2Idx = <idx, idxGlo>
        unsigned long nodeCount = nodeIdx2Idx.size();
        unsigned long nodeStart, nbNodes;
        MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
        int nNodes = nodeStart;
        MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
        nbNodesGlo = nNodes;

        nodeStart -= nodeCount;
        node_start = nodeStart;
        node_count = nodeCount;
        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
        size_t count = 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 nodeIdx = generateNodeIndex(hashValues);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
            if (it != nodeIdx2Idx.end())
            {
              if (dummyMap.count(nodeIdx) == 0)
              {
                dummyMap[nodeIdx].push_back(nodeIdx);
                (it->second)[1] = node_start + count;
                ++count;
              }
            }
          }
        }
        nodeIdxList.resizeAndPreserve(nIdx);
        CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
        dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = 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 nodeIdx1 = generateNodeIndex(hashValues1);
            size_t nodeIdx2 = generateNodeIndex(hashValues2);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
            size_t ownerRank = (itNodeIdx1->second)[0];
            nodeIdxGlo1 = (itNodeIdx1->second)[1];
            nodeIdxGlo2 = (itNodeIdx2->second)[1];

            if (mpiRank == ownerRank)
            {
              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
            }
            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
              edgeHash2Idx[edgeHash].push_back(edgeHash);
              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:
        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>

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

        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;

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

          for (int i = 2; i < (it->second).size();)
          {
            if ((it->second)[i+1] < rankMin)
            {
              rankMin = (it->second)[i+1];
              idx = (it->second)[i];
              (it->second)[i+1] = (it->second)[i-1];
            }
            i += 2;
          }
          if (edgeIdx2Idx.count(idx) == 0)
          {
            if (mpiRank == rankMin)
            {
              edgeIdx2Idx[idx].push_back(rankMin);
              edgeIdx2Idx[idx].push_back(idx);
            }
          }
        }

        unsigned long edgeCount = edgeIdx2Idx.size();
        unsigned long edgeStart, nbEdges;
        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
        int nEdges = edgeStart;
        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
        nbEdgesGlo = nEdges;

        edgeStart -= edgeCount;
        edge_start = edgeStart;
        edge_count = edgeCount;
        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
        count = 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 nodeIdx1 = generateNodeIndex(hashValues1);
            size_t nodeIdx2 = generateNodeIndex(hashValues2);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
            nodeIdxGlo1 = (itNodeIdx1->second)[1];
            nodeIdxGlo2 = (itNodeIdx2->second)[1];

            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
              if (it != edgeIdx2Idx.end())
              {
                if (dummyEdgeMap.count(edgeIdx) == 0)
                {
                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
                  (it->second)[1] = edge_start + count;
                  ++count;
                }
              }
            }
          }
        }
        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = 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 nodeIdx1 = generateNodeIndex(hashValues1);
            size_t nodeIdx2 = generateNodeIndex(hashValues2);
            it1 = nodeIdx2IdxGlo.find(nodeIdx1);
            it2 = nodeIdx2IdxGlo.find(nodeIdx2);
            size_t nodeIdxGlo1 = (it1->second)[1];
            size_t nodeIdxGlo2 = (it2->second)[1];

            if (nodeIdxGlo1 != nodeIdxGlo2)
            {
              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
              int edgeIdxGlo = (itIdx->second)[1];
              size_t faceIdxGlo = nbFacesAccum + nf;

              if (mpiRank == (itIdx->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;
              }
              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);
            size_t myNodeIdx2 = generateNodeIndex(hashValues2);
            if (myNodeIdx1 != myNodeIdx2)
            {
              it1 = nodeIdx2IdxGlo.find(myNodeIdx1);
              it2 = nodeIdx2IdxGlo.find(myNodeIdx2);
              size_t nodeIdxGlo1 = (it1->second)[1];
              size_t nodeIdxGlo2 = (it2->second)[1];
              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
              int edgeIdxGlo = (itIdx->second)[1];

              size_t faceIdxGlo = nbFacesAccum + nf;

              if (mpiRank == (itIdx->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);
                }
              }
              else
              {
                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
                int face1 = it1->second[0];
                int face2 = it1->second[1];
                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
              }
            } // 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 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 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 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
    std::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);

    std::unordered_map <pairDouble, int, boost::hash<pairDouble> > 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
    std::unordered_map <pairInt, int, boost::hash<pairInt> > 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