Changeset 924 for XIOS/trunk/src/node


Ignore:
Timestamp:
09/02/16 09:41:27 (8 years ago)
Author:
oabramkina
Message:

Parallel version of UGRID norms.

The following connectivity parameters have been implemented:

edge_nodes
face_nodes
face_edges
edge_faces
face_faces.

Test with a regular(structured) quadrilateral mesh (test_regular) has been added.

Location:
XIOS/trunk/src/node
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/src/node/domain.cpp

    r911 r924  
    3636      , isRedistributed_(false) 
    3737   { 
    38            //mesh = new CMesh(); 
    39         } 
     38   } 
    4039 
    4140   CDomain::CDomain(const StdString & id) 
    … …  
    4746      , isRedistributed_(false) 
    4847   { 
    49            //mesh = new CMesh(); 
    50         } 
     48         } 
    5149 
    5250   CDomain::~CDomain(void) 
    5351   { 
    54            //delete mesh; 
    5552   } 
    5653 
    5754   ///--------------------------------------------------------------- 
    5855 
    59    void CDomain::assignMesh(const StdString meshName) 
    60    { 
    61      mesh = CMesh::getMesh(meshName); 
     56   void CDomain::assignMesh(const StdString meshName, const int nvertex) 
     57   { 
     58     mesh = CMesh::getMesh(meshName, nvertex); 
    6259   } 
    6360 
    … …  
    18971894 
    18981895    buffer >> lon; 
     1896 
    18991897    if (hasBounds) buffer >> boundslon; 
    19001898 

  • XIOS/trunk/src/node/domain.hpp

    r893 r924  
    6767          
    6868         CMesh* mesh; 
    69          void assignMesh(const StdString); 
     69         void assignMesh(const StdString, const int); 
    7070         
    7171         virtual void parse(xml::CXMLNode & node); 

  • XIOS/trunk/src/node/mesh.cpp

    r901 r924  
    1111/// ////////////////////// Définitions ////////////////////// /// 
    1212 
    13   CMesh::CMesh(void) :  nbFaces{0}, nbNodes{0}, nbEdges{0} 
    14             ,   nvertex{0} 
     13  CMesh::CMesh(void) :  nbNodesGlo{0}, nbEdgesGlo{0} 
     14            ,  node_start{0}, node_count{0} 
     15            ,  edge_start{0}, edge_count{0} 
     16            ,  nbFaces{0}, nbNodes{0}, nbEdges{0} 
    1517            ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false} 
    1618            ,  node_lon(), node_lat() 
    … …  
    1820            ,  face_lon(), face_lat() 
    1921            ,  face_nodes() 
     22            ,  pNodeGlobalIndex{NULL}, pEdgeGlobalIndex{NULL} 
    2023  { 
    2124  } 
    … …  
    2427  CMesh::~CMesh(void) 
    2528  { 
     29    if (pNodeGlobalIndex != NULL) delete pNodeGlobalIndex; 
     30    if (pEdgeGlobalIndex != NULL) delete pEdgeGlobalIndex; 
    2631  } 
    2732 
    2833  std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>(); 
    29   //CMesh* CMesh::getMesh; 
     34  std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >(); 
    3035 
    3136///--------------------------------------------------------------- 
    … …  
    3439 * 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. 
    3540 * \param [in] meshName  The name of a mesh ("name" attribute of a domain). 
     41 * \param [in] nvertex Number of verteces (1 for nodes, 2 for edges, 3 and up for faces). 
    3642 */ 
    37   CMesh* CMesh::getMesh (StdString meshName) 
     43  CMesh* CMesh::getMesh (StdString meshName, int nvertex) 
    3844  { 
     45    CMesh::domainList[meshName].push_back(nvertex); 
     46 
    3947    if ( CMesh::meshList.begin() != CMesh::meshList.end() ) 
    4048    { 
    … …  
    6068 
    6169///---------------------------------------------------------------- 
    62   int hashPair(int first, int second) 
     70  size_t hashPair(size_t first, size_t second) 
    6371  { 
    64     int seed = first + 0x9e3779b9 ; 
    65     seed ^= second + 0x9e3779b9 + (seed << 6) + (seed >> 2); 
     72    HashXIOS<size_t> sizetHash; 
     73    size_t seed = sizetHash(first) + 0x9e3779b9 ; 
     74    seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2); 
     75    return seed ; 
     76  } 
     77 
     78///---------------------------------------------------------------- 
     79  size_t hashPairOrdered(size_t first, size_t second) 
     80  { 
     81    size_t seed; 
     82    HashXIOS<size_t> sizetHash; 
     83    if (first < second) 
     84    { 
     85      seed = sizetHash(first) + 0x9e3779b9 ; 
     86      seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2); 
     87    } 
     88    else 
     89    { 
     90      seed = sizetHash(second) + 0x9e3779b9 ; 
     91      seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2); 
     92    } 
     93    return seed ; 
     94  } 
     95 
     96///---------------------------------------------------------------- 
     97/*! 
     98 * \fn size_t generateEdgeIndex(size_t first, size_t second, int rank) 
     99 * Generates an edge index. 
     100 * If the same edge is generated by two processes, each process will have its own edge index. 
     101 * \param [in] first Edge node. 
     102 * \param [in] second Edge node. 
     103 * \param [in] rank MPI process rank. 
     104 */ 
     105  size_t generateEdgeIndex(size_t first, size_t second, int rank) 
     106  { 
     107    size_t seed = rank ; 
     108    if (first < second) 
     109    { 
     110      seed = hashPair(seed, first); 
     111      seed = hashPair(seed, second); 
     112    } 
     113    else 
     114    { 
     115      seed = hashPair(seed, second); 
     116      seed = hashPair(seed, first); 
     117    } 
     118 
     119    return seed ; 
     120  } 
     121 
     122///---------------------------------------------------------------- 
     123/*! 
     124 * \fn size_t generateNodeIndex(vector<size_t>& valList, int rank) 
     125 * Generates a node index. 
     126 * If the same node is generated by two processes, each process will have its own node index. 
     127 * \param [in] valList Vector storing four node hashes. 
     128 * \param [in] rank MPI process rank. 
     129 */ 
     130  size_t generateNodeIndex(vector<size_t>& valList, int rank) 
     131  { 
     132    // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere 
     133    vector<size_t> vec = valList; 
     134    sort (vec.begin(), vec.end()); 
     135    size_t seed = rank ; 
     136    int it = 0; 
     137    for(; it != vec.size(); ++it) 
     138    { 
     139       seed = hashPair(seed, vec[it]); 
     140    } 
    66141    return seed ; 
    67142  } 
    … …  
    144219///---------------------------------------------------------------- 
    145220/*! 
    146  * \fn CArray<size_t,1>& CMesh::createHashes (double lon, double lat) 
     221 * \fn CArray<size_t,1>& CMesh::createHashes (const double longitude, const double latitude) 
    147222 * Creates two hash values for each dimension, longitude and latitude. 
    148  * \param [in] lon Node longitude in degrees. 
    149  * \param [in] lat Node latitude in degrees ranged from 0 to 360. 
     223 * \param [in] longitude Node longitude in degrees. 
     224 * \param [in] latitude Node latitude in degrees ranged from 0 to 360. 
    150225 */ 
    151226 
    152   vector<size_t> CMesh::createHashes (double lon, double lat) 
     227  vector<size_t> CMesh::createHashes (const double longitude, const double latitude) 
    153228  { 
    154229    double minBoundLon = 0. ; 
    155230    double maxBoundLon = 360. ; 
    156     double minBoundLat = -90 ; 
    157     double maxBoundLat = 90 ; 
     231    double minBoundLat = -90. ; 
     232    double maxBoundLat = 90. ; 
    158233    double prec=1e-11 ; 
    159234    double precLon=prec ; 
    160235    double precLat=prec ; 
     236    double lon = longitude; 
     237    double lat = latitude; 
     238 
     239    if (lon > (360.- prec)) lon = 0.; 
    161240 
    162241    size_t maxsize_t=numeric_limits<size_t>::max() ; 
    … …  
    213292      return std::pair<int,int>(b,a); 
    214293  } 
    215  
    216 ///---------------------------------------------------------------- 
    217 //#include <random> 
    218 //  size_t randNum() 
    219 //  { 
    220 //    typedef mt19937 rng; 
    221 //    size_t max = numeric_limits<size_t>::max(); 
    222 //    uniform_int_distribution<> distr(0, max); 
    223 //    return distr(rng); 
    224 //  } 
    225294 
    226295///---------------------------------------------------------------- 
    … …  
    237306            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat) 
    238307  { 
    239     nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     308    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
    240309 
    241310    if (nvertex == 1) 
    … …  
    357426            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 
    358427                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) : 
    359                         (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);; 
     428                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.); 
    360429            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 
    361430            ++nbEdges; 
    … …  
    418487 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces. 
    419488 * Precision check is implemented with two hash values for each dimension, longitude and latitude. 
     489 * \param [in] comm 
    420490 * \param [in] lonvalue  Array of longitudes. 
    421491 * \param [in] latvalue  Array of latitudes. 
    … …  
    423493 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type. 
    424494 */ 
    425   void CMesh::createMeshEpsilon(const MPI_Comm& comm, const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue, 
    426             const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat) 
     495  void CMesh::createMeshEpsilon(const MPI_Comm& comm, 
     496                                const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue, 
     497                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat) 
    427498  { 
    428499 
    429     nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     500    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     501    int mpiRank, mpiSize; 
     502    MPI_Comm_rank(comm, &mpiRank); 
     503    MPI_Comm_size(comm, &mpiSize); 
    430504 
    431505    if (nvertex == 1) 
    432506    { 
    433507      nbNodes = lonvalue.numElements(); 
    434  
    435 //      CClientClientDHTSizet::Index2VectorInfoTypeMap hash2Idx;      // map <hash, idx> 
    436       vector<size_t> hashValues(4);                                 // temporary vector for storing hashes for each node 
    437       vector<size_t> globalIndexes(nbNodes);                        // Probably a map ? 
    438       CArray<size_t,1> hashList(nbNodes*4);                         // Do I need it? 
    439       CArray<size_t,1> idxList(nbNodes); 
    440  
    441       // Assign a unique index for each node represented by four hashes 
    442       // The first hash will serve as the unique index 
    443       size_t randomIdx; 
    444       int rank, size; 
    445       MPI_Comm_rank(comm, &rank); 
    446       MPI_Comm_size(comm, &size); 
    447       srand((unsigned)time(NULL)+rank*size); 
    448  
    449       for (size_t nn = 0; nn < nbNodes; ++nn) 
    450       { 
    451         hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn)); 
    452         idxList(nn) = hashValues[nn]; 
    453 //        randomIdx = rand() % numeric_limits<size_t>::max(); 
    454 //        idxList(nn) = randomIdx; 
    455 //        for (size_t nh = 0; nh < 4; ++nh) 
    456 //        { 
    457 //          hash2Idx[hashValues[nh]].push_back(randomIdx); 
    458 //          hashList(nn*4 + nh) = hashValues[nh]; 
    459 //        } 
    460       } 
    461  
    462  
    463 //      CClientClientDHTSizet dhtSizet(hash2Idx, comm); 
    464 //      dhtSizet.computeIndexInfoMapping(hashList); 
    465  //     CClientClientDHTSizet::Index2VectorInfoTypeMap& hashIdxList = dhtSizet.getInfoIndexMap(); 
    466  
    467 //      // (2.1) Create a list of indexes for each hush 
    468 //      CClientClientDHTSizet dhtSizet(hash2Idx, comm); 
    469 //      dhtSizet.computeIndexInfoMapping(hashList); 
    470 //      CClientClientDHTSizet::Index2VectorInfoTypeMap& hashIdxList = dhtSizet.getInfoIndexMap(); 
    471 // 
    472 //      // (2.2) If more than one index assigned to a hush, set the larger value to be equal to the smaller 
    473 //      int iidx = 0; 
    474 //      for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = hashIdxList.begin(); it != hashIdxList.end(); ++it, ++iidx) 
    475 //      { 
    476 //        vector<size_t> tmp = it->second; 
    477 //        size_t idxMinValue = it->second[0]; 
    478 //        for (int i = 1; i < it->second.size(); ++i) 
    479 //        { 
    480 //          if (it->second[i] < idxMinValue ) 
    481 //              idxMinValue = it->second[i]; 
    482 //        } 
    483 //        if ( (iidx % 4) == 0 ) 
    484 //          idxList(iidx/4) = idxMinValue; 
    485 //      } 
    486  
    487       // Unique global indexing 
    488 //      CDHTAutoIndexing dhtSizetIdx = CDHTAutoIndexing(idxList, comm); 
    489 //      CClientClientDHTSizet* pDhtSizet = &dhtSizetIdx; 
    490 //      pDhtSizet->computeIndexInfoMapping(idxList); 
    491  
    492       //globalIndexes = dhtSizetIdx.getGlobalIndex(); 
    493  
    494508      node_lon.resize(nbNodes); 
    495509      node_lat.resize(nbNodes); 
    … …  
    497511      node_lat = latvalue; 
    498512 
     513      // Global node indexes 
     514      vector<size_t> hashValues(4); 
     515      CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo; 
     516      for (size_t nn = 0; nn < nbNodes; ++nn) 
     517      { 
     518        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn)); 
     519        for (size_t nh = 0; nh < 4; ++nh) 
     520        { 
     521          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes + nn); 
     522        } 
     523      } 
     524      pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm); 
    499525      nodesAreWritten = true; 
    500  
    501526    } 
    502527 
    … …  
    504529    { 
    505530      nbEdges = bounds_lon.shape()[1]; 
    506  
    507       vector<size_t> hashValues(4); 
    508       for (int ne = 0; ne < nbEdges; ++ne) 
     531      edge_lon.resize(nbEdges); 
     532      edge_lat.resize(nbEdges); 
     533      edge_lon = lonvalue; 
     534      edge_lat = latvalue; 
     535      edge_nodes.resize(nvertex, nbEdges); 
     536      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo; 
     537      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 
     538 
     539      // Case (1): node indexes have been generated by domain "nodes" 
     540      if (nodesAreWritten) 
    509541      { 
    510         for (int nv = 0; nv < nvertex; ++nv) 
    511         { 
    512           hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne)); 
    513 //          for (size_t nh = 0; nh < 4; ++nh) 
    514 //          { 
    515 //            hash2Idx[hashValues[nh]].push_back(randomIdx); 
    516 //            hashList(nn*4 + nh) = hashValues[nh]; 
    517 //          } 
    518  
    519  
     542        vector<size_t> hashValues(4); 
     543        CArray<size_t,1> nodeHashList(nbEdges*nvertex*4); 
     544        for (int ne = 0; ne < nbEdges; ++ne)      // size_t doesn't work with CArray<double, 2> 
     545        { 
     546          for (int nv = 0; nv < nvertex; ++nv) 
     547          { 
     548            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne)); 
     549            for (int nh = 0; nh < 4; ++nh) 
     550            { 
     551              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh]; 
     552            } 
     553          } 
     554        } 
     555 
     556        // Recuperating the node global indexing and writing edge_nodes 
     557        // Creating map edgeHash2IdxGlo <hash, idxGlo> 
     558        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList); 
     559        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 
     560        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it; 
     561        vector <size_t> edgeNodes; 
     562        for (int ne = 0; ne < nbEdges; ++ne) 
     563        { 
     564          for (int nv = 0; nv < nvertex; ++nv) 
     565          { 
     566            int nh = 0; 
     567            it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh)); 
     568            // The loop below is needed in case if a hash generated by domain "edges" differs 
     569            // from that generated by domain "nodes" because of possible precision issues 
     570            while (it == nodeHash2IdxGlo.end()) 
     571            { 
     572              ++nh; 
     573              it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh)); 
     574            } 
     575            edge_nodes(nv,ne) = it->second[0]; 
     576            edgeNodes.push_back(it->second[0]); 
     577          } 
     578          edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne); 
     579        } 
     580      } // nodesAreWritten 
     581 
     582      // Case (2): node indexes have not been generated previously 
     583      else 
     584      { 
     585        // (2.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> > 
     586        vector<size_t> hashValues(4); 
     587        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
     588        CArray<size_t,1> nodeHashList(nbEdges*nvertex*4); 
     589        for (int ne = 0; ne < nbEdges; ++ne) 
     590        { 
     591          for (int nv = 0; nv < nvertex; ++nv) 
     592          { 
     593            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne)); 
     594            for (int nh = 0; nh < 4; ++nh) 
     595            { 
     596              if (nodeHash2Idx[hashValues[nh]].size() == 0) 
     597              { 
     598                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues, mpiRank)); 
     599                nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 
     600              } 
     601              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh]; 
     602            } 
     603          } 
     604        } 
     605 
     606        // (2.2) Generating global node indexes 
     607        // The ownership criterion: priority of the process holding the smaller index 
     608        // Maps generated in this step are: 
     609        // nodeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
     610        // nodeIdx2IdxMin = <idx, idxMin> 
     611        // nodeIdx2IdxGlo = <idxMin, idxMin> 
     612 
     613        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
     614        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 
     615        CArray<size_t,1> nodeIdxMinList(nbEdges*nvertex); 
     616 
     617        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
     618        dhtNodeHash.computeIndexInfoMapping(nodeHashList); 
     619        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
     620        size_t iIdxMin = 0; 
     621 
     622        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
     623        { 
     624          size_t idx = (it->second)[0]; 
     625          size_t idxMin = (it->second)[0]; 
     626          for (int i = 2; i < (it->second).size();) 
     627          { 
     628            if (mpiRank == (it->second)[i+1]) 
     629            { 
     630              idx = (it->second)[i]; 
     631            } 
     632            if ((it->second)[i] < idxMin) 
     633            { 
     634                idxMin = (it->second)[i]; 
     635                (it->second)[i] = (it->second)[i-2]; 
     636            } 
     637            i += 2; 
     638          } 
     639          (it->second)[0] = idxMin; 
     640          if (nodeIdx2IdxMin.count(idx) == 0) 
     641          { 
     642            nodeIdx2IdxMin[idx].push_back(idxMin); 
     643            if (idx == idxMin) 
     644              nodeIdx2IdxGlo[idxMin].push_back(idxMin); 
     645            nodeIdxMinList(iIdxMin) = idxMin; 
     646            ++iIdxMin; 
     647          } 
     648        } 
     649        nodeIdxMinList.resizeAndPreserve(iIdxMin); 
     650        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm); 
     651        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm); 
     652        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList); 
     653        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap(); 
     654        // nodeIdx2IdxGlo holds global indexes only for nodes owned by a process 
     655        // nodeIdxMin2IdxGlo holds global indexes for all nodes generated by a process 
     656 
     657        // (2.3) Saving variables: node_lon, node_lat, edge_nodes 
     658        // Creating map nodeHash2IdxGlo <hash, idxGlo> 
     659        // Creating map edgeHash2IdxGlo <hash, idxGlo> 
     660        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal(); 
     661        node_count = dhtNodeIdxGlo.getIndexCount(); 
     662        node_start = dhtNodeIdxGlo.getIndexStart(); 
     663        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo; 
     664        node_lon.resize(node_count); 
     665        node_lat.resize(node_count); 
     666        vector <size_t> edgeNodes; 
     667        size_t idxGlo = 0; 
     668 
     669        for (int ne = 0; ne < nbEdges; ++ne) 
     670        { 
     671          for (int nv = 0; nv < nvertex; ++nv) 
     672          { 
     673            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne)); 
     674            size_t myIdx = generateNodeIndex(hashValues, mpiRank); 
     675            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxMin.find(myIdx); 
     676            size_t ownerIdx = (itIdx->second)[0]; 
     677 
     678            if (myIdx == ownerIdx) 
     679            { 
     680              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdx2IdxGlo.find(myIdx); 
     681              idxGlo = (itIdxGlo->second)[0]; 
     682//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne)); 
     683              node_lon(idxGlo - node_start) = bounds_lon(nv, ne); 
     684              node_lat(idxGlo - node_start) = bounds_lat(nv, ne); 
     685            } 
     686            else 
     687            { 
     688              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdxMin2IdxGlo.find(ownerIdx); 
     689              idxGlo = (itIdxGlo->second)[0]; 
     690            } 
     691 
     692            edge_nodes(nv,ne) = idxGlo; 
     693            for (int nh = 0; nh < 4; ++nh) 
     694              nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo); 
     695            edgeNodes.push_back(idxGlo); 
     696          } 
     697          edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne); 
     698          edgeNodes.clear(); 
     699        } 
     700        pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm); 
     701      } // !nodesAreWritten 
     702 
     703      pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm); 
     704      edgesAreWritten = true; 
     705    } //nvertex = 2 
     706 
     707    else 
     708    { 
     709      nbFaces = bounds_lon.shape()[1]; 
     710      face_lon.resize(nbFaces); 
     711      face_lat.resize(nbFaces); 
     712      face_lon = lonvalue; 
     713      face_lat = latvalue; 
     714      face_nodes.resize(nvertex, nbFaces); 
     715      face_edges.resize(nvertex, nbFaces); 
     716 
     717      // Case (1): edges have been previously generated 
     718      if (edgesAreWritten) 
     719      { 
     720        // (1.1) Recuperating node global indexing and saving face_nodes 
     721        vector<size_t> hashValues(4); 
     722        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
     723        for (int nf = 0; nf < nbFaces; ++nf) 
     724        { 
     725          for (int nv = 0; nv < nvertex; ++nv) 
     726            { 
     727              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     728              for (int nh = 0; nh < 4; ++nh) 
     729                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh]; 
     730            } 
     731        } 
     732        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList); 
     733        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 
     734        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
     735        CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
     736        for (int nf = 0; nf < nbFaces; ++nf) 
     737        { 
     738          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     739          { 
     740            int nh1 = 0; 
     741            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     742            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     743            while (it1 == nodeHash2IdxGlo.end()) 
     744            { 
     745              ++nh1; 
     746              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     747            } 
     748            int nh2 = 0; 
     749            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     750            while (it2 == nodeHash2IdxGlo.end()) 
     751            { 
     752              ++nh2; 
     753              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     754            } 
     755            face_nodes(nv1,nf) = it1->second[0]; 
     756            edgeHashList(nf*nvertex + nv1) = hashPairOrdered(it1->second[0], it2->second[0]); 
     757 
     758          } 
     759        } 
     760 
     761        // (1.2) Recuperating edge global indexing and saving face_edges 
     762        pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList); 
     763        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap(); 
     764        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash; 
     765        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank; 
     766 
     767        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
     768        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex); 
     769        size_t iIdx = 0; 
     770 
     771        for (int nf = 0; nf < nbFaces; ++nf) 
     772        { 
     773          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     774          { 
     775            int nh1 = 0; 
     776            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     777            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     778            while (it1 == nodeHash2IdxGlo.end()) 
     779            { 
     780              ++nh1; 
     781              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     782            } 
     783            int nh2 = 0; 
     784            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     785            while (it2 == nodeHash2IdxGlo.end()) 
     786            { 
     787              ++nh2; 
     788              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     789            } 
     790            size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     791            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
     792            itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 
     793            size_t edgeIdxGlo = itEdgeHash->second[0]; 
     794            face_edges(nv1,nf) = edgeIdxGlo; 
     795            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
     796            { 
     797              edgeIdxGloList(iIdx) = edgeIdxGlo; 
     798              ++iIdx; 
     799            } 
     800            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
     801            edgeHash2Rank[edgeHash].push_back(itEdgeHash->first); 
     802            edgeHash2Rank[edgeHash].push_back(mpiRank); 
     803          } 
     804        } 
     805        edgeIdxGloList.resizeAndPreserve(iIdx); 
     806 
     807        // (1.3) Saving remaining variables edge_faces and face_faces 
     808 
     809        // Establishing edge ownership 
     810        // The ownership criterion: priority of the process with smaller rank 
     811        CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm); 
     812        dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 
     813        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 
     814        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo;            // map is needed purely for counting 
     815 
     816        // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>> 
     817        // edgeHashMin2IdxGlo edgeHash2Info = <edgeHashMin, idxGlo> 
     818        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
     819        { 
     820          vector <size_t> edgeInfo = it->second; 
     821          if (edgeInfo.size() == 4)                       // two processes generate the same edge 
     822            if (edgeInfo[1] > edgeInfo[3]) 
     823              edgeInfo[1] = edgeInfo[3]; 
     824          if (edgeInfo[1] == mpiRank) 
     825            edgeHashMin2IdxGlo[it->first].push_back(edgeInfo[0]); 
     826        } 
     827 
     828        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeHashMin2IdxGlo, comm); 
     829        edge_count = dhtEdgeIdxGlo.getIndexCount(); 
     830        edge_start = dhtEdgeIdxGlo.getIndexStart(); 
     831 
     832        edge_faces.resize(2, edge_count); 
     833        face_faces.resize(nvertex, nbFaces); 
     834 
     835        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     836        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList); 
     837        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
     838        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2; 
     839 
     840        for (int nf = 0; nf < nbFaces; ++nf) 
     841        { 
     842          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     843          { 
     844            int nh1 = 0; 
     845            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     846            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     847            while (it1 == nodeHash2IdxGlo.end()) 
     848            { 
     849              ++nh1; 
     850              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     851            } 
     852            int nh2 = 0; 
     853            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     854            while (it2 == nodeHash2IdxGlo.end()) 
     855            { 
     856              ++nh2; 
     857              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     858            } 
     859            size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     860            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
     861            itEdgeHash = edgeHash2Info.find(edgeHash); 
     862            size_t edgeOwnerRank = itEdgeHash->second[1]; 
     863            int edgeIdxGlo = itEdgeHash->second[0]; 
     864 
     865            if (mpiRank == edgeOwnerRank) 
     866            { 
     867              itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     868              int face1 = itFace1->second[0]; 
     869              if (itFace1->second.size() == 1) 
     870              { 
     871                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     872                edge_faces(1, edgeIdxGlo - edge_start) = -999; 
     873                face_faces(nv1, nf) = -1; 
     874              } 
     875              else 
     876              { 
     877                int face2 = itFace1->second[1]; 
     878                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     879                edge_faces(1, edgeIdxGlo - edge_start) = face2; 
     880                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     881              } 
     882            } 
     883            else 
     884            { 
     885              itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     886              int face1 = itFace1->second[0]; 
     887              int face2 = itFace1->second[1]; 
     888              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     889            } 
     890          } 
     891        } 
     892      } // edgesAreWritten 
     893 
     894      // Case (2): nodes have been previously generated 
     895      else if (nodesAreWritten) 
     896      { 
     897        // (2.1) Generating nodeHashList 
     898        vector<size_t> hashValues(4); 
     899        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
     900        for (int nf = 0; nf < nbFaces; ++nf) 
     901        { 
     902          for (int nv = 0; nv < nvertex; ++nv) 
     903            { 
     904              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     905              for (int nh = 0; nh < 4; ++nh) 
     906                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh]; 
     907            } 
     908        } 
     909 
     910        // (2.2) Recuperating node global indexing and saving face_nodes 
     911        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList 
     912        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList); 
     913        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 
     914        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 
     915        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
     916        CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
     917        for (int nf = 0; nf < nbFaces; ++nf) 
     918        { 
     919          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     920          { 
     921            int nh1 = 0; 
     922            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     923            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     924            while (it1 == nodeHash2IdxGlo.end()) 
     925            { 
     926              ++nh1; 
     927              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     928            } 
     929            int nh2 = 0; 
     930            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     931            while (it2 == nodeHash2IdxGlo.end()) 
     932            { 
     933              ++nh2; 
     934              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     935            } 
     936            face_nodes(nv1,nf) = it1->second[0]; 
     937            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
     938            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(it1->second[0], it2->second[0], mpiRank)); 
     939            edgeHash2Idx[edgeHash].push_back(mpiRank); 
     940            edgeHashList(nf*nvertex + nv1) = edgeHash; 
     941          } 
     942        } 
     943 
     944        // (2.2) Generating global edge indexes 
     945        // The ownership criterion: priority of the process holding the smaller index 
     946        // Maps generated in this step are: 
     947        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
     948        // edgeIdx2IdxMin = = <idx, idxMin> 
     949        // edgeIdx2IdxGlo = <idxMin, idxGlo> 
     950 
     951        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm); 
     952        dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 
     953        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 
     954 
     955        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 
     956        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 
     957        CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex); 
     958        size_t iIdxMin = 0; 
     959 
     960        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
     961        { 
     962          size_t idxMin = (it->second)[0]; 
     963          size_t idx = (it->second)[0]; 
     964          for (int i = 2; i < (it->second).size();) 
     965          { 
     966            if (mpiRank == (it->second)[i+1]) 
     967            { 
     968              idx = (it->second)[i]; 
     969            } 
     970            if ((it->second)[i] < idxMin) 
     971            { 
     972                idxMin = (it->second)[i]; 
     973                (it->second)[i] = (it->second)[i-2]; 
     974            } 
     975            i += 2; 
     976          } 
     977          (it->second)[0] = idxMin; 
     978          if (edgeIdx2IdxMin.count(idx) == 0) 
     979          { 
     980            edgeIdx2IdxMin[idx].push_back(idxMin); 
     981            if (idx == idxMin) 
     982              edgeIdx2IdxGlo[idxMin].push_back(idxMin); 
     983            edgeIdxMinList(iIdxMin) = idxMin; 
     984            ++iIdxMin; 
     985          } 
     986        } 
     987        edgeIdxMinList.resizeAndPreserve(iIdxMin); 
     988        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm); 
     989        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm); 
     990        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList); 
     991        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap(); 
     992        // edgeIdx2IdxGlo holds global indexes only for edges owned by a process 
     993        // edgeIdxMin2IdxGlo holds global indexes for all edges generated by a process 
     994 
     995        // (2.3) Saving variables: edge_lon, edge_lat, face_edges 
     996        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal(); 
     997        edge_count = dhtEdgeIdxGlo.getIndexCount(); 
     998        edge_start = dhtEdgeIdxGlo.getIndexStart(); 
     999        edge_lon.resize(edge_count); 
     1000        edge_lat.resize(edge_count); 
     1001        edge_nodes.resize(2, edge_count); 
     1002        face_edges.resize(nvertex, nbFaces); 
     1003 
     1004        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
     1005        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex); 
     1006        size_t iIdx = 0; 
     1007 
     1008        for (int nf = 0; nf < nbFaces; ++nf) 
     1009        { 
     1010          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1011          { 
     1012            // Getting global indexes of edge's nodes 
     1013            int nh1 = 0; 
     1014            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1015            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     1016            while (it1 == nodeHash2IdxGlo.end()) 
     1017            { 
     1018              ++nh1; 
     1019              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     1020            } 
     1021            int nh2 = 0; 
     1022            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     1023            while (it2 == nodeHash2IdxGlo.end()) 
     1024            { 
     1025              ++nh2; 
     1026              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     1027            } 
     1028            // Getting edge global index 
     1029            size_t nodeIdxGlo1 = it1->second[0]; 
     1030            size_t nodeIdxGlo2 = it2->second[0]; 
     1031            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
     1032            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
     1033            size_t ownerIdx = (itIdx->second)[0]; 
     1034            int edgeIdxGlo = 0; 
     1035            size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     1036 
     1037            if (myIdx == ownerIdx) 
     1038            { 
     1039              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1040              edgeIdxGlo = (itIdxGlo->second)[0]; 
     1041              edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ? 
     1042                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) : 
     1043                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.); 
     1044              edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 
     1045              edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 
     1046              edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 
     1047            } 
     1048            else 
     1049            { 
     1050              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1051              edgeIdxGlo = (itIdxGlo->second)[0]; 
     1052            } 
     1053            face_edges(nv1,nf) = edgeIdxGlo; 
     1054            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
     1055            { 
     1056              edgeIdxGloList(iIdx) = edgeIdxGlo; 
     1057              ++iIdx; 
     1058            } 
     1059            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
     1060          } 
     1061        } 
     1062        edgeIdxGloList.resizeAndPreserve(iIdx); 
     1063 
     1064        // (2.4) Saving remaining variables edge_faces and face_faces 
     1065        edge_faces.resize(2, edge_count); 
     1066        face_faces.resize(nvertex, nbFaces); 
     1067 
     1068        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     1069        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList); 
     1070        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
     1071        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 
     1072 
     1073        for (int nf = 0; nf < nbFaces; ++nf) 
     1074        { 
     1075          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1076          { 
     1077            // Getting global indexes of edge's nodes 
     1078            int nh1 = 0; 
     1079            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1080            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     1081            while (it1 == nodeHash2IdxGlo.end()) 
     1082            { 
     1083              ++nh1; 
     1084              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     1085            } 
     1086            int nh2 = 0; 
     1087            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     1088            while (it2 == nodeHash2IdxGlo.end()) 
     1089            { 
     1090              ++nh2; 
     1091              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     1092            } 
     1093            size_t nodeIdxGlo1 = it1->second[0]; 
     1094            size_t nodeIdxGlo2 = it2->second[0]; 
     1095 
     1096            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
     1097            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
     1098            size_t ownerIdx = (itIdx->second)[0]; 
     1099            size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     1100 
     1101            if (myIdx == ownerIdx) 
     1102            { 
     1103              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1104              int edgeIdxGlo = (itIdxGlo->second)[0]; 
     1105              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     1106              int face1 = it1->second[0]; 
     1107              if (it1->second.size() == 1) 
     1108              { 
     1109                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     1110                edge_faces(1, edgeIdxGlo - edge_start) = -999; 
     1111                face_faces(nv1, nf) = -1; 
     1112              } 
     1113              else 
     1114              { 
     1115                int face2 = it1->second[1]; 
     1116                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     1117                edge_faces(1, edgeIdxGlo - edge_start) = face2; 
     1118                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     1119              } 
     1120            } 
     1121            else 
     1122            { 
     1123              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1124              size_t edgeIdxGlo = (itIdxGlo->second)[0]; 
     1125              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     1126              int face1 = it1->second[0]; 
     1127              int face2 = it1->second[1]; 
     1128              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     1129            } 
     1130          } 
     1131        } 
     1132      } // nodesAreWritten 
     1133 
     1134      // Case (3): Neither nodes nor edges have been previously generated 
     1135      else 
     1136      { 
     1137        // (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> > 
     1138        vector<size_t> hashValues(4); 
     1139        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
     1140        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
     1141        size_t iHash = 0; 
     1142        for (int nf = 0; nf < nbFaces; ++nf) 
     1143        { 
     1144          for (int nv = 0; nv < nvertex; ++nv) 
     1145          { 
     1146            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     1147            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 
     1148            for (int nh = 0; nh < 4; ++nh) 
     1149            { 
     1150              if (nodeHash2Idx.count(hashValues[nh])==0) 
     1151              { 
     1152                nodeHash2Idx[hashValues[nh]].push_back(nodeIndex); 
     1153                nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 
     1154                nodeHashList(iHash) = hashValues[nh]; 
     1155                ++iHash; 
     1156              } 
     1157            } 
     1158          } 
     1159        } 
     1160        nodeHashList.resizeAndPreserve(iHash); 
     1161 
     1162        // (3.2) Generating global node indexes 
     1163        // The ownership criterion: priority of the process holding the smaller index 
     1164        // Maps generated in this step are: 
     1165        // nodeHash2Info = <hash, [[idxMin, rankMin], [idx1, rank1], [idx3, rank3]..]> 
     1166        // nodeIdx2IdxMin = = <idx, idxMin> 
     1167        // nodeIdx2IdxGlo = <idxMin, idxGlo> 
     1168 
     1169        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
     1170        dhtNodeHash.computeIndexInfoMapping(nodeHashList); 
     1171        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
     1172 
     1173        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
     1174        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 
     1175        CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4); 
     1176        size_t iIdxMin = 0; 
     1177 
     1178        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
     1179        { 
     1180          size_t idxMin = (it->second)[0]; 
     1181          size_t idx = (it->second)[0]; 
     1182          for (int i = 2; i < (it->second).size();) 
     1183          { 
     1184            if (mpiRank == (it->second)[i+1]) 
     1185            { 
     1186              idx = (it->second)[i]; 
     1187            } 
     1188            if ((it->second)[i] < idxMin) 
     1189            { 
     1190                idxMin = (it->second)[i]; 
     1191                (it->second)[i] = (it->second)[i-2]; 
     1192            } 
     1193            i += 2; 
     1194          } 
     1195          (it->second)[0] = idxMin; 
     1196          if (nodeIdx2IdxMin.count(idx) == 0) 
     1197          { 
     1198            nodeIdx2IdxMin[idx].push_back(idxMin); 
     1199            if (idx == idxMin) 
     1200              nodeIdx2IdxGlo[idxMin].push_back(idxMin); 
     1201            nodeIdxMinList(iIdxMin) = idxMin; 
     1202            ++iIdxMin; 
     1203          } 
     1204        } 
     1205 
     1206        nodeIdxMinList.resizeAndPreserve(iIdxMin); 
     1207        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm); 
     1208        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm); 
     1209        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList); 
     1210        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap(); 
     1211 
     1212        // (3.3) Saving node data: node_lon, node_lat, and face_nodes 
     1213        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList 
     1214        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal(); 
     1215        node_count = dhtNodeIdxGlo.getIndexCount(); 
     1216        node_start = dhtNodeIdxGlo.getIndexStart(); 
     1217        node_lon.resize(node_count); 
     1218        node_lat.resize(node_count); 
     1219        size_t nodeIdxGlo1 = 0; 
     1220        size_t nodeIdxGlo2 = 0; 
     1221        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 
     1222        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 
     1223        CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
     1224 
     1225        for (int nf = 0; nf < nbFaces; ++nf) 
     1226        { 
     1227          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1228          { 
     1229            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1230            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
     1231            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
     1232            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
     1233            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
     1234            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxMin.find(myNodeIdx1); 
     1235            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxMin.find(myNodeIdx2); 
     1236            size_t ownerNodeIdx = (itNodeIdx1->second)[0]; 
     1237 
     1238            if (myNodeIdx1 == ownerNodeIdx) 
     1239            { 
     1240              itNodeIdxGlo1 = nodeIdx2IdxGlo.find(myNodeIdx1); 
     1241              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
     1242              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf); 
     1243              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf); 
     1244            } 
     1245            else 
     1246            { 
     1247              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find(ownerNodeIdx); 
     1248              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
     1249            } 
     1250            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]); 
     1251            nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
     1252            size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1253            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank)); 
     1254            edgeHash2Idx[edgeHash].push_back(mpiRank); 
     1255            edgeHashList(nf*nvertex + nv1) = edgeHash; 
     1256            face_nodes(nv1,nf) = nodeIdxGlo1; 
     1257          } 
     1258        } 
     1259 
     1260        // (3.4) Generating global edge indexes 
     1261        // Maps generated in this step are: 
     1262        // edgeIdx2IdxMin = = <idx, idxMin> 
     1263        // edgeIdx2IdxGlo = <idxMin, idxGlo> 
     1264 
     1265        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm); 
     1266        dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 
     1267        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 
     1268        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
     1269 
     1270        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 
     1271        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 
     1272        CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex); 
     1273        iIdxMin = 0; 
     1274 
     1275        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
     1276        { 
     1277          size_t idxMin = (it->second)[0]; 
     1278          size_t idx = (it->second)[0]; 
     1279 
     1280          for (int i = 2; i < (it->second).size();) 
     1281          { 
     1282            if (mpiRank == (it->second)[i+1]) 
     1283            { 
     1284              idx = (it->second)[i]; 
     1285            } 
     1286            if ((it->second)[i] < idxMin) 
     1287            { 
     1288                idxMin = (it->second)[i]; 
     1289                (it->second)[i] = (it->second)[i-2]; 
     1290            } 
     1291            i += 2; 
     1292          } 
     1293          (it->second)[0] = idxMin; 
     1294          if (edgeIdx2IdxMin.count(idx) == 0) 
     1295          { 
     1296            edgeIdx2IdxMin[idx].push_back(idxMin); 
     1297            if (idx == idxMin) 
     1298              edgeIdx2IdxGlo[idxMin].push_back(idxMin); 
     1299            edgeIdxMinList(iIdxMin) = idxMin; 
     1300            ++iIdxMin; 
     1301          } 
     1302        } 
     1303        edgeIdxMinList.resizeAndPreserve(iIdxMin); 
     1304        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm); 
     1305        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm); 
     1306        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList); 
     1307        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap(); 
     1308 
     1309        // (3.5) Saving variables: edge_lon, edge_lat, face_edges 
     1310        // Creating map edgeIdxGlo2Face <idxGlo, face> 
     1311 
     1312        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal(); 
     1313        edge_count = dhtEdgeIdxGlo.getIndexCount(); 
     1314        edge_start = dhtEdgeIdxGlo.getIndexStart(); 
     1315        edge_lon.resize(edge_count); 
     1316        edge_lat.resize(edge_count); 
     1317        edge_nodes.resize(2, edge_count); 
     1318        face_edges.resize(nvertex, nbFaces); 
     1319 
     1320        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
     1321        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
     1322        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex); 
     1323        size_t iIdx = 0; 
     1324 
     1325        for (int nf = 0; nf < nbFaces; ++nf) 
     1326        { 
     1327          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1328          { 
     1329            // Getting global indexes of edge's nodes 
     1330            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1331            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
     1332            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
     1333 
     1334            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
     1335            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
     1336            it1 = nodeIdx2IdxMin.find(myNodeIdx1); 
     1337            it2 = nodeIdx2IdxMin.find(myNodeIdx2); 
     1338            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]); 
     1339            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]); 
     1340            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
     1341            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
     1342 
     1343            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
     1344            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
     1345            size_t ownerIdx = (itIdx->second)[0]; 
     1346 
     1347            int edgeIdxGlo = 0; 
     1348            size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     1349 
     1350            if (myIdx == ownerIdx) 
     1351            { 
     1352              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1353              edgeIdxGlo = (itIdxGlo->second)[0]; 
     1354              edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ? 
     1355                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) : 
     1356                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.); 
     1357              edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 
     1358              edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 
     1359              edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 
     1360            } 
     1361            else 
     1362            { 
     1363              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1364              edgeIdxGlo = (itIdxGlo->second)[0]; 
     1365            } 
     1366            face_edges(nv1,nf) = edgeIdxGlo; 
     1367            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
     1368            { 
     1369              edgeIdxGloList(iIdx) = edgeIdxGlo; 
     1370              ++iIdx; 
     1371            } 
     1372            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
     1373          } 
     1374        } 
     1375        edgeIdxGloList.resizeAndPreserve(iIdx); 
     1376 
     1377        // (3.6) Saving remaining variables edge_faces and face_faces 
     1378        edge_faces.resize(2, edge_count); 
     1379        face_faces.resize(nvertex, nbFaces); 
     1380 
     1381        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     1382        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList); 
     1383        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
     1384 
     1385        for (int nf = 0; nf < nbFaces; ++nf) 
     1386        { 
     1387          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1388          { 
     1389            // Getting global indexes of edge's nodes 
     1390            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1391            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
     1392            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
     1393 
     1394            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
     1395            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
     1396            it1 = nodeIdx2IdxMin.find(myNodeIdx1); 
     1397            it2 = nodeIdx2IdxMin.find(myNodeIdx2); 
     1398            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]); 
     1399            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]); 
     1400            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
     1401            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
     1402 
     1403            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
     1404            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
     1405            size_t ownerIdx = (itIdx->second)[0]; 
     1406            size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     1407 
     1408            if (myIdx == ownerIdx) 
     1409            { 
     1410              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1411              int edgeIdxGlo = (itIdxGlo->second)[0]; 
     1412              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     1413              int face1 = it1->second[0]; 
     1414              if (it1->second.size() == 1) 
     1415              { 
     1416                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     1417                edge_faces(1, edgeIdxGlo - edge_start) = -999; 
     1418                face_faces(nv1, nf) = -1; 
     1419              } 
     1420              else 
     1421              { 
     1422                size_t face2 = it1->second[1]; 
     1423                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     1424                edge_faces(1, edgeIdxGlo - edge_start) = face2; 
     1425                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     1426              } 
     1427            } 
     1428            else 
     1429            { 
     1430              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1431              size_t edgeIdxGlo = (itIdxGlo->second)[0]; 
     1432              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     1433              int face1 = it1->second[0]; 
     1434              int face2 = it1->second[1]; 
     1435              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     1436            } 
     1437          } 
    5201438        } 
    5211439      } 
    522       // Create nodes and edge_node connectivity 
    523  
    524       node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes 
    525       node_lat.resizeAndPreserve(nbEdges*nvertex); 
    526       edge_nodes.resizeAndPreserve(nvertex, nbEdges); 
    527  
    528 //      for (int ne = 0; ne < nbEdges; ++ne) 
    529 //      { 
    530 //        for (int nv = 0; nv < nvertex; ++nv) 
    531 //        { 
    532 //          if ( nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)) == -1) 
    533 //          { 
    534 //            edge_nodes(nv,ne) = nbNodes ; 
    535 //            node_lon(nbNodes) = bounds_lon(nv, ne); 
    536 //            node_lat(nbNodes) = bounds_lat(nv, ne); 
    537 //            ++nbNodes ; 
    538 //          } 
    539 //          else 
    540 //            edge_nodes(nv,ne) = nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)); 
    541 //        } 
    542 //      } 
    543 //      node_lon.resizeAndPreserve(nbNodes); 
    544 //      node_lat.resizeAndPreserve(nbNodes); 
    545  
    546       // Create edges 
    547       edge_lon.resizeAndPreserve(nbEdges); 
    548       edge_lat.resizeAndPreserve(nbEdges); 
    549  
    550       for (int ne = 0; ne < nbEdges; ++ne) 
    551       { 
    552         if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end()) 
    553         { 
    554           map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ; 
    555           edge_lon(ne) = lonvalue(ne); 
    556           edge_lat(ne) = latvalue(ne); 
    557         } 
    558       } 
    559       edgesAreWritten = true; 
    560     } // nvertex = 2 
    561     else 
    562     { 
    563       nbFaces = bounds_lon.shape()[1]; 
    564  
    565       // Create nodes and face_node connectivity 
    566       node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes 
    567       node_lat.resizeAndPreserve(nbFaces*nvertex); 
    568       face_nodes.resize(nvertex, nbFaces); 
    569  
    570       /*for (int nf = 0; nf < nbFaces; ++nf) 
    571       { 
    572         for (int nv = 0; nv < nvertex; ++nv) 
    573         { 
    574           if ( nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)) == -1) 
    575           { 
    576             face_nodes(nv,nf) = nbNodes ; 
    577             node_lon(nbNodes) = bounds_lon(nv, nf); 
    578             node_lat(nbNodes) = bounds_lat(nv ,nf); 
    579             ++nbNodes ; 
    580           } 
    581           else 
    582           { 
    583             face_nodes(nv,nf) = nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)); 
    584           } 
    585         } 
    586       } 
    587       node_lon.resizeAndPreserve(nbNodes); 
    588       node_lat.resizeAndPreserve(nbNodes);*/ 
    589  
    590  
    591       // Create edges and edge_nodes connectivity 
    592 //      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges 
    593 //      edge_lat.resizeAndPreserve(nbFaces*nvertex); 
    594 //      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex); 
    595 //      for (int nf = 0; nf < nbFaces; ++nf) 
    596 //      { 
    597 //        for (int nv1 = 0; nv1 < nvertex; ++nv1) 
    598 //        { 
    599 //          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
    600 //          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end()) 
    601 //          { 
    602 //            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ; 
    603 //            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf); 
    604 //            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 
    605 //                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) : 
    606 //                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);; 
    607 //            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 
    608 //            ++nbEdges; 
    609 //          } 
    610 //        } 
    611 //      } 
    612 //      edge_nodes.resizeAndPreserve(2, nbEdges); 
    613 //      edge_lon.resizeAndPreserve(nbEdges); 
    614 //      edge_lat.resizeAndPreserve(nbEdges); 
    615  
    616       // Create faces 
    617 //      face_lon.resize(nbFaces); 
    618 //      face_lat.resize(nbFaces); 
    619 //      face_lon = lonvalue; 
    620 //      face_lat = latvalue; 
    621 //      facesAreWritten = true; 
    622  
     1440      facesAreWritten = true; 
    6231441    } // nvertex >= 3 
    6241442 
    … …  
    6261444 
    6271445} // namespace xios 
    628  
    629 //  void CMesh::createMeshEpsilon(const MPI_Comm& comm, const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue, 
    630 //            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat) 
    631 //  { 
    632 // 
    633 //    nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
    634 // 
    635 //    if (nvertex == 1) 
    636 //    { 
    637 //      nbNodes = lonvalue.numElements(); 
    638 //      node_lon.resizeAndPreserve(nbNodes); 
    639 //      node_lat.resizeAndPreserve(nbNodes); 
    640 //      for (int nn = 0; nn < nbNodes; ++nn) 
    641 //      { 
    642 //        if ( nodeIndex(lonvalue(nn), latvalue(nn)) == -1 ) 
    643 //        { 
    644 //          node_lon(nn) = lonvalue(nn); 
    645 //          node_lat(nn) = latvalue(nn); 
    646 //        } 
    647 //      } 
    648 // 
    649 //      nodesAreWritten = true; 
    650 //    } 
    651 //    else if (nvertex == 2) 
    652 //    { 
    653 //      nbEdges = bounds_lon.shape()[1]; 
    654 // 
    655 //      // Create nodes and edge_node connectivity 
    656 //      node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes 
    657 //      node_lat.resizeAndPreserve(nbEdges*nvertex); 
    658 //      edge_nodes.resizeAndPreserve(nvertex, nbEdges); 
    659 // 
    660 //      for (int ne = 0; ne < nbEdges; ++ne) 
    661 //      { 
    662 //        for (int nv = 0; nv < nvertex; ++nv) 
    663 //        { 
    664 //          if ( nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)) == -1) 
    665 //          { 
    666 //            edge_nodes(nv,ne) = nbNodes ; 
    667 //            node_lon(nbNodes) = bounds_lon(nv, ne); 
    668 //            node_lat(nbNodes) = bounds_lat(nv, ne); 
    669 //            ++nbNodes ; 
    670 //          } 
    671 //          else 
    672 //            edge_nodes(nv,ne) = nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)); 
    673 //        } 
    674 //      } 
    675 //      node_lon.resizeAndPreserve(nbNodes); 
    676 //      node_lat.resizeAndPreserve(nbNodes); 
    677 // 
    678 //      // Create edges 
    679 //      edge_lon.resizeAndPreserve(nbEdges); 
    680 //      edge_lat.resizeAndPreserve(nbEdges); 
    681 // 
    682 //      for (int ne = 0; ne < nbEdges; ++ne) 
    683 //      { 
    684 //        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end()) 
    685 //        { 
    686 //          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ; 
    687 //          edge_lon(ne) = lonvalue(ne); 
    688 //          edge_lat(ne) = latvalue(ne); 
    689 //        } 
    690 //      } 
    691 //      edgesAreWritten = true; 
    692 //    } // nvertex = 2 
    693 //    else 
    694 //    { 
    695 //      nbFaces = bounds_lon.shape()[1]; 
    696 // 
    697 //      // Create nodes and face_node connectivity 
    698 //      node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes 
    699 //      node_lat.resizeAndPreserve(nbFaces*nvertex); 
    700 //      face_nodes.resize(nvertex, nbFaces); 
    701 // 
    702 //      for (int nf = 0; nf < nbFaces; ++nf) 
    703 //      { 
    704 //        for (int nv = 0; nv < nvertex; ++nv) 
    705 //        { 
    706 //          if ( nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)) == -1) 
    707 //          { 
    708 //            face_nodes(nv,nf) = nbNodes ; 
    709 //            node_lon(nbNodes) = bounds_lon(nv, nf); 
    710 //            node_lat(nbNodes) = bounds_lat(nv ,nf); 
    711 //            ++nbNodes ; 
    712 //          } 
    713 //          else 
    714 //          { 
    715 //            face_nodes(nv,nf) = nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)); 
    716 //          } 
    717 //        } 
    718 //      } 
    719 //      node_lon.resizeAndPreserve(nbNodes); 
    720 //      node_lat.resizeAndPreserve(nbNodes); 
    721 // 
    722 //      // Create edges and edge_nodes connectivity 
    723 //      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges 
    724 //      edge_lat.resizeAndPreserve(nbFaces*nvertex); 
    725 //      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex); 
    726 //      for (int nf = 0; nf < nbFaces; ++nf) 
    727 //      { 
    728 //        for (int nv1 = 0; nv1 < nvertex; ++nv1) 
    729 //        { 
    730 //          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
    731 //          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end()) 
    732 //          { 
    733 //            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ; 
    734 //            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf); 
    735 //            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 
    736 //                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) : 
    737 //                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);; 
    738 //            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 
    739 //            ++nbEdges; 
    740 //          } 
    741 //        } 
    742 //      } 
    743 //      edge_nodes.resizeAndPreserve(2, nbEdges); 
    744 //      edge_lon.resizeAndPreserve(nbEdges); 
    745 //      edge_lat.resizeAndPreserve(nbEdges); 
    746 // 
    747 //      // Create faces 
    748 //      face_lon.resize(nbFaces); 
    749 //      face_lat.resize(nbFaces); 
    750 //      face_lon = lonvalue; 
    751 //      face_lat = latvalue; 
    752 //      facesAreWritten = true; 
    753 //    } // nvertex >= 3 
    754 // 
    755 //  } // createMeshEpsilon 

  • XIOS/trunk/src/node/mesh.hpp

    r900 r924  
    99  
    1010#include "array_new.hpp" 
    11 #include "client_client_dht_template_impl.hpp" 
    1211#include "dht_auto_indexing.hpp" 
    1312 
    … …  
    3231      ~CMesh(void); 
    3332     
    34       int nbNodes; 
    35       int nbEdges; 
    36       int nbFaces; 
    37       int nvertex;   
     33      int nbNodesGlo; 
     34      int nbEdgesGlo; 
     35 
     36      int node_start; 
     37      int node_count; 
     38      int edge_start; 
     39      int edge_count; 
    3840 
    3941      bool nodesAreWritten; 
    … …  
    5658 
    5759      void createMesh(const CArray<double, 1>&, const CArray<double, 1>&,  
    58             const CArray<double, 2>&, const CArray<double, 2>& ); 
     60                      const CArray<double, 2>&, const CArray<double, 2>& ); 
    5961                         
    60       void createMeshEpsilon(const MPI_Comm&, const CArray<double, 1>&, const CArray<double, 1>&, 
    61             const CArray<double, 2>&, const CArray<double, 2>& ); 
     62      void createMeshEpsilon(const MPI_Comm&, 
     63                             const CArray<double, 1>&, const CArray<double, 1>&, 
     64                             const CArray<double, 2>&, const CArray<double, 2>& ); 
    6265             
    63       static CMesh* getMesh(StdString); 
     66      static CMesh* getMesh(StdString, int); 
    6467 
    6568    private: 
    6669 
     70      int nbNodes; 
     71      int nbEdges; 
     72      int nbFaces; 
     73 
    6774      static std::map <StdString, CMesh> meshList; 
    68       vector<size_t> createHashes (double, double); 
     75      static std::map <StdString, vector<int> > domainList; 
     76      CClientClientDHTSizet* pNodeGlobalIndex;                    // pointer to a map <nodeHash, nodeIdxGlo> 
     77      CClientClientDHTSizet* pEdgeGlobalIndex;                    // pointer to a map <edgeHash, edgeIdxGlo> 
     78 
     79      vector<size_t> createHashes (const double, const double); 
    6980 
    7081      size_t nodeIndex (double, double);                           // redundant in parallel version with epsilon precision 
Note: See TracChangeset for help on using the changeset viewer.