Changeset 929


Ignore:
Timestamp:
09/09/16 16:57:43 (8 years ago)
Author:
oabramkina
Message:

Mesh connectivity:

Fixing a bug for UGRID in case meshes containing cells with different number of vertexes.

Function for determining face neighbors of a local domain added.

Location:
XIOS/trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/src/io/nc4_data_output.cpp

    r924 r929  
    617617              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges); 
    618618              SuperClassWriter::addAttribute("start_index", 0, &face_edges); 
     619              SuperClassWriter::addAttribute("_FillValue", 999999, &face_edges); 
    619620              dim0.clear(); 
    620621              dim0.push_back(dimEdge); 
     
    633634              SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces); 
    634635              SuperClassWriter::addAttribute("start_index", 0, &face_faces); 
     636              SuperClassWriter::addAttribute("_FillValue", 999999, &face_faces); 
    635637              SuperClassWriter::addAttribute("flag_values", -1, &face_faces); 
    636638              SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces); 
  • XIOS/trunk/src/node/mesh.cpp

    r924 r929  
    1414            ,  node_start{0}, node_count{0} 
    1515            ,  edge_start{0}, edge_count{0} 
    16             ,  nbFaces{0}, nbNodes{0}, nbEdges{0} 
     16            ,  nbFaces_{0}, nbNodes_{0}, nbEdges_{0} 
    1717            ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false} 
    1818            ,  node_lon(), node_lat() 
     
    310310    if (nvertex == 1) 
    311311    { 
    312       nbNodes = lonvalue.numElements(); 
    313       node_lon.resizeAndPreserve(nbNodes); 
    314       node_lat.resizeAndPreserve(nbNodes); 
    315       for (int nn = 0; nn < nbNodes; ++nn) 
     312      nbNodes_ = lonvalue.numElements(); 
     313      node_lon.resizeAndPreserve(nbNodes_); 
     314      node_lat.resizeAndPreserve(nbNodes_); 
     315      for (int nn = 0; nn < nbNodes_; ++nn) 
    316316      { 
    317317        if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end()) 
     
    325325    else if (nvertex == 2) 
    326326    { 
    327       nbEdges = bounds_lon.shape()[1]; 
     327      nbEdges_ = bounds_lon.shape()[1]; 
    328328 
    329329      // Create nodes and edge_node connectivity 
    330       node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes 
    331       node_lat.resizeAndPreserve(nbEdges*nvertex); 
    332       edge_nodes.resizeAndPreserve(nvertex, nbEdges); 
    333  
    334       for (int ne = 0; ne < nbEdges; ++ne) 
     330      node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes 
     331      node_lat.resizeAndPreserve(nbEdges_*nvertex); 
     332      edge_nodes.resizeAndPreserve(nvertex, nbEdges_); 
     333 
     334      for (int ne = 0; ne < nbEdges_; ++ne) 
    335335      { 
    336336        for (int nv = 0; nv < nvertex; ++nv) 
     
    338338          if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end()) 
    339339          { 
    340             map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes ; 
    341             edge_nodes(nv,ne) = nbNodes ; 
    342             node_lon(nbNodes) = bounds_lon(nv, ne); 
    343             node_lat(nbNodes) = bounds_lat(nv, ne); 
    344             ++nbNodes ; 
     340            map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ; 
     341            edge_nodes(nv,ne) = nbNodes_ ; 
     342            node_lon(nbNodes_) = bounds_lon(nv, ne); 
     343            node_lat(nbNodes_) = bounds_lat(nv, ne); 
     344            ++nbNodes_ ; 
    345345          } 
    346346          else 
     
    348348        } 
    349349      } 
    350       node_lon.resizeAndPreserve(nbNodes); 
    351       node_lat.resizeAndPreserve(nbNodes); 
     350      node_lon.resizeAndPreserve(nbNodes_); 
     351      node_lat.resizeAndPreserve(nbNodes_); 
    352352 
    353353      // Create edges 
    354       edge_lon.resizeAndPreserve(nbEdges); 
    355       edge_lat.resizeAndPreserve(nbEdges); 
    356  
    357       for (int ne = 0; ne < nbEdges; ++ne) 
     354      edge_lon.resizeAndPreserve(nbEdges_); 
     355      edge_lat.resizeAndPreserve(nbEdges_); 
     356 
     357      for (int ne = 0; ne < nbEdges_; ++ne) 
    358358      { 
    359359        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end()) 
     
    369369    else 
    370370    { 
    371       nbFaces = bounds_lon.shape()[1]; 
     371      nbFaces_ = bounds_lon.shape()[1]; 
    372372   
    373373      // Create nodes and face_node connectivity 
    374       node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes 
    375       node_lat.resizeAndPreserve(nbFaces*nvertex); 
    376       face_nodes.resize(nvertex, nbFaces); 
     374      node_lon.resizeAndPreserve(nbFaces_*nvertex);  // Max possible number of nodes 
     375      node_lat.resizeAndPreserve(nbFaces_*nvertex); 
     376      face_nodes.resize(nvertex, nbFaces_); 
    377377   
    378       for (int nf = 0; nf < nbFaces; ++nf) 
     378      for (int nf = 0; nf < nbFaces_; ++nf) 
    379379      { 
    380380        for (int nv = 0; nv < nvertex; ++nv) 
     
    382382          if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end()) 
    383383          { 
    384             map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes ; 
    385             face_nodes(nv,nf) = nbNodes ; 
    386             node_lon(nbNodes) = bounds_lon(nv, nf); 
    387             node_lat(nbNodes) = bounds_lat(nv ,nf); 
    388             ++nbNodes ; 
     384            map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ; 
     385            face_nodes(nv,nf) = nbNodes_ ; 
     386            node_lon(nbNodes_) = bounds_lon(nv, nf); 
     387            node_lat(nbNodes_) = bounds_lat(nv ,nf); 
     388            ++nbNodes_ ; 
    389389          } 
    390390          else 
     
    394394        } 
    395395      } 
    396       node_lon.resizeAndPreserve(nbNodes); 
    397       node_lat.resizeAndPreserve(nbNodes); 
     396      node_lon.resizeAndPreserve(nbNodes_); 
     397      node_lat.resizeAndPreserve(nbNodes_); 
    398398   
    399399      // Create edges and edge_nodes connectivity 
    400       edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges 
    401       edge_lat.resizeAndPreserve(nbFaces*nvertex); 
    402       edge_nodes.resizeAndPreserve(2, nbFaces*nvertex); 
    403       edge_faces.resize(2, nbFaces*nvertex); 
    404       face_edges.resize(nvertex, nbFaces); 
    405       face_faces.resize(nvertex, nbFaces); 
    406  
    407       vector<int> countEdges(nbFaces*nvertex);   // needed in case if edges have been already generated 
    408       vector<int> countFaces(nbFaces); 
    409       countEdges.assign(nbFaces*nvertex, 0); 
    410       countFaces.assign(nbFaces, 0); 
     400      edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges 
     401      edge_lat.resizeAndPreserve(nbFaces_*nvertex); 
     402      edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex); 
     403      edge_faces.resize(2, nbFaces_*nvertex); 
     404      face_edges.resize(nvertex, nbFaces_); 
     405      face_faces.resize(nvertex, nbFaces_); 
     406 
     407      vector<int> countEdges(nbFaces_*nvertex);   // needed in case if edges have been already generated 
     408      vector<int> countFaces(nbFaces_); 
     409      countEdges.assign(nbFaces_*nvertex, 0); 
     410      countFaces.assign(nbFaces_, 0); 
    411411      int edge; 
    412       for (int nf = 0; nf < nbFaces; ++nf) 
     412      for (int nf = 0; nf < nbFaces_; ++nf) 
    413413      { 
    414414        for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    418418          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end()) 
    419419          { 
    420             map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ; 
     420            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ; 
    421421            face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))]; 
    422             edge_faces(0,nbEdges) = nf; 
    423             edge_faces(1,nbEdges) = -999; 
    424             face_faces(nv1,nf) = -1; 
    425             edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf); 
    426             edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 
     422            edge_faces(0,nbEdges_) = nf; 
     423            edge_faces(1,nbEdges_) = -999; 
     424            face_faces(nv1,nf) = 999999; 
     425            edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf); 
     426            edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 
    427427                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) : 
    428428                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.); 
    429             edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 
    430             ++nbEdges; 
     429            edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 
     430            ++nbEdges_; 
    431431          } 
    432432          else 
     
    439439              if (countEdges[edge]==0) 
    440440              { 
    441                 face_faces(nv1,nf) = -1; 
     441                face_faces(nv1,nf) = 999999; 
    442442              } 
    443443              else 
     
    465465        } 
    466466      } 
    467       edge_nodes.resizeAndPreserve(2, nbEdges); 
    468       edge_faces.resizeAndPreserve(2, nbEdges); 
    469       edge_lon.resizeAndPreserve(nbEdges); 
    470       edge_lat.resizeAndPreserve(nbEdges); 
     467      edge_nodes.resizeAndPreserve(2, nbEdges_); 
     468      edge_faces.resizeAndPreserve(2, nbEdges_); 
     469      edge_lon.resizeAndPreserve(nbEdges_); 
     470      edge_lat.resizeAndPreserve(nbEdges_); 
    471471 
    472472      // Create faces 
    473       face_lon.resize(nbFaces); 
    474       face_lat.resize(nbFaces); 
     473      face_lon.resize(nbFaces_); 
     474      face_lat.resize(nbFaces_); 
    475475      face_lon = lonvalue; 
    476476      face_lat = latvalue; 
     
    502502    MPI_Comm_rank(comm, &mpiRank); 
    503503    MPI_Comm_size(comm, &mpiSize); 
     504    double prec = 1e-11;  // used in calculations of edge_lon/lat 
    504505 
    505506    if (nvertex == 1) 
    506507    { 
    507       nbNodes = lonvalue.numElements(); 
    508       node_lon.resize(nbNodes); 
    509       node_lat.resize(nbNodes); 
     508      nbNodes_ = lonvalue.numElements(); 
     509      node_lon.resize(nbNodes_); 
     510      node_lat.resize(nbNodes_); 
    510511      node_lon = lonvalue; 
    511512      node_lat = latvalue; 
     
    514515      vector<size_t> hashValues(4); 
    515516      CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo; 
    516       for (size_t nn = 0; nn < nbNodes; ++nn) 
     517      for (size_t nn = 0; nn < nbNodes_; ++nn) 
    517518      { 
    518519        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn)); 
    519520        for (size_t nh = 0; nh < 4; ++nh) 
    520521        { 
    521           nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes + nn); 
     522          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes_ + nn); 
    522523        } 
    523524      } 
     
    528529    else if (nvertex == 2) 
    529530    { 
    530       nbEdges = bounds_lon.shape()[1]; 
    531       edge_lon.resize(nbEdges); 
    532       edge_lat.resize(nbEdges); 
     531      nbEdges_ = bounds_lon.shape()[1]; 
     532      edge_lon.resize(nbEdges_); 
     533      edge_lat.resize(nbEdges_); 
    533534      edge_lon = lonvalue; 
    534535      edge_lat = latvalue; 
    535       edge_nodes.resize(nvertex, nbEdges); 
     536      edge_nodes.resize(nvertex, nbEdges_); 
     537 
     538      // For determining the global edge index 
     539      size_t nbEdgesOnProc = nbEdges_; 
     540      size_t nbEdgesAccum; 
     541      MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     542      nbEdgesAccum -= nbEdges_; 
     543 
    536544      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo; 
    537545      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 
     
    541549      { 
    542550        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> 
     551        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4); 
     552        for (int ne = 0; ne < nbEdges_; ++ne)      // size_t doesn't work with CArray<double, 2> 
    545553        { 
    546554          for (int nv = 0; nv < nvertex; ++nv) 
     
    559567        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 
    560568        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it; 
    561         vector <size_t> edgeNodes; 
    562         for (int ne = 0; ne < nbEdges; ++ne) 
     569        size_t nodeIdxGlo1, nodeIdxGlo2; 
     570        for (int ne = 0; ne < nbEdges_; ++ne) 
    563571        { 
    564572          for (int nv = 0; nv < nvertex; ++nv) 
     
    574582            } 
    575583            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); 
     584            if (nv ==0) 
     585              nodeIdxGlo1 = it->second[0]; 
     586            else 
     587              nodeIdxGlo2 = it->second[0]; 
     588          } 
     589          size_t edgeIdxGlo = nbEdgesAccum + ne; 
     590          edgeHash2IdxGlo[ hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2) ].push_back(edgeIdxGlo); 
    579591        } 
    580592      } // nodesAreWritten 
     593 
    581594 
    582595      // Case (2): node indexes have not been generated previously 
     
    586599        vector<size_t> hashValues(4); 
    587600        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
    588         CArray<size_t,1> nodeHashList(nbEdges*nvertex*4); 
    589         for (int ne = 0; ne < nbEdges; ++ne) 
     601        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4); 
     602        for (int ne = 0; ne < nbEdges_; ++ne) 
    590603        { 
    591604          for (int nv = 0; nv < nvertex; ++nv) 
     
    613626        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
    614627        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 
    615         CArray<size_t,1> nodeIdxMinList(nbEdges*nvertex); 
     628        CArray<size_t,1> nodeIdxMinList(nbEdges_*nvertex); 
    616629 
    617630        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
     
    667680        size_t idxGlo = 0; 
    668681 
    669         for (int ne = 0; ne < nbEdges; ++ne) 
     682        for (int ne = 0; ne < nbEdges_; ++ne) 
    670683        { 
    671684          for (int nv = 0; nv < nvertex; ++nv) 
     
    695708            edgeNodes.push_back(idxGlo); 
    696709          } 
    697           edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne); 
     710          if (edgeNodes[0] != edgeNodes[1]) 
     711          { 
     712            size_t edgeIdxGlo = nbEdgesAccum + ne; 
     713            edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(edgeIdxGlo); 
     714          } 
    698715          edgeNodes.clear(); 
    699716        } 
     
    707724    else 
    708725    { 
    709       nbFaces = bounds_lon.shape()[1]; 
    710       face_lon.resize(nbFaces); 
    711       face_lat.resize(nbFaces); 
     726      nbFaces_ = bounds_lon.shape()[1]; 
     727      face_lon.resize(nbFaces_); 
     728      face_lat.resize(nbFaces_); 
    712729      face_lon = lonvalue; 
    713730      face_lat = latvalue; 
    714       face_nodes.resize(nvertex, nbFaces); 
    715       face_edges.resize(nvertex, nbFaces); 
     731      face_nodes.resize(nvertex, nbFaces_); 
     732      face_edges.resize(nvertex, nbFaces_); 
     733 
     734      // For determining the global face index 
     735      size_t nbFacesOnProc = nbFaces_; 
     736      size_t nbFacesAccum; 
     737      MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     738      nbFacesAccum -= nbFaces_; 
    716739 
    717740      // Case (1): edges have been previously generated 
     
    720743        // (1.1) Recuperating node global indexing and saving face_nodes 
    721744        vector<size_t> hashValues(4); 
    722         CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
    723         for (int nf = 0; nf < nbFaces; ++nf) 
     745        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4); 
     746        for (int nf = 0; nf < nbFaces_; ++nf) 
    724747        { 
    725748          for (int nv = 0; nv < nvertex; ++nv) 
     
    733756        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap(); 
    734757        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
    735         CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
    736         for (int nf = 0; nf < nbFaces; ++nf) 
     758        CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 
     759        size_t nEdge = 0; 
     760        for (int nf = 0; nf < nbFaces_; ++nf) 
    737761        { 
    738762          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    754778            } 
    755779            face_nodes(nv1,nf) = it1->second[0]; 
    756             edgeHashList(nf*nvertex + nv1) = hashPairOrdered(it1->second[0], it2->second[0]); 
    757  
    758           } 
    759         } 
     780            if (it1->second[0] != it2->second[0]) 
     781            { 
     782              edgeHashList(nEdge) = hashPairOrdered(it1->second[0], it2->second[0]); 
     783              ++nEdge; 
     784            } 
     785          } 
     786        } 
     787        edgeHashList.resizeAndPreserve(nEdge); 
    760788 
    761789        // (1.2) Recuperating edge global indexing and saving face_edges 
     
    764792        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash; 
    765793        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank; 
    766  
    767794        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
    768         CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex); 
     795        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 
    769796        size_t iIdx = 0; 
    770797 
    771         for (int nf = 0; nf < nbFaces; ++nf) 
     798        for (int nf = 0; nf < nbFaces_; ++nf) 
    772799        { 
    773800          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    788815              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
    789816            } 
    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); 
     817            if (it1->second[0] != it2->second[0]) 
     818            { 
     819              size_t faceIdxGlo = nbFacesAccum + nf; 
     820              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
     821              itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 
     822              size_t edgeIdxGlo = itEdgeHash->second[0]; 
     823              face_edges(nv1,nf) = edgeIdxGlo; 
     824              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
     825              { 
     826                edgeIdxGloList(iIdx) = edgeIdxGlo; 
     827                ++iIdx; 
     828              } 
     829              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
     830              edgeHash2Rank[edgeHash].push_back(itEdgeHash->first); 
     831              edgeHash2Rank[edgeHash].push_back(mpiRank); 
     832            } 
     833            else 
     834            { 
     835              face_edges(nv1,nf) = 999999; 
     836            } 
    803837          } 
    804838        } 
     
    812846        dhtEdgeHash.computeIndexInfoMapping(edgeHashList); 
    813847        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap(); 
    814         CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo;            // map is needed purely for counting 
     848        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo; 
    815849 
    816850        // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>> 
    817         // edgeHashMin2IdxGlo edgeHash2Info = <edgeHashMin, idxGlo> 
     851        // edgeHashMin2IdxGlo = <edgeHashMin, idxGlo> 
    818852        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it) 
    819853        { 
     
    831865 
    832866        edge_faces.resize(2, edge_count); 
    833         face_faces.resize(nvertex, nbFaces); 
     867        face_faces.resize(nvertex, nbFaces_); 
    834868 
    835869        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     
    838872        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2; 
    839873 
    840         for (int nf = 0; nf < nbFaces; ++nf) 
     874        for (int nf = 0; nf < nbFaces_; ++nf) 
    841875        { 
    842876          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    857891              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2)); 
    858892            } 
    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) 
     893 
     894            if (it1->second[0] != it2->second[0]) 
     895            { 
     896              size_t faceIdxGlo = nbFacesAccum + nf; 
     897              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]); 
     898              itEdgeHash = edgeHash2IdxGlo.find(edgeHash); 
     899              if (itEdgeHash != edgeHashMin2IdxGlo.end()) 
    870900              { 
    871                 edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    872                 edge_faces(1, edgeIdxGlo - edge_start) = -999; 
    873                 face_faces(nv1, nf) = -1; 
    874               } 
     901                int edgeIdxGlo = itEdgeHash->second[0]; 
     902                itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     903                int face1 = itFace1->second[0]; 
     904                if (itFace1->second.size() == 1) 
     905                { 
     906                  edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     907                  edge_faces(1, edgeIdxGlo - edge_start) = -999; 
     908                  face_faces(nv1, nf) = 999999; 
     909                } 
     910                else 
     911                { 
     912                  int face2 = itFace1->second[1]; 
     913                  edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     914                  edge_faces(1, edgeIdxGlo - edge_start) = face2; 
     915                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     916                } 
     917              } // edge owner 
    875918              else 
    876919              { 
     920                itEdgeHash = edgeHashMin2IdxGlo.find(edgeHash); 
     921                size_t edgeIdxGlo = itEdgeHash->second[0]; 
     922                itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     923                int face1 = itFace1->second[0]; 
    877924                int face2 = itFace1->second[1]; 
    878                 edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    879                 edge_faces(1, edgeIdxGlo - edge_start) = face2; 
    880925                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
    881               } 
    882             } 
     926              } // not an edge owner 
     927            } // node1 != node2 
    883928            else 
    884929            { 
    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); 
     930              face_faces(nv1, nf) = 999999; 
    889931            } 
    890932          } 
     
    897939        // (2.1) Generating nodeHashList 
    898940        vector<size_t> hashValues(4); 
    899         CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
    900         for (int nf = 0; nf < nbFaces; ++nf) 
     941        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4); 
     942        for (int nf = 0; nf < nbFaces_; ++nf) 
    901943        { 
    902944          for (int nv = 0; nv < nvertex; ++nv) 
     
    914956        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 
    915957        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
    916         CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
    917         for (int nf = 0; nf < nbFaces; ++nf) 
     958        CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 
     959        for (int nf = 0; nf < nbFaces_; ++nf) 
    918960        { 
    919961          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    955997        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 
    956998        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 
    957         CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex); 
     999        CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex); 
    9581000        size_t iIdxMin = 0; 
    9591001 
     
    10001042        edge_lat.resize(edge_count); 
    10011043        edge_nodes.resize(2, edge_count); 
    1002         face_edges.resize(nvertex, nbFaces); 
     1044        face_edges.resize(nvertex, nbFaces_); 
    10031045 
    10041046        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face; 
    1005         CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex); 
     1047        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 
    10061048        size_t iIdx = 0; 
    10071049 
    1008         for (int nf = 0; nf < nbFaces; ++nf) 
     1050        for (int nf = 0; nf < nbFaces_; ++nf) 
    10091051        { 
    10101052          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    10301072            size_t nodeIdxGlo2 = it2->second[0]; 
    10311073            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             } 
     1074            if (nodeIdxGlo1 != nodeIdxGlo2) 
     1075            { 
     1076              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
     1077              size_t ownerIdx = (itIdx->second)[0]; 
     1078              int edgeIdxGlo = 0; 
     1079              size_t faceIdxGlo = nbFacesAccum + nf; 
     1080 
     1081              if (myIdx == ownerIdx) 
     1082              { 
     1083                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1084                edgeIdxGlo = (itIdxGlo->second)[0]; 
     1085                double edgeLon; 
     1086                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf)); 
     1087                if (diffLon < (180.- prec)) 
     1088                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5; 
     1089                else if (diffLon > (180.+ prec)) 
     1090                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.; 
     1091                else 
     1092                  edgeLon = 0.; 
     1093                edge_lon(edgeIdxGlo - edge_start) = edgeLon; 
     1094                edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 
     1095                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 
     1096                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 
     1097              } 
     1098              else 
     1099              { 
     1100                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1101                edgeIdxGlo = (itIdxGlo->second)[0]; 
     1102              } 
     1103              face_edges(nv1,nf) = edgeIdxGlo; 
     1104              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
     1105              { 
     1106                edgeIdxGloList(iIdx) = edgeIdxGlo; 
     1107                ++iIdx; 
     1108              } 
     1109              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
     1110            } // nodeIdxGlo1 != nodeIdxGlo2 
    10481111            else 
    10491112            { 
    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); 
     1113              face_edges(nv1,nf) = 999999; 
     1114            } 
    10601115          } 
    10611116        } 
     
    10641119        // (2.4) Saving remaining variables edge_faces and face_faces 
    10651120        edge_faces.resize(2, edge_count); 
    1066         face_faces.resize(nvertex, nbFaces); 
     1121        face_faces.resize(nvertex, nbFaces_); 
    10671122 
    10681123        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     
    10711126        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 
    10721127 
    1073         for (int nf = 0; nf < nbFaces; ++nf) 
     1128        for (int nf = 0; nf < nbFaces_; ++nf) 
    10741129        { 
    10751130          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    10971152            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx); 
    10981153            size_t ownerIdx = (itIdx->second)[0]; 
    1099             size_t faceIdxGlo = mpiRank * nbFaces + nf; 
     1154            size_t faceIdxGlo = nbFacesAccum + nf; 
    11001155 
    11011156            if (myIdx == ownerIdx) 
     
    11091164                edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    11101165                edge_faces(1, edgeIdxGlo - edge_start) = -999; 
    1111                 face_faces(nv1, nf) = -1; 
     1166                face_faces(nv1, nf) = 999999; 
    11121167              } 
    11131168              else 
     
    11381193        vector<size_t> hashValues(4); 
    11391194        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
    1140         CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
     1195        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4); 
    11411196        size_t iHash = 0; 
    1142         for (int nf = 0; nf < nbFaces; ++nf) 
     1197        for (int nf = 0; nf < nbFaces_; ++nf) 
    11431198        { 
    11441199          for (int nv = 0; nv < nvertex; ++nv) 
     
    11731228        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
    11741229        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo; 
    1175         CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4); 
     1230        CArray<size_t,1> nodeIdxMinList(nbFaces_*nvertex*4); 
    11761231        size_t iIdxMin = 0; 
    11771232 
     
    12211276        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx; 
    12221277        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2; 
    1223         CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
    1224  
    1225         for (int nf = 0; nf < nbFaces; ++nf) 
     1278        CArray<size_t,1> edgeHashList(nbFaces_*nvertex); 
     1279        size_t nEdgeHash = 0; 
     1280 
     1281        for (int nf = 0; nf < nbFaces_; ++nf) 
    12261282        { 
    12271283          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    12501306            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]); 
    12511307            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; 
     1308            if (nodeIdxGlo1 != nodeIdxGlo2) 
     1309            { 
     1310              size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2); 
     1311              edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank)); 
     1312              edgeHash2Idx[edgeHash].push_back(mpiRank); 
     1313              edgeHashList(nEdgeHash) = edgeHash; 
     1314              ++nEdgeHash; 
     1315            } 
    12561316            face_nodes(nv1,nf) = nodeIdxGlo1; 
    12571317          } 
    12581318        } 
     1319        edgeHashList.resizeAndPreserve(nEdgeHash); 
    12591320 
    12601321        // (3.4) Generating global edge indexes 
     
    12701331        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin; 
    12711332        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo; 
    1272         CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex); 
     1333        CArray<size_t,1> edgeIdxMinList(nbFaces_*nvertex); 
    12731334        iIdxMin = 0; 
    12741335 
     
    13161377        edge_lat.resize(edge_count); 
    13171378        edge_nodes.resize(2, edge_count); 
    1318         face_edges.resize(nvertex, nbFaces); 
     1379        face_edges.resize(nvertex, nbFaces_); 
    13191380 
    13201381        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2; 
    13211382        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) 
     1383        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex); 
     1384        size_t nEdge = 0; 
     1385 
     1386        for (int nf = 0; nf < nbFaces_; ++nf) 
    13261387        { 
    13271388          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    13411402            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0]; 
    13421403 
    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             } 
     1404            if (nodeIdxGlo1 != nodeIdxGlo2) 
     1405            { 
     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; 
     1410              size_t faceIdxGlo = nbFacesAccum + nf; 
     1411 
     1412              if (myIdx == ownerIdx) 
     1413              { 
     1414                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1415                edgeIdxGlo = (itIdxGlo->second)[0]; 
     1416                double edgeLon; 
     1417                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf)); 
     1418                if (diffLon < (180.- prec)) 
     1419                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5; 
     1420                else if (diffLon > (180.+ prec)) 
     1421                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.; 
     1422                else 
     1423                  edgeLon = 0.; 
     1424                edge_lon(edgeIdxGlo - edge_start) = edgeLon; 
     1425                edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5; 
     1426                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1; 
     1427                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2; 
     1428              } 
     1429              else 
     1430              { 
     1431                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1432                edgeIdxGlo = (itIdxGlo->second)[0]; 
     1433              } 
     1434              face_edges(nv1,nf) = edgeIdxGlo; 
     1435              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0) 
     1436              { 
     1437                edgeIdxGloList(nEdge) = edgeIdxGlo; 
     1438                ++nEdge; 
     1439              } 
     1440              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo); 
     1441            } // nodeIdxGlo1 != nodeIdxGlo2 
    13611442            else 
    13621443            { 
    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); 
     1444              face_edges(nv1,nf) = 999999; 
     1445            } 
     1446          } 
     1447        } 
     1448        edgeIdxGloList.resizeAndPreserve(nEdge); 
    13761449 
    13771450        // (3.6) Saving remaining variables edge_faces and face_faces 
    13781451        edge_faces.resize(2, edge_count); 
    1379         face_faces.resize(nvertex, nbFaces); 
     1452        face_faces.resize(nvertex, nbFaces_); 
    13801453 
    13811454        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm); 
     
    13831456        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap(); 
    13841457 
    1385         for (int nf = 0; nf < nbFaces; ++nf) 
     1458        for (int nf = 0; nf < nbFaces_; ++nf) 
    13861459        { 
    13871460          for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     
    13941467            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
    13951468            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) 
     1469            if (myNodeIdx1 != myNodeIdx2) 
     1470            { 
     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]; 
     1481              size_t faceIdxGlo = nbFacesAccum + nf; 
     1482 
     1483              if (myIdx == ownerIdx) 
    14151484              { 
    1416                 edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    1417                 edge_faces(1, edgeIdxGlo - edge_start) = -999; 
    1418                 face_faces(nv1, nf) = -1; 
    1419               } 
     1485                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx); 
     1486                int edgeIdxGlo = (itIdxGlo->second)[0]; 
     1487                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     1488                int face1 = it1->second[0]; 
     1489                if (it1->second.size() == 1) 
     1490                { 
     1491                  edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     1492                  edge_faces(1, edgeIdxGlo - edge_start) = -999; 
     1493                  face_faces(nv1, nf) = 999999; 
     1494                } 
     1495                else 
     1496                { 
     1497                  size_t face2 = it1->second[1]; 
     1498                  edge_faces(0, edgeIdxGlo - edge_start) = face1; 
     1499                  edge_faces(1, edgeIdxGlo - edge_start) = face2; 
     1500                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
     1501                } 
     1502              } // myIdx == ownerIdx 
    14201503              else 
    14211504              { 
    1422                 size_t face2 = it1->second[1]; 
    1423                 edge_faces(0, edgeIdxGlo - edge_start) = face1; 
    1424                 edge_faces(1, edgeIdxGlo - edge_start) = face2; 
     1505                CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx); 
     1506                size_t edgeIdxGlo = (itIdxGlo->second)[0]; 
     1507                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo); 
     1508                int face1 = it1->second[0]; 
     1509                int face2 = it1->second[1]; 
    14251510                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1); 
    1426               } 
    1427             } 
     1511              } // myIdx != ownerIdx 
     1512            } // myNodeIdx1 != myNodeIdx2 
    14281513            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             } 
     1514              face_faces(nv1, nf) = 999999; 
    14371515          } 
    14381516        } 
     
    14431521  } // createMeshEpsilon 
    14441522 
     1523  ///---------------------------------------------------------------- 
     1524  /*! 
     1525   * \fn  void CMesh::getNghbFacesNodeType(const MPI_Comm& comm, 
     1526                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1527                               CArray<int, 2>& nghbFaces) 
     1528   * Finds neighboring cells of a local domain for node-type of neighbors. 
     1529   * \param [in] comm 
     1530   * \param [in] bounds_lon Array of boundary longitudes. 
     1531   * \param [in] bounds_lat Array of boundary latitudes. 
     1532   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs. 
     1533   */ 
     1534 
     1535  void CMesh::getNghbFacesNodeType(const MPI_Comm& comm, 
     1536                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1537                               CArray<int, 2>& nghbFaces) 
     1538  { 
     1539    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     1540    int nbFaces = bounds_lon.shape()[1]; 
     1541    nghbFaces.resize(2, nbFaces*10);    // some estimate on max number of neighbouring cells 
     1542 
     1543    int mpiRank, mpiSize; 
     1544    MPI_Comm_rank(comm, &mpiRank); 
     1545    MPI_Comm_size(comm, &mpiSize); 
     1546 
     1547    // For determining the global face index 
     1548    size_t nbFacesOnProc = nbFaces; 
     1549    size_t nbFacesAccum; 
     1550    MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     1551    nbFacesAccum -= nbFaces; 
     1552 
     1553    // (1) Generating unique node indexes 
     1554    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> > 
     1555    vector<size_t> hashValues(4); 
     1556    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
     1557    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
     1558    size_t iIdx = 0; 
     1559    for (int nf = 0; nf < nbFaces; ++nf) 
     1560    { 
     1561      for (int nv = 0; nv < nvertex; ++nv) 
     1562      { 
     1563        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     1564        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 
     1565        for (int nh = 0; nh < 4; ++nh) 
     1566        { 
     1567          if (nodeHash2Idx.count(hashValues[nh])==0) 
     1568         { 
     1569            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex); 
     1570            nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 
     1571           nodeHashList(iIdx) = hashValues[nh]; 
     1572           ++iIdx; 
     1573         } 
     1574        } 
     1575      } 
     1576    } 
     1577    nodeHashList.resizeAndPreserve(iIdx); 
     1578 
     1579    // (1.2) Generating node indexes 
     1580    // The ownership criterion: priority of the process holding the smaller index 
     1581    // Maps generated in this step are: 
     1582    // nodeHash2Info = <hash, idx1, idx2, idx3....> 
     1583    // nodeIdx2IdxMin = <idx, idxMin> 
     1584    // idxMin is a unique node identifier 
     1585 
     1586    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
     1587    dhtNodeHash.computeIndexInfoMapping(nodeHashList); 
     1588    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
     1589 
     1590    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
     1591 
     1592    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
     1593    { 
     1594      size_t idxMin = (it->second)[0]; 
     1595      size_t idx = (it->second)[0]; 
     1596      for (int i = 2; i < (it->second).size();) 
     1597      { 
     1598        if (mpiRank == (it->second)[i+1]) 
     1599        { 
     1600          idx = (it->second)[i]; 
     1601        } 
     1602        if ((it->second)[i] < idxMin) 
     1603        { 
     1604          idxMin = (it->second)[i]; 
     1605          (it->second)[i] = (it->second)[i-2]; 
     1606          (it->second)[i+1] = (it->second)[i-1]; 
     1607        } 
     1608        i += 2; 
     1609      } 
     1610      (it->second)[0] = idxMin; 
     1611      if (nodeIdx2IdxMin.count(idx) == 0) 
     1612      { 
     1613        nodeIdx2IdxMin[idx].push_back(idxMin); 
     1614      } 
     1615    } 
     1616 
     1617    // (2) Creating maps nodeIdxMin2Face = <nodeIdxMin, [face1, face2, ...]> 
     1618    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it; 
     1619    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdxMin2Face; 
     1620    CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4); 
     1621 
     1622    size_t nNode = 0; 
     1623 
     1624    for (int nf = 0; nf < nbFaces; ++nf) 
     1625    { 
     1626      for (int nv = 0; nv < nvertex; ++nv) 
     1627      { 
     1628        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     1629        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank); 
     1630        it = nodeIdx2IdxMin.find(myNodeIdx); 
     1631        size_t nodeIdxMin = (it->second)[0]; 
     1632        size_t faceIdx = nbFacesAccum + nf; 
     1633        if (nodeIdxMin2Face.count(nodeIdxMin) == 0) 
     1634        { 
     1635          nodeIdxMinList(nNode) = nodeIdxMin; 
     1636          ++nNode; 
     1637        } 
     1638        nodeIdxMin2Face[nodeIdxMin].push_back(faceIdx); 
     1639        nodeIdxMin2Face[nodeIdxMin].push_back(mpiRank); 
     1640      } 
     1641    } 
     1642    nodeIdxMinList.resizeAndPreserve(nNode); 
     1643 
     1644    // (3) Face_face connectivity 
     1645 
     1646    // nodeIdxMin2Info = <nodeIdxMin, [face1, face2,...]> 
     1647    CClientClientDHTSizet dhtNode2Face (nodeIdxMin2Face, comm); 
     1648    dhtNode2Face.computeIndexInfoMapping(nodeIdxMinList); 
     1649    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2Info = dhtNode2Face.getInfoIndexMap(); 
     1650    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map 
     1651 
     1652    int nbNghb = 0; 
     1653    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode; 
     1654 
     1655    for (int nf = 0; nf < nbFaces; ++nf) 
     1656    { 
     1657      for (int nv = 0; nv < nvertex; ++nv) 
     1658      { 
     1659        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     1660        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank); 
     1661        itNode = nodeIdx2IdxMin.find(myNodeIdx); 
     1662        size_t nodeIdxMin = (itNode->second)[0]; 
     1663 
     1664        itNode = nodeIdxMin2Info.find(nodeIdxMin); 
     1665        for (int i = 0; i < itNode->second.size();) 
     1666        { 
     1667          size_t face = itNode->second[i]; 
     1668          size_t rank = itNode->second[i+1]; 
     1669          if (rank != mpiRank) 
     1670            if (mapFaces.count(face) == 0) 
     1671            { 
     1672              nghbFaces(0, nbNghb) = face; 
     1673              nghbFaces(1, nbNghb) = rank; 
     1674              ++nbNghb; 
     1675              mapFaces[face].push_back(face); 
     1676            } 
     1677          i += 2; 
     1678        } 
     1679      } 
     1680    } 
     1681    nghbFaces.resizeAndPreserve(2, nbNghb); 
     1682  } // getNghbFacesNodeType 
     1683 
     1684  ///---------------------------------------------------------------- 
     1685  /*! 
     1686   * \fn  void CMesh::getNghbFacesEdgeType(const MPI_Comm& comm, 
     1687                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1688                               CArray<int, 2>& nghbFaces) 
     1689   * Finds neighboring cells of a local domain for edge-type of neighbors. 
     1690   * \param [in] comm 
     1691   * \param [in] bounds_lon Array of boundary longitudes. 
     1692   * \param [in] bounds_lat Array of boundary latitudes. 
     1693   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs. 
     1694   */ 
     1695 
     1696  void CMesh::getNghbFacesEdgeType(const MPI_Comm& comm, 
     1697                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1698                               CArray<int, 2>& nghbFaces) 
     1699  { 
     1700    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     1701    int nbFaces = bounds_lon.shape()[1]; 
     1702    nghbFaces.resize(2, nbFaces*10);    // estimate of max number of neighbouring cells 
     1703 
     1704    int mpiRank, mpiSize; 
     1705    MPI_Comm_rank(comm, &mpiRank); 
     1706    MPI_Comm_size(comm, &mpiSize); 
     1707 
     1708    // For determining the global face index 
     1709    size_t nbFacesOnProc = nbFaces; 
     1710    size_t nbFacesAccum; 
     1711    MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
     1712    nbFacesAccum -= nbFaces; 
     1713 
     1714    // (1) Generating unique node indexes 
     1715    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> > 
     1716    vector<size_t> hashValues(4); 
     1717    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx; 
     1718    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4); 
     1719    size_t iIdx = 0; 
     1720    for (int nf = 0; nf < nbFaces; ++nf) 
     1721    { 
     1722      for (int nv = 0; nv < nvertex; ++nv) 
     1723      { 
     1724        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf)); 
     1725        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank); 
     1726        for (int nh = 0; nh < 4; ++nh) 
     1727        { 
     1728          if (nodeHash2Idx.count(hashValues[nh])==0) 
     1729          { 
     1730            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex); 
     1731            nodeHash2Idx[hashValues[nh]].push_back(mpiRank); 
     1732            nodeHashList(iIdx) = hashValues[nh]; 
     1733            ++iIdx; 
     1734          } 
     1735        } 
     1736      } 
     1737    } 
     1738    nodeHashList.resizeAndPreserve(iIdx); 
     1739 
     1740    // (1.2) Generating node indexes 
     1741    // The ownership criterion: priority of the process holding the smaller index 
     1742    // Maps generated in this step are: 
     1743    // nodeHash2Info = <hash, idx1, idx2, idx3....> 
     1744    // nodeIdx2IdxMin = <idx, idxMin> 
     1745    // idxMin is a unique node identifier 
     1746 
     1747    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm); 
     1748    dhtNodeHash.computeIndexInfoMapping(nodeHashList); 
     1749    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap(); 
     1750 
     1751    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin; 
     1752 
     1753    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it) 
     1754    { 
     1755      size_t idxMin = (it->second)[0]; 
     1756      size_t idx = (it->second)[0]; 
     1757      for (int i = 2; i < (it->second).size();) 
     1758      { 
     1759        if (mpiRank == (it->second)[i+1]) 
     1760        { 
     1761          idx = (it->second)[i]; 
     1762        } 
     1763        if ((it->second)[i] < idxMin) 
     1764        { 
     1765          idxMin = (it->second)[i]; 
     1766          (it->second)[i] = (it->second)[i-2]; 
     1767          (it->second)[i+1] = (it->second)[i-1]; 
     1768        } 
     1769        i += 2; 
     1770      } 
     1771      (it->second)[0] = idxMin; 
     1772      if (nodeIdx2IdxMin.count(idx) == 0) 
     1773      { 
     1774        nodeIdx2IdxMin[idx].push_back(idxMin); 
     1775      } 
     1776    } 
     1777 
     1778    // (2) Creating map edgeHash2Face = <edgeHash, [[face1, rank1], [face2, rank2]]>, where rank1 = rank2 = ... 
     1779 
     1780    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2, it; 
     1781    CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Face; 
     1782    CArray<size_t,1> edgeHashList(nbFaces*nvertex); 
     1783 
     1784    size_t nEdge = 0; 
     1785 
     1786    for (int nf = 0; nf < nbFaces; ++nf) 
     1787    { 
     1788      for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1789      { 
     1790        // Getting indexes of edge's nodes 
     1791        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1792        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
     1793        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
     1794        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
     1795        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
     1796        it1 = nodeIdx2IdxMin.find(myNodeIdx1); 
     1797        it2 = nodeIdx2IdxMin.find(myNodeIdx2); 
     1798        size_t nodeIdxMin1 = (it1->second)[0]; 
     1799        size_t nodeIdxMin2 = (it2->second)[0]; 
     1800        size_t faceIdx = nbFacesAccum + nf; 
     1801 
     1802        if (nodeIdxMin1 != nodeIdxMin2) 
     1803        { 
     1804          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2); 
     1805          if (edgeHash2Face.count(edgeHash) == 0) 
     1806          { 
     1807            edgeHashList(nEdge) = edgeHash; 
     1808            ++nEdge; 
     1809          } 
     1810          edgeHash2Face[edgeHash].push_back(faceIdx); 
     1811          edgeHash2Face[edgeHash].push_back(mpiRank); 
     1812        } // nodeIdxMin1 != nodeIdxMin2 
     1813      } 
     1814    } 
     1815    edgeHashList.resizeAndPreserve(nEdge); 
     1816 
     1817    // (3) Face_face connectivity 
     1818 
     1819    int nbNghb = 0; 
     1820    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode1, itNode2; 
     1821 
     1822    // edgeHash2Info = <edgeHash, [[face1, rank1], [face2, rank2]]> 
     1823    CClientClientDHTSizet dhtEdge2Face (edgeHash2Face, comm); 
     1824    dhtEdge2Face.computeIndexInfoMapping(edgeHashList); 
     1825    CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdge2Face.getInfoIndexMap(); 
     1826    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map 
     1827 
     1828    for (int nf = 0; nf < nbFaces; ++nf) 
     1829    { 
     1830      nbNghb = 0; 
     1831      for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     1832      { 
     1833        // Getting indexes of edge's nodes 
     1834        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     1835        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf)); 
     1836        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf)); 
     1837 
     1838        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank); 
     1839        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank); 
     1840        itNode1 = nodeIdx2IdxMin.find(myNodeIdx1); 
     1841        itNode2 = nodeIdx2IdxMin.find(myNodeIdx2); 
     1842        size_t nodeIdxMin1 = (itNode1->second)[0]; 
     1843        size_t nodeIdxMin2 = (itNode2->second)[0]; 
     1844 
     1845        if (nodeIdxMin1 != nodeIdxMin2) 
     1846        { 
     1847          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2); 
     1848          it1 = edgeHash2Info.find(edgeHash); 
     1849 
     1850          for (int i = 0; i < it1->second.size();) 
     1851          { 
     1852            size_t face = it1->second[i]; 
     1853            size_t rank = it1->second[i+1]; 
     1854            if (rank != mpiRank) 
     1855              if (mapFaces.count(face) == 0) 
     1856              { 
     1857                nghbFaces(0, nbNghb) = face; 
     1858                nghbFaces(1, nbNghb) = rank; 
     1859                ++nbNghb; 
     1860                mapFaces[face].push_back(face); 
     1861              } 
     1862            i += 2; 
     1863          } 
     1864        } // nodeIdxMin1 != nodeIdxMin2 
     1865      } 
     1866    } 
     1867    nghbFaces.resizeAndPreserve(2, nbNghb); 
     1868  } // getNghbFacesEdgeType 
     1869 
     1870  ///---------------------------------------------------------------- 
     1871  /*! 
     1872   * \fn void getDomainExtended(const int nghbType, const MPI_Comm& comm, 
     1873                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1874                                CArray<size_t, 1>& nghbFaces) 
     1875   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges. 
     1876   * \param [in] comm 
     1877   * \param [in] bounds_lon Array of boundary longitudes. 
     1878   * \param [in] bounds_lat Array of boundary latitudes. 
     1879   * \param [out] nghbFaces Array containing neighboring faces and their ranks. The ordering is face1, rank1,.... 
     1880   */ 
     1881 
     1882  void CMesh::getDomainExtended(const int nghbType, const MPI_Comm& comm, 
     1883      const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1884      CArray<int, 2>& nghbFaces) 
     1885  { 
     1886    if (nghbType == 0) 
     1887      getNghbFacesNodeType(comm, bounds_lon, bounds_lat, nghbFaces); 
     1888    else 
     1889      getNghbFacesEdgeType(comm, bounds_lon, bounds_lat, nghbFaces); 
     1890  } // getDomainExtended 
     1891 
    14451892} // namespace xios 
  • XIOS/trunk/src/node/mesh.hpp

    r924 r929  
    5757      CArray<int, 2> face_faces; 
    5858 
    59       void createMesh(const CArray<double, 1>&, const CArray<double, 1>&,  
     59      void createMesh(const CArray<double, 1>&, const CArray<double, 1>&, 
    6060                      const CArray<double, 2>&, const CArray<double, 2>& ); 
    6161                         
     
    6363                             const CArray<double, 1>&, const CArray<double, 1>&, 
    6464                             const CArray<double, 2>&, const CArray<double, 2>& ); 
     65 
     66      void getDomainExtended(const int, const MPI_Comm&, 
     67                              const CArray<double, 2>&, const CArray<double, 2>&, 
     68                              CArray<int, 2>&); 
    6569             
    6670      static CMesh* getMesh(StdString, int); 
     
    6872    private: 
    6973 
    70       int nbNodes; 
    71       int nbEdges; 
    72       int nbFaces; 
     74      int nbNodes_; 
     75      int nbEdges_; 
     76      int nbFaces_; 
    7377 
    7478      static std::map <StdString, CMesh> meshList; 
     
    7680      CClientClientDHTSizet* pNodeGlobalIndex;                    // pointer to a map <nodeHash, nodeIdxGlo> 
    7781      CClientClientDHTSizet* pEdgeGlobalIndex;                    // pointer to a map <edgeHash, edgeIdxGlo> 
     82      void getNghbFacesNodeType(const MPI_Comm&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 
     83      void getNghbFacesEdgeType(const MPI_Comm&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 
    7884 
    7985      vector<size_t> createHashes (const double, const double); 
Note: See TracChangeset for help on using the changeset viewer.