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

Last change on this file since 924 was 924, checked in by oabramkina, 8 years ago

Parallel version of UGRID norms.

The following connectivity parameters have been implemented:

edge_nodes
face_nodes
face_edges
edge_faces
face_faces.

Test with a regular(structured) quadrilateral mesh (test_regular) has been added.
File size: 57.3 KB
Line 
1/*!
2  \file mesh.cpp
3  \author Olga Abramkina
4  \brief Definition of class CMesh.
5*/
6
7#include "mesh.hpp"
8
9namespace xios {
10
11/// ////////////////////// Définitions ////////////////////// ///
12
13  CMesh::CMesh(void) :  nbNodesGlo{0}, nbEdgesGlo{0}
14            ,  node_start{0}, node_count{0}
15            ,  edge_start{0}, edge_count{0}
16            ,  nbFaces{0}, nbNodes{0}, nbEdges{0}
17            ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false}
18            ,  node_lon(), node_lat()
19            ,  edge_lon(), edge_lat(), edge_nodes()
20            ,  face_lon(), face_lat()
21            ,  face_nodes()
22            ,  pNodeGlobalIndex{NULL}, pEdgeGlobalIndex{NULL}
23  {
24  }
25
26
27  CMesh::~CMesh(void)
28  {
29    if (pNodeGlobalIndex != NULL) delete pNodeGlobalIndex;
30    if (pEdgeGlobalIndex != NULL) delete pEdgeGlobalIndex;
31  }
32
33  std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>();
34  std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >();
35
36///---------------------------------------------------------------
37/*!
38 * \fn bool CMesh::getMesh (StdString meshName)
39 * Returns a pointer to a mesh. If a mesh has not been created, creates it and adds its name to the list of meshes meshList.
40 * \param [in] meshName  The name of a mesh ("name" attribute of a domain).
41 * \param [in] nvertex Number of verteces (1 for nodes, 2 for edges, 3 and up for faces).
42 */
43  CMesh* CMesh::getMesh (StdString meshName, int nvertex)
44  {
45    CMesh::domainList[meshName].push_back(nvertex);
46
47    if ( CMesh::meshList.begin() != CMesh::meshList.end() )
48    {
49      for (std::map<StdString, CMesh>::iterator it=CMesh::meshList.begin(); it!=CMesh::meshList.end(); ++it)
50      {
51        if (it->first == meshName)
52          return &meshList[meshName];
53        else
54        {
55          CMesh newMesh;
56          CMesh::meshList.insert( make_pair(meshName, newMesh) );
57          return &meshList[meshName];
58        }
59      }
60    }
61    else
62    {
63      CMesh newMesh;
64      CMesh::meshList.insert( make_pair(meshName, newMesh) );
65      return &meshList[meshName];
66    }
67  }
68
69///----------------------------------------------------------------
70  size_t hashPair(size_t first, size_t second)
71  {
72    HashXIOS<size_t> sizetHash;
73    size_t seed = sizetHash(first) + 0x9e3779b9 ;
74    seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
75    return seed ;
76  }
77
78///----------------------------------------------------------------
79  size_t hashPairOrdered(size_t first, size_t second)
80  {
81    size_t seed;
82    HashXIOS<size_t> sizetHash;
83    if (first < second)
84    {
85      seed = sizetHash(first) + 0x9e3779b9 ;
86      seed ^= sizetHash(second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
87    }
88    else
89    {
90      seed = sizetHash(second) + 0x9e3779b9 ;
91      seed ^= sizetHash(first) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
92    }
93    return seed ;
94  }
95
96///----------------------------------------------------------------
97/*!
98 * \fn size_t generateEdgeIndex(size_t first, size_t second, int rank)
99 * Generates an edge index.
100 * If the same edge is generated by two processes, each process will have its own edge index.
101 * \param [in] first Edge node.
102 * \param [in] second Edge node.
103 * \param [in] rank MPI process rank.
104 */
105  size_t generateEdgeIndex(size_t first, size_t second, int rank)
106  {
107    size_t seed = rank ;
108    if (first < second)
109    {
110      seed = hashPair(seed, first);
111      seed = hashPair(seed, second);
112    }
113    else
114    {
115      seed = hashPair(seed, second);
116      seed = hashPair(seed, first);
117    }
118
119    return seed ;
120  }
121
122///----------------------------------------------------------------
123/*!
124 * \fn size_t generateNodeIndex(vector<size_t>& valList, int rank)
125 * Generates a node index.
126 * If the same node is generated by two processes, each process will have its own node index.
127 * \param [in] valList Vector storing four node hashes.
128 * \param [in] rank MPI process rank.
129 */
130  size_t generateNodeIndex(vector<size_t>& valList, int rank)
131  {
132    // Sort is needed to avoid problems for nodes with lon = 0 generated by faces in east and west semisphere
133    vector<size_t> vec = valList;
134    sort (vec.begin(), vec.end());
135    size_t seed = rank ;
136    int it = 0;
137    for(; it != vec.size(); ++it)
138    {
139       seed = hashPair(seed, vec[it]);
140    }
141    return seed ;
142  }
143
144///----------------------------------------------------------------
145/*!
146 * \fn size_t CMesh::nodeIndex (double lon, double lat)
147 * Returns its index if a node exists; otherwise adds the node and returns -1.
148 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
149 * \param [in] lon Node longitude in degrees.
150 * \param [in] lat Node latitude in degrees ranged from 0 to 360.
151 * \return node index if a node exists; -1 otherwise
152 */
153  size_t CMesh::nodeIndex (double lon, double lat)
154  {
155    double minBoundLon = 0. ;
156    double maxBoundLon = 360. ;
157    double minBoundLat = -90 ;
158    double maxBoundLat = 90 ;
159    double prec=1e-11 ;
160    double precLon=prec ;
161    double precLat=prec ;
162
163    size_t maxsize_t=numeric_limits<size_t>::max() ;
164    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
165    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
166
167    size_t iMinLon=0 ;
168    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
169    size_t iMinLat=0 ;
170    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
171
172    size_t hash0,hash1,hash2,hash3 ;
173    size_t lon0,lon1,lat0,lat1 ;
174
175    lon0=(lon-minBoundLon)/precLon ;
176    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
177    {
178      if (lon0==iMinLon) lon1=iMaxLon ;
179      else lon1=lon0-1 ;
180    }
181    else
182    {
183      if (lon0==iMaxLon) lon1=iMinLon ;
184      else lon1=lon0+1 ;
185    }
186
187    lat0=(lat-minBoundLat)/precLat ;
188    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
189    {
190      if (lat0==iMinLat) lat1=lat0 ;
191      else lat1=lat0-1 ;
192    }
193    else
194    {
195      if (lat0==iMaxLat) lat1=lat0 ;
196      else lat1=lat0+1 ;
197    }
198
199    hash0=hashPair(lon0,lat0) ;
200    hash1=hashPair(lon0,lat1) ;
201    hash2=hashPair(lon1,lat0) ;
202    hash3=hashPair(lon1,lat1) ;
203
204    boost::unordered_map<size_t, size_t>::iterator end = hashed_map_nodes.end() ;
205    size_t mapSize = hashed_map_nodes.size();
206    if (hashed_map_nodes.find(hash0)==end && hashed_map_nodes.find(hash1)==end && hashed_map_nodes.find(hash2)==end && hashed_map_nodes.find(hash3)==end)
207    {
208      hashed_map_nodes[hash0] = mapSize ;
209      hashed_map_nodes[hash1] = mapSize + 1;
210      hashed_map_nodes[hash2] = mapSize + 2;
211      hashed_map_nodes[hash3] = mapSize + 3;
212      return -1;
213    }
214    else
215      return ( (hashed_map_nodes[hash0]+1) / 4 );
216
217  } // nodeIndex()
218
219///----------------------------------------------------------------
220/*!
221 * \fn CArray<size_t,1>& CMesh::createHashes (const double longitude, const double latitude)
222 * Creates two hash values for each dimension, longitude and latitude.
223 * \param [in] longitude Node longitude in degrees.
224 * \param [in] latitude Node latitude in degrees ranged from 0 to 360.
225 */
226
227  vector<size_t> CMesh::createHashes (const double longitude, const double latitude)
228  {
229    double minBoundLon = 0. ;
230    double maxBoundLon = 360. ;
231    double minBoundLat = -90. ;
232    double maxBoundLat = 90. ;
233    double prec=1e-11 ;
234    double precLon=prec ;
235    double precLat=prec ;
236    double lon = longitude;
237    double lat = latitude;
238
239    if (lon > (360.- prec)) lon = 0.;
240
241    size_t maxsize_t=numeric_limits<size_t>::max() ;
242    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
243    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
244
245    size_t iMinLon=0 ;
246    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
247    size_t iMinLat=0 ;
248    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
249
250    vector<size_t> hash(4);
251    size_t lon0,lon1,lat0,lat1 ;
252
253    lon0=(lon-minBoundLon)/precLon ;
254    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
255    {
256      if (lon0==iMinLon) lon1=iMaxLon ;
257      else lon1=lon0-1 ;
258    }
259    else
260    {
261      if (lon0==iMaxLon) lon1=iMinLon ;
262      else lon1=lon0+1 ;
263    }
264
265    lat0=(lat-minBoundLat)/precLat ;
266    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
267    {
268      if (lat0==iMinLat) lat1=lat0 ;
269      else lat1=lat0-1 ;
270    }
271    else
272    {
273      if (lat0==iMaxLat) lat1=lat0 ;
274      else lat1=lat0+1 ;
275    }
276
277    hash[0] = hashPair(lon0,lat0) ;
278    hash[1] = hashPair(lon0,lat1) ;
279    hash[2] = hashPair(lon1,lat0) ;
280    hash[3] = hashPair(lon1,lat1) ;
281
282    return hash;
283
284  } // createHashes
285
286///----------------------------------------------------------------
287  std::pair<int,int> make_ordered_pair(int a, int b)
288  {
289    if ( a < b )
290      return std::pair<int,int>(a,b);
291    else
292      return std::pair<int,int>(b,a);
293  }
294
295///----------------------------------------------------------------
296/*!
297 * \fn void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
298            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
299 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
300 * \param [in] lonvalue  Array of longitudes.
301 * \param [in] latvalue  Array of latitudes.
302 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
303 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
304 */
305  void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
306            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
307  {
308    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
309
310    if (nvertex == 1)
311    {
312      nbNodes = lonvalue.numElements();
313      node_lon.resizeAndPreserve(nbNodes);
314      node_lat.resizeAndPreserve(nbNodes);
315      for (int nn = 0; nn < nbNodes; ++nn)
316      {
317        if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end())
318        {
319          map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ;
320          node_lon(nn) = lonvalue(nn);
321          node_lat(nn) = latvalue(nn);
322        }
323      }
324    }
325    else if (nvertex == 2)
326    {
327      nbEdges = bounds_lon.shape()[1];
328
329      // Create nodes and edge_node connectivity
330      node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes
331      node_lat.resizeAndPreserve(nbEdges*nvertex);
332      edge_nodes.resizeAndPreserve(nvertex, nbEdges);
333
334      for (int ne = 0; ne < nbEdges; ++ne)
335      {
336        for (int nv = 0; nv < nvertex; ++nv)
337        {
338          if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end())
339          {
340            map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes ;
341            edge_nodes(nv,ne) = nbNodes ;
342            node_lon(nbNodes) = bounds_lon(nv, ne);
343            node_lat(nbNodes) = bounds_lat(nv, ne);
344            ++nbNodes ;
345          }
346          else
347            edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))];
348        }
349      }
350      node_lon.resizeAndPreserve(nbNodes);
351      node_lat.resizeAndPreserve(nbNodes);
352
353      // Create edges
354      edge_lon.resizeAndPreserve(nbEdges);
355      edge_lat.resizeAndPreserve(nbEdges);
356
357      for (int ne = 0; ne < nbEdges; ++ne)
358      {
359        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
360        {
361          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
362          edge_lon(ne) = lonvalue(ne);
363          edge_lat(ne) = latvalue(ne);
364        }
365
366      }
367      edgesAreWritten = true;
368    }
369    else
370    {
371      nbFaces = bounds_lon.shape()[1];
372 
373      // Create nodes and face_node connectivity
374      node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes
375      node_lat.resizeAndPreserve(nbFaces*nvertex);
376      face_nodes.resize(nvertex, nbFaces);
377 
378      for (int nf = 0; nf < nbFaces; ++nf)
379      {
380        for (int nv = 0; nv < nvertex; ++nv)
381        {
382          if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end())
383          {
384            map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes ;
385            face_nodes(nv,nf) = nbNodes ;
386            node_lon(nbNodes) = bounds_lon(nv, nf);
387            node_lat(nbNodes) = bounds_lat(nv ,nf);
388            ++nbNodes ;
389          }
390          else
391          {
392            face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
393          }
394        }
395      }
396      node_lon.resizeAndPreserve(nbNodes);
397      node_lat.resizeAndPreserve(nbNodes);
398 
399      // Create edges and edge_nodes connectivity
400      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges
401      edge_lat.resizeAndPreserve(nbFaces*nvertex);
402      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex);
403      edge_faces.resize(2, nbFaces*nvertex);
404      face_edges.resize(nvertex, nbFaces);
405      face_faces.resize(nvertex, nbFaces);
406
407      vector<int> countEdges(nbFaces*nvertex);   // needed in case if edges have been already generated
408      vector<int> countFaces(nbFaces);
409      countEdges.assign(nbFaces*nvertex, 0);
410      countFaces.assign(nbFaces, 0);
411      int edge;
412      for (int nf = 0; nf < nbFaces; ++nf)
413      {
414        for (int nv1 = 0; nv1 < nvertex; ++nv1)
415        {
416          int nv = 0;
417          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
418          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
419          {
420            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ;
421            face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
422            edge_faces(0,nbEdges) = nf;
423            edge_faces(1,nbEdges) = -999;
424            face_faces(nv1,nf) = -1;
425            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf);
426            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
427                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
428                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);
429            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
430            ++nbEdges;
431          }
432          else
433          {
434            edge = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
435            face_edges(nv1,nf) = edge;
436            if (edgesAreWritten)
437            {
438              edge_faces(countEdges[edge], edge) = nf;
439              if (countEdges[edge]==0)
440              {
441                face_faces(nv1,nf) = -1;
442              }
443              else
444              {
445                int face1 = nf; // = edge_faces(1,edge)
446                int face2 = edge_faces(0,edge);
447                face_faces(countFaces[face1], face1) =  face2;
448                face_faces(countFaces[face2], face2) =  face1;
449                ++(countFaces[face1]);
450                ++(countFaces[face2]);
451              }
452            }
453            else
454            {
455              edge_faces(1,edge) = nf;
456              int face1 = nf; // = edge_faces(1,edge)
457              int face2 = edge_faces(0,edge);
458              face_faces(countFaces[face1], face1) =  face2;
459              face_faces(countFaces[face2], face2) =  face1;
460              ++(countFaces[face1]);
461              ++(countFaces[face2]);
462            }
463            ++(countEdges[edge]);
464          }
465        }
466      }
467      edge_nodes.resizeAndPreserve(2, nbEdges);
468      edge_faces.resizeAndPreserve(2, nbEdges);
469      edge_lon.resizeAndPreserve(nbEdges);
470      edge_lat.resizeAndPreserve(nbEdges);
471
472      // Create faces
473      face_lon.resize(nbFaces);
474      face_lat.resize(nbFaces);
475      face_lon = lonvalue;
476      face_lat = latvalue;
477      facesAreWritten = true;
478
479    } // nvertex > 2
480   
481  } // createMesh()
482
483///----------------------------------------------------------------
484/*!
485 * \fn void CMesh::createMeshEpsilon(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
486            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
487 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
488 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
489 * \param [in] comm
490 * \param [in] lonvalue  Array of longitudes.
491 * \param [in] latvalue  Array of latitudes.
492 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
493 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
494 */
495  void CMesh::createMeshEpsilon(const MPI_Comm& comm,
496                                const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
497                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
498  {
499
500    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
501    int mpiRank, mpiSize;
502    MPI_Comm_rank(comm, &mpiRank);
503    MPI_Comm_size(comm, &mpiSize);
504
505    if (nvertex == 1)
506    {
507      nbNodes = lonvalue.numElements();
508      node_lon.resize(nbNodes);
509      node_lat.resize(nbNodes);
510      node_lon = lonvalue;
511      node_lat = latvalue;
512
513      // Global node indexes
514      vector<size_t> hashValues(4);
515      CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
516      for (size_t nn = 0; nn < nbNodes; ++nn)
517      {
518        hashValues = CMesh::createHashes(lonvalue(nn), latvalue(nn));
519        for (size_t nh = 0; nh < 4; ++nh)
520        {
521          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes + nn);
522        }
523      }
524      pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
525      nodesAreWritten = true;
526    }
527
528    else if (nvertex == 2)
529    {
530      nbEdges = bounds_lon.shape()[1];
531      edge_lon.resize(nbEdges);
532      edge_lat.resize(nbEdges);
533      edge_lon = lonvalue;
534      edge_lat = latvalue;
535      edge_nodes.resize(nvertex, nbEdges);
536      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo;
537      CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
538
539      // Case (1): node indexes have been generated by domain "nodes"
540      if (nodesAreWritten)
541      {
542        vector<size_t> hashValues(4);
543        CArray<size_t,1> nodeHashList(nbEdges*nvertex*4);
544        for (int ne = 0; ne < nbEdges; ++ne)      // size_t doesn't work with CArray<double, 2>
545        {
546          for (int nv = 0; nv < nvertex; ++nv)
547          {
548            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
549            for (int nh = 0; nh < 4; ++nh)
550            {
551              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
552            }
553          }
554        }
555
556        // Recuperating the node global indexing and writing edge_nodes
557        // Creating map edgeHash2IdxGlo <hash, idxGlo>
558        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
559        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
560        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
561        vector <size_t> edgeNodes;
562        for (int ne = 0; ne < nbEdges; ++ne)
563        {
564          for (int nv = 0; nv < nvertex; ++nv)
565          {
566            int nh = 0;
567            it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
568            // The loop below is needed in case if a hash generated by domain "edges" differs
569            // from that generated by domain "nodes" because of possible precision issues
570            while (it == nodeHash2IdxGlo.end())
571            {
572              ++nh;
573              it = nodeHash2IdxGlo.find(nodeHashList((ne*nvertex + nv)*4 + nh));
574            }
575            edge_nodes(nv,ne) = it->second[0];
576            edgeNodes.push_back(it->second[0]);
577          }
578          edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne);
579        }
580      } // nodesAreWritten
581
582      // Case (2): node indexes have not been generated previously
583      else
584      {
585        // (2.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
586        vector<size_t> hashValues(4);
587        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
588        CArray<size_t,1> nodeHashList(nbEdges*nvertex*4);
589        for (int ne = 0; ne < nbEdges; ++ne)
590        {
591          for (int nv = 0; nv < nvertex; ++nv)
592          {
593            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
594            for (int nh = 0; nh < 4; ++nh)
595            {
596              if (nodeHash2Idx[hashValues[nh]].size() == 0)
597              {
598                nodeHash2Idx[hashValues[nh]].push_back(generateNodeIndex(hashValues, mpiRank));
599                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
600              }
601              nodeHashList((ne*nvertex + nv)*4 + nh) = hashValues[nh];
602            }
603          }
604        }
605
606        // (2.2) Generating global node indexes
607        // The ownership criterion: priority of the process holding the smaller index
608        // Maps generated in this step are:
609        // nodeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
610        // nodeIdx2IdxMin = <idx, idxMin>
611        // nodeIdx2IdxGlo = <idxMin, idxMin>
612
613        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
614        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
615        CArray<size_t,1> nodeIdxMinList(nbEdges*nvertex);
616
617        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
618        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
619        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
620        size_t iIdxMin = 0;
621
622        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
623        {
624          size_t idx = (it->second)[0];
625          size_t idxMin = (it->second)[0];
626          for (int i = 2; i < (it->second).size();)
627          {
628            if (mpiRank == (it->second)[i+1])
629            {
630              idx = (it->second)[i];
631            }
632            if ((it->second)[i] < idxMin)
633            {
634                idxMin = (it->second)[i];
635                (it->second)[i] = (it->second)[i-2];
636            }
637            i += 2;
638          }
639          (it->second)[0] = idxMin;
640          if (nodeIdx2IdxMin.count(idx) == 0)
641          {
642            nodeIdx2IdxMin[idx].push_back(idxMin);
643            if (idx == idxMin)
644              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
645            nodeIdxMinList(iIdxMin) = idxMin;
646            ++iIdxMin;
647          }
648        }
649        nodeIdxMinList.resizeAndPreserve(iIdxMin);
650        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
651        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
652        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
653        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
654        // nodeIdx2IdxGlo holds global indexes only for nodes owned by a process
655        // nodeIdxMin2IdxGlo holds global indexes for all nodes generated by a process
656
657        // (2.3) Saving variables: node_lon, node_lat, edge_nodes
658        // Creating map nodeHash2IdxGlo <hash, idxGlo>
659        // Creating map edgeHash2IdxGlo <hash, idxGlo>
660        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
661        node_count = dhtNodeIdxGlo.getIndexCount();
662        node_start = dhtNodeIdxGlo.getIndexStart();
663        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2IdxGlo;
664        node_lon.resize(node_count);
665        node_lat.resize(node_count);
666        vector <size_t> edgeNodes;
667        size_t idxGlo = 0;
668
669        for (int ne = 0; ne < nbEdges; ++ne)
670        {
671          for (int nv = 0; nv < nvertex; ++nv)
672          {
673            hashValues = CMesh::createHashes(bounds_lon(nv, ne), bounds_lat(nv, ne));
674            size_t myIdx = generateNodeIndex(hashValues, mpiRank);
675            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = nodeIdx2IdxMin.find(myIdx);
676            size_t ownerIdx = (itIdx->second)[0];
677
678            if (myIdx == ownerIdx)
679            {
680              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdx2IdxGlo.find(myIdx);
681              idxGlo = (itIdxGlo->second)[0];
682//              node_lon(idxGlo - node_start) = (bounds_lon(nv, ne) == 360.) ? (0.) : (bounds_lon(nv, ne));
683              node_lon(idxGlo - node_start) = bounds_lon(nv, ne);
684              node_lat(idxGlo - node_start) = bounds_lat(nv, ne);
685            }
686            else
687            {
688              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = nodeIdxMin2IdxGlo.find(ownerIdx);
689              idxGlo = (itIdxGlo->second)[0];
690            }
691
692            edge_nodes(nv,ne) = idxGlo;
693            for (int nh = 0; nh < 4; ++nh)
694              nodeHash2IdxGlo[hashValues[nh]].push_back(idxGlo);
695            edgeNodes.push_back(idxGlo);
696          }
697          edgeHash2IdxGlo[ hashPairOrdered(edgeNodes[0], edgeNodes[1]) ].push_back(nbEdges*mpiRank + ne);
698          edgeNodes.clear();
699        }
700        pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
701      } // !nodesAreWritten
702
703      pEdgeGlobalIndex = new CClientClientDHTSizet (edgeHash2IdxGlo, comm);
704      edgesAreWritten = true;
705    } //nvertex = 2
706
707    else
708    {
709      nbFaces = bounds_lon.shape()[1];
710      face_lon.resize(nbFaces);
711      face_lat.resize(nbFaces);
712      face_lon = lonvalue;
713      face_lat = latvalue;
714      face_nodes.resize(nvertex, nbFaces);
715      face_edges.resize(nvertex, nbFaces);
716
717      // Case (1): edges have been previously generated
718      if (edgesAreWritten)
719      {
720        // (1.1) Recuperating node global indexing and saving face_nodes
721        vector<size_t> hashValues(4);
722        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
723        for (int nf = 0; nf < nbFaces; ++nf)
724        {
725          for (int nv = 0; nv < nvertex; ++nv)
726            {
727              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
728              for (int nh = 0; nh < 4; ++nh)
729                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
730            }
731        }
732        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
733        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
734        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
735        CArray<size_t,1> edgeHashList(nbFaces*nvertex);
736        for (int nf = 0; nf < nbFaces; ++nf)
737        {
738          for (int nv1 = 0; nv1 < nvertex; ++nv1)
739          {
740            int nh1 = 0;
741            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
742            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
743            while (it1 == nodeHash2IdxGlo.end())
744            {
745              ++nh1;
746              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
747            }
748            int nh2 = 0;
749            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
750            while (it2 == nodeHash2IdxGlo.end())
751            {
752              ++nh2;
753              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
754            }
755            face_nodes(nv1,nf) = it1->second[0];
756            edgeHashList(nf*nvertex + nv1) = hashPairOrdered(it1->second[0], it2->second[0]);
757
758          }
759        }
760
761        // (1.2) Recuperating edge global indexing and saving face_edges
762        pEdgeGlobalIndex->computeIndexInfoMapping(edgeHashList);
763        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2IdxGlo = pEdgeGlobalIndex->getInfoIndexMap();
764        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itEdgeHash;
765        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Rank;
766
767        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
768        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex);
769        size_t iIdx = 0;
770
771        for (int nf = 0; nf < nbFaces; ++nf)
772        {
773          for (int nv1 = 0; nv1 < nvertex; ++nv1)
774          {
775            int nh1 = 0;
776            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
777            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
778            while (it1 == nodeHash2IdxGlo.end())
779            {
780              ++nh1;
781              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
782            }
783            int nh2 = 0;
784            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
785            while (it2 == nodeHash2IdxGlo.end())
786            {
787              ++nh2;
788              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
789            }
790            size_t faceIdxGlo = mpiRank * nbFaces + nf;
791            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
792            itEdgeHash = edgeHash2IdxGlo.find(edgeHash);
793            size_t edgeIdxGlo = itEdgeHash->second[0];
794            face_edges(nv1,nf) = edgeIdxGlo;
795            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
796            {
797              edgeIdxGloList(iIdx) = edgeIdxGlo;
798              ++iIdx;
799            }
800            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
801            edgeHash2Rank[edgeHash].push_back(itEdgeHash->first);
802            edgeHash2Rank[edgeHash].push_back(mpiRank);
803          }
804        }
805        edgeIdxGloList.resizeAndPreserve(iIdx);
806
807        // (1.3) Saving remaining variables edge_faces and face_faces
808
809        // Establishing edge ownership
810        // The ownership criterion: priority of the process with smaller rank
811        CClientClientDHTSizet dhtEdgeHash (edgeHash2Rank, comm);
812        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
813        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
814        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHashMin2IdxGlo;            // map is needed purely for counting
815
816        // edgeHash2Info = <edgeHash, <idxGlo, rank1, idxGlo, rank2>>
817        // edgeHashMin2IdxGlo edgeHash2Info = <edgeHashMin, idxGlo>
818        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
819        {
820          vector <size_t> edgeInfo = it->second;
821          if (edgeInfo.size() == 4)                       // two processes generate the same edge
822            if (edgeInfo[1] > edgeInfo[3])
823              edgeInfo[1] = edgeInfo[3];
824          if (edgeInfo[1] == mpiRank)
825            edgeHashMin2IdxGlo[it->first].push_back(edgeInfo[0]);
826        }
827
828        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeHashMin2IdxGlo, comm);
829        edge_count = dhtEdgeIdxGlo.getIndexCount();
830        edge_start = dhtEdgeIdxGlo.getIndexStart();
831
832        edge_faces.resize(2, edge_count);
833        face_faces.resize(nvertex, nbFaces);
834
835        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
836        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
837        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
838        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itFace1, itFace2;
839
840        for (int nf = 0; nf < nbFaces; ++nf)
841        {
842          for (int nv1 = 0; nv1 < nvertex; ++nv1)
843          {
844            int nh1 = 0;
845            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
846            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
847            while (it1 == nodeHash2IdxGlo.end())
848            {
849              ++nh1;
850              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
851            }
852            int nh2 = 0;
853            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
854            while (it2 == nodeHash2IdxGlo.end())
855            {
856              ++nh2;
857              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
858            }
859            size_t faceIdxGlo = mpiRank * nbFaces + nf;
860            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
861            itEdgeHash = edgeHash2Info.find(edgeHash);
862            size_t edgeOwnerRank = itEdgeHash->second[1];
863            int edgeIdxGlo = itEdgeHash->second[0];
864
865            if (mpiRank == edgeOwnerRank)
866            {
867              itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
868              int face1 = itFace1->second[0];
869              if (itFace1->second.size() == 1)
870              {
871                edge_faces(0, edgeIdxGlo - edge_start) = face1;
872                edge_faces(1, edgeIdxGlo - edge_start) = -999;
873                face_faces(nv1, nf) = -1;
874              }
875              else
876              {
877                int face2 = itFace1->second[1];
878                edge_faces(0, edgeIdxGlo - edge_start) = face1;
879                edge_faces(1, edgeIdxGlo - edge_start) = face2;
880                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
881              }
882            }
883            else
884            {
885              itFace1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
886              int face1 = itFace1->second[0];
887              int face2 = itFace1->second[1];
888              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
889            }
890          }
891        }
892      } // edgesAreWritten
893
894      // Case (2): nodes have been previously generated
895      else if (nodesAreWritten)
896      {
897        // (2.1) Generating nodeHashList
898        vector<size_t> hashValues(4);
899        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
900        for (int nf = 0; nf < nbFaces; ++nf)
901        {
902          for (int nv = 0; nv < nvertex; ++nv)
903            {
904              hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
905              for (int nh = 0; nh < 4; ++nh)
906                nodeHashList((nf*nvertex + nv)*4 + nh) = hashValues[nh];
907            }
908        }
909
910        // (2.2) Recuperating node global indexing and saving face_nodes
911        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
912        pNodeGlobalIndex->computeIndexInfoMapping(nodeHashList);
913        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2IdxGlo = pNodeGlobalIndex->getInfoIndexMap();
914        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
915        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
916        CArray<size_t,1> edgeHashList(nbFaces*nvertex);
917        for (int nf = 0; nf < nbFaces; ++nf)
918        {
919          for (int nv1 = 0; nv1 < nvertex; ++nv1)
920          {
921            int nh1 = 0;
922            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
923            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
924            while (it1 == nodeHash2IdxGlo.end())
925            {
926              ++nh1;
927              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
928            }
929            int nh2 = 0;
930            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
931            while (it2 == nodeHash2IdxGlo.end())
932            {
933              ++nh2;
934              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
935            }
936            face_nodes(nv1,nf) = it1->second[0];
937            size_t edgeHash = hashPairOrdered(it1->second[0], it2->second[0]);
938            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(it1->second[0], it2->second[0], mpiRank));
939            edgeHash2Idx[edgeHash].push_back(mpiRank);
940            edgeHashList(nf*nvertex + nv1) = edgeHash;
941          }
942        }
943
944        // (2.2) Generating global edge indexes
945        // The ownership criterion: priority of the process holding the smaller index
946        // Maps generated in this step are:
947        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
948        // edgeIdx2IdxMin = = <idx, idxMin>
949        // edgeIdx2IdxGlo = <idxMin, idxGlo>
950
951        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
952        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
953        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
954
955        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
956        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
957        CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex);
958        size_t iIdxMin = 0;
959
960        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
961        {
962          size_t idxMin = (it->second)[0];
963          size_t idx = (it->second)[0];
964          for (int i = 2; i < (it->second).size();)
965          {
966            if (mpiRank == (it->second)[i+1])
967            {
968              idx = (it->second)[i];
969            }
970            if ((it->second)[i] < idxMin)
971            {
972                idxMin = (it->second)[i];
973                (it->second)[i] = (it->second)[i-2];
974            }
975            i += 2;
976          }
977          (it->second)[0] = idxMin;
978          if (edgeIdx2IdxMin.count(idx) == 0)
979          {
980            edgeIdx2IdxMin[idx].push_back(idxMin);
981            if (idx == idxMin)
982              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
983            edgeIdxMinList(iIdxMin) = idxMin;
984            ++iIdxMin;
985          }
986        }
987        edgeIdxMinList.resizeAndPreserve(iIdxMin);
988        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
989        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
990        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
991        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
992        // edgeIdx2IdxGlo holds global indexes only for edges owned by a process
993        // edgeIdxMin2IdxGlo holds global indexes for all edges generated by a process
994
995        // (2.3) Saving variables: edge_lon, edge_lat, face_edges
996        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
997        edge_count = dhtEdgeIdxGlo.getIndexCount();
998        edge_start = dhtEdgeIdxGlo.getIndexStart();
999        edge_lon.resize(edge_count);
1000        edge_lat.resize(edge_count);
1001        edge_nodes.resize(2, edge_count);
1002        face_edges.resize(nvertex, nbFaces);
1003
1004        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
1005        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex);
1006        size_t iIdx = 0;
1007
1008        for (int nf = 0; nf < nbFaces; ++nf)
1009        {
1010          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1011          {
1012            // Getting global indexes of edge's nodes
1013            int nh1 = 0;
1014            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1015            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1016            while (it1 == nodeHash2IdxGlo.end())
1017            {
1018              ++nh1;
1019              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1020            }
1021            int nh2 = 0;
1022            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1023            while (it2 == nodeHash2IdxGlo.end())
1024            {
1025              ++nh2;
1026              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1027            }
1028            // Getting edge global index
1029            size_t nodeIdxGlo1 = it1->second[0];
1030            size_t nodeIdxGlo2 = it2->second[0];
1031            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
1032            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
1033            size_t ownerIdx = (itIdx->second)[0];
1034            int edgeIdxGlo = 0;
1035            size_t faceIdxGlo = mpiRank * nbFaces + nf;
1036
1037            if (myIdx == ownerIdx)
1038            {
1039              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
1040              edgeIdxGlo = (itIdxGlo->second)[0];
1041              edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ?
1042                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) :
1043                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.);
1044              edge_lat(edgeIdxGlo - edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
1045              edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
1046              edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
1047            }
1048            else
1049            {
1050              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
1051              edgeIdxGlo = (itIdxGlo->second)[0];
1052            }
1053            face_edges(nv1,nf) = edgeIdxGlo;
1054            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
1055            {
1056              edgeIdxGloList(iIdx) = edgeIdxGlo;
1057              ++iIdx;
1058            }
1059            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
1060          }
1061        }
1062        edgeIdxGloList.resizeAndPreserve(iIdx);
1063
1064        // (2.4) Saving remaining variables edge_faces and face_faces
1065        edge_faces.resize(2, edge_count);
1066        face_faces.resize(nvertex, nbFaces);
1067
1068        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
1069        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
1070        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
1071        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
1072
1073        for (int nf = 0; nf < nbFaces; ++nf)
1074        {
1075          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1076          {
1077            // Getting global indexes of edge's nodes
1078            int nh1 = 0;
1079            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1080            it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1081            while (it1 == nodeHash2IdxGlo.end())
1082            {
1083              ++nh1;
1084              it1 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh1));
1085            }
1086            int nh2 = 0;
1087            it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv2)*4 + nh2));
1088            while (it2 == nodeHash2IdxGlo.end())
1089            {
1090              ++nh2;
1091              it2 = nodeHash2IdxGlo.find(nodeHashList((nf*nvertex + nv1)*4 + nh2));
1092            }
1093            size_t nodeIdxGlo1 = it1->second[0];
1094            size_t nodeIdxGlo2 = it2->second[0];
1095
1096            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
1097            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
1098            size_t ownerIdx = (itIdx->second)[0];
1099            size_t faceIdxGlo = mpiRank * nbFaces + nf;
1100
1101            if (myIdx == ownerIdx)
1102            {
1103              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
1104              int edgeIdxGlo = (itIdxGlo->second)[0];
1105              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1106              int face1 = it1->second[0];
1107              if (it1->second.size() == 1)
1108              {
1109                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1110                edge_faces(1, edgeIdxGlo - edge_start) = -999;
1111                face_faces(nv1, nf) = -1;
1112              }
1113              else
1114              {
1115                int face2 = it1->second[1];
1116                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1117                edge_faces(1, edgeIdxGlo - edge_start) = face2;
1118                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1119              }
1120            }
1121            else
1122            {
1123              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
1124              size_t edgeIdxGlo = (itIdxGlo->second)[0];
1125              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1126              int face1 = it1->second[0];
1127              int face2 = it1->second[1];
1128              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1129            }
1130          }
1131        }
1132      } // nodesAreWritten
1133
1134      // Case (3): Neither nodes nor edges have been previously generated
1135      else
1136      {
1137        // (3.1) Creating a list of hashes for each node and a map nodeHash2Idx <hash, <idx,rank> >
1138        vector<size_t> hashValues(4);
1139        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeHash2Idx;
1140        CArray<size_t,1> nodeHashList(nbFaces*nvertex*4);
1141        size_t iHash = 0;
1142        for (int nf = 0; nf < nbFaces; ++nf)
1143        {
1144          for (int nv = 0; nv < nvertex; ++nv)
1145          {
1146            hashValues = CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv, nf));
1147            size_t nodeIndex = generateNodeIndex(hashValues, mpiRank);
1148            for (int nh = 0; nh < 4; ++nh)
1149            {
1150              if (nodeHash2Idx.count(hashValues[nh])==0)
1151              {
1152                nodeHash2Idx[hashValues[nh]].push_back(nodeIndex);
1153                nodeHash2Idx[hashValues[nh]].push_back(mpiRank);
1154                nodeHashList(iHash) = hashValues[nh];
1155                ++iHash;
1156              }
1157            }
1158          }
1159        }
1160        nodeHashList.resizeAndPreserve(iHash);
1161
1162        // (3.2) Generating global node indexes
1163        // The ownership criterion: priority of the process holding the smaller index
1164        // Maps generated in this step are:
1165        // nodeHash2Info = <hash, [[idxMin, rankMin], [idx1, rank1], [idx3, rank3]..]>
1166        // nodeIdx2IdxMin = = <idx, idxMin>
1167        // nodeIdx2IdxGlo = <idxMin, idxGlo>
1168
1169        CClientClientDHTSizet dhtNodeHash(nodeHash2Idx, comm);
1170        dhtNodeHash.computeIndexInfoMapping(nodeHashList);
1171        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeHash2Info = dhtNodeHash.getInfoIndexMap();
1172
1173        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxMin;
1174        CClientClientDHTSizet::Index2VectorInfoTypeMap nodeIdx2IdxGlo;
1175        CArray<size_t,1> nodeIdxMinList(nbFaces*nvertex*4);
1176        size_t iIdxMin = 0;
1177
1178        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = nodeHash2Info.begin(); it != nodeHash2Info.end(); ++it)
1179        {
1180          size_t idxMin = (it->second)[0];
1181          size_t idx = (it->second)[0];
1182          for (int i = 2; i < (it->second).size();)
1183          {
1184            if (mpiRank == (it->second)[i+1])
1185            {
1186              idx = (it->second)[i];
1187            }
1188            if ((it->second)[i] < idxMin)
1189            {
1190                idxMin = (it->second)[i];
1191                (it->second)[i] = (it->second)[i-2];
1192            }
1193            i += 2;
1194          }
1195          (it->second)[0] = idxMin;
1196          if (nodeIdx2IdxMin.count(idx) == 0)
1197          {
1198            nodeIdx2IdxMin[idx].push_back(idxMin);
1199            if (idx == idxMin)
1200              nodeIdx2IdxGlo[idxMin].push_back(idxMin);
1201            nodeIdxMinList(iIdxMin) = idxMin;
1202            ++iIdxMin;
1203          }
1204        }
1205
1206        nodeIdxMinList.resizeAndPreserve(iIdxMin);
1207        CDHTAutoIndexing dhtNodeIdxGlo = CDHTAutoIndexing(nodeIdx2IdxGlo, comm);
1208        CClientClientDHTSizet dhtNodeIdx(nodeIdx2IdxGlo, comm);
1209        dhtNodeIdx.computeIndexInfoMapping(nodeIdxMinList);
1210        CClientClientDHTSizet::Index2VectorInfoTypeMap& nodeIdxMin2IdxGlo = dhtNodeIdx.getInfoIndexMap();
1211
1212        // (3.3) Saving node data: node_lon, node_lat, and face_nodes
1213        // Generating edgeHash2Info = <hash, <idx, rank>> and edgeHashList
1214        nbNodesGlo = dhtNodeIdxGlo.getNbIndexesGlobal();
1215        node_count = dhtNodeIdxGlo.getIndexCount();
1216        node_start = dhtNodeIdxGlo.getIndexStart();
1217        node_lon.resize(node_count);
1218        node_lat.resize(node_count);
1219        size_t nodeIdxGlo1 = 0;
1220        size_t nodeIdxGlo2 = 0;
1221        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2Idx;
1222        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdxGlo1, itNodeIdxGlo2;
1223        CArray<size_t,1> edgeHashList(nbFaces*nvertex);
1224
1225        for (int nf = 0; nf < nbFaces; ++nf)
1226        {
1227          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1228          {
1229            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1230            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1231            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1232            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1233            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1234            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx1 = nodeIdx2IdxMin.find(myNodeIdx1);
1235            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itNodeIdx2 = nodeIdx2IdxMin.find(myNodeIdx2);
1236            size_t ownerNodeIdx = (itNodeIdx1->second)[0];
1237
1238            if (myNodeIdx1 == ownerNodeIdx)
1239            {
1240              itNodeIdxGlo1 = nodeIdx2IdxGlo.find(myNodeIdx1);
1241              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
1242              node_lon(nodeIdxGlo1 - node_start) = bounds_lon(nv1, nf);
1243              node_lat(nodeIdxGlo1 - node_start) = bounds_lat(nv1, nf);
1244            }
1245            else
1246            {
1247              itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find(ownerNodeIdx);
1248              nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
1249            }
1250            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((itNodeIdx2->second)[0]);
1251            nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
1252            size_t edgeHash = hashPairOrdered(nodeIdxGlo1, nodeIdxGlo2);
1253            edgeHash2Idx[edgeHash].push_back(generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank));
1254            edgeHash2Idx[edgeHash].push_back(mpiRank);
1255            edgeHashList(nf*nvertex + nv1) = edgeHash;
1256            face_nodes(nv1,nf) = nodeIdxGlo1;
1257          }
1258        }
1259
1260        // (3.4) Generating global edge indexes
1261        // Maps generated in this step are:
1262        // edgeIdx2IdxMin = = <idx, idxMin>
1263        // edgeIdx2IdxGlo = <idxMin, idxGlo>
1264
1265        CClientClientDHTSizet dhtEdgeHash(edgeHash2Idx, comm);
1266        dhtEdgeHash.computeIndexInfoMapping(edgeHashList);
1267        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeHash2Info = dhtEdgeHash.getInfoIndexMap();
1268        // edgeHash2Info = <hash, [[idx1, rank1], [idx2, rank2], [idx3, rank3]..]>
1269
1270        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxMin;
1271        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdx2IdxGlo;
1272        CArray<size_t,1> edgeIdxMinList(nbFaces*nvertex);
1273        iIdxMin = 0;
1274
1275        for (CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it = edgeHash2Info.begin(); it != edgeHash2Info.end(); ++it)
1276        {
1277          size_t idxMin = (it->second)[0];
1278          size_t idx = (it->second)[0];
1279
1280          for (int i = 2; i < (it->second).size();)
1281          {
1282            if (mpiRank == (it->second)[i+1])
1283            {
1284              idx = (it->second)[i];
1285            }
1286            if ((it->second)[i] < idxMin)
1287            {
1288                idxMin = (it->second)[i];
1289                (it->second)[i] = (it->second)[i-2];
1290            }
1291            i += 2;
1292          }
1293          (it->second)[0] = idxMin;
1294          if (edgeIdx2IdxMin.count(idx) == 0)
1295          {
1296            edgeIdx2IdxMin[idx].push_back(idxMin);
1297            if (idx == idxMin)
1298              edgeIdx2IdxGlo[idxMin].push_back(idxMin);
1299            edgeIdxMinList(iIdxMin) = idxMin;
1300            ++iIdxMin;
1301          }
1302        }
1303        edgeIdxMinList.resizeAndPreserve(iIdxMin);
1304        CDHTAutoIndexing dhtEdgeIdxGlo = CDHTAutoIndexing(edgeIdx2IdxGlo, comm);
1305        CClientClientDHTSizet dhtEdgeIdx(edgeIdx2IdxGlo, comm);
1306        dhtEdgeIdx.computeIndexInfoMapping(edgeIdxMinList);
1307        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxMin2IdxGlo = dhtEdgeIdx.getInfoIndexMap();
1308
1309        // (3.5) Saving variables: edge_lon, edge_lat, face_edges
1310        // Creating map edgeIdxGlo2Face <idxGlo, face>
1311
1312        nbEdgesGlo = dhtEdgeIdxGlo.getNbIndexesGlobal();
1313        edge_count = dhtEdgeIdxGlo.getIndexCount();
1314        edge_start = dhtEdgeIdxGlo.getIndexStart();
1315        edge_lon.resize(edge_count);
1316        edge_lat.resize(edge_count);
1317        edge_nodes.resize(2, edge_count);
1318        face_edges.resize(nvertex, nbFaces);
1319
1320        CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it1, it2;
1321        CClientClientDHTSizet::Index2VectorInfoTypeMap edgeIdxGlo2Face;
1322        CArray<size_t,1> edgeIdxGloList(nbFaces*nvertex);
1323        size_t iIdx = 0;
1324
1325        for (int nf = 0; nf < nbFaces; ++nf)
1326        {
1327          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1328          {
1329            // Getting global indexes of edge's nodes
1330            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1331            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1332            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1333
1334            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1335            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1336            it1 = nodeIdx2IdxMin.find(myNodeIdx1);
1337            it2 = nodeIdx2IdxMin.find(myNodeIdx2);
1338            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
1339            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
1340            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
1341            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
1342
1343            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
1344            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
1345            size_t ownerIdx = (itIdx->second)[0];
1346
1347            int edgeIdxGlo = 0;
1348            size_t faceIdxGlo = mpiRank * nbFaces + nf;
1349
1350            if (myIdx == ownerIdx)
1351            {
1352              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
1353              edgeIdxGlo = (itIdxGlo->second)[0];
1354              edge_lon(edgeIdxGlo - edge_start) = ( abs( bounds_lon(nv1, nf) - bounds_lon(nv2, nf)) < 180.) ?
1355                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5) :
1356                          (( bounds_lon(nv1, nf) + bounds_lon(nv2, nf)) * 0.5 -180.);
1357              edge_lat(edgeIdxGlo-edge_start) = ( bounds_lat(nv1, nf) + bounds_lat(nv2, nf) ) * 0.5;
1358              edge_nodes(0, edgeIdxGlo - edge_start) = nodeIdxGlo1;
1359              edge_nodes(1, edgeIdxGlo - edge_start) = nodeIdxGlo2;
1360            }
1361            else
1362            {
1363              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
1364              edgeIdxGlo = (itIdxGlo->second)[0];
1365            }
1366            face_edges(nv1,nf) = edgeIdxGlo;
1367            if (edgeIdxGlo2Face.count(edgeIdxGlo) == 0)
1368            {
1369              edgeIdxGloList(iIdx) = edgeIdxGlo;
1370              ++iIdx;
1371            }
1372            edgeIdxGlo2Face[edgeIdxGlo].push_back(faceIdxGlo);
1373          }
1374        }
1375        edgeIdxGloList.resizeAndPreserve(iIdx);
1376
1377        // (3.6) Saving remaining variables edge_faces and face_faces
1378        edge_faces.resize(2, edge_count);
1379        face_faces.resize(nvertex, nbFaces);
1380
1381        CClientClientDHTSizet dhtEdge2Face (edgeIdxGlo2Face, comm);
1382        dhtEdge2Face.computeIndexInfoMapping(edgeIdxGloList);
1383        CClientClientDHTSizet::Index2VectorInfoTypeMap& edgeIdxGlo2FaceIdx = dhtEdge2Face.getInfoIndexMap();
1384
1385        for (int nf = 0; nf < nbFaces; ++nf)
1386        {
1387          for (int nv1 = 0; nv1 < nvertex; ++nv1)
1388          {
1389            // Getting global indexes of edge's nodes
1390            int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
1391            vector<size_t> hashValues1 = CMesh::createHashes(bounds_lon(nv1, nf), bounds_lat(nv1, nf));
1392            vector<size_t> hashValues2 = CMesh::createHashes(bounds_lon(nv2, nf), bounds_lat(nv2, nf));
1393
1394            size_t myNodeIdx1 = generateNodeIndex(hashValues1, mpiRank);
1395            size_t myNodeIdx2 = generateNodeIndex(hashValues2, mpiRank);
1396            it1 = nodeIdx2IdxMin.find(myNodeIdx1);
1397            it2 = nodeIdx2IdxMin.find(myNodeIdx2);
1398            itNodeIdxGlo1 = nodeIdxMin2IdxGlo.find((it1->second)[0]);
1399            itNodeIdxGlo2 = nodeIdxMin2IdxGlo.find((it2->second)[0]);
1400            size_t nodeIdxGlo1 = (itNodeIdxGlo1->second)[0];
1401            size_t nodeIdxGlo2 = (itNodeIdxGlo2->second)[0];
1402
1403            size_t myIdx = generateEdgeIndex(nodeIdxGlo1, nodeIdxGlo2, mpiRank);
1404            CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdx = edgeIdx2IdxMin.find(myIdx);
1405            size_t ownerIdx = (itIdx->second)[0];
1406            size_t faceIdxGlo = mpiRank * nbFaces + nf;
1407
1408            if (myIdx == ownerIdx)
1409            {
1410              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdx2IdxGlo.find(myIdx);
1411              int edgeIdxGlo = (itIdxGlo->second)[0];
1412              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1413              int face1 = it1->second[0];
1414              if (it1->second.size() == 1)
1415              {
1416                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1417                edge_faces(1, edgeIdxGlo - edge_start) = -999;
1418                face_faces(nv1, nf) = -1;
1419              }
1420              else
1421              {
1422                size_t face2 = it1->second[1];
1423                edge_faces(0, edgeIdxGlo - edge_start) = face1;
1424                edge_faces(1, edgeIdxGlo - edge_start) = face2;
1425                face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1426              }
1427            }
1428            else
1429            {
1430              CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator itIdxGlo = edgeIdxMin2IdxGlo.find(ownerIdx);
1431              size_t edgeIdxGlo = (itIdxGlo->second)[0];
1432              it1 = edgeIdxGlo2FaceIdx.find(edgeIdxGlo);
1433              int face1 = it1->second[0];
1434              int face2 = it1->second[1];
1435              face_faces(nv1, nf) = (faceIdxGlo == face1 ? face2 : face1);
1436            }
1437          }
1438        }
1439      }
1440      facesAreWritten = true;
1441    } // nvertex >= 3
1442
1443  } // createMeshEpsilon
1444
1445} // namespace xios
Note: See TracBrowser for help on using the repository browser.