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

Last change on this file since 2230 was 1936, checked in by ymipsl, 4 years ago

XIOS coupling branch :
Blitz array resizeAndPreserve crash when size=0 => use resize instead

YM
  • Property svn:executable set to *
File size: 83.0 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 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    MPI_Comm_rank(comm, &mpiRank);
424    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      MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_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        if (nbHash==0) nodeHashList.resize(nbHash);
542        else nodeHashList.resizeAndPreserve(nbHash);
543
544        // (2.2) Generating global node indexes
545        // The ownership criterion: priority of the process of smaller index
546        // Maps generated in this step are:
547        // Maps generated in this step are:
548         // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
549         // nodeIdx2Idx = <idx, <rankOwner, idx>>
550
551         CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
552         dhtNodeHash.computeIndexInfoMapping(nodeHashList);
553         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
554
555
556         CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
557         CArray<size_t,1> nodeIdxList(nbEdges_*nvertex*4);
558         size_t nIdx = 0;
559
560         for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
561         {
562           size_t rankMin = (it->second)[1];
563           size_t idx = (it->second)[0];
564           for (int i = 2; i < (it->second).size();)
565           {
566             if ( (it->second)[i+1] < rankMin)
567             {
568               idx = (it->second)[i];
569               rankMin = (it->second)[i+1];
570               (it->second)[i+1] = (it->second)[i-1];
571             }
572             i += 2;
573           }
574           if (nodeIdx2Idx.count(idx) == 0)
575           {
576             if (mpiRank == rankMin)
577             {
578               nodeIdx2Idx[idx].push_back(rankMin);
579               nodeIdx2Idx[idx].push_back(idx);
580             }
581             nodeIdxList(nIdx) = idx;
582             ++nIdx;
583           }
584         }
585         nodeIdxList.resizeAndPreserve(nIdx);
586
587         // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
588         // Solution: global node indexing by hand.
589         // Maps modified in this step:
590         // nodeIdx2Idx = <idx, idxGlo>
591         unsigned long nodeCount = nodeIdx2Idx.size();
592         unsigned long nodeStart, nbNodes;
593         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
594         int nNodes = nodeStart;
595         MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
596         nbNodesGlo = nNodes;
597
598         nodeStart -= nodeCount;
599         node_start = nodeStart;
600         node_count = nodeCount;
601         CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
602         size_t count = 0;
603
604         for (int ne = 0; ne < nbEdges_; ++ne)
605         {
606           for (int nv = 0; nv < nvertex; ++nv)
607           {
608             vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
609             size_t nodeIdx = generateNodeIndex(hashValues);
610             CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
611             if (it != nodeIdx2Idx.end())
612             {
613               if (dummyMap.count(nodeIdx) == 0)
614               {
615                 dummyMap[nodeIdx].push_back(nodeIdx);
616                 (it->second)[1] = node_start + count;
617                 ++count;
618               }
619             }
620           }
621         }
622
623         CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
624         dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
625         CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
626
627        // (2.3) Saving variables: node_lon, node_lat, edge_nodes
628        // Creating map nodeHash2IdxGlo <hash, idxGlo>
629        // Creating map edgeHash2IdxGlo <hash, idxGlo>
630//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
631//        node_count = dhtNodeIdxGlo.getIndexCount();
632//        node_start = dhtNodeIdxGlo.getIndexStart();
633        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
634        node_lon.resize(node_count);
635        node_lat.resize(node_count);
636        vector <size_t> edgeNodes;
637        size_t idxGlo = 0;
638
639        for (int ne = 0; ne < nbEdges_; ++ne)
640        {
641          for (int nv = 0; nv < nvertex; ++nv)
642          {
643            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
644            size_t myIdx = generateNodeIndex(hashValues);
645            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxGlo.find(myIdx);
646            idxGlo = (itIdx->second)[1];
647
648            if (mpiRank == (itIdx->second)[0])
649            {
650//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
651              node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
652              node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
653            }
654            edge_nodes(nv,ne) = idxGlo;
655            for (int nh = 0; nh < 4; ++nh)
656              nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo);
657            edgeNodes.push_back(idxGlo);
658          }
659          if (edgeNodes[0] != edgeNodes[1])
660          {
661            size_t edgeIdxGlo = nbEdgesAccum + ne;
662            edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(edgeIdxGlo);
663          }
664          edgeNodes.clear();
665        }
666        pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
667      } // !nodesAreWritten
668
669      pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm);
670      edgesAreWritten = true;
671    } //nvertex = 2
672
673    else
674    {
675      nbFaces_ = bounds_lon.shape()[1];
676      face_lon.resize(nbFaces_);
677      face_lat.resize(nbFaces_);
678      face_lon = lonvalue;
679      face_lat = latvalue;
680      face_nodes.resize(nvertex, nbFaces_);
681      face_edges.resize(nvertex, nbFaces_);
682
683      // For determining the global face index
684      unsigned long nbFacesOnProc = nbFaces_;
685      unsigned long nbFacesAccum;
686      MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
687      nbFacesAccum -= nbFaces_;
688
689      // Case (1): edges have been previously generated
690      if (edgesAreWritten)
691      {
692        // (1.1) Recuperating node global indexing and saving face_nodes
693        vector<size_t> hashValues(4);
694        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
695        for (int nf = 0; nf < nbFaces_; ++nf)
696        {
697          for (int nv = 0; nv < nvertex; ++nv)
698            {
699              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
700              for (int nh = 0; nh < 4; ++nh)
701                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
702            }
703        }
704        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
705        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
706        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
707        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
708        size_t nEdge = 0;
709        for (int nf = 0; nf < nbFaces_; ++nf)
710        {
711          for (int nv1 = 0; nv1 < nvertex; ++nv1)
712          {
713            int nh1 = 0;
714            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
715            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
716            while (it1 == nodeHash2IdxGlo.end())
717            {
718              ++nh1;
719              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
720            }
721            int nh2 = 0;
722            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
723            while (it2 == nodeHash2IdxGlo.end())
724            {
725              ++nh2;
726              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
727            }
728            face_nodes(nv1,nf) = it1->second[0];
729            if (it1->second[0] != it2->second[0])
730            {
731              edgeHashList(nEdge) = hashPairOrdered(it1->second[0], it2->second[0]);
732              ++nEdge;
733            }
734          }
735        }
736        edgeHashList.resizeAndPreserve(nEdge);
737
738        // (1.2) Recuperating edge global indexing and saving face_edges
739        pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList);
740        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap();
741        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash;
742        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
743        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
744        CArray<size_t,1> edgeIdxList(nbFaces_*nvertex);
745        size_t iIdx = 0;
746
747        for (int nf = 0; nf < nbFaces_; ++nf)
748        {
749          for (int nv1 = 0; nv1 < nvertex; ++nv1)
750          {
751            int nh1 = 0;
752            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
753            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
754            while (it1 == nodeHash2IdxGlo.end())
755            {
756              ++nh1;
757              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
758            }
759            int nh2 = 0;
760            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
761            while (it2 == nodeHash2IdxGlo.end())
762            {
763              ++nh2;
764              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
765            }
766            if (it1->second[0] != it2->second[0])
767            {
768              size_t faceIdxGlo = nbFacesAccum + nf;
769              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
770              itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
771              size_t edgeIdxGlo = itEdgeHash->second[0];
772              face_edges(nv1,nf) = edgeIdxGlo;
773              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
774              {
775                edgeIdxList(iIdx) = edgeIdxGlo;
776                ++iIdx;
777              }
778              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
779              edgeHash2Rank[edgeHash].push_back(mpiRank);
780              edgeHash2Rank[edgeHash].push_back(itEdgeHash->second[0]);
781            }
782            else
783            {
784              face_edges(nv1,nf) = 999999;
785            }
786          }
787        }
788        edgeIdxList.resizeAndPreserve(iIdx);
789
790        // (1.3) Saving remaining variables edge_faces and face_faces
791
792        // Establishing edge ownership
793        // The ownership criterion: priority of the process with smaller rank
794        CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm);
795        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
796        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
797
798        // edgeHash2Info = <edgeHash, < rank1, idxGlo, rank2, idxGlo>>
799        int edgeCount = 0;
800        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
801        {
802          vector <size_t> edgeInfo = it->second;
803          if (edgeInfo[0] == mpiRank)
804          {
805            ++edgeCount;
806          }
807        }
808
809        unsigned long edgeStart, nbEdges;
810        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
811        int nEdges = edgeStart;
812        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
813        nbEdgesGlo = nEdges;
814
815        // edges to be splitted equally between procs
816        if ( (nbEdgesGlo % mpiSize) == 0)
817        {
818          edge_count = nbEdgesGlo/mpiSize;
819          edge_start = mpiRank*edge_count;
820        }
821        else
822        {
823          if (mpiRank == (mpiSize - 1) )
824          {
825            edge_count = nbEdgesGlo/mpiSize;
826            edge_start = mpiRank*(nbEdgesGlo/mpiSize + 1);
827          }
828          else
829          {
830            edge_count = nbEdgesGlo/mpiSize + 1;
831            edge_start = mpiRank*edge_count;
832          }
833        }
834        CArray<size_t,1> edgeIdxGloList(edge_count);
835        for (int i = 0; i < edge_count; ++i)
836        {
837          edgeIdxGloList(i) = i + edge_start;
838        }
839
840        CClientClientDHTSizet dhtEdgeIdxGlo2Face (edgeIdxGlo2Face, comm);
841        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
842        dhtEdgeIdxGlo2Face.computeIndexInfoMapping(edgeIdxGloList);
843        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdgeIdxGlo2Face.getInfoIndexMap();
844        dhtEdge2Face.computeIndexInfoMapping(edgeIdxList);
845        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2FaceIdx = dhtEdge2Face.getInfoIndexMap();
846
847
848        edge_faces.resize(2, edge_count);
849        for (int i = 0; i < edge_count; ++i)
850        {
851          CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdxGlo2FaceIdx.find(i + edge_start);
852          int indexGlo = it->first;
853          vector<size_t> faces = it->second;
854          int face1 = faces[0];
855          edge_faces(0, indexGlo - edge_start) = face1;
856          if (faces.size() == 2)
857          {
858            int face2 = faces[1];
859            edge_faces(1, indexGlo - edge_start) = face2;
860          }
861          else
862          {
863            edge_faces(1, indexGlo - edge_start) = -999;
864          }
865        }
866
867        size_t tmp;
868        vector<size_t> tmpVec;
869        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2FaceIdx.begin(); it != edgeIdx2FaceIdx.end(); it++)
870        {
871          tmp = it->first;
872          tmpVec = it->second;
873          tmp++;
874        }
875
876        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2, itIndex;
877        face_faces.resize(nvertex, nbFaces_);
878        for (int nf = 0; nf < nbFaces_; ++nf)
879        {
880          for (int nv1 = 0; nv1 < nvertex; ++nv1)
881          {
882            int nh1 = 0;
883            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
884            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
885            while (it1 == nodeHash2IdxGlo.end())
886            {
887              ++nh1;
888              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
889            }
890            int nh2 = 0;
891            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
892            while (it2 == nodeHash2IdxGlo.end())
893            {
894              ++nh2;
895              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
896            }
897
898            if (it1->second[0] != it2->second[0])
899            {
900              size_t faceIdxGlo = nbFacesAccum + nf;
901              size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
902              itEdgeHash = edgeHash2Info.find(edgeHash);
903              int edgeIdxGlo = (itEdgeHash->second)[1];
904
905              if ( (itEdgeHash->second)[0] == mpiRank)
906              {
907                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
908                int face1 = itFace1->second[0];
909                if (itFace1->second.size() == 1)
910                {
911                  face_faces(nv1, nf) = 999999;
912                }
913                else
914                {
915                  int face2 = itFace1->second[1];
916                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
917                }
918              } // edge owner
919              else
920              {
921                itFace1 = edgeIdx2FaceIdx.find(edgeIdxGlo);
922                int face1 = itFace1->second[0];
923                int face2 = itFace1->second[1];
924                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
925              } // not an edge owner
926            } // node1 != node2
927            else
928            {
929              face_faces(nv1, nf) = 999999;
930            }
931          }
932        }
933      } // edgesAreWritten
934
935      // Case (2): nodes have been previously generated
936      else if (nodesAreWritten)
937      {
938        // (2.1) Generating nodeHashList
939        vector<size_t> hashValues(4);
940        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
941        for (int nf = 0; nf < nbFaces_; ++nf)
942        {
943          for (int nv = 0; nv < nvertex; ++nv)
944            {
945              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
946              for (int nh = 0; nh < 4; ++nh)
947                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
948            }
949        }
950
951        // (2.2) Recuperating node global indexing and saving face_nodes
952        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
953        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
954        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
955        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
956        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
957        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
958        int nEdgeHash = 0;
959        for (int nf = 0; nf < nbFaces_; ++nf)
960        {
961          for (int nv1 = 0; nv1 < nvertex; ++nv1)
962          {
963            int nh1 = 0;
964            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
965            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
966            while (it1 == nodeHash2IdxGlo.end())
967            {
968              ++nh1;
969              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
970            }
971            int nh2 = 0;
972            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
973            while (it2 == nodeHash2IdxGlo.end())
974            {
975              ++nh2;
976              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
977            }
978            face_nodes(nv1,nf) = it1->second[0];
979            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
980            if (edgeHash2Idx.count(edgeHash) == 0)
981            {
982              edgeHash2Idx[edgeHash].push_back(edgeHash);
983              edgeHash2Idx[edgeHash].push_back(mpiRank);
984              edgeHashList(nEdgeHash) = edgeHash;
985              ++nEdgeHash;
986            }
987          }
988        }
989        if (nEdgeHash==0) edgeHashList.resize(nEdgeHash);
990        else edgeHashList.resizeAndPreserve(nEdgeHash);
991
992        // (2.3) Generating global edge indexes
993        // The ownership criterion: priority of the process with smaller rank
994        // Maps generated in this step are:
995        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
996        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
997
998        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
999        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
1000        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
1001        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
1002
1003        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
1004
1005        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
1006        {
1007          size_t rankMin = (it->second)[1];
1008          size_t idx = (it->second)[0];
1009
1010          for (int i = 2; i < (it->second).size();)
1011          {
1012            if ((it->second)[i+1] < rankMin)
1013            {
1014              rankMin = (it->second)[i+1];
1015              idx = (it->second)[i];
1016              (it->second)[i+1] = (it->second)[i-1];
1017            }
1018            i += 2;
1019          }
1020          if (edgeIdx2Idx.count(idx) == 0)
1021          {
1022            if (mpiRank == rankMin)
1023            {
1024              edgeIdx2Idx[idx].push_back(rankMin);
1025              edgeIdx2Idx[idx].push_back(idx);
1026            }
1027          }
1028        }
1029
1030        unsigned long edgeCount = edgeIdx2Idx.size();
1031        unsigned long edgeStart, nbEdges;
1032        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
1033        int nEdges = edgeStart;
1034        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
1035        nbEdgesGlo = nEdges;
1036
1037        edgeStart -= edgeCount;
1038        edge_start = edgeStart;
1039        edge_count = edgeCount;
1040        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
1041        int count = 0;
1042
1043        for (int nf = 0; nf < nbFaces_; ++nf)
1044        {
1045          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1046          {
1047            // Getting global indexes of edge's nodes
1048            int nh1 = 0;
1049            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1050            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1051            while (it1 == nodeHash2IdxGlo.end())
1052            {
1053              ++nh1;
1054              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1055            }
1056            int nh2 = 0;
1057            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1058            while (it2 == nodeHash2IdxGlo.end())
1059            {
1060              ++nh2;
1061              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1062            }
1063            size_t nodeIdxGlo1 = it1->second[0];
1064            size_t nodeIdxGlo2 = it2->second[0];
1065
1066            if (nodeIdxGlo1 != nodeIdxGlo2)
1067            {
1068              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1069              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
1070              if (it != edgeIdx2Idx.end())
1071              {
1072                if (dummyEdgeMap.count(edgeIdx) == 0)
1073                {
1074                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
1075                  (it->second)[1] = edge_start + count;
1076                  ++count;
1077                }
1078              }
1079            }
1080          }
1081        }
1082
1083        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
1084        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
1085        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
1086
1087        // (2.4) Saving variables: edge_lon, edge_lat, face_edges
1088        edge_lon.resize(edge_count);
1089        edge_lat.resize(edge_count);
1090        edge_nodes.resize(2, edge_count);
1091        face_edges.resize(nvertex, nbFaces_);
1092
1093        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
1094        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
1095        size_t iIdx = 0;
1096
1097        for (int nf = 0; nf < nbFaces_; ++nf)
1098        {
1099          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1100          {
1101            // Getting global indexes of edge's nodes
1102            int nh1 = 0;
1103            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1104            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1105            while (it1 == nodeHash2IdxGlo.end())
1106            {
1107              ++nh1;
1108              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1109            }
1110            int nh2 = 0;
1111            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1112            while (it2 == nodeHash2IdxGlo.end())
1113            {
1114              ++nh2;
1115              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1116            }
1117            // Getting edge global index
1118            size_t nodeIdxGlo1 = it1->second[0];
1119            size_t nodeIdxGlo2 = it2->second[0];
1120            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1121            if (nodeIdxGlo1 != nodeIdxGlo2)
1122            {
1123              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
1124              int edgeIdxGlo = (itIdx->second)[1];
1125              size_t faceIdxGlo = nbFacesAccum + nf;
1126
1127              if (mpiRank == (itIdx->second)[0])
1128              {
1129                double edgeLon;
1130                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
1131                if (diffLon < (180.- prec))
1132                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
1133                else if (diffLon > (180.+ prec))
1134                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
1135                else
1136                  edgeLon = 0.;
1137                edge_lon(edgeIdxGlo - edge_start) = edgeLon;
1138                edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
1139                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
1140                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
1141              }
1142              face_edges(nv1,nf) = edgeIdxGlo;
1143              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
1144              {
1145                edgeIdxGloList(iIdx) = edgeIdxGlo;
1146                ++iIdx;
1147              }
1148              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
1149            } // nodeIdxGlo1 != nodeIdxGlo2
1150            else
1151            {
1152              face_edges(nv1,nf) = 999999;
1153            }
1154          }
1155        }
1156        edgeIdxGloList.resizeAndPreserve(iIdx);
1157
1158        // (2.5) Saving remaining variables edge_faces and face_faces
1159        edge_faces.resize(2, edge_count);
1160        face_faces.resize(nvertex, nbFaces_);
1161
1162        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
1163        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
1164        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
1165        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
1166        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx;
1167
1168        for (int nf = 0; nf < nbFaces_; ++nf)
1169        {
1170          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1171          {
1172            // Getting global indexes of edge's nodes
1173            int nh1 = 0;
1174            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1175            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1176            while (it1 == nodeHash2IdxGlo.end())
1177            {
1178              ++nh1;
1179              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1180            }
1181            int nh2 = 0;
1182            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1183            while (it2 == nodeHash2IdxGlo.end())
1184            {
1185              ++nh2;
1186              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1187            }
1188            size_t nodeIdxGlo1 = it1->second[0];
1189            size_t nodeIdxGlo2 = it2->second[0];
1190
1191            size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1192            itIdx = edgeIdx2IdxGlo.find(myIdx);
1193            size_t faceIdxGlo = nbFacesAccum + nf;
1194            int edgeIdxGlo = (itIdx->second)[1];
1195
1196            if (mpiRank == (itIdx->second)[0])
1197            {
1198              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1199              int face1 = it1->second[0];
1200              if (it1->second.size() == 1)
1201              {
1202                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1203                edge_faces(1, edgeIdxGlo - edge_start) = -999;
1204                face_faces(nv1, nf) = 999999;
1205              }
1206              else
1207              {
1208                int face2 = it1->second[1];
1209                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1210                edge_faces(1, edgeIdxGlo - edge_start) = face2;
1211                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1212              }
1213            }
1214            else
1215            {
1216              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1217              int face1 = it1->second[0];
1218              int face2 = it1->second[1];
1219              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1220            }
1221          }
1222        }
1223      } // nodesAreWritten
1224
1225      // Case (3): Neither nodes nor edges have been previously generated
1226      else
1227      {
1228        // (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1229        vector<size_t> hashValues(4);
1230        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1231        CArray<size_t,1> nodeHashList(nbFaces_*nvertex*4);
1232        size_t iHash = 0;
1233        for (int nf = 0; nf < nbFaces_; ++nf)
1234        {
1235          for (int nv = 0; nv < nvertex; ++nv)
1236          {
1237            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1238//            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1239            size_t nodeIndex = generateNodeIndex(hashValues);
1240            for (int nh = 0; nh < 4; ++nh)
1241            {
1242              if (nodeHash2Idx.count(hashValues[nh])==0)
1243              {
1244                nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1245                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1246                nodeHashList(iHash) = hashValues[nh];
1247                ++iHash;
1248              }
1249            }
1250          }
1251        }
1252        nodeHashList.resizeAndPreserve(iHash);
1253
1254        // (3.2) Generating global node indexes
1255        // The ownership criterion: priority of the process with smaller rank.
1256        // With any other criterion it is not possible to have consistent node indexing for different number of procs.
1257        // Maps generated in this step are:
1258        // nodeHash2Info = <hash, [[idx, rankMin], [idx, rank1], [idx, rank3]..]>
1259        // nodeIdx2Idx = <idx, <rankOwner, idx>>
1260
1261        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1262        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1263        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1264
1265        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2Idx;
1266        CArray<size_t,1> nodeIdxList(nbFaces_*nvertex*4);
1267        size_t nIdx = 0;
1268
1269        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1270        {
1271          size_t rankMin = (it->second)[1];
1272          size_t idx = (it->second)[0];
1273          for (int i = 2; i < (it->second).size();)
1274          {
1275            if ( (it->second)[i+1] < rankMin)
1276            {
1277              idx = (it->second)[i];
1278              rankMin = (it->second)[i+1];
1279              (it->second)[i+1] = (it->second)[i-1];
1280            }
1281            i += 2;
1282          }
1283          if (nodeIdx2Idx.count(idx) == 0)
1284          {
1285            if (mpiRank == rankMin)
1286            {
1287              nodeIdx2Idx[idx].push_back(rankMin);
1288              nodeIdx2Idx[idx].push_back(idx);
1289            }
1290            nodeIdxList(nIdx) = idx;
1291            ++nIdx;
1292          }
1293        }
1294
1295//        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2Idx, comm);
1296        // CDHTAutoIndexing will not give consistent node numbering for varying number of procs. =>
1297        // Solution: global node indexing by hand.
1298        // Maps modified in this step:
1299        // nodeIdx2Idx = <idx, idxGlo>
1300        unsigned long nodeCount = nodeIdx2Idx.size();
1301        unsigned long nodeStart, nbNodes;
1302        MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
1303        int nNodes = nodeStart;
1304        MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
1305        nbNodesGlo = nNodes;
1306
1307        nodeStart -= nodeCount;
1308        node_start = nodeStart;
1309        node_count = nodeCount;
1310        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyMap; // just a dummy map used to ensure that each node is numbered only once
1311        size_t count = 0;
1312
1313        for (int nf = 0; nf < nbFaces_; ++nf)
1314        {
1315          for (int nv = 0; nv < nvertex; ++nv)
1316          {
1317            vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1318            size_t nodeIdx = generateNodeIndex(hashValues);
1319            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeIdx2Idx.find(nodeIdx);
1320            if (it != nodeIdx2Idx.end())
1321            {
1322              if (dummyMap.count(nodeIdx) == 0)
1323              {
1324                dummyMap[nodeIdx].push_back(nodeIdx);
1325                (it->second)[1] = node_start + count;
1326                ++count;
1327              }
1328            }
1329          }
1330        }
1331        if (nIdx==0) nodeIdxList.resize(nIdx);
1332        else nodeIdxList.resizeAndPreserve(nIdx);
1333        CClientClientDHTSizet dhtNodeIdx(nodeIdx2Idx, comm);
1334        dhtNodeIdx.computeIndexInfoMapping(nodeIdxList);
1335        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdx2IdxGlo = dhtNodeIdx.getInfoIndexMap();
1336
1337        // (3.3) Saving node data: node_lon, node_lat, and face_nodes
1338        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
1339//        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
1340//        node_count = dhtNodeIdxGlo.getIndexCount();
1341//        node_start = dhtNodeIdxGlo.getIndexStart();
1342        node_lon.resize(node_count);
1343        node_lat.resize(node_count);
1344        size_t nodeIdxGlo1 = 0;
1345        size_t nodeIdxGlo2 = 0;
1346        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
1347        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
1348        CArray<size_t,1> edgeHashList(nbFaces_*nvertex);
1349        size_t nEdgeHash = 0;
1350
1351        for (int nf = 0; nf < nbFaces_; ++nf)
1352        {
1353          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1354          {
1355            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1356            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1357            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1358            size_t nodeIdx1 = generateNodeIndex(hashValues1);
1359            size_t nodeIdx2 = generateNodeIndex(hashValues2);
1360            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
1361            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
1362            size_t ownerRank = (itNodeIdx1->second)[0];
1363            nodeIdxGlo1 = (itNodeIdx1->second)[1];
1364            nodeIdxGlo2 = (itNodeIdx2->second)[1];
1365
1366            if (mpiRank == ownerRank)
1367            {
1368              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
1369              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
1370            }
1371            if (nodeIdxGlo1 != nodeIdxGlo2)
1372            {
1373              size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1374              edgeHash2Idx[edgeHash].push_back(edgeHash);
1375              edgeHash2Idx[edgeHash].push_back(mpiRank);
1376              edgeHashList(nEdgeHash) = edgeHash;
1377              ++nEdgeHash;
1378            }
1379            face_nodes(nv1,nf) = nodeIdxGlo1;
1380          }
1381        }
1382        if (nEdgeHash==0) edgeHashList.resize(nEdgeHash);
1383        else edgeHashList.resizeAndPreserve(nEdgeHash);
1384
1385        // (3.4) Generating global edge indexes
1386        // Maps generated in this step are:
1387        // edgeIdx2Idx = = <idx, <rankOwner, idx>>
1388        // edgeIdx2IdxGlo = <idxMin, <rankOwner, idxGlo>>
1389
1390        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
1391        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
1392        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
1393        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
1394
1395        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2Idx;
1396
1397        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
1398        {
1399          size_t rankMin = (it->second)[1];
1400          size_t idx = (it->second)[0];
1401
1402          for (int i = 2; i < (it->second).size();)
1403          {
1404            if ((it->second)[i+1] < rankMin)
1405            {
1406              rankMin = (it->second)[i+1];
1407              idx = (it->second)[i];
1408              (it->second)[i+1] = (it->second)[i-1];
1409            }
1410            i += 2;
1411          }
1412          if (edgeIdx2Idx.count(idx) == 0)
1413          {
1414            if (mpiRank == rankMin)
1415            {
1416              edgeIdx2Idx[idx].push_back(rankMin);
1417              edgeIdx2Idx[idx].push_back(idx);
1418            }
1419          }
1420        }
1421
1422        unsigned long edgeCount = edgeIdx2Idx.size();
1423        unsigned long edgeStart, nbEdges;
1424        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
1425        int nEdges = edgeStart;
1426        MPI_Bcast(&nEdges, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
1427        nbEdgesGlo = nEdges;
1428
1429        edgeStart -= edgeCount;
1430        edge_start = edgeStart;
1431        edge_count = edgeCount;
1432        CClientClientDHTSizet::Index2VectorInfoTypeMap dummyEdgeMap;
1433        count = 0;
1434
1435        for (int nf = 0; nf < nbFaces_; ++nf)
1436        {
1437          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1438          {
1439            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1440            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1441            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1442            size_t nodeIdx1 = generateNodeIndex(hashValues1);
1443            size_t nodeIdx2 = generateNodeIndex(hashValues2);
1444            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxGlo.find(nodeIdx1);
1445            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxGlo.find(nodeIdx2);
1446            nodeIdxGlo1 = (itNodeIdx1->second)[1];
1447            nodeIdxGlo2 = (itNodeIdx2->second)[1];
1448
1449            if (nodeIdxGlo1 != nodeIdxGlo2)
1450            {
1451              size_t edgeIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1452              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeIdx2Idx.find(edgeIdx);
1453              if (it != edgeIdx2Idx.end())
1454              {
1455                if (dummyEdgeMap.count(edgeIdx) == 0)
1456                {
1457                  dummyEdgeMap[edgeIdx].push_back(edgeIdx);
1458                  (it->second)[1] = edge_start + count;
1459                  ++count;
1460                }
1461              }
1462            }
1463          }
1464        }
1465        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2Idx, comm);
1466        dhtEdgeIdx.computeIndexInfoMapping(edgeHashList);
1467        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdx2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
1468
1469        // (3.5) Saving variables: edge_lon, edge_lat, face_edges
1470        // Creating map edgeIdxGlo2Face <idxGlo, face>
1471//        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
1472//        edge_count = dhtEdgeIdxGlo.getIndexCount();
1473//        edge_start = dhtEdgeIdxGlo.getIndexStart();
1474
1475        edge_lon.resize(edge_count);
1476        edge_lat.resize(edge_count);
1477        edge_nodes.resize(2, edge_count);
1478        face_edges.resize(nvertex, nbFaces_);
1479
1480        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
1481        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
1482        CArray<size_t,1> edgeIdxGloList(nbFaces_*nvertex);
1483        size_t nEdge = 0;
1484
1485        for (int nf = 0; nf < nbFaces_; ++nf)
1486        {
1487          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1488          {
1489            // Getting global indexes of edge's nodes
1490            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1491            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1492            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1493
1494            size_t nodeIdx1 = generateNodeIndex(hashValues1);
1495            size_t nodeIdx2 = generateNodeIndex(hashValues2);
1496            it1 = nodeIdx2IdxGlo.find(nodeIdx1);
1497            it2 = nodeIdx2IdxGlo.find(nodeIdx2);
1498            size_t nodeIdxGlo1 = (it1->second)[1];
1499            size_t nodeIdxGlo2 = (it2->second)[1];
1500
1501            if (nodeIdxGlo1 != nodeIdxGlo2)
1502            {
1503              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1504              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
1505              int edgeIdxGlo = (itIdx->second)[1];
1506              size_t faceIdxGlo = nbFacesAccum + nf;
1507
1508              if (mpiRank == (itIdx->second)[0])
1509              {
1510                double edgeLon;
1511                double diffLon = abs(bounds_lon(nv1, nf) - bounds_lon(nv2, nf));
1512                if (diffLon < (180.- prec))
1513                  edgeLon = ( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5;
1514                else if (diffLon > (180.+ prec))
1515                  edgeLon = (bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.;
1516                else
1517                  edgeLon = 0.;
1518                edge_lon(edgeIdxGlo - edge_start) = edgeLon;
1519                edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
1520                edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
1521                edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
1522              }
1523              face_edges(nv1,nf) = edgeIdxGlo;
1524              if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
1525              {
1526                edgeIdxGloList(nEdge) = edgeIdxGlo;
1527                ++nEdge;
1528              }
1529              edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
1530            } // nodeIdxGlo1 != nodeIdxGlo2
1531            else
1532            {
1533              face_edges(nv1,nf) = 999999;
1534            }
1535          }
1536        }
1537        edgeIdxGloList.resizeAndPreserve(nEdge);
1538
1539        // (3.6) Saving remaining variables edge_faces and face_faces
1540        edge_faces.resize(2, edge_count);
1541        face_faces.resize(nvertex, nbFaces_);
1542
1543        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
1544        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
1545        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
1546
1547        for (int nf = 0; nf < nbFaces_; ++nf)
1548        {
1549          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1550          {
1551            // Getting global indexes of edge's nodes
1552            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1553            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1554            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1555
1556            size_t myNodeIdx1 = generateNodeIndex(hashValues1);
1557            size_t myNodeIdx2 = generateNodeIndex(hashValues2);
1558            if (myNodeIdx1 != myNodeIdx2)
1559            {
1560              it1 = nodeIdx2IdxGlo.find(myNodeIdx1);
1561              it2 = nodeIdx2IdxGlo.find(myNodeIdx2);
1562              size_t nodeIdxGlo1 = (it1->second)[1];
1563              size_t nodeIdxGlo2 = (it2->second)[1];
1564              size_t myIdx = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1565              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxGlo.find(myIdx);
1566              int edgeIdxGlo = (itIdx->second)[1];
1567
1568              size_t faceIdxGlo = nbFacesAccum + nf;
1569
1570              if (mpiRank == (itIdx->second)[0])
1571              {
1572                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1573                int face1 = it1->second[0];
1574                if (it1->second.size() == 1)
1575                {
1576                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
1577                  edge_faces(1, edgeIdxGlo - edge_start) = -999;
1578                  face_faces(nv1, nf) = 999999;
1579                }
1580                else
1581                {
1582                  size_t face2 = it1->second[1];
1583                  edge_faces(0, edgeIdxGlo - edge_start) = face1;
1584                  edge_faces(1, edgeIdxGlo - edge_start) = face2;
1585                  face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1586                }
1587              }
1588              else
1589              {
1590                it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1591                int face1 = it1->second[0];
1592                int face2 = it1->second[1];
1593                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1594              }
1595            } // myNodeIdx1 != myNodeIdx2
1596            else
1597              face_faces(nv1, nf) = 999999;
1598          }
1599        }
1600
1601      }
1602      facesAreWritten = true;
1603    } // nvertex >= 3
1604
1605  } // createMeshEpsilon
1606
1607  ///----------------------------------------------------------------
1608  /*!
1609   * \fn  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1610                                              const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1611                                              CArray<int, 2>& nghbFaces)
1612   * Finds neighboring cells of a local domain for node-type of neighbors.
1613   * \param [in] comm
1614   * \param [in] face_idx Array with global indexes.
1615   * \param [in] bounds_lon Array of boundary longitudes.
1616   * \param [in] bounds_lat Array of boundary latitudes.
1617   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
1618   */
1619
1620  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1621                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1622                               CArray<int, 2>& nghbFaces)
1623  {
1624    int nvertex = bounds_lon.rows();
1625    int nbFaces = bounds_lon.shape()[1];
1626    nghbFaces.resize(2, nbFaces*10);    // some estimate on max number of neighbouring cells
1627
1628    int mpiRank, mpiSize;
1629    MPI_Comm_rank(comm, &mpiRank);
1630    MPI_Comm_size(comm, &mpiSize);
1631
1632    // (1) Generating unique node indexes
1633    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1634    vector<size_t> hashValues(4);
1635    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1636    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
1637    size_t iIdx = 0;
1638    for (int nf = 0; nf < nbFaces; ++nf)
1639    {
1640      for (int nv = 0; nv < nvertex; ++nv)
1641      {
1642        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1643        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1644        for (int nh = 0; nh < 4; ++nh)
1645        {
1646          if (nodeHash2Idx.count(hashValues[nh])==0)
1647         {
1648            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1649            nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1650           nodeHashList(iIdx) = hashValues[nh];
1651           ++iIdx;
1652         }
1653        }
1654      }
1655    }
1656    nodeHashList.resizeAndPreserve(iIdx);
1657
1658    // (1.2) Generating node indexes
1659    // The ownership criterion: priority of the process holding the smaller index
1660    // Maps generated in this step are:
1661    // nodeHash2Info = <hash, idx1, idx2, idx3....>
1662    // nodeIdx2IdxMin = <idx, idxMin>
1663    // idxMin is a unique node identifier
1664
1665    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1666    dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1667    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1668
1669    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
1670
1671    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1672    {
1673      size_t idxMin = (it->second)[0];
1674      size_t idx = (it->second)[0];
1675      for (int i = 2; i < (it->second).size();)
1676      {
1677        if (mpiRank == (it->second)[i+1])
1678        {
1679          idx = (it->second)[i];
1680        }
1681        if ((it->second)[i] < idxMin)
1682        {
1683          idxMin = (it->second)[i];
1684          (it->second)[i] = (it->second)[i-2];
1685          (it->second)[i+1] = (it->second)[i-1];
1686        }
1687        i += 2;
1688      }
1689      (it->second)[0] = idxMin;
1690      if (nodeIdx2IdxMin.count(idx) == 0)
1691      {
1692        nodeIdx2IdxMin[idx].push_back(idxMin);
1693      }
1694    }
1695
1696    // (2) Creating maps nodeIdxMin2Face = <nodeIdxMin, [face1, face2, ...]>
1697    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
1698    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdxMin2Face;
1699    CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4);
1700
1701    size_t nNode = 0;
1702
1703    for (int nf = 0; nf < nbFaces; ++nf)
1704    {
1705      for (int nv = 0; nv < nvertex; ++nv)
1706      {
1707        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1708        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
1709        it = nodeIdx2IdxMin.find(myNodeIdx);
1710        size_t nodeIdxMin = (it->second)[0];
1711        size_t faceIdx = face_idx(nf);
1712        if (nodeIdxMin2Face.count(nodeIdxMin) == 0)
1713        {
1714          nodeIdxMinList(nNode) = nodeIdxMin;
1715          ++nNode;
1716        }
1717        nodeIdxMin2Face[nodeIdxMin].push_back(faceIdx);
1718        nodeIdxMin2Face[nodeIdxMin].push_back(mpiRank);
1719      }
1720    }
1721    nodeIdxMinList.resizeAndPreserve(nNode);
1722
1723    // (3) Face_face connectivity
1724
1725    // nodeIdxMin2Info = <nodeIdxMin, [face1, face2,...]>
1726    CClientClientDHTSizet dhtNode2Face (nodeIdxMin2Face, comm);
1727    dhtNode2Face.computeIndexInfoMapping(nodeIdxMinList);
1728    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2Info = dhtNode2Face.getInfoIndexMap();
1729    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map
1730
1731    int nbNghb = 0;
1732    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode;
1733
1734    for (int nf = 0; nf < nbFaces; ++nf)
1735    {
1736      for (int nv = 0; nv < nvertex; ++nv)
1737      {
1738        vector<size_t> hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1739        size_t myNodeIdx = generateNodeIndex(hashValues, mpiRank);
1740        itNode = nodeIdx2IdxMin.find(myNodeIdx);
1741        size_t nodeIdxMin = (itNode->second)[0];
1742
1743        itNode = nodeIdxMin2Info.find(nodeIdxMin);
1744        for (int i = 0; i < itNode->second.size();)
1745        {
1746          size_t face = itNode->second[i];
1747          size_t rank = itNode->second[i+1];
1748          if (rank != mpiRank)
1749            if (mapFaces.count(face) == 0)
1750            {
1751              nghbFaces(0, nbNghb) = face;
1752              nghbFaces(1, nbNghb) = rank;
1753              ++nbNghb;
1754              mapFaces[face].push_back(face);
1755            }
1756          i += 2;
1757        }
1758      }
1759    }
1760    if (nbNghb==0) nghbFaces.resize(2, nbNghb);
1761    else nghbFaces.resizeAndPreserve(2, nbNghb);
1762  } // getGloNghbFacesNodeType
1763
1764  ///----------------------------------------------------------------
1765  /*!
1766   * \fn  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1767                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1768                               CArray<int, 2>& nghbFaces)
1769   * Finds neighboring cells of a local domain for edge-type of neighbors.
1770   * \param [in] comm
1771   * \param [in] face_idx Array with global indexes.
1772   * \param [in] bounds_lon Array of boundary longitudes.
1773   * \param [in] bounds_lat Array of boundary latitudes.
1774   * \param [out] nghbFaces 2D array of storing global indexes of neighboring cells and their owner procs.
1775   */
1776
1777  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1778                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1779                               CArray<int, 2>& nghbFaces)
1780  {
1781    int nvertex = bounds_lon.rows();
1782    int nbFaces = bounds_lon.shape()[1];
1783    nghbFaces.resize(2, nbFaces*10);    // estimate of max number of neighbouring cells
1784
1785    int mpiRank, mpiSize;
1786    MPI_Comm_rank(comm, &mpiRank);
1787    MPI_Comm_size(comm, &mpiSize);
1788
1789    // (1) Generating unique node indexes
1790    // (1.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1791    vector<size_t> hashValues(4);
1792    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1793    CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
1794    size_t iIdx = 0;
1795    for (int nf = 0; nf < nbFaces; ++nf)
1796    {
1797      for (int nv = 0; nv < nvertex; ++nv)
1798      {
1799        hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1800        size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1801        for (int nh = 0; nh < 4; ++nh)
1802        {
1803          if (nodeHash2Idx.count(hashValues[nh])==0)
1804          {
1805            nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1806            nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1807            nodeHashList(iIdx) = hashValues[nh];
1808            ++iIdx;
1809          }
1810        }
1811      }
1812    }
1813    if (iIdx==0) nodeHashList.resize(iIdx);
1814    else nodeHashList.resizeAndPreserve(iIdx);
1815
1816    // (1.2) Generating node indexes
1817    // The ownership criterion: priority of the process holding the smaller index
1818    // Maps generated in this step are:
1819    // nodeHash2Info = <hash, idx1, idx2, idx3....>
1820    // nodeIdx2IdxMin = <idx, idxMin>
1821    // idxMin is a unique node identifier
1822
1823    CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1824    dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1825    CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1826
1827    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
1828
1829    for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1830    {
1831      size_t idxMin = (it->second)[0];
1832      size_t idx = (it->second)[0];
1833      for (int i = 2; i < (it->second).size();)
1834      {
1835        if (mpiRank == (it->second)[i+1])
1836        {
1837          idx = (it->second)[i];
1838        }
1839        if ((it->second)[i] < idxMin)
1840        {
1841          idxMin = (it->second)[i];
1842          (it->second)[i] = (it->second)[i-2];
1843          (it->second)[i+1] = (it->second)[i-1];
1844        }
1845        i += 2;
1846      }
1847      (it->second)[0] = idxMin;
1848      if (nodeIdx2IdxMin.count(idx) == 0)
1849      {
1850        nodeIdx2IdxMin[idx].push_back(idxMin);
1851      }
1852    }
1853
1854    // (2) Creating map edgeHash2Face = <edgeHash, [[face1, rank1], [face2, rank2]]>, where rank1 = rank2 = ...
1855
1856    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2, it;
1857    CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Face;
1858    CArray<size_t,1> edgeHashList(nbFaces*nvertex);
1859
1860    size_t nEdge = 0;
1861
1862    for (int nf = 0; nf < nbFaces; ++nf)
1863    {
1864      for (int nv1 = 0; nv1 < nvertex; ++nv1)
1865      {
1866        // Getting indexes of edge's nodes
1867        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1868        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1869        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1870        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1871        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1872        it1 = nodeIdx2IdxMin.find(myNodeIdx1);
1873        it2 = nodeIdx2IdxMin.find(myNodeIdx2);
1874        size_t nodeIdxMin1 = (it1->second)[0];
1875        size_t nodeIdxMin2 = (it2->second)[0];
1876        size_t faceIdx = face_idx(nf);
1877
1878        if (nodeIdxMin1 != nodeIdxMin2)
1879        {
1880          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
1881          if (edgeHash2Face.count(edgeHash) == 0)
1882          {
1883            edgeHashList(nEdge) = edgeHash;
1884            ++nEdge;
1885          }
1886          edgeHash2Face[edgeHash].push_back(faceIdx);
1887          edgeHash2Face[edgeHash].push_back(mpiRank);
1888        } // nodeIdxMin1 != nodeIdxMin2
1889      }
1890    }
1891    edgeHashList.resizeAndPreserve(nEdge);
1892
1893    // (3) Face_face connectivity
1894
1895    int nbNghb = 0;
1896    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNode1, itNode2;
1897
1898    // edgeHash2Info = <edgeHash, [[face1, rank1], [face2, rank2]]>
1899    CClientClientDHTSizet dhtEdge2Face (edgeHash2Face, comm);
1900    dhtEdge2Face.computeIndexInfoMapping(edgeHashList);
1901    CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdge2Face.getInfoIndexMap();
1902    CClientClientDHTSizet::Index2VectorInfoTypeMap mapFaces;  // auxiliar map
1903
1904    for (int nf = 0; nf < nbFaces; ++nf)
1905    {
1906      for (int nv1 = 0; nv1 < nvertex; ++nv1)
1907      {
1908        // Getting indexes of edge's nodes
1909        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1910        vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1911        vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1912
1913        size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1914        size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1915        itNode1 = nodeIdx2IdxMin.find(myNodeIdx1);
1916        itNode2 = nodeIdx2IdxMin.find(myNodeIdx2);
1917        size_t nodeIdxMin1 = (itNode1->second)[0];
1918        size_t nodeIdxMin2 = (itNode2->second)[0];
1919
1920        if (nodeIdxMin1 != nodeIdxMin2)
1921        {
1922          size_t edgeHash = hashPairOrdered(nodeIdxMin1, nodeIdxMin2);
1923          it1 = edgeHash2Info.find(edgeHash);
1924
1925          for (int i = 0; i < it1->second.size();)
1926          {
1927            size_t face = it1->second[i];
1928            size_t rank = it1->second[i+1];
1929            if (rank != mpiRank)
1930              if (mapFaces.count(face) == 0)
1931              {
1932                nghbFaces(0, nbNghb) = face;
1933                nghbFaces(1, nbNghb) = rank;
1934                ++nbNghb;
1935                mapFaces[face].push_back(face);
1936              }
1937            i += 2;
1938          }
1939        } // nodeIdxMin1 != nodeIdxMin2
1940      }
1941    }
1942    if (nbNghb==0) nghbFaces.resize(2, nbNghb);
1943    else nghbFaces.resizeAndPreserve(2, nbNghb);
1944  } // getGloNghbFacesEdgeType
1945
1946  ///----------------------------------------------------------------
1947  /*!
1948   * \fn void getGlobalNghbFaces (const int nghbType, const MPI_Comm& comm, const CArray<int, 1>& face_idx,
1949                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1950                                  CArray<size_t, 1>& nghbFaces)
1951   * Finds neighboring faces owned by other procs.
1952   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
1953   * \param [in] comm
1954   * \param [in] face_idx Array with global indexes.
1955   * \param [in] bounds_lon Array of boundary longitudes.
1956   * \param [in] bounds_lat Array of boundary latitudes.
1957   * \param [out] nghbFaces 2D array containing neighboring faces and owner ranks.
1958   */
1959
1960  void CMesh::getGlobalNghbFaces(const int nghbType, const MPI_Comm& comm,
1961                                 const CArray<int, 1>& face_idx,
1962                                 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1963                                 CArray<int, 2>& nghbFaces)
1964  {
1965    if (nghbType == 0)
1966      getGloNghbFacesNodeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
1967    else
1968      getGloNghbFacesEdgeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces);
1969  } // getGlobalNghbFaces
1970
1971  ///----------------------------------------------------------------
1972  /*!
1973   * \fn void getLocalNghbFaces (const int nghbType,
1974                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1975                                  CArray<size_t, 1>& nghbFaces)
1976   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges.
1977   * \param [in] bounds_lon Array of boundary longitudes.
1978   * \param [in] bounds_lat Array of boundary latitudes.
1979   * \param [out] nghbFaces 1D array containing neighboring faces.
1980   */
1981
1982  void CMesh::getLocalNghbFaces(const int nghbType, const CArray<int, 1>& face_idx,
1983                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1984                                CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
1985  {
1986    if (nghbType == 0)
1987      getLocNghbFacesNodeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
1988    else
1989      getLocNghbFacesEdgeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces);
1990  } // getLocalNghbFaces
1991
1992  ///----------------------------------------------------------------
1993  /*!
1994   * \fn void getLocNghbFacesNodeType (const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
1995                                       CArray<int, 2>& nghbFaces)
1996   * \param [in] face_idx Array with local face indexing.
1997   * \param [in] bounds_lon Array of boundary longitudes.
1998   * \param [in] bounds_lat Array of boundary latitudes.
1999   * \param [out] nghbFaces 2D array containing neighboring faces.
2000   * \param [out] nbNghbFaces Array containing number of neighboring faces.
2001   */
2002
2003  void CMesh::getLocNghbFacesNodeType (const CArray<int, 1>& face_idx,
2004                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2005                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
2006  {
2007    int nvertex = bounds_lon.rows();
2008    int nbFaces = bounds_lon.shape()[1];
2009    int nbNodes = 0;
2010    nbNghbFaces.resize(nbFaces);
2011    nbNghbFaces = 0;
2012
2013    // nodeToFaces connectivity
2014    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces;
2015    for (int nf = 0; nf < nbFaces; ++nf)
2016      for (int nv = 0; nv < nvertex; ++nv)
2017      {
2018        size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0];
2019        nodeToFaces[nodeHash].push_back(face_idx(nf));
2020      }
2021
2022    // faceToFaces connectivity
2023    std::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
2024    int maxNb = 20;                            // some assumption on the max possible number of neighboring cells
2025    faceToFaces.resize(maxNb, nbFaces);
2026    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
2027    for (it = nodeToFaces.begin(); it != nodeToFaces.end(); ++it)
2028    {
2029      int size = it->second.size();
2030      for (int i = 0; i < (size-1); ++i)
2031      {
2032        int face1 = it->second[i];
2033        for (int j = i+1; j < size; ++j)
2034        {
2035          int face2 = it->second[j];
2036          if (face2 != face1)
2037          {
2038            int hashFace = hashPairOrdered(face1, face2);
2039            if (mapFaces.count(hashFace) == 0)
2040            {
2041              faceToFaces(nbNghbFaces(face1), face1) = face2;
2042              faceToFaces(nbNghbFaces(face2), face2) = face1;
2043              ++nbNghbFaces(face1);
2044              ++nbNghbFaces(face2);
2045              mapFaces[hashFace] = hashFace;
2046            }
2047          }
2048        }
2049      }
2050    }
2051  } //getLocNghbFacesNodeType
2052
2053
2054  ///----------------------------------------------------------------
2055  /*!
2056   * \fn void getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
2057   *                                   const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2058   *                                   CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces)
2059   * \param [in] face_idx Array with local face indexing.
2060   * \param [in] bounds_lon Array of boundary longitudes.
2061   * \param [in] bounds_lat Array of boundary latitudes.
2062   * \param [out] nghbFaces 2D array containing neighboring faces.
2063   * \param [out] nbNghbFaces Array containing number of neighboring faces.
2064   */
2065
2066  void CMesh::getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx,
2067                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat,
2068                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces)
2069  {
2070    int nvertex = bounds_lon.rows();
2071    int nbFaces = bounds_lon.shape()[1];
2072    int nbNodes = 0;
2073    int nbEdges = 0;
2074    nbNghbFaces.resize(nbFaces);
2075    nbNghbFaces = 0;
2076
2077    // faceToNodes connectivity
2078    CArray<double, 2> faceToNodes (nvertex, nbFaces);
2079
2080    std::unordered_map <pairDouble, int, boost::hash<pairDouble> > mapNodes;
2081
2082    for (int nf = 0; nf < nbFaces; ++nf)
2083      for (int nv = 0; nv < nvertex; ++nv)
2084      {
2085        if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end())
2086        {
2087          mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes;
2088          faceToNodes(nv,nf) = nbNodes ;
2089          ++nbNodes ;
2090        }
2091        else
2092          faceToNodes(nv,nf) = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
2093      }
2094
2095    // faceToFaces connectivity
2096    std::unordered_map <pairInt, int, boost::hash<pairInt> > mapEdges;
2097    faceToFaces.resize(nvertex, nbFaces);
2098    CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible
2099
2100    for (int nf = 0; nf < nbFaces; ++nf)
2101    {
2102      for (int nv1 = 0; nv1 < nvertex; ++nv1)
2103      {
2104        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
2105        int face = face_idx(nf);
2106        int node1 = faceToNodes(nv1,face);
2107        int node2 = faceToNodes(nv2,face);
2108        if (node1 != node2)
2109        {
2110          if (mapEdges.find(make_ordered_pair (node1, node2)) == mapEdges.end())
2111          {
2112            mapEdges[make_ordered_pair (node1, node2)] = nbEdges;
2113            edgeToFaces(0,nbEdges) = face;
2114            ++nbEdges;
2115          }
2116          else
2117          {
2118            int edge = mapEdges[make_ordered_pair (node1, node2)];
2119            edgeToFaces(1, edge) = face;
2120            int face1 = face;
2121            int face2 = edgeToFaces(0,edge);
2122            faceToFaces(nbNghbFaces(face1), face1) =  face2;
2123            faceToFaces(nbNghbFaces(face2), face2) =  face1;
2124            ++nbNghbFaces(face1);
2125            ++nbNghbFaces(face2);
2126          }
2127        } // node1 != node2
2128      } // nv
2129    } // nf
2130
2131  } //getLocNghbFacesEdgeType
2132
2133
2134} // namespace xios
Note: See TracBrowser for help on using the repository browser.