source: XIOS2/trunk/src/node/mesh.cpp @ 2515

Last change on this file since 2515 was 2515, checked in by jderouillat, 13 months ago

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