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

Last change on this file since 2250 was 2250, checked in by ymipsl, 2 years ago

Fix issue for Ugrid convention output.

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