Changeset 1002 for XIOS/trunk


Ignore:
Timestamp:
11/23/16 16:22:48 (7 years ago)
Author:
oabramkina
Message:

UGRID: consistent node and edge numbering for varying number of procs.

Class CDHTAutoIndexing deprecated.

Location:
XIOS/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/arch/arch-X64_CURIE.fcm

    r591 r1002  
    99%BASE_CFLAGS    -diag-disable 1125 -diag-disable 279 
    1010%PROD_CFLAGS    -O3 -D BOOST_DISABLE_ASSERTS 
    11 %DEV_CFLAGS     -g -traceback 
     11#%DEV_CFLAGS     -g -traceback 
     12%DEV_CFLAGS     -g 
    1213%DEBUG_CFLAGS   -DBZ_DEBUG -g -traceback -fno-inline 
    1314 
    1415%BASE_FFLAGS    -D__NONE__  
    1516%PROD_FFLAGS    -O3 
    16 %DEV_FFLAGS     -g -O2 -traceback 
     17#%DEV_FFLAGS     -g -O2 -traceback 
     18%DEV_FFLAGS     -g -O2  
    1719%DEBUG_FFLAGS   -g -traceback 
    1820 

  • XIOS/trunk/src/dht_auto_indexing.cpp

    r924 r1002  
    7979    for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx) 
    8080    { 
    81       (itIdx->second)[0] = beginIndexOnProc_ + idx; 
     81//      (itIdx->second)[0] = beginIndexOnProc_ + idx; 
     82      (itIdx->second)[1] = beginIndexOnProc_ + idx; 
    8283      globalIndex_[idx] = beginIndexOnProc_ + idx; 
    8384      ++idx ; 

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

    r993 r1002  
    9191      seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2); 
    9292    } 
    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  
    11993    return seed ; 
    12094  } 
    … …  
    141115    return seed ; 
    142116  } 
     117 
     118  ///---------------------------------------------------------------- 
     119  /*! 
     120   * \fn size_t generateNodeIndex(vector<size_t>& valList) 
     121   * Generates a node index unique for all processes. 
     122   * \param [in] valList Vector storing four node hashes. 
     123   */ 
     124    size_t generateNodeIndex(vector<size_t>& valList) 
     125    { 
     126      // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere 
     127      vector<size_t> vec = valList; 
     128      sort (vec.begin(), vec.end()); 
     129      size_t seed = vec[0] ; 
     130      int it = 1; 
     131      for(; it != vec.size(); ++it) 
     132      { 
     133         seed = hashPair(seed, vec[it]); 
     134      } 
     135      return seed ; 
     136    } 
     137 
    143138 
    144139///---------------------------------------------------------------- 
    … …  
    600595        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
    601596        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4); 
     597        int nbHash = 0; 
    602598        for (int ne = 0; ne < nbEdges_; ++ne) 
    603599        { 
    … …  
    609605              if (nodeHash2Idx[hashValues[nh]].size() == 0) 
    610606              { 
    611                 nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues, mpiRank)); 
     607                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues)); 
    612608                nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 
     609                nodeHashList(nbHash) = hashValues[nh]; 
     610                ++nbHash; 
    613611              } 
    614               nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh]; 
    615             } 
    616           } 
    617         } 
     612            } 
     613          } 
     614        } 
     615        nodeHashList.resizeAndPreserve(nbHash); 
    618616 
    619617        // (2.2) Generating global node indexes 
    620         // The ownership criterion: priority of the process holding the smaller index 
     618        // The ownership criterion: priority of the process of smaller index 
    621619        // Maps generated in this step are: 
    622         // nodeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
    623         // nodeIdx2IdxMin = <idx, idxMin> 
    624         // nodeIdx2IdxGlo = <idxMin, idxMin> 
    625  
    626         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
    627         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 
    628         CArray<size_t,1> nodeIdxMinList(nbEdges_*nvertex); 
    629  
    630         CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
    631         dhtNodeHash.computeIndexInfoMapping(nodeHashList); 
    632         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
    633         size_t iIdxMin = 0; 
    634  
    635         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
    636         { 
    637           size_t idx = (it->second)[0]; 
    638           size_t idxMin = (it->second)[0]; 
    639           for (int i = 2; i < (it->second).size();) 
    640           { 
    641             if (mpiRank == (it->second)[i+1]) 
    642             { 
    643               idx = (it->second)[i]; 
    644             } 
    645             if ((it->second)[i] < idxMin) 
    646             { 
    647                 idxMin = (it->second)[i]; 
    648                 (it->second)[i] = (it->second)[i-2]; 
    649             } 
    650             i += 2; 
    651           } 
    652           (it->second)[0] = idxMin; 
    653           if (nodeIdx2IdxMin.count(idx) == 0) 
    654           { 
    655             nodeIdx2IdxMin[idx].push_back(idxMin); 
    656             if (idx == idxMin) 
    657               nodeIdx2IdxGlo[idxMin].push_back(idxMin); 
    658             nodeIdxMinList(iIdxMin) = idxMin; 
    659             ++iIdxMin; 
    660           } 
    661         } 
    662         nodeIdxMinList.resizeAndPreserve(iIdxMin); 
    663         CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm); 
    664         CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm); 
    665         dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList); 
    666         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap(); 
    667         // nodeIdx2IdxGlo holds global indexes only for nodes owned by a process 
    668         // nodeIdxMin2IdxGlo holds global indexes for all nodes generated by a process 
     620        // Maps generated in this step are: 
     621         // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]> 
     622         // nodeIdx2Idx = <idx, <rankOwner, idx>> 
     623 
     624         CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
     625         dhtNodeHash.computeIndexInfoMapping(nodeHashList); 
     626         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
     627 
     628 
     629         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx; 
     630         CArray<size_t,1> nodeIdxList(nbEdges_*nvertex*4); 
     631         size_t nIdx = 0; 
     632 
     633         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
     634         { 
     635           size_t rankMin = (it->second)[1]; 
     636           size_t idx = (it->second)[0]; 
     637           for (int i = 2; i < (it->second).size();) 
     638           { 
     639             if ( (it->second)[i+1] < rankMin) 
     640             { 
     641               idx = (it->second)[i]; 
     642               rankMin = (it->second)[i+1]; 
     643               (it->second)[i+1] = (it->second)[i-1]; 
     644             } 
     645             i += 2; 
     646           } 
     647           if (nodeIdx2Idx.count(idx) == 0) 
     648           { 
     649             if (mpiRank == rankMin) 
     650             { 
     651               nodeIdx2Idx[idx].push_back(rankMin); 
     652               nodeIdx2Idx[idx].push_back(idx); 
     653             } 
     654             nodeIdxList(nIdx) = idx; 
     655             ++nIdx; 
     656           } 
     657         } 
     658         nodeIdxList.resizeAndPreserve(nIdx); 
     659 
     660         // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. => 
     661         // Solution: global node indexing by hand. 
     662         // Maps modified in this step: 
     663         // nodeIdx2Idx = <idx, idxGlo> 
     664         int nodeCount = nodeIdx2Idx.size(); 
     665         int nodeStart, nbNodes; 
     666         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     667         int nNodes = nodeStart; 
     668         MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm); 
     669         nbNodesGlo = nNodes; 
     670 
     671         nodeStart -= nodeCount; 
     672         node_start = nodeStart; 
     673         node_count = nodeCount; 
     674         CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once 
     675         size_t count = 0; 
     676 
     677         for (int ne = 0; ne < nbEdges_; ++ne) 
     678         { 
     679           for (int nv = 0; nv < nvertex; ++nv) 
     680           { 
     681             vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne)); 
     682             size_t nodeIdx = generateNodeIndex(hashValues); 
     683             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx); 
     684             if (it != nodeIdx2Idx.end()) 
     685             { 
     686               if (dummyMap.count(nodeIdx) == 0) 
     687               { 
     688                 dummyMap[nodeIdx].push_back(nodeIdx); 
     689                 (it->second)[1] = node_start + count; 
     690                 ++count; 
     691               } 
     692             } 
     693           } 
     694         } 
     695 
     696         CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm); 
     697         dhtNodeIdx.computeIndexInfoMapping(nodeIdxList); 
     698         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap(); 
    669699 
    670700        // (2.3) Saving variables: node_lon, node_lat, edge_nodes 
    671701        // Creating map nodeHash2IdxGlo <hash, idxGlo> 
    672702        // Creating map edgeHash2IdxGlo <hash, idxGlo> 
    673         nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal(); 
    674         node_count = dhtNodeIdxGlo.getIndexCount(); 
    675         node_start = dhtNodeIdxGlo.getIndexStart(); 
     703//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal(); 
     704//        node_count = dhtNodeIdxGlo.getIndexCount(); 
     705//        node_start = dhtNodeIdxGlo.getIndexStart(); 
    676706        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo; 
    677707        node_lon.resize(node_count); 
    … …  
    685715          { 
    686716            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne)); 
    687             size_t myIdx = generateNodeIndex(hashValues, mpiRank); 
    688             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxMin.find(myIdx); 
    689             size_t ownerIdx = (itIdx->second)[0]; 
    690  
    691             if (myIdx == ownerIdx) 
    692             { 
    693               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdx2IdxGlo.find(myIdx); 
    694               idxGlo = (itIdxGlo->second)[0]; 
     717            size_t myIdx = generateNodeIndex(hashValues); 
     718            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxGlo.find(myIdx); 
     719            idxGlo = (itIdx->second)[1]; 
     720 
     721            if (mpiRank == (itIdx->second)[0]) 
     722            { 
    695723//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne)); 
    696724              node_lon(idxGlo - node_start) = bounds_lon(nv, ne); 
    697725              node_lat(idxGlo - node_start) = bounds_lat(nv, ne); 
    698726            } 
    699             else 
    700             { 
    701               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdxMin2IdxGlo.find(ownerIdx); 
    702               idxGlo = (itIdxGlo->second)[0]; 
    703             } 
    704  
    705727            edge_nodes(nv,ne) = idxGlo; 
    706728            for (int nh = 0; nh < 4; ++nh) 
    … …  
    793815        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank; 
    794816        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
    795         CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 
     817        CArray<size_t,1> edgeIdxList(nbFaces_*nvertex); 
    796818        size_t iIdx = 0; 
    797819 
    … …  
    824846              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
    825847              { 
    826                 edgeIdxGloList(iIdx) = edgeIdxGlo; 
     848                edgeIdxList(iIdx) = edgeIdxGlo; 
    827849                ++iIdx; 
    828850              } 
    829851              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
    830               edgeHash2Rank[edgeHash].push_back(itEdgeHash->first); 
    831852              edgeHash2Rank[edgeHash].push_back(mpiRank); 
     853              edgeHash2Rank[edgeHash].push_back(itEdgeHash->second[0]); 
    832854            } 
    833855            else 
    … …  
    837859          } 
    838860        } 
    839         edgeIdxGloList.resizeAndPreserve(iIdx); 
     861        edgeIdxList.resizeAndPreserve(iIdx); 
    840862 
    841863        // (1.3) Saving remaining variables edge_faces and face_faces 
    … …  
    846868        dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 
    847869        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 
    848         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo; 
    849  
    850         // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>> 
    851         // edgeHashMin2IdxGlo = <edgeHashMin, idxGlo> 
     870 
     871        // edgeHash2Info = <edgeHash, < rank1, idxGlo, rank2, idxGlo>> 
     872        int edgeCount = 0; 
    852873        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
    853874        { 
    854875          vector <size_t> edgeInfo = it->second; 
    855           if (edgeInfo.size() == 4)                       // two processes generate the same edge 
    856             if (edgeInfo[1] > edgeInfo[3]) 
    857               edgeInfo[1] = edgeInfo[3]; 
    858           if (edgeInfo[1] == mpiRank) 
    859             edgeHashMin2IdxGlo[it->first].push_back(edgeInfo[0]); 
    860         } 
    861  
    862         CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeHashMin2IdxGlo, comm); 
    863         edge_count = dhtEdgeIdxGlo.getIndexCount(); 
    864         edge_start = dhtEdgeIdxGlo.getIndexStart(); 
     876          if (edgeInfo[0] == mpiRank) 
     877          { 
     878            ++edgeCount; 
     879          } 
     880        } 
     881 
     882        int edgeStart, nbEdges; 
     883        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     884        int nEdges = edgeStart; 
     885        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm); 
     886        nbEdgesGlo = nEdges; 
     887 
     888        // edges to be splitted equally between procs 
     889        if ( (nbEdgesGlo % mpiSize) == 0) 
     890        { 
     891          edge_count = nbEdgesGlo/mpiSize; 
     892          edge_start = mpiRank*edge_count; 
     893        } 
     894        else 
     895        { 
     896          if (mpiRank == (mpiSize - 1) ) 
     897          { 
     898            edge_count = nbEdgesGlo/mpiSize; 
     899            edge_start = mpiRank*(nbEdgesGlo/mpiSize + 1); 
     900          } 
     901          else 
     902          { 
     903            edge_count = nbEdgesGlo/mpiSize + 1; 
     904            edge_start = mpiRank*edge_count; 
     905          } 
     906        } 
     907        CArray<size_t,1> edgeIdxGloList(edge_count); 
     908        for (int i = 0; i < edge_count; ++i) 
     909        { 
     910          edgeIdxGloList(i) = i + edge_start; 
     911        } 
     912 
     913        CClientClientDHTSizet dhtEdgeIdxGlo2Face (edgeIdxGlo2Face, comm); 
     914        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     915        dhtEdgeIdxGlo2Face.computeIndexInfoMapping(edgeIdxGloList); 
     916        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdgeIdxGlo2Face.getInfoIndexMap(); 
     917        dhtEdge2Face.computeIndexInfoMapping(edgeIdxList); 
     918        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
     919 
    865920 
    866921        edge_faces.resize(2, edge_count); 
     922        for (int i = 0; i < edge_count; ++i) 
     923        { 
     924          CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdxGlo2FaceIdx.find(i + edge_start); 
     925          int indexGlo = it->first; 
     926          vector<size_t> faces = it->second; 
     927          int face1 = faces[0]; 
     928          edge_faces(0, indexGlo - edge_start) = face1; 
     929          if (faces.size() == 2) 
     930          { 
     931            int face2 = faces[1]; 
     932            edge_faces(1, indexGlo - edge_start) = face2; 
     933          } 
     934          else 
     935          { 
     936            edge_faces(1, indexGlo - edge_start) = -999; 
     937          } 
     938        } 
     939 
     940        size_t tmp; 
     941        vector<size_t> tmpVec; 
     942        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2FaceIdx.begin(); it != edgeIdx2FaceIdx.end(); it++) 
     943        { 
     944          tmp = it->first; 
     945          tmpVec = it->second; 
     946          tmp++; 
     947        } 
     948 
     949        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2, itIndex; 
    867950        face_faces.resize(nvertex, nbFaces_); 
    868  
    869         CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
    870         dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList); 
    871         CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
    872         CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2; 
    873  
    874951        for (int nf = 0; nf < nbFaces_; ++nf) 
    875952        { 
    … …  
    896973              size_t faceIdxGlo = nbFacesAccum + nf; 
    897974              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
    898               itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 
    899               if (itEdgeHash != edgeHashMin2IdxGlo.end()) 
     975              itEdgeHash = edgeHash2Info.find(edgeHash); 
     976              int edgeIdxGlo = (itEdgeHash->second)[1]; 
     977 
     978              if ( (itEdgeHash->second)[0] == mpiRank) 
    900979              { 
    901                 int edgeIdxGlo = itEdgeHash->second[0]; 
    902                 itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     980                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo); 
    903981                int face1 = itFace1->second[0]; 
    904982                if (itFace1->second.size() == 1) 
    905983                { 
    906                   edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    907                   edge_faces(1, edgeIdxGlo - edge_start) = -999; 
    908984                  face_faces(nv1, nf) = 999999; 
    909985                } 
    … …  
    911987                { 
    912988                  int face2 = itFace1->second[1]; 
    913                   edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    914                   edge_faces(1, edgeIdxGlo - edge_start) = face2; 
    915989                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
    916990                } 
    … …  
    918992              else 
    919993              { 
    920                 itEdgeHash = edgeHashMin2IdxGlo.find(edgeHash); 
    921                 size_t edgeIdxGlo = itEdgeHash->second[0]; 
    922                 itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     994                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo); 
    923995                int face1 = itFace1->second[0]; 
    924996                int face2 = itFace1->second[1]; 
    … …  
    9571029        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
    9581030        CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 
     1031        int nEdgeHash = 0; 
    9591032        for (int nf = 0; nf < nbFaces_; ++nf) 
    9601033        { 
    … …  
    9781051            face_nodes(nv1,nf) = it1->second[0]; 
    9791052            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
    980             edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(it1->second[0], it2->second[0], mpiRank)); 
    981             edgeHash2Idx[edgeHash].push_back(mpiRank); 
    982             edgeHashList(nf*nvertex + nv1) = edgeHash; 
    983           } 
    984         } 
    985  
    986         // (2.2) Generating global edge indexes 
    987         // The ownership criterion: priority of the process holding the smaller index 
     1053            if (edgeHash2Idx.count(edgeHash) == 0) 
     1054            { 
     1055              edgeHash2Idx[edgeHash].push_back(edgeHash); 
     1056              edgeHash2Idx[edgeHash].push_back(mpiRank); 
     1057              edgeHashList(nEdgeHash) = edgeHash; 
     1058              ++nEdgeHash; 
     1059            } 
     1060          } 
     1061        } 
     1062        edgeHashList.resizeAndPreserve(nEdgeHash); 
     1063 
     1064        // (2.3) Generating global edge indexes 
     1065        // The ownership criterion: priority of the process with smaller rank 
    9881066        // Maps generated in this step are: 
    989         // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
    990         // edgeIdx2IdxMin = = <idx, idxMin> 
    991         // edgeIdx2IdxGlo = <idxMin, idxGlo> 
     1067        // edgeIdx2Idx = = <idx, <rankOwner, idx>> 
     1068        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>> 
    9921069 
    9931070        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm); 
    9941071        dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 
    9951072        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 
    996  
    997         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 
    998         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 
    999         CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex); 
    1000         size_t iIdxMin = 0; 
     1073        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
     1074 
     1075        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx; 
    10011076 
    10021077        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
    10031078        { 
    1004           size_t idxMin = (it->second)[0]; 
     1079          size_t rankMin = (it->second)[1]; 
    10051080          size_t idx = (it->second)[0]; 
     1081 
    10061082          for (int i = 2; i < (it->second).size();) 
    10071083          { 
    1008             if (mpiRank == (it->second)[i+1]) 
    1009             { 
     1084            if ((it->second)[i+1] < rankMin) 
     1085            { 
     1086              rankMin = (it->second)[i+1]; 
    10101087              idx = (it->second)[i]; 
    1011             } 
    1012             if ((it->second)[i] < idxMin) 
    1013             { 
    1014                 idxMin = (it->second)[i]; 
    1015                 (it->second)[i] = (it->second)[i-2]; 
     1088              (it->second)[i+1] = (it->second)[i-1]; 
    10161089            } 
    10171090            i += 2; 
    10181091          } 
    1019           (it->second)[0] = idxMin; 
    1020           if (edgeIdx2IdxMin.count(idx) == 0) 
    1021           { 
    1022             edgeIdx2IdxMin[idx].push_back(idxMin); 
    1023             if (idx == idxMin) 
    1024               edgeIdx2IdxGlo[idxMin].push_back(idxMin); 
    1025             edgeIdxMinList(iIdxMin) = idxMin; 
    1026             ++iIdxMin; 
    1027           } 
    1028         } 
    1029         edgeIdxMinList.resizeAndPreserve(iIdxMin); 
    1030         CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm); 
    1031         CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm); 
    1032         dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList); 
    1033         CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap(); 
    1034         // edgeIdx2IdxGlo holds global indexes only for edges owned by a process 
    1035         // edgeIdxMin2IdxGlo holds global indexes for all edges generated by a process 
    1036  
    1037         // (2.3) Saving variables: edge_lon, edge_lat, face_edges 
    1038         nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal(); 
    1039         edge_count = dhtEdgeIdxGlo.getIndexCount(); 
    1040         edge_start = dhtEdgeIdxGlo.getIndexStart(); 
    1041         edge_lon.resize(edge_count); 
    1042         edge_lat.resize(edge_count); 
    1043         edge_nodes.resize(2, edge_count); 
    1044         face_edges.resize(nvertex, nbFaces_); 
    1045  
    1046         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
    1047         CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 
    1048         size_t iIdx = 0; 
     1092          if (edgeIdx2Idx.count(idx) == 0) 
     1093          { 
     1094            if (mpiRank == rankMin) 
     1095            { 
     1096              edgeIdx2Idx[idx].push_back(rankMin); 
     1097              edgeIdx2Idx[idx].push_back(idx); 
     1098            } 
     1099          } 
     1100        } 
     1101 
     1102        int edgeCount = edgeIdx2Idx.size(); 
     1103        int edgeStart, nbEdges; 
     1104        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     1105        int nEdges = edgeStart; 
     1106        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm); 
     1107        nbEdgesGlo = nEdges; 
     1108 
     1109        edgeStart -= edgeCount; 
     1110        edge_start = edgeStart; 
     1111        edge_count = edgeCount; 
     1112        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap; 
     1113        int count = 0; 
    10491114 
    10501115        for (int nf = 0; nf < nbFaces_; ++nf) 
    … …  
    10681133              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
    10691134            } 
     1135            size_t nodeIdxGlo1 = it1->second[0]; 
     1136            size_t nodeIdxGlo2 = it2->second[0]; 
     1137 
     1138            if (nodeIdxGlo1 != nodeIdxGlo2) 
     1139            { 
     1140              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1141              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx); 
     1142              if (it != edgeIdx2Idx.end()) 
     1143              { 
     1144                if (dummyEdgeMap.count(edgeIdx) == 0) 
     1145                { 
     1146                  dummyEdgeMap[edgeIdx].push_back(edgeIdx); 
     1147                  (it->second)[1] = edge_start + count; 
     1148                  ++count; 
     1149                } 
     1150              } 
     1151            } 
     1152          } 
     1153        } 
     1154 
     1155        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm); 
     1156        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList); 
     1157        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap(); 
     1158 
     1159        // (2.4) Saving variables: edge_lon, edge_lat, face_edges 
     1160        edge_lon.resize(edge_count); 
     1161        edge_lat.resize(edge_count); 
     1162        edge_nodes.resize(2, edge_count); 
     1163        face_edges.resize(nvertex, nbFaces_); 
     1164 
     1165        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
     1166        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 
     1167        size_t iIdx = 0; 
     1168 
     1169        for (int nf = 0; nf < nbFaces_; ++nf) 
     1170        { 
     1171          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1172          { 
     1173            // Getting global indexes of edge's nodes 
     1174            int nh1 = 0; 
     1175            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1176            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     1177            while (it1 == nodeHash2IdxGlo.end()) 
     1178            { 
     1179              ++nh1; 
     1180              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1)); 
     1181            } 
     1182            int nh2 = 0; 
     1183            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2)); 
     1184            while (it2 == nodeHash2IdxGlo.end()) 
     1185            { 
     1186              ++nh2; 
     1187              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
     1188            } 
    10701189            // Getting edge global index 
    10711190            size_t nodeIdxGlo1 = it1->second[0]; 
    10721191            size_t nodeIdxGlo2 = it2->second[0]; 
    1073             size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
     1192            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
    10741193            if (nodeIdxGlo1 != nodeIdxGlo2) 
    10751194            { 
    1076               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
    1077               size_t ownerIdx = (itIdx->second)[0]; 
    1078               int edgeIdxGlo = 0; 
     1195              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx); 
     1196              int edgeIdxGlo = (itIdx->second)[1]; 
    10791197              size_t faceIdxGlo = nbFacesAccum + nf; 
    10801198 
    1081               if (myIdx == ownerIdx) 
     1199              if (mpiRank == (itIdx->second)[0]) 
    10821200              { 
    1083                 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
    1084                 edgeIdxGlo = (itIdxGlo->second)[0]; 
    10851201                double edgeLon; 
    10861202                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf)); 
    … …  
    10961212                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 
    10971213              } 
    1098               else 
    1099               { 
    1100                 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
    1101                 edgeIdxGlo = (itIdxGlo->second)[0]; 
    1102               } 
    11031214              face_edges(nv1,nf) = edgeIdxGlo; 
    11041215              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
    … …  
    11171228        edgeIdxGloList.resizeAndPreserve(iIdx); 
    11181229 
    1119         // (2.4) Saving remaining variables edge_faces and face_faces 
     1230        // (2.5) Saving remaining variables edge_faces and face_faces 
    11201231        edge_faces.resize(2, edge_count); 
    11211232        face_faces.resize(nvertex, nbFaces_); 
    … …  
    11251236        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
    11261237        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 
     1238        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx; 
    11271239 
    11281240        for (int nf = 0; nf < nbFaces_; ++nf) 
    … …  
    11491261            size_t nodeIdxGlo2 = it2->second[0]; 
    11501262 
    1151             size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
    1152             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
    1153             size_t ownerIdx = (itIdx->second)[0]; 
     1263            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1264            itIdx = edgeIdx2IdxGlo.find(myIdx); 
    11541265            size_t faceIdxGlo = nbFacesAccum + nf; 
    1155  
    1156             if (myIdx == ownerIdx) 
    1157             { 
    1158               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
    1159               int edgeIdxGlo = (itIdxGlo->second)[0]; 
     1266            int edgeIdxGlo = (itIdx->second)[1]; 
     1267 
     1268            if (mpiRank == (itIdx->second)[0]) 
     1269            { 
    11601270              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
    11611271              int face1 = it1->second[0]; 
    … …  
    11761286            else 
    11771287            { 
    1178               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
    1179               size_t edgeIdxGlo = (itIdxGlo->second)[0]; 
    11801288              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
    11811289              int face1 = it1->second[0]; 
    … …  
    12001308          { 
    12011309            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
    1202             size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 
     1310//            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 
     1311            size_t nodeIndex = generateNodeIndex(hashValues); 
    12031312            for (int nh = 0; nh < 4; ++nh) 
    12041313            { 
    … …  
    12161325 
    12171326        // (3.2) Generating global node indexes 
    1218         // The ownership criterion: priority of the process holding the smaller index 
     1327        // The ownership criterion: priority of the process with smaller rank. 
     1328        // With any other criterion it is not possible to have consistent node indexing for different number of procs. 
    12191329        // Maps generated in this step are: 
    1220         // nodeHash2Info = <hash, [[idxMin, rankMin], [idx1, rank1], [idx3, rank3]..]> 
    1221         // nodeIdx2IdxMin = = <idx, idxMin> 
    1222         // nodeIdx2IdxGlo = <idxMin, idxGlo> 
     1330        // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]> 
     1331        // nodeIdx2Idx = <idx, <rankOwner, idx>> 
    12231332 
    12241333        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
    … …  
    12261335        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
    12271336 
    1228         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
    1229         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 
    1230         CArray<size_t,1> nodeIdxMinList(nbFaces_*nvertex*4); 
    1231         size_t iIdxMin = 0; 
     1337        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx; 
     1338        CArray<size_t,1> nodeIdxList(nbFaces_*nvertex*4); 
     1339        size_t nIdx = 0; 
    12321340 
    12331341        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
    12341342        { 
    1235           size_t idxMin = (it->second)[0]; 
     1343          size_t rankMin = (it->second)[1]; 
    12361344          size_t idx = (it->second)[0]; 
    12371345          for (int i = 2; i < (it->second).size();) 
    12381346          { 
    1239             if (mpiRank == (it->second)[i+1]) 
     1347            if ( (it->second)[i+1] < rankMin) 
    12401348            { 
    12411349              idx = (it->second)[i]; 
    1242             } 
    1243             if ((it->second)[i] < idxMin) 
    1244             { 
    1245                 idxMin = (it->second)[i]; 
    1246                 (it->second)[i] = (it->second)[i-2]; 
     1350              rankMin = (it->second)[i+1]; 
     1351              (it->second)[i+1] = (it->second)[i-1]; 
    12471352            } 
    12481353            i += 2; 
    12491354          } 
    1250           (it->second)[0] = idxMin; 
    1251           if (nodeIdx2IdxMin.count(idx) == 0) 
    1252           { 
    1253             nodeIdx2IdxMin[idx].push_back(idxMin); 
    1254             if (idx == idxMin) 
    1255               nodeIdx2IdxGlo[idxMin].push_back(idxMin); 
    1256             nodeIdxMinList(iIdxMin) = idxMin; 
    1257             ++iIdxMin; 
    1258           } 
    1259         } 
    1260  
    1261         nodeIdxMinList.resizeAndPreserve(iIdxMin); 
    1262         CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm); 
    1263         CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm); 
    1264         dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList); 
    1265         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap(); 
     1355          if (nodeIdx2Idx.count(idx) == 0) 
     1356          { 
     1357            if (mpiRank == rankMin) 
     1358            { 
     1359              nodeIdx2Idx[idx].push_back(rankMin); 
     1360              nodeIdx2Idx[idx].push_back(idx); 
     1361            } 
     1362            nodeIdxList(nIdx) = idx; 
     1363            ++nIdx; 
     1364          } 
     1365        } 
     1366 
     1367//        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2Idx, comm); 
     1368        // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. => 
     1369        // Solution: global node indexing by hand. 
     1370        // Maps modified in this step: 
     1371        // nodeIdx2Idx = <idx, idxGlo> 
     1372        int nodeCount = nodeIdx2Idx.size(); 
     1373        int nodeStart, nbNodes; 
     1374        MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     1375        int nNodes = nodeStart; 
     1376        MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm); 
     1377        nbNodesGlo = nNodes; 
     1378 
     1379        nodeStart -= nodeCount; 
     1380        node_start = nodeStart; 
     1381        node_count = nodeCount; 
     1382        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once 
     1383        size_t count = 0; 
     1384 
     1385        for (int nf = 0; nf < nbFaces_; ++nf) 
     1386        { 
     1387          for (int nv = 0; nv < nvertex; ++nv) 
     1388          { 
     1389            vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     1390            size_t nodeIdx = generateNodeIndex(hashValues); 
     1391            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx); 
     1392            if (it != nodeIdx2Idx.end()) 
     1393            { 
     1394              if (dummyMap.count(nodeIdx) == 0) 
     1395              { 
     1396                dummyMap[nodeIdx].push_back(nodeIdx); 
     1397                (it->second)[1] = node_start + count; 
     1398                ++count; 
     1399              } 
     1400            } 
     1401          } 
     1402        } 
     1403        nodeIdxList.resizeAndPreserve(nIdx); 
     1404        CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm); 
     1405        dhtNodeIdx.computeIndexInfoMapping(nodeIdxList); 
     1406        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap(); 
    12661407 
    12671408        // (3.3) Saving node data: node_lon, node_lat, and face_nodes 
    12681409        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList 
    1269         nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal(); 
    1270         node_count = dhtNodeIdxGlo.getIndexCount(); 
    1271         node_start = dhtNodeIdxGlo.getIndexStart(); 
     1410//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal(); 
     1411//        node_count = dhtNodeIdxGlo.getIndexCount(); 
     1412//        node_start = dhtNodeIdxGlo.getIndexStart(); 
    12721413        node_lon.resize(node_count); 
    12731414        node_lat.resize(node_count); 
    … …  
    12861427            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
    12871428            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
    1288             size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
    1289             size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
    1290             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxMin.find(myNodeIdx1); 
    1291             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxMin.find(myNodeIdx2); 
    1292             size_t ownerNodeIdx = (itNodeIdx1->second)[0]; 
    1293  
    1294             if (myNodeIdx1 == ownerNodeIdx) 
    1295             { 
    1296               itNodeIdxGlo1 = nodeIdx2IdxGlo.find(myNodeIdx1); 
    1297               nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
     1429            size_t nodeIdx1 = generateNodeIndex(hashValues1); 
     1430            size_t nodeIdx2 = generateNodeIndex(hashValues2); 
     1431            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1); 
     1432            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2); 
     1433            size_t ownerRank = (itNodeIdx1->second)[0]; 
     1434            nodeIdxGlo1 = (itNodeIdx1->second)[1]; 
     1435            nodeIdxGlo2 = (itNodeIdx2->second)[1]; 
     1436 
     1437            if (mpiRank == ownerRank) 
     1438            { 
    12981439              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf); 
    12991440              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf); 
    13001441            } 
    1301             else 
    1302             { 
    1303               itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find(ownerNodeIdx); 
    1304               nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
    1305             } 
    1306             itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]); 
    1307             nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
    13081442            if (nodeIdxGlo1 != nodeIdxGlo2) 
    13091443            { 
    13101444              size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
    1311               edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank)); 
     1445              edgeHash2Idx[edgeHash].push_back(edgeHash); 
    13121446              edgeHash2Idx[edgeHash].push_back(mpiRank); 
    13131447              edgeHashList(nEdgeHash) = edgeHash; 
    … …  
    13211455        // (3.4) Generating global edge indexes 
    13221456        // Maps generated in this step are: 
    1323         // edgeIdx2IdxMin = = <idx, idxMin> 
    1324         // edgeIdx2IdxGlo = <idxMin, idxGlo> 
     1457        // edgeIdx2Idx = = <idx, <rankOwner, idx>> 
     1458        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>> 
    13251459 
    13261460        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm); 
    … …  
    13291463        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]> 
    13301464 
    1331         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 
    1332         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 
    1333         CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex); 
    1334         iIdxMin = 0; 
     1465        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx; 
    13351466 
    13361467        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
    13371468        { 
    1338           size_t idxMin = (it->second)[0]; 
     1469          size_t rankMin = (it->second)[1]; 
    13391470          size_t idx = (it->second)[0]; 
    13401471 
    13411472          for (int i = 2; i < (it->second).size();) 
    13421473          { 
    1343             if (mpiRank == (it->second)[i+1]) 
    1344             { 
     1474            if ((it->second)[i+1] < rankMin) 
     1475            { 
     1476              rankMin = (it->second)[i+1]; 
    13451477              idx = (it->second)[i]; 
    1346             } 
    1347             if ((it->second)[i] < idxMin) 
    1348             { 
    1349                 idxMin = (it->second)[i]; 
    1350                 (it->second)[i] = (it->second)[i-2]; 
     1478              (it->second)[i+1] = (it->second)[i-1]; 
    13511479            } 
    13521480            i += 2; 
    13531481          } 
    1354           (it->second)[0] = idxMin; 
    1355           if (edgeIdx2IdxMin.count(idx) == 0) 
    1356           { 
    1357             edgeIdx2IdxMin[idx].push_back(idxMin); 
    1358             if (idx == idxMin) 
    1359               edgeIdx2IdxGlo[idxMin].push_back(idxMin); 
    1360             edgeIdxMinList(iIdxMin) = idxMin; 
    1361             ++iIdxMin; 
    1362           } 
    1363         } 
    1364         edgeIdxMinList.resizeAndPreserve(iIdxMin); 
    1365         CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm); 
    1366         CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm); 
    1367         dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList); 
    1368         CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap(); 
     1482          if (edgeIdx2Idx.count(idx) == 0) 
     1483          { 
     1484            if (mpiRank == rankMin) 
     1485            { 
     1486              edgeIdx2Idx[idx].push_back(rankMin); 
     1487              edgeIdx2Idx[idx].push_back(idx); 
     1488            } 
     1489          } 
     1490        } 
     1491 
     1492        int edgeCount = edgeIdx2Idx.size(); 
     1493        int edgeStart, nbEdges; 
     1494        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     1495        int nEdges = edgeStart; 
     1496        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm); 
     1497        nbEdgesGlo = nEdges; 
     1498 
     1499        edgeStart -= edgeCount; 
     1500        edge_start = edgeStart; 
     1501        edge_count = edgeCount; 
     1502        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap; 
     1503        count = 0; 
     1504 
     1505        for (int nf = 0; nf < nbFaces_; ++nf) 
     1506        { 
     1507          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1508          { 
     1509            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1510            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
     1511            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
     1512            size_t nodeIdx1 = generateNodeIndex(hashValues1); 
     1513            size_t nodeIdx2 = generateNodeIndex(hashValues2); 
     1514            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1); 
     1515            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2); 
     1516            nodeIdxGlo1 = (itNodeIdx1->second)[1]; 
     1517            nodeIdxGlo2 = (itNodeIdx2->second)[1]; 
     1518 
     1519            if (nodeIdxGlo1 != nodeIdxGlo2) 
     1520            { 
     1521              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1522              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx); 
     1523              if (it != edgeIdx2Idx.end()) 
     1524              { 
     1525                if (dummyEdgeMap.count(edgeIdx) == 0) 
     1526                { 
     1527                  dummyEdgeMap[edgeIdx].push_back(edgeIdx); 
     1528                  (it->second)[1] = edge_start + count; 
     1529                  ++count; 
     1530                } 
     1531              } 
     1532            } 
     1533          } 
     1534        } 
     1535        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm); 
     1536        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList); 
     1537        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap(); 
    13691538 
    13701539        // (3.5) Saving variables: edge_lon, edge_lat, face_edges 
    13711540        // Creating map edgeIdxGlo2Face <idxGlo, face> 
    1372  
    1373         nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal(); 
    1374         edge_count = dhtEdgeIdxGlo.getIndexCount(); 
    1375         edge_start = dhtEdgeIdxGlo.getIndexStart(); 
     1541//        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal(); 
     1542//        edge_count = dhtEdgeIdxGlo.getIndexCount(); 
     1543//        edge_start = dhtEdgeIdxGlo.getIndexStart(); 
     1544 
    13761545        edge_lon.resize(edge_count); 
    13771546        edge_lat.resize(edge_count); 
    … …  
    13931562            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
    13941563 
    1395             size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
    1396             size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
    1397             it1 = nodeIdx2IdxMin.find(myNodeIdx1); 
    1398             it2 = nodeIdx2IdxMin.find(myNodeIdx2); 
    1399             itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]); 
    1400             itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]); 
    1401             size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
    1402             size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
     1564            size_t nodeIdx1 = generateNodeIndex(hashValues1); 
     1565            size_t nodeIdx2 = generateNodeIndex(hashValues2); 
     1566            it1 = nodeIdx2IdxGlo.find(nodeIdx1); 
     1567            it2 = nodeIdx2IdxGlo.find(nodeIdx2); 
     1568            size_t nodeIdxGlo1 = (it1->second)[1]; 
     1569            size_t nodeIdxGlo2 = (it2->second)[1]; 
    14031570 
    14041571            if (nodeIdxGlo1 != nodeIdxGlo2) 
    14051572            { 
    1406               size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
    1407               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
    1408               size_t ownerIdx = (itIdx->second)[0]; 
    1409               int edgeIdxGlo; 
     1573              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1574              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx); 
     1575              int edgeIdxGlo = (itIdx->second)[1]; 
    14101576              size_t faceIdxGlo = nbFacesAccum + nf; 
    14111577 
    1412               if (myIdx == ownerIdx) 
     1578              if (mpiRank == (itIdx->second)[0]) 
    14131579              { 
    1414                 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
    1415                 edgeIdxGlo = (itIdxGlo->second)[0]; 
    14161580                double edgeLon; 
    14171581                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf)); 
    … …  
    14271591                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 
    14281592              } 
    1429               else 
    1430               { 
    1431                 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
    1432                 edgeIdxGlo = (itIdxGlo->second)[0]; 
    1433               } 
    14341593              face_edges(nv1,nf) = edgeIdxGlo; 
    14351594              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
    … …  
    14651624            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
    14661625 
    1467             size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
    1468             size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
     1626            size_t myNodeIdx1 = generateNodeIndex(hashValues1); 
     1627            size_t myNodeIdx2 = generateNodeIndex(hashValues2); 
    14691628            if (myNodeIdx1 != myNodeIdx2) 
    14701629            { 
    1471               it1 = nodeIdx2IdxMin.find(myNodeIdx1); 
    1472               it2 = nodeIdx2IdxMin.find(myNodeIdx2); 
    1473               itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]); 
    1474               itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]); 
    1475               size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0]; 
    1476               size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
    1477  
    1478               size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank); 
    1479               CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
    1480               size_t ownerIdx = (itIdx->second)[0]; 
     1630              it1 = nodeIdx2IdxGlo.find(myNodeIdx1); 
     1631              it2 = nodeIdx2IdxGlo.find(myNodeIdx2); 
     1632              size_t nodeIdxGlo1 = (it1->second)[1]; 
     1633              size_t nodeIdxGlo2 = (it2->second)[1]; 
     1634              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1635              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx); 
     1636              int edgeIdxGlo = (itIdx->second)[1]; 
     1637 
    14811638              size_t faceIdxGlo = nbFacesAccum + nf; 
    14821639 
    1483               if (myIdx == ownerIdx) 
     1640              if (mpiRank == (itIdx->second)[0]) 
    14841641              { 
    1485                 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
    1486                 int edgeIdxGlo = (itIdxGlo->second)[0]; 
    14871642                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
    14881643                int face1 = it1->second[0]; 
    … …  
    15001655                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
    15011656                } 
    1502               } // myIdx == ownerIdx 
     1657              } 
    15031658              else 
    15041659              { 
    1505                 CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
    1506                 size_t edgeIdxGlo = (itIdxGlo->second)[0]; 
    15071660                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
    15081661                int face1 = it1->second[0]; 
    15091662                int face2 = it1->second[1]; 
    15101663                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
    1511               } // myIdx != ownerIdx 
     1664              } 
    15121665            } // myNodeIdx1 != myNodeIdx2 
    15131666            else 
    … …  
    15151668          } 
    15161669        } 
     1670 
    15171671      } 
    15181672      facesAreWritten = true; 

  • XIOS/trunk/src/test/test_regular.f90

    r993 r1002  
    2020  INTEGER :: sizeComm, rank    ! SIZE is a fortran function 
    2121 
    22   INTEGER :: nlon = 5 !100  
    23   INTEGER :: nlat = 5 !100 
     22  INTEGER :: nlon = 100  
     23  INTEGER :: nlat = 100 
    2424  INTEGER :: ncell  
    2525  INTEGER :: ilat, ilon, ind 
Note: See TracChangeset for help on using the changeset viewer.