source: XIOS/dev/branch_yushan_merged/src/node/mesh.cpp @ 1134

Last change on this file since 1134 was 1134, checked in by yushan, 7 years ago

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