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

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

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