source: XIOS/dev/branch_openmp/src/node/mesh.cpp @ 1328

Last change on this file since 1328 was 1328, checked in by yushan, 6 years ago

dev_omp
  • Property svn:executable set to *
File size: 85.3 KB
Line 
1/*!
2  \file mesh.cpp
3  \author Olga Abramkina
4  \brief Definition of class CMesh.
5*/
6
7#include "mesh.hpp"
8using namespace ep_lib;
9
10namespace xios {
11
12/// ////////////////////// Définitions ////////////////////// ///
13
14  CMesh::CMesh(void) :  nbNodesGlo(0), nbEdgesGlo(0)
15            ,  node_start(0), node_count(0)
16            ,  edge_start(0), edge_count(0)
17            ,  nbFaces_(0), nbNodes_(0), nbEdges_(0)
18            ,  nodesAreWritten(false), edgesAreWritten(false), facesAreWritten(false)
19            ,  node_lon(), node_lat()
20            ,  edge_lon(), edge_lat(), edge_nodes()
21            ,  face_lon(), face_lat()
22            ,  face_nodes()
23            ,  pNodeGlobalIndex(NULL), pEdgeGlobalIndex(NULL)
24  {
25  }
26
27
28  CMesh::~CMesh(void)
29  {
30    if (pNodeGlobalIndex != NULL) delete pNodeGlobalIndex;
31    if (pEdgeGlobalIndex != NULL) delete pEdgeGlobalIndex;
32  }
33
34  //std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>();
35  //std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >();
36
37  std::map <StdString, CMesh> *CMesh::meshList_ptr = 0;
38  std::map <StdString, vector<int> > *CMesh::domainList_ptr = 0;
39
40///---------------------------------------------------------------
41/*!
42 * \fn bool CMesh::getMesh (StdString meshName)
43 * 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.
44 * \param [in] meshName  The name of a mesh ("name" attribute of a domain).
45 * \param [in] nvertex Number of verteces (1 for nodes, 2 for edges, 3 and up for faces).
46 */
47  CMesh* CMesh::getMesh (StdString meshName, int nvertex)
48  {
49    if(CMesh::domainList_ptr == NULL) CMesh::domainList_ptr = new std::map <StdString, vector<int> >();
50    if(CMesh::meshList_ptr == NULL)   CMesh::meshList_ptr   = new std::map <StdString, CMesh>();
51
52    //CMesh::domainList[meshName].push_back(nvertex);
53    CMesh::domainList_ptr->at(meshName).push_back(nvertex);
54
55    //if ( CMesh::meshList.begin() != CMesh::meshList.end() )
56    if ( CMesh::meshList_ptr->begin() != CMesh::meshList_ptr->end() )
57    {
58      //for (std::map<StdString, CMesh>::iterator it=CMesh::meshList.begin(); it!=CMesh::meshList.end(); ++it)
59      for (std::map<StdString, CMesh>::iterator it=CMesh::meshList_ptr->begin(); it!=CMesh::meshList_ptr->end(); ++it)
60      {
61        if (it->first == meshName)
62          //return &meshList[meshName];
63          return &meshList_ptr->at(meshName);
64        else
65        {
66          CMesh newMesh;
67          //CMesh::meshList.insert( make_pair(meshName, newMesh) );
68          CMesh::meshList_ptr->insert( make_pair(meshName, newMesh) );
69          //return &meshList[meshName];
70          return &meshList_ptr->at(meshName);
71        }
72      }
73    }
74    else
75    {
76      CMesh newMesh;
77      //CMesh::meshList.insert( make_pair(meshName, newMesh) );
78      CMesh::meshList_ptr->insert( make_pair(meshName, newMesh) );
79      //return &meshList[meshName];
80      return &meshList_ptr->at(meshName);
81    }
82  }
83
84///----------------------------------------------------------------
85  size_t hashPair(size_t first, size_t second)
86  {
87    HashXIOS<size_t> sizetHash;
88    size_t seed = sizetHash(first) + 0x9e3779b9 ;
89    seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
90    return seed ;
91  }
92
93///----------------------------------------------------------------
94  size_t hashPairOrdered(size_t first, size_t second)
95  {
96    size_t seed;
97    HashXIOS<size_t> sizetHash;
98    if (first < second)
99    {
100      seed = sizetHash(first) + 0x9e3779b9 ;
101      seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
102    }
103    else
104    {
105      seed = sizetHash(second) + 0x9e3779b9 ;
106      seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
107    }
108    return seed ;
109  }
110
111///----------------------------------------------------------------
112/*!
113 * \fn size_t generateNodeIndex(vector<size_t>& valList, int rank)
114 * Generates a node index.
115 * If the same node is generated by two processes, each process will have its own node index.
116 * \param [in] valList Vector storing four node hashes.
117 * \param [in] rank MPI process rank.
118 */
119  size_t generateNodeIndex(vector<size_t>& valList, int rank)
120  {
121    // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
122    vector<size_t> vec = valList;
123    sort (vec.begin(), vec.end());
124    size_t seed = rank ;
125    int it = 0;
126    for(; it != vec.size(); ++it)
127    {
128       seed = hashPair(seed, vec[it]);
129    }
130    return seed ;
131  }
132
133  ///----------------------------------------------------------------
134  /*!
135   * \fn size_t generateNodeIndex(vector<size_t>& valList)
136   * Generates a node index unique for all processes.
137   * \param [in] valList Vector storing four node hashes.
138   */
139    size_t generateNodeIndex(vector<size_t>& valList)
140    {
141      // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
142      vector<size_t> vec = valList;
143      sort (vec.begin(), vec.end());
144      size_t seed = vec[0] ;
145      int it = 1;
146      for(; it != vec.size(); ++it)
147      {
148         seed = hashPair(seed, vec[it]);
149      }
150      return seed ;
151    }
152
153
154///----------------------------------------------------------------
155/*!
156 * \fn size_t CMesh::nodeIndex (double lon, double lat)
157 * Returns its index if a node exists; otherwise adds the node and returns -1.
158 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
159 * \param [in] lon Node longitude in degrees.
160 * \param [in] lat Node latitude in degrees ranged from 0 to 360.
161 * \return node index if a node exists; -1 otherwise
162 */
163  size_t CMesh::nodeIndex (double lon, double lat)
164  {
165    double minBoundLon = 0. ;
166    double maxBoundLon = 360. ;
167    double minBoundLat = -90 ;
168    double maxBoundLat = 90 ;
169    double prec=1e-11 ;
170    double precLon=prec ;
171    double precLat=prec ;
172
173    size_t maxsize_t=numeric_limits<size_t>::max() ;
174    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
175    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
176
177    size_t iMinLon=0 ;
178    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
179    size_t iMinLat=0 ;
180    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
181
182    size_t hash0,hash1,hash2,hash3 ;
183    size_t lon0,lon1,lat0,lat1 ;
184
185    lon0=(lon-minBoundLon)/precLon ;
186    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
187    {
188      if (lon0==iMinLon) lon1=iMaxLon ;
189      else lon1=lon0-1 ;
190    }
191    else
192    {
193      if (lon0==iMaxLon) lon1=iMinLon ;
194      else lon1=lon0+1 ;
195    }
196
197    lat0=(lat-minBoundLat)/precLat ;
198    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
199    {
200      if (lat0==iMinLat) lat1=lat0 ;
201      else lat1=lat0-1 ;
202    }
203    else
204    {
205      if (lat0==iMaxLat) lat1=lat0 ;
206      else lat1=lat0+1 ;
207    }
208
209    hash0=hashPair(lon0,lat0) ;
210    hash1=hashPair(lon0,lat1) ;
211    hash2=hashPair(lon1,lat0) ;
212    hash3=hashPair(lon1,lat1) ;
213
214    boost::unordered_map<size_t, size_t>::iterator end = hashed_map_nodes.end() ;
215    size_t mapSize = hashed_map_nodes.size();
216    if (hashed_map_nodes.find(hash0)==end && hashed_map_nodes.find(hash1)==end && hashed_map_nodes.find(hash2)==end && hashed_map_nodes.find(hash3)==end)
217    {
218      hashed_map_nodes[hash0] = mapSize ;
219      hashed_map_nodes[hash1] = mapSize + 1;
220      hashed_map_nodes[hash2] = mapSize + 2;
221      hashed_map_nodes[hash3] = mapSize + 3;
222      return -1;
223    }
224    else
225      return ( (hashed_map_nodes[hash0]+1) / 4 );
226
227  } // nodeIndex()
228
229///----------------------------------------------------------------
230/*!
231 * \fn CArray<size_t,1>& CMesh::createHashes (const double longitude, const double latitude)
232 * Creates two hash values for each dimension, longitude and latitude.
233 * \param [in] longitude Node longitude in degrees.
234 * \param [in] latitude Node latitude in degrees ranged from 0 to 360.
235 */
236
237  vector<size_t> CMesh::createHashes (const double longitude, const double latitude)
238  {
239    double minBoundLon = 0. ;
240    double maxBoundLon = 360. ;
241    double minBoundLat = -90. ;
242    double maxBoundLat = 90. ;
243    double prec=1e-11 ;
244    double precLon=prec ;
245    double precLat=prec ;
246    double lon = longitude;
247    double lat = latitude;
248
249    if (lon > (360.- prec)) lon = 0.;
250
251    size_t maxsize_t=numeric_limits<size_t>::max() ;
252    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
253    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
254
255    size_t iMinLon=0 ;
256    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
257    size_t iMinLat=0 ;
258    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
259
260    vector<size_t> hash(4);
261    size_t lon0,lon1,lat0,lat1 ;
262
263    lon0=(lon-minBoundLon)/precLon ;
264    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
265    {
266      if (lon0==iMinLon) lon1=iMaxLon ;
267      else lon1=lon0-1 ;
268    }
269    else
270    {
271      if (lon0==iMaxLon) lon1=iMinLon ;
272      else lon1=lon0+1 ;
273    }
274
275    lat0=(lat-minBoundLat)/precLat ;
276    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
277    {
278      if (lat0==iMinLat) lat1=lat0 ;
279      else lat1=lat0-1 ;
280    }
281    else
282    {
283      if (lat0==iMaxLat) lat1=lat0 ;
284      else lat1=lat0+1 ;
285    }
286
287    hash[0] = hashPair(lon0,lat0) ;
288    hash[1] = hashPair(lon0,lat1) ;
289    hash[2] = hashPair(lon1,lat0) ;
290    hash[3] = hashPair(lon1,lat1) ;
291
292    return hash;
293
294  } // createHashes
295
296///----------------------------------------------------------------
297  std::pair<int,int> make_ordered_pair(int a, int b)
298  {
299    if ( a < b )
300      return std::pair<int,int>(a,b);
301    else
302      return std::pair<int,int>(b,a);
303  }
304
305///----------------------------------------------------------------
306/*!
307 * \fn void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
308            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
309 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
310 * \param [in] lonvalue  Array of longitudes.
311 * \param [in] latvalue  Array of latitudes.
312 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
313 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
314 */
315  void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
316            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
317  {
318    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
319
320    if (nvertex == 1)
321    {
322      nbNodes_ = lonvalue.numElements();
323      node_lon.resizeAndPreserve(nbNodes_);
324      node_lat.resizeAndPreserve(nbNodes_);
325      for (int nn = 0; nn < nbNodes_; ++nn)
326      {
327        if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end())
328        {
329          map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ;
330          node_lon(nn) = lonvalue(nn);
331          node_lat(nn) = latvalue(nn);
332        }
333      }
334    }
335    else if (nvertex == 2)
336    {
337      nbEdges_ = bounds_lon.shape()[1];
338
339      // Create nodes and edge_node connectivity
340      node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes
341      node_lat.resizeAndPreserve(nbEdges_*nvertex);
342      edge_nodes.resizeAndPreserve(nvertex, nbEdges_);
343
344      for (int ne = 0; ne < nbEdges_; ++ne)
345      {
346        for (int nv = 0; nv < nvertex; ++nv)
347        {
348          if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end())
349          {
350            map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ;
351            edge_nodes(nv,ne) = nbNodes_ ;
352            node_lon(nbNodes_) = bounds_lon(nv, ne);
353            node_lat(nbNodes_) = bounds_lat(nv, ne);
354            ++nbNodes_ ;
355          }
356          else
357            edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))];
358        }
359      }
360      node_lon.resizeAndPreserve(nbNodes_);
361      node_lat.resizeAndPreserve(nbNodes_);
362
363      // Create edges
364      edge_lon.resizeAndPreserve(nbEdges_);
365      edge_lat.resizeAndPreserve(nbEdges_);
366
367      for (int ne = 0; ne < nbEdges_; ++ne)
368      {
369        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
370        {
371          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
372          edge_lon(ne) = lonvalue(ne);
373          edge_lat(ne) = latvalue(ne);
374        }
375
376      }
377      edgesAreWritten = true;
378    }
379    else
380    {
381      nbFaces_ = bounds_lon.shape()[1];
382 
383      // Create nodes and face_node connectivity
384      node_lon.resizeAndPreserve(nbFaces_*nvertex);  // Max possible number of nodes
385      node_lat.resizeAndPreserve(nbFaces_*nvertex);
386      face_nodes.resize(nvertex, nbFaces_);
387 
388      for (int nf = 0; nf < nbFaces_; ++nf)
389      {
390        for (int nv = 0; nv < nvertex; ++nv)
391        {
392          if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end())
393          {
394            map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ;
395            face_nodes(nv,nf) = nbNodes_ ;
396            node_lon(nbNodes_) = bounds_lon(nv, nf);
397            node_lat(nbNodes_) = bounds_lat(nv ,nf);
398            ++nbNodes_ ;
399          }
400          else
401          {
402            face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
403          }
404        }
405      }
406      node_lon.resizeAndPreserve(nbNodes_);
407      node_lat.resizeAndPreserve(nbNodes_);
408 
409      // Create edges and edge_nodes connectivity
410      edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges
411      edge_lat.resizeAndPreserve(nbFaces_*nvertex);
412      edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex);
413      edge_faces.resize(2, nbFaces_*nvertex);
414      face_edges.resize(nvertex, nbFaces_);
415      face_faces.resize(nvertex, nbFaces_);
416
417      vector<int> countEdges(nbFaces_*nvertex);   // needed in case if edges have been already generated
418      vector<int> countFaces(nbFaces_);
419      countEdges.assign(nbFaces_*nvertex, 0);
420      countFaces.assign(nbFaces_, 0);
421      int edge;
422      for (int nf = 0; nf < nbFaces_; ++nf)
423      {
424        for (int nv1 = 0; nv1 < nvertex; ++nv1)
425        {
426          int nv = 0;
427          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
428          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
429          {
430            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ;
431            face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
432            edge_faces(0,nbEdges_) = nf;
433            edge_faces(1,nbEdges_) = -999;
434            face_faces(nv1,nf) = 999999;
435            edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf);
436            edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
437                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
438                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);
439            edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
440            ++nbEdges_;
441          }
442          else
443          {
444            edge = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
445            face_edges(nv1,nf) = edge;
446            if (edgesAreWritten)
447            {
448              edge_faces(countEdges[edge], edge) = nf;
449              if (countEdges[edge]==0)
450              {
451                face_faces(nv1,nf) = 999999;
452              }
453              else
454              {
455                int face1 = nf; // = edge_faces(1,edge)
456                int face2 = edge_faces(0,edge);
457                face_faces(countFaces[face1], face1) =  face2;
458                face_faces(countFaces[face2], face2) =  face1;
459                ++(countFaces[face1]);
460                ++(countFaces[face2]);
461              }
462            }
463            else
464            {
465              edge_faces(1,edge) = nf;
466              int face1 = nf; // = edge_faces(1,edge)
467              int face2 = edge_faces(0,edge);
468              face_faces(countFaces[face1], face1) =  face2;
469              face_faces(countFaces[face2], face2) =  face1;
470              ++(countFaces[face1]);
471              ++(countFaces[face2]);
472            }
473            ++(countEdges[edge]);
474          }
475        }
476      }
477      edge_nodes.resizeAndPreserve(2, nbEdges_);
478      edge_faces.resizeAndPreserve(2, nbEdges_);
479      edge_lon.resizeAndPreserve(nbEdges_);
480      edge_lat.resizeAndPreserve(nbEdges_);
481
482      // Create faces
483      face_lon.resize(nbFaces_);
484      face_lat.resize(nbFaces_);
485      face_lon = lonvalue;
486      face_lat = latvalue;
487      facesAreWritten = true;
488
489    } // nvertex > 2
490   
491  } // createMesh()
492
493///----------------------------------------------------------------
494/*!
495 * \fn void CMesh::createMeshEpsilon(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
496            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
497 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
498 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
499 * \param [in] comm
500 * \param [in] lonvalue  Array of longitudes.
501 * \param [in] latvalue  Array of latitudes.
502 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
503 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
504 */
505  void CMesh::createMeshEpsilon(const MPI_Comm& comm,
506                                const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
507                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
508  {
509
510    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
511    int mpiRank, mpiSize;
512    MPI_Comm_rank(comm, &mpiRank);
513    MPI_Comm_size(comm, &mpiSize);
514    double prec = 1e-11;  // used in calculations of edge_lon/lat
515
516    if (nvertex == 1)
517    {
518      nbNodes_ = lonvalue.numElements();
519      node_lon.resize(nbNodes_);
520      node_lat.resize(nbNodes_);
521      node_lon = lonvalue;
522      node_lat = latvalue;
523
524      // Global node indexes
525      vector<size_t> hashValues(4);
526      CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
527      for (size_t nn = 0; nn < nbNodes_; ++nn)
528      {
529        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn));
530        for (size_t nh = 0; nh < 4; ++nh)
531        {
532          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes_ + nn);
533        }
534      }
535      pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
536      nodesAreWritten = true;
537    }
538
539    else if (nvertex == 2)
540    {
541      nbEdges_ = bounds_lon.shape()[1];
542      edge_lon.resize(nbEdges_);
543      edge_lat.resize(nbEdges_);
544      edge_lon = lonvalue;
545      edge_lat = latvalue;
546      edge_nodes.resize(nvertex, nbEdges_);
547
548      // For determining the global edge index
549      size_t nbEdgesOnProc = nbEdges_;
550      size_t nbEdgesAccum;
551      MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
552      nbEdgesAccum -= nbEdges_;
553
554      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo;
555      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
556
557      // Case (1): node indexes have been generated by domain "nodes"
558      if (nodesAreWritten)
559      {
560        vector<size_t> hashValues(4);
561        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4);
562        for (int ne = 0; ne < nbEdges_; ++ne)      // size_t doesn't work with CArray<double, 2>
563        {
564          for (int nv = 0; nv < nvertex; ++nv)
565          {
566            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
567            for (int nh = 0; nh < 4; ++nh)
568            {
569              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
570            }
571          }
572        }
573
574        // Recuperating the node global indexing and writing edge_nodes
575        // Creating map edgeHash2IdxGlo <hash, idxGlo>
576        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
577        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
578        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
579        size_t nodeIdxGlo1, nodeIdxGlo2;
580        for (int ne = 0; ne < nbEdges_; ++ne)
581        {
582          for (int nv = 0; nv < nvertex; ++nv)
583          {
584            int nh = 0;
585            it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
586            // The loop below is needed in case if a hash generated by domain "edges" differs
587            // from that generated by domain "nodes" because of possible precision issues
588            while (it == nodeHash2IdxGlo.end())
589            {
590              ++nh;
591              it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
592            }
593            edge_nodes(nv,ne) = it->second[0];
594            if (nv ==0)
595              nodeIdxGlo1 = it->second[0];
596            else
597              nodeIdxGlo2 = it->second[0];
598          }
599          size_t edgeIdxGlo = nbEdgesAccum + ne;
600          edgeHash2IdxGlo[ hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2) ].push_back(edgeIdxGlo);
601        }
602      } // nodesAreWritten
603
604
605      // Case (2): node indexes have not been generated previously
606      else
607      {
608        // (2.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
609        vector<size_t> hashValues(4);
610        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
611        CArray<size_t,1> nodeHashList(nbEdges_*nvertex*4);
612        int nbHash = 0;
613        for (int ne = 0; ne < nbEdges_; ++ne)
614        {
615          for (int nv = 0; nv < nvertex; ++nv)
616          {
617            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
618            for (int nh = 0; nh < 4; ++nh)
619            {
620              if (nodeHash2Idx[hashValues[nh]].size() == 0)
621              {
622                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues));
623                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
624                nodeHashList(nbHash) = hashValues[nh];
625                ++nbHash;
626              }
627            }
628          }
629        }
630        nodeHashList.resizeAndPreserve(nbHash);
631
632        // (2.2) Generating global node indexes
633        // The ownership criterion: priority of the process of smaller index
634        // Maps generated in this step are:
635        // Maps generated in this step are:
636         // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
637         // nodeIdx2Idx = <idx, <rankOwner, idx>>
638
639         CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
640         dhtNodeHash.computeIndexInfoMapping(nodeHashList);
641         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
642
643
644         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
645         CArray<size_t,1> nodeIdxList(nbEdges_*nvertex*4);
646         size_t nIdx = 0;
647
648         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
649         {
650           size_t rankMin = (it->second)[1];
651           size_t idx = (it->second)[0];
652           for (int i = 2; i < (it->second).size();)
653           {
654             if ( (it->second)[i+1] < rankMin)
655             {
656               idx = (it->second)[i];
657               rankMin = (it->second)[i+1];
658               (it->second)[i+1] = (it->second)[i-1];
659             }
660             i += 2;
661           }
662           if (nodeIdx2Idx.count(idx) == 0)
663           {
664             if (mpiRank == rankMin)
665             {
666               nodeIdx2Idx[idx].push_back(rankMin);
667               nodeIdx2Idx[idx].push_back(idx);
668             }
669             nodeIdxList(nIdx) = idx;
670             ++nIdx;
671           }
672         }
673         nodeIdxList.resizeAndPreserve(nIdx);
674
675         // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
676         // Solution: global node indexing by hand.
677         // Maps modified in this step:
678         // nodeIdx2Idx = <idx, idxGlo>
679         int nodeCount = nodeIdx2Idx.size();
680         int nodeStart, nbNodes;
681         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
682         int nNodes = nodeStart;
683         MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
684         nbNodesGlo = nNodes;
685
686         nodeStart -= nodeCount;
687         node_start = nodeStart;
688         node_count = nodeCount;
689         CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
690         size_t count = 0;
691
692         for (int ne = 0; ne < nbEdges_; ++ne)
693         {
694           for (int nv = 0; nv < nvertex; ++nv)
695           {
696             vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
697             size_t nodeIdx = generateNodeIndex(hashValues);
698             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
699             if (it != nodeIdx2Idx.end())
700             {
701               if (dummyMap.count(nodeIdx) == 0)
702               {
703                 dummyMap[nodeIdx].push_back(nodeIdx);
704                 (it->second)[1] = node_start + count;
705                 ++count;
706               }
707             }
708           }
709         }
710
711         CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
712         dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
713         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
714
715        // (2.3) Saving variables: node_lon, node_lat, edge_nodes
716        // Creating map nodeHash2IdxGlo <hash, idxGlo>
717        // Creating map edgeHash2IdxGlo <hash, idxGlo>
718//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
719//        node_count = dhtNodeIdxGlo.getIndexCount();
720//        node_start = dhtNodeIdxGlo.getIndexStart();
721        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
722        node_lon.resize(node_count);
723        node_lat.resize(node_count);
724        vector <size_t> edgeNodes;
725        size_t idxGlo = 0;
726
727        for (int ne = 0; ne < nbEdges_; ++ne)
728        {
729          for (int nv = 0; nv < nvertex; ++nv)
730          {
731            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
732            size_t myIdx = generateNodeIndex(hashValues);
733            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxGlo.find(myIdx);
734            idxGlo = (itIdx->second)[1];
735
736            if (mpiRank == (itIdx->second)[0])
737            {
738//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
739              node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
740              node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
741            }
742            edge_nodes(nv,ne) = idxGlo;
743            for (int nh = 0; nh < 4; ++nh)
744              nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo);
745            edgeNodes.push_back(idxGlo);
746          }
747          if (edgeNodes[0] != edgeNodes[1])
748          {
749            size_t edgeIdxGlo = nbEdgesAccum + ne;
750            edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(edgeIdxGlo);
751          }
752          edgeNodes.clear();
753        }
754        pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
755      } // !nodesAreWritten
756
757      pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm);
758      edgesAreWritten = true;
759    } //nvertex = 2
760
761    else
762    {
763      nbFaces_ = bounds_lon.shape()[1];
764      face_lon.resize(nbFaces_);
765      face_lat.resize(nbFaces_);
766      face_lon = lonvalue;
767      face_lat = latvalue;
768      face_nodes.resize(nvertex, nbFaces_);
769      face_edges.resize(nvertex, nbFaces_);
770
771      // For determining the global face index
772      size_t nbFacesOnProc = nbFaces_;
773      size_t nbFacesAccum;
774      MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
775      nbFacesAccum -= nbFaces_;
776
777      // Case (1): edges have been previously generated
778      if (edgesAreWritten)
779      {
780        // (1.1) Recuperating node global indexing and saving face_nodes
781        vector<size_t> hashValues(4);
782        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
783        for (int nf = 0; nf < nbFaces_; ++nf)
784        {
785          for (int nv = 0; nv < nvertex; ++nv)
786            {
787              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
788              for (int nh = 0; nh < 4; ++nh)
789                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
790            }
791        }
792        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
793        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
794        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
795        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
796        size_t nEdge = 0;
797        for (int nf = 0; nf < nbFaces_; ++nf)
798        {
799          for (int nv1 = 0; nv1 < nvertex; ++nv1)
800          {
801            int nh1 = 0;
802            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
803            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
804            while (it1 == nodeHash2IdxGlo.end())
805            {
806              ++nh1;
807              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
808            }
809            int nh2 = 0;
810            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
811            while (it2 == nodeHash2IdxGlo.end())
812            {
813              ++nh2;
814              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
815            }
816            face_nodes(nv1,nf) = it1->second[0];
817            if (it1->second[0] != it2->second[0])
818            {
819              edgeHashList(nEdge) = hashPairOrdered(it1->second[0], it2->second[0]);
820              ++nEdge;
821            }
822          }
823        }
824        edgeHashList.resizeAndPreserve(nEdge);
825
826        // (1.2) Recuperating edge global indexing and saving face_edges
827        pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList);
828        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap();
829        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash;
830        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
831        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
832        CArray<size_t,1> edgeIdxList(nbFaces_*nvertex);
833        size_t iIdx = 0;
834
835        for (int nf = 0; nf < nbFaces_; ++nf)
836        {
837          for (int nv1 = 0; nv1 < nvertex; ++nv1)
838          {
839            int nh1 = 0;
840            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
841            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
842            while (it1 == nodeHash2IdxGlo.end())
843            {
844              ++nh1;
845              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
846            }
847            int nh2 = 0;
848            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
849            while (it2 == nodeHash2IdxGlo.end())
850            {
851              ++nh2;
852              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
853            }
854            if (it1->second[0] != it2->second[0])
855            {
856              size_t faceIdxGlo = nbFacesAccum + nf;
857              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
858              itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
859              size_t edgeIdxGlo = itEdgeHash->second[0];
860              face_edges(nv1,nf) = edgeIdxGlo;
861              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
862              {
863                edgeIdxList(iIdx) = edgeIdxGlo;
864                ++iIdx;
865              }
866              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
867              edgeHash2Rank[edgeHash].push_back(mpiRank);
868              edgeHash2Rank[edgeHash].push_back(itEdgeHash->second[0]);
869            }
870            else
871            {
872              face_edges(nv1,nf) = 999999;
873            }
874          }
875        }
876        edgeIdxList.resizeAndPreserve(iIdx);
877
878        // (1.3) Saving remaining variables edge_faces and face_faces
879
880        // Establishing edge ownership
881        // The ownership criterion: priority of the process with smaller rank
882        CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm);
883        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
884        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
885
886        // edgeHash2Info = <edgeHash, < rank1, idxGlo, rank2, idxGlo>>
887        int edgeCount = 0;
888        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
889        {
890          vector <size_t> edgeInfo = it->second;
891          if (edgeInfo[0] == mpiRank)
892          {
893            ++edgeCount;
894          }
895        }
896
897        int edgeStart, nbEdges;
898        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
899        int nEdges = edgeStart;
900        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
901        nbEdgesGlo = nEdges;
902
903        // edges to be splitted equally between procs
904        if ( (nbEdgesGlo % mpiSize) == 0)
905        {
906          edge_count = nbEdgesGlo/mpiSize;
907          edge_start = mpiRank*edge_count;
908        }
909        else
910        {
911          if (mpiRank == (mpiSize - 1) )
912          {
913            edge_count = nbEdgesGlo/mpiSize;
914            edge_start = mpiRank*(nbEdgesGlo/mpiSize + 1);
915          }
916          else
917          {
918            edge_count = nbEdgesGlo/mpiSize + 1;
919            edge_start = mpiRank*edge_count;
920          }
921        }
922        CArray<size_t,1> edgeIdxGloList(edge_count);
923        for (int i = 0; i < edge_count; ++i)
924        {
925          edgeIdxGloList(i) = i + edge_start;
926        }
927
928        CClientClientDHTSizet dhtEdgeIdxGlo2Face (edgeIdxGlo2Face, comm);
929        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
930        dhtEdgeIdxGlo2Face.computeIndexInfoMapping(edgeIdxGloList);
931        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdgeIdxGlo2Face.getInfoIndexMap();
932        dhtEdge2Face.computeIndexInfoMapping(edgeIdxList);
933        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2FaceIdx = dhtEdge2Face.getInfoIndexMap();
934
935
936        edge_faces.resize(2, edge_count);
937        for (int i = 0; i < edge_count; ++i)
938        {
939          CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdxGlo2FaceIdx.find(i + edge_start);
940          int indexGlo = it->first;
941          vector<size_t> faces = it->second;
942          int face1 = faces[0];
943          edge_faces(0, indexGlo - edge_start) = face1;
944          if (faces.size() == 2)
945          {
946            int face2 = faces[1];
947            edge_faces(1, indexGlo - edge_start) = face2;
948          }
949          else
950          {
951            edge_faces(1, indexGlo - edge_start) = -999;
952          }
953        }
954
955        size_t tmp;
956        vector<size_t> tmpVec;
957        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2FaceIdx.begin(); it != edgeIdx2FaceIdx.end(); it++)
958        {
959          tmp = it->first;
960          tmpVec = it->second;
961          tmp++;
962        }
963
964        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2, itIndex;
965        face_faces.resize(nvertex, nbFaces_);
966        for (int nf = 0; nf < nbFaces_; ++nf)
967        {
968          for (int nv1 = 0; nv1 < nvertex; ++nv1)
969          {
970            int nh1 = 0;
971            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
972            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
973            while (it1 == nodeHash2IdxGlo.end())
974            {
975              ++nh1;
976              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
977            }
978            int nh2 = 0;
979            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
980            while (it2 == nodeHash2IdxGlo.end())
981            {
982              ++nh2;
983              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
984            }
985
986            if (it1->second[0] != it2->second[0])
987            {
988              size_t faceIdxGlo = nbFacesAccum + nf;
989              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
990              itEdgeHash = edgeHash2Info.find(edgeHash);
991              int edgeIdxGlo = (itEdgeHash->second)[1];
992
993              if ( (itEdgeHash->second)[0] == mpiRank)
994              {
995                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
996                int face1 = itFace1->second[0];
997                if (itFace1->second.size() == 1)
998                {
999                  face_faces(nv1, nf) = 999999;
1000                }
1001                else
1002                {
1003                  int face2 = itFace1->second[1];
1004                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1005                }
1006              } // edge owner
1007              else
1008              {
1009                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
1010                int face1 = itFace1->second[0];
1011                int face2 = itFace1->second[1];
1012                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1013              } // not an edge owner
1014            } // node1 != node2
1015            else
1016            {
1017              face_faces(nv1, nf) = 999999;
1018            }
1019          }
1020        }
1021      } // edgesAreWritten
1022
1023      // Case (2): nodes have been previously generated
1024      else if (nodesAreWritten)
1025      {
1026        // (2.1) Generating nodeHashList
1027        vector<size_t> hashValues(4);
1028        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
1029        for (int nf = 0; nf < nbFaces_; ++nf)
1030        {
1031          for (int nv = 0; nv < nvertex; ++nv)
1032            {
1033              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1034              for (int nh = 0; nh < 4; ++nh)
1035                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
1036            }
1037        }
1038
1039        // (2.2) Recuperating node global indexing and saving face_nodes
1040        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
1041        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
1042        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
1043        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
1044        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
1045        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
1046        int nEdgeHash = 0;
1047        for (int nf = 0; nf < nbFaces_; ++nf)
1048        {
1049          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1050          {
1051            int nh1 = 0;
1052            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1053            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1054            while (it1 == nodeHash2IdxGlo.end())
1055            {
1056              ++nh1;
1057              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1058            }
1059            int nh2 = 0;
1060            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1061            while (it2 == nodeHash2IdxGlo.end())
1062            {
1063              ++nh2;
1064              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1065            }
1066            face_nodes(nv1,nf) = it1->second[0];
1067            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
1068            if (edgeHash2Idx.count(edgeHash) == 0)
1069            {
1070              edgeHash2Idx[edgeHash].push_back(edgeHash);
1071              edgeHash2Idx[edgeHash].push_back(mpiRank);
1072              edgeHashList(nEdgeHash) = edgeHash;
1073              ++nEdgeHash;
1074            }
1075          }
1076        }
1077        edgeHashList.resizeAndPreserve(nEdgeHash);
1078
1079        // (2.3) Generating global edge indexes
1080        // The ownership criterion: priority of the process with smaller rank
1081        // Maps generated in this step are:
1082        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
1083        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
1084
1085        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
1086        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
1087        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
1088        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
1089
1090        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
1091
1092        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
1093        {
1094          size_t rankMin = (it->second)[1];
1095          size_t idx = (it->second)[0];
1096
1097          for (int i = 2; i < (it->second).size();)
1098          {
1099            if ((it->second)[i+1] < rankMin)
1100            {
1101              rankMin = (it->second)[i+1];
1102              idx = (it->second)[i];
1103              (it->second)[i+1] = (it->second)[i-1];
1104            }
1105            i += 2;
1106          }
1107          if (edgeIdx2Idx.count(idx) == 0)
1108          {
1109            if (mpiRank == rankMin)
1110            {
1111              edgeIdx2Idx[idx].push_back(rankMin);
1112              edgeIdx2Idx[idx].push_back(idx);
1113            }
1114          }
1115        }
1116
1117        int edgeCount = edgeIdx2Idx.size();
1118        int edgeStart, nbEdges;
1119        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
1120        int nEdges = edgeStart;
1121        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
1122        nbEdgesGlo = nEdges;
1123
1124        edgeStart -= edgeCount;
1125        edge_start = edgeStart;
1126        edge_count = edgeCount;
1127        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
1128        int count = 0;
1129
1130        for (int nf = 0; nf < nbFaces_; ++nf)
1131        {
1132          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1133          {
1134            // Getting global indexes of edge's nodes
1135            int nh1 = 0;
1136            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1137            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1138            while (it1 == nodeHash2IdxGlo.end())
1139            {
1140              ++nh1;
1141              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1142            }
1143            int nh2 = 0;
1144            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1145            while (it2 == nodeHash2IdxGlo.end())
1146            {
1147              ++nh2;
1148              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1149            }
1150            size_t nodeIdxGlo1 = it1->second[0];
1151            size_t nodeIdxGlo2 = it2->second[0];
1152
1153            if (nodeIdxGlo1 != nodeIdxGlo2)
1154            {
1155              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1156              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
1157              if (it != edgeIdx2Idx.end())
1158              {
1159                if (dummyEdgeMap.count(edgeIdx) == 0)
1160                {
1161                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
1162                  (it->second)[1] = edge_start + count;
1163                  ++count;
1164                }
1165              }
1166            }
1167          }
1168        }
1169
1170        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
1171        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
1172        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
1173
1174        // (2.4) Saving variables: edge_lon, edge_lat, face_edges
1175        edge_lon.resize(edge_count);
1176        edge_lat.resize(edge_count);
1177        edge_nodes.resize(2, edge_count);
1178        face_edges.resize(nvertex, nbFaces_);
1179
1180        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
1181        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
1182        size_t iIdx = 0;
1183
1184        for (int nf = 0; nf < nbFaces_; ++nf)
1185        {
1186          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1187          {
1188            // Getting global indexes of edge's nodes
1189            int nh1 = 0;
1190            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1191            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1192            while (it1 == nodeHash2IdxGlo.end())
1193            {
1194              ++nh1;
1195              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1196            }
1197            int nh2 = 0;
1198            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1199            while (it2 == nodeHash2IdxGlo.end())
1200            {
1201              ++nh2;
1202              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1203            }
1204            // Getting edge global index
1205            size_t nodeIdxGlo1 = it1->second[0];
1206            size_t nodeIdxGlo2 = it2->second[0];
1207            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1208            if (nodeIdxGlo1 != nodeIdxGlo2)
1209            {
1210              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
1211              int edgeIdxGlo = (itIdx->second)[1];
1212              size_t faceIdxGlo = nbFacesAccum + nf;
1213
1214              if (mpiRank == (itIdx->second)[0])
1215              {
1216                double edgeLon;
1217                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
1218                if (diffLon < (180.- prec))
1219                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
1220                else if (diffLon > (180.+ prec))
1221                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
1222                else
1223                  edgeLon = 0.;
1224                edge_lon(edgeIdxGlo - edge_start) = edgeLon;
1225                edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
1226                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
1227                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
1228              }
1229              face_edges(nv1,nf) = edgeIdxGlo;
1230              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
1231              {
1232                edgeIdxGloList(iIdx) = edgeIdxGlo;
1233                ++iIdx;
1234              }
1235              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
1236            } // nodeIdxGlo1 != nodeIdxGlo2
1237            else
1238            {
1239              face_edges(nv1,nf) = 999999;
1240            }
1241          }
1242        }
1243        edgeIdxGloList.resizeAndPreserve(iIdx);
1244
1245        // (2.5) Saving remaining variables edge_faces and face_faces
1246        edge_faces.resize(2, edge_count);
1247        face_faces.resize(nvertex, nbFaces_);
1248
1249        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
1250        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
1251        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
1252        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
1253        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx;
1254
1255        for (int nf = 0; nf < nbFaces_; ++nf)
1256        {
1257          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1258          {
1259            // Getting global indexes of edge's nodes
1260            int nh1 = 0;
1261            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1262            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1263            while (it1 == nodeHash2IdxGlo.end())
1264            {
1265              ++nh1;
1266              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1267            }
1268            int nh2 = 0;
1269            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1270            while (it2 == nodeHash2IdxGlo.end())
1271            {
1272              ++nh2;
1273              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1274            }
1275            size_t nodeIdxGlo1 = it1->second[0];
1276            size_t nodeIdxGlo2 = it2->second[0];
1277
1278            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1279            itIdx = edgeIdx2IdxGlo.find(myIdx);
1280            size_t faceIdxGlo = nbFacesAccum + nf;
1281            int edgeIdxGlo = (itIdx->second)[1];
1282
1283            if (mpiRank == (itIdx->second)[0])
1284            {
1285              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1286              int face1 = it1->second[0];
1287              if (it1->second.size() == 1)
1288              {
1289                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1290                edge_faces(1, edgeIdxGlo - edge_start) = -999;
1291                face_faces(nv1, nf) = 999999;
1292              }
1293              else
1294              {
1295                int face2 = it1->second[1];
1296                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1297                edge_faces(1, edgeIdxGlo - edge_start) = face2;
1298                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1299              }
1300            }
1301            else
1302            {
1303              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1304              int face1 = it1->second[0];
1305              int face2 = it1->second[1];
1306              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1307            }
1308          }
1309        }
1310      } // nodesAreWritten
1311
1312      // Case (3): Neither nodes nor edges have been previously generated
1313      else
1314      {
1315        // (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1316        vector<size_t> hashValues(4);
1317        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1318        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
1319        size_t iHash = 0;
1320        for (int nf = 0; nf < nbFaces_; ++nf)
1321        {
1322          for (int nv = 0; nv < nvertex; ++nv)
1323          {
1324            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1325//            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1326            size_t nodeIndex = generateNodeIndex(hashValues);
1327            for (int nh = 0; nh < 4; ++nh)
1328            {
1329              if (nodeHash2Idx.count(hashValues[nh])==0)
1330              {
1331                nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1332                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1333                nodeHashList(iHash) = hashValues[nh];
1334                ++iHash;
1335              }
1336            }
1337          }
1338        }
1339        nodeHashList.resizeAndPreserve(iHash);
1340
1341        // (3.2) Generating global node indexes
1342        // The ownership criterion: priority of the process with smaller rank.
1343        // With any other criterion it is not possible to have consistent node indexing for different number of procs.
1344        // Maps generated in this step are:
1345        // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
1346        // nodeIdx2Idx = <idx, <rankOwner, idx>>
1347
1348        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1349        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1350        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1351
1352        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
1353        CArray<size_t,1> nodeIdxList(nbFaces_*nvertex*4);
1354        size_t nIdx = 0;
1355
1356        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1357        {
1358          size_t rankMin = (it->second)[1];
1359          size_t idx = (it->second)[0];
1360          for (int i = 2; i < (it->second).size();)
1361          {
1362            if ( (it->second)[i+1] < rankMin)
1363            {
1364              idx = (it->second)[i];
1365              rankMin = (it->second)[i+1];
1366              (it->second)[i+1] = (it->second)[i-1];
1367            }
1368            i += 2;
1369          }
1370          if (nodeIdx2Idx.count(idx) == 0)
1371          {
1372            if (mpiRank == rankMin)
1373            {
1374              nodeIdx2Idx[idx].push_back(rankMin);
1375              nodeIdx2Idx[idx].push_back(idx);
1376            }
1377            nodeIdxList(nIdx) = idx;
1378            ++nIdx;
1379          }
1380        }
1381
1382//        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2Idx, comm);
1383        // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
1384        // Solution: global node indexing by hand.
1385        // Maps modified in this step:
1386        // nodeIdx2Idx = <idx, idxGlo>
1387        int nodeCount = nodeIdx2Idx.size();
1388        int nodeStart, nbNodes;
1389        MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
1390        int nNodes = nodeStart;
1391        MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
1392        nbNodesGlo = nNodes;
1393
1394        nodeStart -= nodeCount;
1395        node_start = nodeStart;
1396        node_count = nodeCount;
1397        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
1398        size_t count = 0;
1399
1400        for (int nf = 0; nf < nbFaces_; ++nf)
1401        {
1402          for (int nv = 0; nv < nvertex; ++nv)
1403          {
1404            vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1405            size_t nodeIdx = generateNodeIndex(hashValues);
1406            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
1407            if (it != nodeIdx2Idx.end())
1408            {
1409              if (dummyMap.count(nodeIdx) == 0)
1410              {
1411                dummyMap[nodeIdx].push_back(nodeIdx);
1412                (it->second)[1] = node_start + count;
1413                ++count;
1414              }
1415            }
1416          }
1417        }
1418        nodeIdxList.resizeAndPreserve(nIdx);
1419        CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
1420        dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
1421        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
1422
1423        // (3.3) Saving node data: node_lon, node_lat, and face_nodes
1424        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
1425//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
1426//        node_count = dhtNodeIdxGlo.getIndexCount();
1427//        node_start = dhtNodeIdxGlo.getIndexStart();
1428        node_lon.resize(node_count);
1429        node_lat.resize(node_count);
1430        size_t nodeIdxGlo1 = 0;
1431        size_t nodeIdxGlo2 = 0;
1432        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
1433        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
1434        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
1435        size_t nEdgeHash = 0;
1436
1437        for (int nf = 0; nf < nbFaces_; ++nf)
1438        {
1439          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1440          {
1441            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1442            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1443            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1444            size_t nodeIdx1 = generateNodeIndex(hashValues1);
1445            size_t nodeIdx2 = generateNodeIndex(hashValues2);
1446            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
1447            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
1448            size_t ownerRank = (itNodeIdx1->second)[0];
1449            nodeIdxGlo1 = (itNodeIdx1->second)[1];
1450            nodeIdxGlo2 = (itNodeIdx2->second)[1];
1451
1452            if (mpiRank == ownerRank)
1453            {
1454              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
1455              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
1456            }
1457            if (nodeIdxGlo1 != nodeIdxGlo2)
1458            {
1459              size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1460              edgeHash2Idx[edgeHash].push_back(edgeHash);
1461              edgeHash2Idx[edgeHash].push_back(mpiRank);
1462              edgeHashList(nEdgeHash) = edgeHash;
1463              ++nEdgeHash;
1464            }
1465            face_nodes(nv1,nf) = nodeIdxGlo1;
1466          }
1467        }
1468        edgeHashList.resizeAndPreserve(nEdgeHash);
1469
1470        // (3.4) Generating global edge indexes
1471        // Maps generated in this step are:
1472        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
1473        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
1474
1475        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
1476        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
1477        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
1478        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
1479
1480        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
1481
1482        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
1483        {
1484          size_t rankMin = (it->second)[1];
1485          size_t idx = (it->second)[0];
1486
1487          for (int i = 2; i < (it->second).size();)
1488          {
1489            if ((it->second)[i+1] < rankMin)
1490            {
1491              rankMin = (it->second)[i+1];
1492              idx = (it->second)[i];
1493              (it->second)[i+1] = (it->second)[i-1];
1494            }
1495            i += 2;
1496          }
1497          if (edgeIdx2Idx.count(idx) == 0)
1498          {
1499            if (mpiRank == rankMin)
1500            {
1501              edgeIdx2Idx[idx].push_back(rankMin);
1502              edgeIdx2Idx[idx].push_back(idx);
1503            }
1504          }
1505        }
1506
1507        int edgeCount = edgeIdx2Idx.size();
1508        int edgeStart, nbEdges;
1509        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
1510        int nEdges = edgeStart;
1511        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
1512        nbEdgesGlo = nEdges;
1513
1514        edgeStart -= edgeCount;
1515        edge_start = edgeStart;
1516        edge_count = edgeCount;
1517        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
1518        count = 0;
1519
1520        for (int nf = 0; nf < nbFaces_; ++nf)
1521        {
1522          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1523          {
1524            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1525            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1526            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1527            size_t nodeIdx1 = generateNodeIndex(hashValues1);
1528            size_t nodeIdx2 = generateNodeIndex(hashValues2);
1529            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
1530            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
1531            nodeIdxGlo1 = (itNodeIdx1->second)[1];
1532            nodeIdxGlo2 = (itNodeIdx2->second)[1];
1533
1534            if (nodeIdxGlo1 != nodeIdxGlo2)
1535            {
1536              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1537              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
1538              if (it != edgeIdx2Idx.end())
1539              {
1540                if (dummyEdgeMap.count(edgeIdx) == 0)
1541                {
1542                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
1543                  (it->second)[1] = edge_start + count;
1544                  ++count;
1545                }
1546              }
1547            }
1548          }
1549        }
1550        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
1551        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
1552        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
1553
1554        // (3.5) Saving variables: edge_lon, edge_lat, face_edges
1555        // Creating map edgeIdxGlo2Face <idxGlo, face>
1556//        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
1557//        edge_count = dhtEdgeIdxGlo.getIndexCount();
1558//        edge_start = dhtEdgeIdxGlo.getIndexStart();
1559
1560        edge_lon.resize(edge_count);
1561        edge_lat.resize(edge_count);
1562        edge_nodes.resize(2, edge_count);
1563        face_edges.resize(nvertex, nbFaces_);
1564
1565        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
1566        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
1567        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
1568        size_t nEdge = 0;
1569
1570        for (int nf = 0; nf < nbFaces_; ++nf)
1571        {
1572          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1573          {
1574            // Getting global indexes of edge's nodes
1575            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1576            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1577            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1578
1579            size_t nodeIdx1 = generateNodeIndex(hashValues1);
1580            size_t nodeIdx2 = generateNodeIndex(hashValues2);
1581            it1 = nodeIdx2IdxGlo.find(nodeIdx1);
1582            it2 = nodeIdx2IdxGlo.find(nodeIdx2);
1583            size_t nodeIdxGlo1 = (it1->second)[1];
1584            size_t nodeIdxGlo2 = (it2->second)[1];
1585
1586            if (nodeIdxGlo1 != nodeIdxGlo2)
1587            {
1588              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1589              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
1590              int edgeIdxGlo = (itIdx->second)[1];
1591              size_t faceIdxGlo = nbFacesAccum + nf;
1592
1593              if (mpiRank == (itIdx->second)[0])
1594              {
1595                double edgeLon;
1596                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
1597                if (diffLon < (180.- prec))
1598                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
1599                else if (diffLon > (180.+ prec))
1600                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
1601                else
1602                  edgeLon = 0.;
1603                edge_lon(edgeIdxGlo - edge_start) = edgeLon;
1604                edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
1605                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
1606                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
1607              }
1608              face_edges(nv1,nf) = edgeIdxGlo;
1609              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
1610              {
1611                edgeIdxGloList(nEdge) = edgeIdxGlo;
1612                ++nEdge;
1613              }
1614              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
1615            } // nodeIdxGlo1 != nodeIdxGlo2
1616            else
1617            {
1618              face_edges(nv1,nf) = 999999;
1619            }
1620          }
1621        }
1622        edgeIdxGloList.resizeAndPreserve(nEdge);
1623
1624        // (3.6) Saving remaining variables edge_faces and face_faces
1625        edge_faces.resize(2, edge_count);
1626        face_faces.resize(nvertex, nbFaces_);
1627
1628        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
1629        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
1630        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
1631
1632        for (int nf = 0; nf < nbFaces_; ++nf)
1633        {
1634          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1635          {
1636            // Getting global indexes of edge's nodes
1637            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1638            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1639            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1640
1641            size_t myNodeIdx1 = generateNodeIndex(hashValues1);
1642            size_t myNodeIdx2 = generateNodeIndex(hashValues2);
1643            if (myNodeIdx1 != myNodeIdx2)
1644            {
1645              it1 = nodeIdx2IdxGlo.find(myNodeIdx1);
1646              it2 = nodeIdx2IdxGlo.find(myNodeIdx2);
1647              size_t nodeIdxGlo1 = (it1->second)[1];
1648              size_t nodeIdxGlo2 = (it2->second)[1];
1649              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1650              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
1651              int edgeIdxGlo = (itIdx->second)[1];
1652
1653              size_t faceIdxGlo = nbFacesAccum + nf;
1654
1655              if (mpiRank == (itIdx->second)[0])
1656              {
1657                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1658                int face1 = it1->second[0];
1659                if (it1->second.size() == 1)
1660                {
1661                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
1662                  edge_faces(1, edgeIdxGlo - edge_start) = -999;
1663                  face_faces(nv1, nf) = 999999;
1664                }
1665                else
1666                {
1667                  size_t face2 = it1->second[1];
1668                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
1669                  edge_faces(1, edgeIdxGlo - edge_start) = face2;
1670                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1671                }
1672              }
1673              else
1674              {
1675                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1676                int face1 = it1->second[0];
1677                int face2 = it1->second[1];
1678                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1679              }
1680            } // myNodeIdx1 != myNodeIdx2
1681            else
1682              face_faces(nv1, nf) = 999999;
1683          }
1684        }
1685
1686      }
1687      facesAreWritten = true;
1688    } // nvertex >= 3
1689
1690  } // createMeshEpsilon
1691
1692  ///----------------------------------------------------------------
1693  /*!
1694   * \fn  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1695                                              const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1696                                              CArray<int, 2>& nghbFaces)
1697   * Finds neighboring cells of a local domain for node-type of neighbors.
1698   * \param [in] comm
1699   * \param [in] face_idx Array with global indexes.
1700   * \param [in] bounds_lon Array of boundary longitudes.
1701   * \param [in] bounds_lat Array of boundary latitudes.
1702   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
1703   */
1704
1705  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1706                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1707                               CArray<int, 2>& nghbFaces)
1708  {
1709    int nvertex = bounds_lon.rows();
1710    int nbFaces = bounds_lon.shape()[1];
1711    nghbFaces.resize(2, nbFaces*10);    // some estimate on max number of neighbouring cells
1712
1713    int mpiRank, mpiSize;
1714    MPI_Comm_rank(comm, &mpiRank);
1715    MPI_Comm_size(comm, &mpiSize);
1716
1717    // (1) Generating unique node indexes
1718    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1719    vector<size_t> hashValues(4);
1720    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1721    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
1722    size_t iIdx = 0;
1723    for (int nf = 0; nf < nbFaces; ++nf)
1724    {
1725      for (int nv = 0; nv < nvertex; ++nv)
1726      {
1727        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1728        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1729        for (int nh = 0; nh < 4; ++nh)
1730        {
1731          if (nodeHash2Idx.count(hashValues[nh])==0)
1732         {
1733            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1734            nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1735           nodeHashList(iIdx) = hashValues[nh];
1736           ++iIdx;
1737         }
1738        }
1739      }
1740    }
1741    nodeHashList.resizeAndPreserve(iIdx);
1742
1743    // (1.2) Generating node indexes
1744    // The ownership criterion: priority of the process holding the smaller index
1745    // Maps generated in this step are:
1746    // nodeHash2Info = <hash, idx1, idx2, idx3....>
1747    // nodeIdx2IdxMin = <idx, idxMin>
1748    // idxMin is a unique node identifier
1749
1750    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1751    dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1752    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1753
1754    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
1755
1756    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1757    {
1758      size_t idxMin = (it->second)[0];
1759      size_t idx = (it->second)[0];
1760      for (int i = 2; i < (it->second).size();)
1761      {
1762        if (mpiRank == (it->second)[i+1])
1763        {
1764          idx = (it->second)[i];
1765        }
1766        if ((it->second)[i] < idxMin)
1767        {
1768          idxMin = (it->second)[i];
1769          (it->second)[i] = (it->second)[i-2];
1770          (it->second)[i+1] = (it->second)[i-1];
1771        }
1772        i += 2;
1773      }
1774      (it->second)[0] = idxMin;
1775      if (nodeIdx2IdxMin.count(idx) == 0)
1776      {
1777        nodeIdx2IdxMin[idx].push_back(idxMin);
1778      }
1779    }
1780
1781    // (2) Creating maps nodeIdxMin2Face = <nodeIdxMin, [face1, face2, ...]>
1782    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
1783    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdxMin2Face;
1784    CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4);
1785
1786    size_t nNode = 0;
1787
1788    for (int nf = 0; nf < nbFaces; ++nf)
1789    {
1790      for (int nv = 0; nv < nvertex; ++nv)
1791      {
1792        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1793        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
1794        it = nodeIdx2IdxMin.find(myNodeIdx);
1795        size_t nodeIdxMin = (it->second)[0];
1796        size_t faceIdx = face_idx(nf);
1797        if (nodeIdxMin2Face.count(nodeIdxMin) == 0)
1798        {
1799          nodeIdxMinList(nNode) = nodeIdxMin;
1800          ++nNode;
1801        }
1802        nodeIdxMin2Face[nodeIdxMin].push_back(faceIdx);
1803        nodeIdxMin2Face[nodeIdxMin].push_back(mpiRank);
1804      }
1805    }
1806    nodeIdxMinList.resizeAndPreserve(nNode);
1807
1808    // (3) Face_face connectivity
1809
1810    // nodeIdxMin2Info = <nodeIdxMin, [face1, face2,...]>
1811    CClientClientDHTSizet dhtNode2Face (nodeIdxMin2Face, comm);
1812    dhtNode2Face.computeIndexInfoMapping(nodeIdxMinList);
1813    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2Info = dhtNode2Face.getInfoIndexMap();
1814    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map
1815
1816    int nbNghb = 0;
1817    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode;
1818
1819    for (int nf = 0; nf < nbFaces; ++nf)
1820    {
1821      for (int nv = 0; nv < nvertex; ++nv)
1822      {
1823        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1824        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
1825        itNode = nodeIdx2IdxMin.find(myNodeIdx);
1826        size_t nodeIdxMin = (itNode->second)[0];
1827
1828        itNode = nodeIdxMin2Info.find(nodeIdxMin);
1829        for (int i = 0; i < itNode->second.size();)
1830        {
1831          size_t face = itNode->second[i];
1832          size_t rank = itNode->second[i+1];
1833          if (rank != mpiRank)
1834            if (mapFaces.count(face) == 0)
1835            {
1836              nghbFaces(0, nbNghb) = face;
1837              nghbFaces(1, nbNghb) = rank;
1838              ++nbNghb;
1839              mapFaces[face].push_back(face);
1840            }
1841          i += 2;
1842        }
1843      }
1844    }
1845    nghbFaces.resizeAndPreserve(2, nbNghb);
1846  } // getGloNghbFacesNodeType
1847
1848  ///----------------------------------------------------------------
1849  /*!
1850   * \fn  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1851                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1852                               CArray<int, 2>& nghbFaces)
1853   * Finds neighboring cells of a local domain for edge-type of neighbors.
1854   * \param [in] comm
1855   * \param [in] face_idx Array with global indexes.
1856   * \param [in] bounds_lon Array of boundary longitudes.
1857   * \param [in] bounds_lat Array of boundary latitudes.
1858   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
1859   */
1860
1861  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1862                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1863                               CArray<int, 2>& nghbFaces)
1864  {
1865    int nvertex = bounds_lon.rows();
1866    int nbFaces = bounds_lon.shape()[1];
1867    nghbFaces.resize(2, nbFaces*10);    // estimate of max number of neighbouring cells
1868
1869    int mpiRank, mpiSize;
1870    MPI_Comm_rank(comm, &mpiRank);
1871    MPI_Comm_size(comm, &mpiSize);
1872
1873    // (1) Generating unique node indexes
1874    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1875    vector<size_t> hashValues(4);
1876    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1877    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
1878    size_t iIdx = 0;
1879    for (int nf = 0; nf < nbFaces; ++nf)
1880    {
1881      for (int nv = 0; nv < nvertex; ++nv)
1882      {
1883        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1884        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1885        for (int nh = 0; nh < 4; ++nh)
1886        {
1887          if (nodeHash2Idx.count(hashValues[nh])==0)
1888          {
1889            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1890            nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1891            nodeHashList(iIdx) = hashValues[nh];
1892            ++iIdx;
1893          }
1894        }
1895      }
1896    }
1897    nodeHashList.resizeAndPreserve(iIdx);
1898
1899    // (1.2) Generating node indexes
1900    // The ownership criterion: priority of the process holding the smaller index
1901    // Maps generated in this step are:
1902    // nodeHash2Info = <hash, idx1, idx2, idx3....>
1903    // nodeIdx2IdxMin = <idx, idxMin>
1904    // idxMin is a unique node identifier
1905
1906    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1907    dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1908    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1909
1910    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
1911
1912    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1913    {
1914      size_t idxMin = (it->second)[0];
1915      size_t idx = (it->second)[0];
1916      for (int i = 2; i < (it->second).size();)
1917      {
1918        if (mpiRank == (it->second)[i+1])
1919        {
1920          idx = (it->second)[i];
1921        }
1922        if ((it->second)[i] < idxMin)
1923        {
1924          idxMin = (it->second)[i];
1925          (it->second)[i] = (it->second)[i-2];
1926          (it->second)[i+1] = (it->second)[i-1];
1927        }
1928        i += 2;
1929      }
1930      (it->second)[0] = idxMin;
1931      if (nodeIdx2IdxMin.count(idx) == 0)
1932      {
1933        nodeIdx2IdxMin[idx].push_back(idxMin);
1934      }
1935    }
1936
1937    // (2) Creating map edgeHash2Face = <edgeHash, [[face1, rank1], [face2, rank2]]>, where rank1 = rank2 = ...
1938
1939    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2, it;
1940    CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Face;
1941    CArray<size_t,1> edgeHashList(nbFaces*nvertex);
1942
1943    size_t nEdge = 0;
1944
1945    for (int nf = 0; nf < nbFaces; ++nf)
1946    {
1947      for (int nv1 = 0; nv1 < nvertex; ++nv1)
1948      {
1949        // Getting indexes of edge's nodes
1950        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1951        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1952        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1953        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1954        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1955        it1 = nodeIdx2IdxMin.find(myNodeIdx1);
1956        it2 = nodeIdx2IdxMin.find(myNodeIdx2);
1957        size_t nodeIdxMin1 = (it1->second)[0];
1958        size_t nodeIdxMin2 = (it2->second)[0];
1959        size_t faceIdx = face_idx(nf);
1960
1961        if (nodeIdxMin1 != nodeIdxMin2)
1962        {
1963          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
1964          if (edgeHash2Face.count(edgeHash) == 0)
1965          {
1966            edgeHashList(nEdge) = edgeHash;
1967            ++nEdge;
1968          }
1969          edgeHash2Face[edgeHash].push_back(faceIdx);
1970          edgeHash2Face[edgeHash].push_back(mpiRank);
1971        } // nodeIdxMin1 != nodeIdxMin2
1972      }
1973    }
1974    edgeHashList.resizeAndPreserve(nEdge);
1975
1976    // (3) Face_face connectivity
1977
1978    int nbNghb = 0;
1979    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode1, itNode2;
1980
1981    // edgeHash2Info = <edgeHash, [[face1, rank1], [face2, rank2]]>
1982    CClientClientDHTSizet dhtEdge2Face (edgeHash2Face, comm);
1983    dhtEdge2Face.computeIndexInfoMapping(edgeHashList);
1984    CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdge2Face.getInfoIndexMap();
1985    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map
1986
1987    for (int nf = 0; nf < nbFaces; ++nf)
1988    {
1989      for (int nv1 = 0; nv1 < nvertex; ++nv1)
1990      {
1991        // Getting indexes of edge's nodes
1992        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1993        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1994        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1995
1996        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1997        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1998        itNode1 = nodeIdx2IdxMin.find(myNodeIdx1);
1999        itNode2 = nodeIdx2IdxMin.find(myNodeIdx2);
2000        size_t nodeIdxMin1 = (itNode1->second)[0];
2001        size_t nodeIdxMin2 = (itNode2->second)[0];
2002
2003        if (nodeIdxMin1 != nodeIdxMin2)
2004        {
2005          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
2006          it1 = edgeHash2Info.find(edgeHash);
2007
2008          for (int i = 0; i < it1->second.size();)
2009          {
2010            size_t face = it1->second[i];
2011            size_t rank = it1->second[i+1];
2012            if (rank != mpiRank)
2013              if (mapFaces.count(face) == 0)
2014              {
2015                nghbFaces(0, nbNghb) = face;
2016                nghbFaces(1, nbNghb) = rank;
2017                ++nbNghb;
2018                mapFaces[face].push_back(face);
2019              }
2020            i += 2;
2021          }
2022        } // nodeIdxMin1 != nodeIdxMin2
2023      }
2024    }
2025    nghbFaces.resizeAndPreserve(2, nbNghb);
2026  } // getGloNghbFacesEdgeType
2027
2028  ///----------------------------------------------------------------
2029  /*!
2030   * \fn void getGlobalNghbFaces (const int nghbType, const MPI_Comm& comm, const CArray<int, 1>& face_idx,
2031                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2032                                  CArray<size_t, 1>& nghbFaces)
2033   * Finds neighboring faces owned by other procs.
2034   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
2035   * \param [in] comm
2036   * \param [in] face_idx Array with global indexes.
2037   * \param [in] bounds_lon Array of boundary longitudes.
2038   * \param [in] bounds_lat Array of boundary latitudes.
2039   * \param [out] nghbFaces 2D array containing neighboring faces and owner ranks.
2040   */
2041
2042  void CMesh::getGlobalNghbFaces(const int nghbType, const MPI_Comm& comm,
2043                                 const CArray<int, 1>& face_idx,
2044                                 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2045                                 CArray<int, 2>& nghbFaces)
2046  {
2047    if (nghbType == 0)
2048      getGloNghbFacesNodeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
2049    else
2050      getGloNghbFacesEdgeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
2051  } // getGlobalNghbFaces
2052
2053  ///----------------------------------------------------------------
2054  /*!
2055   * \fn void getLocalNghbFaces (const int nghbType,
2056                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2057                                  CArray<size_t, 1>& nghbFaces)
2058   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
2059   * \param [in] bounds_lon Array of boundary longitudes.
2060   * \param [in] bounds_lat Array of boundary latitudes.
2061   * \param [out] nghbFaces 1D array containing neighboring faces.
2062   */
2063
2064  void CMesh::getLocalNghbFaces(const int nghbType, const CArray<int, 1>& face_idx,
2065                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2066                                CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
2067  {
2068    if (nghbType == 0)
2069      getLocNghbFacesNodeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
2070    else
2071      getLocNghbFacesEdgeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
2072  } // getLocalNghbFaces
2073
2074  ///----------------------------------------------------------------
2075  /*!
2076   * \fn void getLocNghbFacesNodeType (const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2077                                       CArray<int, 2>& nghbFaces)
2078   * \param [in] face_idx Array with local face indexing.
2079   * \param [in] bounds_lon Array of boundary longitudes.
2080   * \param [in] bounds_lat Array of boundary latitudes.
2081   * \param [out] nghbFaces 2D array containing neighboring faces.
2082   * \param [out] nbNghbFaces Array containing number of neighboring faces.
2083   */
2084
2085  void CMesh::getLocNghbFacesNodeType (const CArray<int, 1>& face_idx,
2086                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2087                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
2088  {
2089    int nvertex = bounds_lon.rows();
2090    int nbFaces = bounds_lon.shape()[1];
2091    int nbNodes = 0;
2092    nbNghbFaces.resize(nbFaces);
2093    nbNghbFaces = 0;
2094
2095    // nodeToFaces connectivity
2096    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces;
2097    for (int nf = 0; nf < nbFaces; ++nf)
2098      for (int nv = 0; nv < nvertex; ++nv)
2099      {
2100        size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0];
2101        nodeToFaces[nodeHash].push_back(face_idx(nf));
2102      }
2103
2104    // faceToFaces connectivity
2105    boost::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
2106    int maxNb = 20;                            // some assumption on the max possible number of neighboring cells
2107    faceToFaces.resize(maxNb, nbFaces);
2108    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
2109    for (it = nodeToFaces.begin(); it != nodeToFaces.end(); ++it)
2110    {
2111      int size = it->second.size();
2112      for (int i = 0; i < (size-1); ++i)
2113      {
2114        int face1 = it->second[i];
2115        for (int j = i+1; j < size; ++j)
2116        {
2117          int face2 = it->second[j];
2118          if (face2 != face1)
2119          {
2120            int hashFace = hashPairOrdered(face1, face2);
2121            if (mapFaces.count(hashFace) == 0)
2122            {
2123              faceToFaces(nbNghbFaces(face1), face1) = face2;
2124              faceToFaces(nbNghbFaces(face2), face2) = face1;
2125              ++nbNghbFaces(face1);
2126              ++nbNghbFaces(face2);
2127              mapFaces[hashFace] = hashFace;
2128            }
2129          }
2130        }
2131      }
2132    }
2133  } //getLocNghbFacesNodeType
2134
2135
2136  ///----------------------------------------------------------------
2137  /*!
2138   * \fn void getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
2139   *                                   const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2140   *                                   CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
2141   * \param [in] face_idx Array with local face indexing.
2142   * \param [in] bounds_lon Array of boundary longitudes.
2143   * \param [in] bounds_lat Array of boundary latitudes.
2144   * \param [out] nghbFaces 2D array containing neighboring faces.
2145   * \param [out] nbNghbFaces Array containing number of neighboring faces.
2146   */
2147
2148  void CMesh::getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
2149                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2150                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
2151  {
2152    int nvertex = bounds_lon.rows();
2153    int nbFaces = bounds_lon.shape()[1];
2154    int nbNodes = 0;
2155    int nbEdges = 0;
2156    nbNghbFaces.resize(nbFaces);
2157    nbNghbFaces = 0;
2158
2159    // faceToNodes connectivity
2160    CArray<double, 2> faceToNodes (nvertex, nbFaces);
2161
2162    boost::unordered_map <pair<double,double>, int> mapNodes;
2163
2164    for (int nf = 0; nf < nbFaces; ++nf)
2165      for (int nv = 0; nv < nvertex; ++nv)
2166      {
2167        if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end())
2168        {
2169          mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes;
2170          faceToNodes(nv,nf) = nbNodes ;
2171          ++nbNodes ;
2172        }
2173        else
2174          faceToNodes(nv,nf) = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
2175      }
2176
2177    // faceToFaces connectivity
2178    boost::unordered_map <pair<int,int>, int> mapEdges;
2179    faceToFaces.resize(nvertex, nbFaces);
2180    CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible
2181
2182    for (int nf = 0; nf < nbFaces; ++nf)
2183    {
2184      for (int nv1 = 0; nv1 < nvertex; ++nv1)
2185      {
2186        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
2187        int face = face_idx(nf);
2188        int node1 = faceToNodes(nv1,face);
2189        int node2 = faceToNodes(nv2,face);
2190        if (node1 != node2)
2191        {
2192          if (mapEdges.find(make_ordered_pair (node1, node2)) == mapEdges.end())
2193          {
2194            mapEdges[make_ordered_pair (node1, node2)] = nbEdges;
2195            edgeToFaces(0,nbEdges) = face;
2196            ++nbEdges;
2197          }
2198          else
2199          {
2200            int edge = mapEdges[make_ordered_pair (node1, node2)];
2201            edgeToFaces(1, edge) = face;
2202            int face1 = face;
2203            int face2 = edgeToFaces(0,edge);
2204            faceToFaces(nbNghbFaces(face1), face1) =  face2;
2205            faceToFaces(nbNghbFaces(face2), face2) =  face1;
2206            ++nbNghbFaces(face1);
2207            ++nbNghbFaces(face2);
2208          }
2209        } // node1 != node2
2210      } // nv
2211    } // nf
2212
2213  } //getLocNghbFacesEdgeType
2214
2215
2216} // namespace xios
Note: See TracBrowser for help on using the repository browser.