source: XIOS/dev/dev_ym/XIOS_COUPLING/src/node/mesh.cpp @ 2284

Last change on this file since 2284 was 2284, checked in by jderouillat, 2 years ago

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