source: XIOS/dev/branch_yushan/src/node/mesh.cpp @ 1109

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

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