source: XIOS3/branches/xios-3.0-beta/src/node/mesh.cpp @ 2417

Last change on this file since 2417 was 2417, checked in by jderouillat, 21 months ago

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