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

Last change on this file since 1638 was 1638, checked in by yushan, 6 years ago

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