source: XIOS/trunk/src/node/mesh.cpp @ 929

Last change on this file since 929 was 929, checked in by oabramkina, 5 years ago

Mesh connectivity:

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

Function for determining face neighbors of a local domain added.

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