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

Last change on this file since 2200 was 2200, checked in by ymipsl, 3 years ago

Bugs fix for UGRID convention output

  • global indexation was not taking into account
  • coherence problem in connectivity for node, edge and face mesh
  • add UGRID testcase based on sphere cube

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