Changeset 924 for XIOS


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

Parallel version of UGRID norms.

The following connectivity parameters have been implemented:

edge_nodes
face_nodes
face_edges
edge_faces
face_faces.

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

Location:
XIOS/trunk
Files:
4 added
10 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/src/client_client_dht_decl.cpp

    r832 r924  
    1313 
    1414template class CClientClientDHTTemplate<int>; 
     15template class CClientClientDHTTemplate<size_t>; 
    1516template class CClientClientDHTTemplate<PairIntInt>; 
    1617 

  • XIOS/trunk/src/client_client_dht_template.hpp

    r892 r924  
    125125 
    126126typedef CClientClientDHTTemplate<int> CClientClientDHTInt; 
     127typedef CClientClientDHTTemplate<size_t> CClientClientDHTSizet; 
    127128typedef CClientClientDHTTemplate<PairIntInt> CClientClientDHTPairIntInt; 
    128129 

  • XIOS/trunk/src/dht_auto_indexing.cpp

    r892 r924  
    1212{ 
    1313 
    14 /*! 
    15   Constructor with initial distribution information and the corresponding index 
    16   Each process holds a piece of hash value, which has no information (global index) associated. 
    17 The constructor tries to assign this information to each process basing on their rank. 
    18 The global index (information) are assigned across processes in the incremental manner. 
    19 The process has lower rank will have low global index. 
    20   \param [in] hashValue Carray contains hash value 
    21   \param [in] clientIntraComm communicator of clients 
    22 */ 
     14  CDHTAutoIndexing::~CDHTAutoIndexing() 
     15  { 
     16  } 
    2317 
    24 CDHTAutoIndexing::CDHTAutoIndexing(const CArray<size_t,1>& hashValue, 
    25                                    const MPI_Comm& clientIntraComm) 
    26   : CClientClientDHTTemplate<size_t>(clientIntraComm) 
    27 { 
    28   Index2VectorInfoTypeMap indexToVecInfoMap; 
    29   indexToVecInfoMap.rehash(std::ceil(hashValue.size()/indexToVecInfoMap.max_load_factor())); 
    30   CArray<size_t,1>::const_iterator it = hashValue.begin(), ite = hashValue.end(); 
     18  /*! 
     19    \param [in] 
     20    \param [in] 
     21  */ 
    3122 
    32   // Just use a dummy value to initialize 
    33   int nbLvl = this->getNbLevel(); 
    34   for (; it != ite; ++it) 
     23  CDHTAutoIndexing::CDHTAutoIndexing(const CArray<size_t,1>& hashValue, 
     24                                     const MPI_Comm& clientIntraComm) 
     25    : CClientClientDHTTemplate<size_t>(clientIntraComm) 
    3526  { 
    36     indexToVecInfoMap[*it].resize(1); 
    37     indexToVecInfoMap[*it][0] = 0; 
     27 
     28    nbIndexOnProc_ = hashValue.size(); 
     29    size_t nbIndexAccum; 
     30    MPI_Scan(&nbIndexOnProc_, &nbIndexAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, clientIntraComm); 
     31 
     32    // Broadcasting the total number of indexes 
     33    int rank, size; 
     34    MPI_Comm_rank(clientIntraComm, &rank); 
     35    MPI_Comm_size(clientIntraComm, &size); 
     36    if (rank == (size-1)) nbIndexesGlobal_ = nbIndexAccum; 
     37    MPI_Bcast(&nbIndexesGlobal_, 1, MPI_UNSIGNED_LONG, size-1, clientIntraComm); 
     38 
     39    CArray<size_t,1>::const_iterator itbIdx = hashValue.begin(), itIdx, 
     40                                     iteIdx = hashValue.end(); 
     41 
     42    size_t idx = 0; 
     43    beginIndexOnProc_ = nbIndexAccum - nbIndexOnProc_; 
     44    globalIndex_.resize(nbIndexOnProc_); 
     45    for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx) 
     46    { 
     47      globalIndex_[idx] = beginIndexOnProc_ + idx; 
     48      ++idx ; 
     49    } 
    3850  } 
    39   computeDistributedIndex(indexToVecInfoMap, clientIntraComm, nbLvl-1); 
    4051 
    41   // Find out number of index on other proc 
    42   size_t nbIndexOnProc = index2InfoMapping_.size(); 
    43   size_t nbIndexAccum; 
     52  /*! 
     53   * \fn void CDHTAutoIndexing::CDHTAutoIndexing(Index2VectorInfoTypeMap& hashInitMap, const MPI_Comm& clientIntraComm) 
     54   * Assigns a global index to unique input indexes. 
     55   * The returned map has unique indexes as a key and global indexes as a mapped value. 
     56   * \param [in] hashInitMap map<size_t, vector<size_t>> is a map of unique indexes. 
     57   * \param [in] clientIntraComm 
     58  */ 
     59  CDHTAutoIndexing::CDHTAutoIndexing(Index2VectorInfoTypeMap& hashInitMap, 
     60                                     const MPI_Comm& clientIntraComm) 
     61    : CClientClientDHTTemplate<size_t>(clientIntraComm) 
     62  { 
    4463 
    45   MPI_Scan(&nbIndexOnProc, &nbIndexAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, clientIntraComm); 
     64    nbIndexOnProc_ = hashInitMap.size(); 
     65    size_t nbIndexAccum; 
     66    MPI_Scan(&nbIndexOnProc_, &nbIndexAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, clientIntraComm); 
    4667 
    47   Index2VectorInfoTypeMap::iterator itbIdx = index2InfoMapping_.begin(), itIdx, 
    48                                     iteIdx = index2InfoMapping_.end(); 
    49   size_t idx = 0; 
    50   size_t idxBegin = nbIndexAccum - nbIndexOnProc; 
    51   for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx) 
     68    int rank, size; 
     69    MPI_Comm_rank(clientIntraComm, &rank); 
     70    MPI_Comm_size(clientIntraComm, &size); 
     71    if (rank == (size-1)) nbIndexesGlobal_ = nbIndexAccum; 
     72    MPI_Bcast(&nbIndexesGlobal_, 1, MPI_UNSIGNED_LONG, size-1, clientIntraComm); 
     73 
     74    Index2VectorInfoTypeMap::iterator itbIdx = hashInitMap.begin(), itIdx, 
     75                                      iteIdx = hashInitMap.end(); 
     76    size_t idx = 0; 
     77    beginIndexOnProc_ = nbIndexAccum - nbIndexOnProc_; 
     78    globalIndex_.resize(nbIndexOnProc_); 
     79    for (itIdx = itbIdx; itIdx != iteIdx; ++itIdx) 
     80    { 
     81      (itIdx->second)[0] = beginIndexOnProc_ + idx; 
     82      globalIndex_[idx] = beginIndexOnProc_ + idx; 
     83      ++idx ; 
     84    } 
     85  } 
     86 
     87  /*! 
     88   * \fn size_t CDHTAutoIndexing::getNbIndexesGlobal() const 
     89   * Returns the total number of global indexes. 
     90  */ 
     91  size_t CDHTAutoIndexing::getNbIndexesGlobal() const 
    5292  { 
    53     (itIdx->second)[0] = idxBegin + idx; 
    54     ++idx ; 
     93    return nbIndexesGlobal_; 
    5594  } 
    56 } 
     95 
     96  /*! 
     97   * \fn size_t CDHTAutoIndexing::getIndexStart() const 
     98   * Returns the starting global index for a proc. 
     99  */ 
     100  size_t CDHTAutoIndexing::getIndexStart() const 
     101  { 
     102    return beginIndexOnProc_; 
     103  } 
     104 
     105  /*! 
     106   * \fn size_t CDHTAutoIndexing::getIndexCount() const 
     107   * Returns the number of global indexes on a proc. 
     108  */ 
     109  size_t CDHTAutoIndexing::getIndexCount() const 
     110  { 
     111    return nbIndexOnProc_; 
     112  } 
    57113 
    58114} 

  • XIOS/trunk/src/dht_auto_indexing.hpp

    r892 r924  
    2222class CDHTAutoIndexing: public CClientClientDHTTemplate<size_t> 
    2323{ 
    24 public: 
    25   CDHTAutoIndexing(const CArray<size_t,1>& hashValue, 
    26                    const MPI_Comm& clientIntraComm); 
     24  public: 
     25 
     26    CDHTAutoIndexing(const CArray<size_t,1>& hashValue, 
     27                     const MPI_Comm& clientIntraComm); 
     28 
     29    CDHTAutoIndexing(Index2VectorInfoTypeMap& hashInitMap, 
     30                     const MPI_Comm& clientIntraComm); 
     31 
     32    size_t getNbIndexesGlobal() const; 
     33    size_t getIndexStart() const; 
     34    size_t getIndexCount() const; 
    2735 
    2836    /** Default destructor */ 
    2937    virtual ~CDHTAutoIndexing(); 
     38 
     39  protected: 
     40    std::vector<size_t> globalIndex_; 
     41    size_t nbIndexesGlobal_; 
     42    size_t nbIndexOnProc_ ; 
     43    size_t beginIndexOnProc_ ; 
    3044}; 
    3145 

  • XIOS/trunk/src/io/nc4_data_output.cpp

    r902 r924  
    2626      CNc4DataOutput::CNc4DataOutput 
    2727         (const StdString & filename, bool exist, bool useClassicFormat, bool useCFConvention, 
    28           MPI_Comm comm_file,bool multifile, bool isCollective, const StdString& timeCounterName) 
     28          MPI_Comm comm_file, bool multifile, bool isCollective, const StdString& timeCounterName) 
    2929            : SuperClass() 
    3030            , SuperClassWriter(filename, exist, useClassicFormat, useCFConvention, &comm_file, multifile, timeCounterName) 
    … …  
    449449      StdString domid = domain->getDomainOutputName(); 
    450450      StdString domainName = domain->name; 
    451       domain->assignMesh(domainName); 
    452       domain->mesh->createMesh(domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv); 
    453       //domain->mesh->createMeshEpsilon(server->intraComm, domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv); 
     451      domain->assignMesh(domainName, domain->nvertex); 
     452      //domain->mesh->createMesh(domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv); 
     453      domain->mesh->createMeshEpsilon(server->intraComm, domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv); 
    454454 
    455455      StdString node_x = domainName + "_node_x"; 
    … …  
    479479        SuperClassWriter::addVariable(domainName, NC_INT, dim0); 
    480480        SuperClassWriter::addAttribute("cf_role", StdString("mesh_topology"), &domainName); 
     481        SuperClassWriter::addAttribute("long_name", StdString("Topology data of 2D unstructured mesh"), &domainName); 
     482        SuperClassWriter::addAttribute("topology_dimension", 2, &domainName); 
    481483        SuperClassWriter::addAttribute("node_coordinates", node_x + " " + node_y, &domainName); 
    482484      } 
    … …  
    498500                SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0); 
    499501                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x); 
    500                 SuperClassWriter::addAttribute("longname_name", StdString("Longitude of mesh nodes."), &node_x); 
     502                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x); 
    501503                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x); 
    502504                SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0); 
    503505                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y); 
    504                 SuperClassWriter::addAttribute("longname_name", StdString("Latitude of mesh nodes."), &node_y); 
     506                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y); 
    505507                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y); 
    506508              } 
    … …  
    512514              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y)) 
    513515              { 
    514                 SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodes); 
     516                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodesGlo); 
    515517                dim0.clear(); 
    516518                dim0.push_back(dimNode); 
    517519                SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0); 
    518520                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x); 
    519                 SuperClassWriter::addAttribute("longname_name", StdString("Longitude of mesh nodes."), &node_x); 
     521                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x); 
    520522                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x); 
    521523                SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0); 
    522524                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y); 
    523                 SuperClassWriter::addAttribute("longname_name", StdString("Latitude of mesh nodes."), &node_y); 
     525                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y); 
     526                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y); 
     527              } 
     528              SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName); 
     529              SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName); 
     530              SuperClassWriter::addDimension(dimEdge, domain->ni_glo); 
     531              dim0.clear(); 
     532              dim0.push_back(dimEdge); 
     533              SuperClassWriter::addVariable(edge_x, NC_FLOAT, dim0); 
     534              SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x); 
     535              SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh edges."), &edge_x); 
     536              SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x); 
     537              SuperClassWriter::addVariable(edge_y, NC_FLOAT, dim0); 
     538              SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y); 
     539              SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh edges."), &edge_y); 
     540              SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y); 
     541              dim0.clear(); 
     542              dim0.push_back(dimEdge); 
     543              dim0.push_back(dimTwo); 
     544              SuperClassWriter::addVariable(edge_nodes, NC_INT, dim0); 
     545              SuperClassWriter::addAttribute("cf_role", StdString("edge_node_connectivity"), &edge_nodes); 
     546              SuperClassWriter::addAttribute("long_name", StdString("Maps every edge/link to two nodes that it connects."), &edge_nodes); 
     547              SuperClassWriter::addAttribute("start_index", 0, &edge_nodes); 
     548            } // domain->nvertex == 2 
     549 
     550            // Adding faces, edges, and nodes, if edges and nodes have not been defined previously 
     551            if (domain->nvertex > 2) 
     552            { 
     553              // Nodes 
     554              if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y)) 
     555              { 
     556                SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodesGlo); 
     557                dim0.clear(); 
     558                dim0.push_back(dimNode); 
     559                SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0); 
     560                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x); 
     561                SuperClassWriter::addAttribute("long_name", StdString("Longitude of mesh nodes."), &node_x); 
     562                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x); 
     563                SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0); 
     564                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y); 
     565                SuperClassWriter::addAttribute("long_name", StdString("Latitude of mesh nodes."), &node_y); 
    524566                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y); 
    525567              } 
    526568              if (!SuperClassWriter::varExist(edge_x) || !SuperClassWriter::varExist(edge_y)) 
    527569              { 
     570                SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName); 
    528571                SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName); 
    529                 SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName); 
    530                 SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdges); 
     572                SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdgesGlo); 
    531573                dim0.clear(); 
    532574                dim0.push_back(dimEdge); 
    533575                SuperClassWriter::addVariable(edge_x, NC_FLOAT, dim0); 
    534576                SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x); 
    535                 SuperClassWriter::addAttribute("longname_name", StdString("Characteristic longitude of mesh edges."), &edge_x); 
     577                SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh edges."), &edge_x); 
    536578                SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x); 
    537579                SuperClassWriter::addVariable(edge_y, NC_FLOAT, dim0); 
    538580                SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y); 
    539                 SuperClassWriter::addAttribute("longname_name", StdString("Characteristic latitude of mesh edges."), &edge_y); 
     581                SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh edges."), &edge_y); 
    540582                SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y); 
    541583                dim0.clear(); 
    … …  
    547589                SuperClassWriter::addAttribute("start_index", 0, &edge_nodes); 
    548590              } 
    549             } // domain->nvertex == 2 
    550  
    551             // Adding faces, edges, and nodes, if edges and nodes have not been defined previously 
    552             if (domain->nvertex > 2) 
    553             { 
    554               // Nodes 
    555               if (!SuperClassWriter::varExist(node_x) || !SuperClassWriter::varExist(node_y)) 
    556               { 
    557                 SuperClassWriter::addDimension(dimNode, domain->mesh->nbNodes); 
    558                 dim0.clear(); 
    559                 dim0.push_back(dimNode); 
    560                 SuperClassWriter::addVariable(node_x, NC_FLOAT, dim0); 
    561                 SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &node_x); 
    562                 SuperClassWriter::addAttribute("longname_name", StdString("Longitude of mesh nodes."), &node_x); 
    563                 SuperClassWriter::addAttribute("units", StdString("degrees_east"), &node_x); 
    564                 SuperClassWriter::addVariable(node_y, NC_FLOAT, dim0); 
    565                 SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &node_y); 
    566                 SuperClassWriter::addAttribute("longname_name", StdString("Latitude of mesh nodes."), &node_y); 
    567                 SuperClassWriter::addAttribute("units", StdString("degrees_north"), &node_y); 
    568               } 
    569               if (!SuperClassWriter::varExist(edge_x) || !SuperClassWriter::varExist(edge_y)) 
    570               { 
    571                 SuperClassWriter::addAttribute("edge_coordinates", edge_x + " " + edge_y, &domainName); 
    572                 SuperClassWriter::addAttribute("edge_node_connectivity", edge_nodes, &domainName); 
    573                 SuperClassWriter::addDimension(dimEdge, domain->mesh->nbEdges); 
    574                 dim0.clear(); 
    575                 dim0.push_back(dimEdge); 
    576                 SuperClassWriter::addVariable(edge_x, NC_FLOAT, dim0); 
    577                 SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &edge_x); 
    578                 SuperClassWriter::addAttribute("longname_name", StdString("Characteristic longitude of mesh edges."), &edge_x); 
    579                 SuperClassWriter::addAttribute("units", StdString("degrees_east"), &edge_x); 
    580                 SuperClassWriter::addVariable(edge_y, NC_FLOAT, dim0); 
    581                 SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &edge_y); 
    582                 SuperClassWriter::addAttribute("longname_name", StdString("Characteristic latitude of mesh edges."), &edge_y); 
    583                 SuperClassWriter::addAttribute("units", StdString("degrees_north"), &edge_y); 
    584                 dim0.clear(); 
    585                 dim0.push_back(dimEdge); 
    586                 dim0.push_back(dimTwo); 
    587                 SuperClassWriter::addVariable(edge_nodes, NC_INT, dim0); 
    588                 SuperClassWriter::addAttribute("cf_role", StdString("edge_node_connectivity"), &edge_nodes); 
    589                 SuperClassWriter::addAttribute("long_name", StdString("Maps every edge/link to two nodes that it connects."), &edge_nodes); 
    590                 SuperClassWriter::addAttribute("start_index", 0, &edge_nodes); 
    591               } 
    592                 SuperClassWriter::addAttribute("face_coordinates", face_x + " " + face_y, &domainName); 
    593                 SuperClassWriter::addAttribute("face_node_connectivity", face_nodes, &domainName); 
    594                 SuperClassWriter::addDimension(dimFace, domain->mesh->nbFaces); 
    595                 SuperClassWriter::addDimension(dimVertex, domain->mesh->nvertex); 
    596                 dim0.clear(); 
    597                 dim0.push_back(dimFace); 
    598                 SuperClassWriter::addVariable(face_x, NC_FLOAT, dim0); 
    599                 SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &face_x); 
    600                 SuperClassWriter::addAttribute("longname_name", StdString("Characteristic longitude of mesh faces."), &face_x); 
    601                 SuperClassWriter::addAttribute("units", StdString("degrees_east"), &face_x); 
    602                 SuperClassWriter::addVariable(face_y, NC_FLOAT, dim0); 
    603                 SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &face_y); 
    604                 SuperClassWriter::addAttribute("longname_name", StdString("Characteristic latitude of mesh faces."), &face_y); 
    605                 SuperClassWriter::addAttribute("units", StdString("degrees_north"), &face_y); 
    606                 dim0.clear(); 
    607                 dim0.push_back(dimFace); 
    608                 dim0.push_back(dimVertex); 
    609                 SuperClassWriter::addVariable(face_nodes, NC_INT, dim0); 
    610                 SuperClassWriter::addAttribute("cf_role", StdString("face_node_connectivity"), &face_nodes); 
    611                 SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its corner nodes."), &face_nodes); 
    612                 SuperClassWriter::addAttribute("start_index", 0, &face_nodes); 
    613                 dim0.clear(); 
    614                 dim0.push_back(dimFace); 
    615                 dim0.push_back(dimVertex); 
    616                 SuperClassWriter::addVariable(face_edges, NC_INT, dim0); 
    617                 SuperClassWriter::addAttribute("cf_role", StdString("face_edge_connectivity"), &face_edges); 
    618                 SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges); 
    619                 SuperClassWriter::addAttribute("start_index", 0, &face_edges); 
    620                 dim0.clear(); 
    621                 dim0.push_back(dimEdge); 
    622                 dim0.push_back(dimTwo); 
    623                 SuperClassWriter::addVariable(edge_faces, NC_INT, dim0); 
    624                 SuperClassWriter::addAttribute("cf_role", StdString("edge_face connectivity"), &edge_faces); 
    625                 SuperClassWriter::addAttribute("long_name", StdString("neighbor faces for edges"), &edge_faces); 
    626                 SuperClassWriter::addAttribute("start_index", 0, &edge_faces); 
    627                 SuperClassWriter::addAttribute("_FillValue", -999, &edge_faces); 
    628                 SuperClassWriter::addAttribute("comment", StdString("missing neighbor faces are indicated using _FillValue"), &edge_faces); 
    629                 dim0.clear(); 
    630                 dim0.push_back(dimFace); 
    631                 dim0.push_back(dimVertex); 
    632                 SuperClassWriter::addVariable(face_faces, NC_INT, dim0); 
    633                 SuperClassWriter::addAttribute("cf_role", StdString("face_face connectivity"), &face_faces); 
    634                 SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces); 
    635                 SuperClassWriter::addAttribute("start_index", 0, &face_faces); 
    636                 SuperClassWriter::addAttribute("flag_values", -1, &face_faces); 
    637                 SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces); 
     591              SuperClassWriter::addAttribute("face_coordinates", face_x + " " + face_y, &domainName); 
     592              SuperClassWriter::addAttribute("face_node_connectivity", face_nodes, &domainName); 
     593              SuperClassWriter::addDimension(dimFace, domain->ni_glo); 
     594              SuperClassWriter::addDimension(dimVertex, domain->nvertex); 
     595              dim0.clear(); 
     596              dim0.push_back(dimFace); 
     597              SuperClassWriter::addVariable(face_x, NC_FLOAT, dim0); 
     598              SuperClassWriter::addAttribute("standard_name", StdString("longitude"), &face_x); 
     599              SuperClassWriter::addAttribute("long_name", StdString("Characteristic longitude of mesh faces."), &face_x); 
     600              SuperClassWriter::addAttribute("units", StdString("degrees_east"), &face_x); 
     601              SuperClassWriter::addVariable(face_y, NC_FLOAT, dim0); 
     602              SuperClassWriter::addAttribute("standard_name", StdString("latitude"), &face_y); 
     603              SuperClassWriter::addAttribute("long_name", StdString("Characteristic latitude of mesh faces."), &face_y); 
     604              SuperClassWriter::addAttribute("units", StdString("degrees_north"), &face_y); 
     605              dim0.clear(); 
     606              dim0.push_back(dimFace); 
     607              dim0.push_back(dimVertex); 
     608              SuperClassWriter::addVariable(face_nodes, NC_INT, dim0); 
     609              SuperClassWriter::addAttribute("cf_role", StdString("face_node_connectivity"), &face_nodes); 
     610              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its corner nodes."), &face_nodes); 
     611              SuperClassWriter::addAttribute("start_index", 0, &face_nodes); 
     612              dim0.clear(); 
     613              dim0.push_back(dimFace); 
     614              dim0.push_back(dimVertex); 
     615              SuperClassWriter::addVariable(face_edges, NC_INT, dim0); 
     616              SuperClassWriter::addAttribute("cf_role", StdString("face_edge_connectivity"), &face_edges); 
     617              SuperClassWriter::addAttribute("long_name", StdString("Maps every face to its edges."), &face_edges); 
     618              SuperClassWriter::addAttribute("start_index", 0, &face_edges); 
     619              dim0.clear(); 
     620              dim0.push_back(dimEdge); 
     621              dim0.push_back(dimTwo); 
     622              SuperClassWriter::addVariable(edge_faces, NC_INT, dim0); 
     623              SuperClassWriter::addAttribute("cf_role", StdString("edge_face connectivity"), &edge_faces); 
     624              SuperClassWriter::addAttribute("long_name", StdString("neighbor faces for edges"), &edge_faces); 
     625              SuperClassWriter::addAttribute("start_index", 0, &edge_faces); 
     626              SuperClassWriter::addAttribute("_FillValue", -999, &edge_faces); 
     627              SuperClassWriter::addAttribute("comment", StdString("missing neighbor faces are indicated using _FillValue"), &edge_faces); 
     628              dim0.clear(); 
     629              dim0.push_back(dimFace); 
     630              dim0.push_back(dimVertex); 
     631              SuperClassWriter::addVariable(face_faces, NC_INT, dim0); 
     632              SuperClassWriter::addAttribute("cf_role", StdString("face_face connectivity"), &face_faces); 
     633              SuperClassWriter::addAttribute("long_name", StdString("Indicates which other faces neighbor each face"), &face_faces); 
     634              SuperClassWriter::addAttribute("start_index", 0, &face_faces); 
     635              SuperClassWriter::addAttribute("flag_values", -1, &face_faces); 
     636              SuperClassWriter::addAttribute("flag_meanings", StdString("out_of_mesh"), &face_faces); 
    638637            } // domain->nvertex > 2 
    639638 
    640639            SuperClassWriter::definition_end(); 
    641640 
    642             std::vector<StdSize> start(1) ; 
    643             std::vector<StdSize> count(1) ; 
    644             if (domain->isEmpty()) 
    645              { 
    646                start[0]=0 ; 
    647                count[0]=0 ; 
    648              } 
    649              else 
    650              { 
    651                start[0]=domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
    652                count[0]=domain->zoom_ni_srv ; 
    653              } 
     641            std::vector<StdSize> startEdges(1) ; 
     642            std::vector<StdSize> countEdges(1) ; 
     643            std::vector<StdSize> startNodes(1) ; 
     644            std::vector<StdSize> countNodes(1) ; 
     645            std::vector<StdSize> startFaces(1) ; 
     646            std::vector<StdSize> countFaces(1) ; 
     647            std::vector<StdSize> startEdgeNodes(2) ; 
     648            std::vector<StdSize> countEdgeNodes(2) ; 
     649            std::vector<StdSize> startEdgeFaces(2) ; 
     650            std::vector<StdSize> countEdgeFaces(2) ; 
     651            std::vector<StdSize> startFaceConctv(2) ; 
     652            std::vector<StdSize> countFaceConctv(2) ; 
    654653 
    655654            if (!isWrittenDomain(domid)) 
    … …  
    657656              if (domain->nvertex == 1) 
    658657              { 
    659                 SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &start, &count); 
    660                 SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &start, &count); 
     658                if (domain->isEmpty()) 
     659                 { 
     660                   startNodes[0]=0 ; 
     661                   countNodes[0]=0 ; 
     662                 } 
     663                 else 
     664                 { 
     665                   startNodes[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     666                   countNodes[0] = domain->zoom_ni_srv ; 
     667                 } 
     668 
     669                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes); 
     670                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes); 
    661671              } 
    662672              else if (domain->nvertex == 2) 
    663                             { 
    664                 SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0); 
    665                 SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0); 
    666                 SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0); 
    667                 SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0); 
    668                 SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes); 
    669                             } 
     673              { 
     674                if (domain->isEmpty()) 
     675                 { 
     676                  startEdges[0]=0 ; 
     677                  countEdges[0]=0 ; 
     678                  startNodes[0]=0 ; 
     679                  countNodes[0]=0 ; 
     680                  startEdgeNodes[0]=0; 
     681                  startEdgeNodes[1]=0; 
     682                  countEdgeNodes[0]=0; 
     683                  countEdgeNodes[1]=0; 
     684 
     685                 } 
     686                 else 
     687                 { 
     688                   startEdges[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     689                   countEdges[0] = domain->zoom_ni_srv ; 
     690                   startNodes[0] = domain->mesh->node_start; 
     691                   countNodes[0] = domain->mesh->node_count; 
     692                   startEdgeNodes[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     693                   startEdgeNodes[1] = 0; 
     694                   countEdgeNodes[0] = domain->zoom_ni_srv; 
     695                   countEdgeNodes[1]= 2; 
     696                 } 
     697                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes); 
     698                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes); 
     699                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges); 
     700                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges); 
     701                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes); 
     702              } 
    670703              else 
    671704              { 
    672                 SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0); 
    673                 SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0); 
    674                 SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0); 
    675                 SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0); 
    676                 SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes); 
    677                 SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0); 
    678                 SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0); 
    679                 SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes); 
    680                 SuperClassWriter::writeData(domain->mesh->face_edges, face_edges); 
    681                 SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces); 
    682                 SuperClassWriter::writeData(domain->mesh->face_faces, face_faces); 
     705                if (domain->isEmpty()) 
     706                 { 
     707                   startFaces[0] = 0 ; 
     708                   countFaces[0] = 0 ; 
     709                   startNodes[0] = 0; 
     710                   countNodes[0] = 0; 
     711                   startEdges[0] = 0; 
     712                   countEdges[0] = 0; 
     713                   startEdgeFaces[0] = 0; 
     714                   startEdgeFaces[1] = 0; 
     715                   countEdgeFaces[0] = 0; 
     716                   countEdgeFaces[1] = 0; 
     717                   startFaceConctv[0] = 0; 
     718                   startFaceConctv[1] = 0; 
     719                   countFaceConctv[0] = 0; 
     720                   countFaceConctv[1] = 0; 
     721                 } 
     722                 else 
     723                 { 
     724                   startFaces[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     725                   countFaces[0] = domain->zoom_ni_srv ; 
     726                   startNodes[0] = domain->mesh->node_start; 
     727                   countNodes[0] = domain->mesh->node_count; 
     728                   startEdges[0] = domain->mesh->edge_start; 
     729                   countEdges[0] = domain->mesh->edge_count; 
     730                   startEdgeNodes[0] = domain->mesh->edge_start; 
     731                   startEdgeNodes[1] = 0; 
     732                   countEdgeNodes[0] = domain->mesh->edge_count; 
     733                   countEdgeNodes[1]= 2; 
     734                   startEdgeFaces[0] = domain->mesh->edge_start; 
     735                   startEdgeFaces[1]= 0; 
     736                   countEdgeFaces[0] = domain->mesh->edge_count; 
     737                   countEdgeFaces[1]= 2; 
     738                   startFaceConctv[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     739                   startFaceConctv[1] = 0; 
     740                   countFaceConctv[0] = domain->zoom_ni_srv; 
     741                   countFaceConctv[1] = domain->nvertex; 
     742                 } 
     743                SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0, &startNodes, &countNodes); 
     744                SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0, &startNodes, &countNodes); 
     745                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges); 
     746                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges); 
     747                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes); 
     748                SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0, &startFaces, &countFaces); 
     749                SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0, &startFaces, &countFaces); 
     750                SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes, isCollective, 0, &startFaceConctv, &countFaceConctv); 
     751                SuperClassWriter::writeData(domain->mesh->face_edges, face_edges, isCollective, 0, &startFaceConctv, &countFaceConctv); 
     752                SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces, isCollective, 0, &startEdgeFaces, &countEdgeFaces); 
     753                SuperClassWriter::writeData(domain->mesh->face_faces, face_faces, isCollective, 0, &startFaceConctv, &countFaceConctv); 
    683754              } 
    684755              setWrittenDomain(domid); 
    … …  
    686757            else 
    687758            { 
    688               if (domain->nvertex == 1) 
    689               { 
    690                 if ( (!domain->mesh->edgesAreWritten) && (!domain->mesh->facesAreWritten) ) 
    691                 { 
    692                   SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0); 
    693                   SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0); 
    694                 } 
    695               } 
    696759              if (domain->nvertex == 2) 
    697760              { 
    698                 if (!domain->mesh->facesAreWritten) 
    699                 { 
    700                   if (!domain->mesh->nodesAreWritten) 
    701                   { 
    702                     SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0); 
    703                     SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0); 
    704                   } 
    705                   SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0); 
    706                   SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0); 
    707                   SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes); 
    708                 } 
     761                startEdges[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     762                countEdges[0] = domain->zoom_ni_srv ; 
     763                startEdgeNodes[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     764                startEdgeNodes[1] = 0; 
     765                countEdgeNodes[0] = domain->zoom_ni_srv; 
     766                countEdgeNodes[1]= 2; 
     767                SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges); 
     768                SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges); 
     769                SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes); 
    709770              }   
    710771              if (domain->nvertex > 2) 
    711772              { 
     773 
    712774                if (!domain->mesh->edgesAreWritten) 
    713775                { 
    714                   if (!domain->mesh->nodesAreWritten) 
    715                   { 
    716                     SuperClassWriter::writeData(domain->mesh->node_lat, node_y, isCollective, 0); 
    717                     SuperClassWriter::writeData(domain->mesh->node_lon, node_x, isCollective, 0); 
    718                   } 
    719                   SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0); 
    720                   SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0); 
    721                   SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes); 
     776                  startEdges[0] = domain->mesh->edge_start; 
     777                  countEdges[0] = domain->mesh->edge_count; 
     778                  startEdgeNodes[0] = domain->mesh->edge_start; 
     779                  startEdgeNodes[1] = 0; 
     780                  countEdgeNodes[0] = domain->mesh->edge_count; 
     781                  countEdgeNodes[1]= 2; 
     782                  SuperClassWriter::writeData(domain->mesh->edge_lat, edge_y, isCollective, 0, &startEdges, &countEdges); 
     783                  SuperClassWriter::writeData(domain->mesh->edge_lon, edge_x, isCollective, 0, &startEdges, &countEdges); 
     784                  SuperClassWriter::writeData(domain->mesh->edge_nodes, edge_nodes, isCollective, 0, &startEdgeNodes, &countEdgeNodes); 
    722785                } 
    723                 SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0); 
    724                 SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0); 
    725                 SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes); 
    726                 SuperClassWriter::writeData(domain->mesh->face_edges, face_edges); 
    727                 SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces); 
    728                 SuperClassWriter::writeData(domain->mesh->face_faces, face_faces); 
     786                startFaces[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     787                countFaces[0] = domain->zoom_ni_srv; 
     788                startEdgeFaces[0] = domain->mesh->edge_start; 
     789                startEdgeFaces[1]= 0; 
     790                countEdgeFaces[0] = domain->mesh->edge_count; 
     791                countEdgeFaces[1]= 2; 
     792                startFaceConctv[0] = domain->zoom_ibegin_srv-domain->global_zoom_ibegin; 
     793                startFaceConctv[1] = 0; 
     794                countFaceConctv[0] = domain->zoom_ni_srv; 
     795                countFaceConctv[1]= domain->nvertex; 
     796                SuperClassWriter::writeData(domain->mesh->face_lat, face_y, isCollective, 0, &startFaces, &countFaces); 
     797                SuperClassWriter::writeData(domain->mesh->face_lon, face_x, isCollective, 0, &startFaces, &countFaces); 
     798                SuperClassWriter::writeData(domain->mesh->face_nodes, face_nodes, isCollective, 0, &startFaceConctv, &countFaceConctv); 
     799                SuperClassWriter::writeData(domain->mesh->face_edges, face_edges, isCollective, 0, &startFaceConctv, &countFaceConctv); 
     800                SuperClassWriter::writeData(domain->mesh->edge_faces, edge_faces, isCollective, 0, &startEdgeFaces, &countEdgeFaces); 
     801                SuperClassWriter::writeData(domain->mesh->face_faces, face_faces, isCollective, 0, &startFaceConctv, &countFaceConctv); 
    729802              } 
    730803            }// isWrittenDomain 
    … …  
    754827          msg.append(context->getId()); msg.append("\n"); 
    755828          msg.append(e.what()); 
    756           ERROR("CNc4DataOutput::writeUnstructuredDomain(CDomain* domain)", << msg); 
     829          ERROR("CNc4DataOutput::writeUnstructuredDomainUgrid(CDomain* domain)", << msg); 
    757830        } 
    758831 

  • XIOS/trunk/src/io/onetcdf4_decl.cpp

    r685 r924  
    1010  
    1111  macro(int, 1) 
     12  macro(int, 2) 
    1213  macro(double, 1) 
    1314  macro(double, 2) 

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

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

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

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

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

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

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

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