source: XIOS/dev/dev_trunk_omp/src/node/mesh.cpp @ 1646

Last change on this file since 1646 was 1646, checked in by yushan, 5 years ago

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