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

Last change on this file since 1507 was 1507, checked in by oabramkina, 6 years ago

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