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

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

Corrected submit for UGRID.
File size: 16.4 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) :  nbFaces{0}, nbNodes{0}, nbEdges{0}
14            ,   nvertex{0}
15            ,  nodesAreWritten{false}, edgesAreWritten{false}, facesAreWritten{false}
16            ,  node_lon(), node_lat()
17            ,  edge_lon(), edge_lat(), edge_nodes()
18            ,  face_lon(), face_lat()
19            ,  face_nodes()
20  {
21  }
22
23
24  CMesh::~CMesh(void)
25  {
26  }
27
28  std::map <StdString, CMesh*> CMesh::meshList = std::map <StdString, CMesh*>();
29  CMesh* CMesh::getMesh;
30
31///---------------------------------------------------------------
32/*!
33 * \fn bool CMesh::isWritten (StdString meshName)
34 * Checks if a mesh has been written, updates the list of meshes stored in meshList
35 * \param [in] meshName  The name of a mesh ("name" attribute of a domain).
36 */
37  bool CMesh::isWritten (StdString meshName)
38  {
39    if ( CMesh::meshList.begin() != CMesh::meshList.end() )
40    {
41      if ( CMesh::meshList.find(meshName) != CMesh::meshList.end() )
42      {
43        CMesh::getMesh = meshList[meshName];
44        return true;
45      }
46      else
47      {
48        CMesh::meshList.insert( make_pair(meshName, this) );
49        return false;
50      }
51    }
52    else
53    {
54      CMesh::meshList.insert( make_pair(meshName, this) );
55      return false;
56    }
57  }
58
59///----------------------------------------------------------------
60  int hashPair(int first, int second)
61  {
62    int seed = first + 0x9e3779b9 ;
63    seed ^= second + 0x9e3779b9 + (seed << 6) + (seed >> 2);
64    return seed ;
65  }
66
67///----------------------------------------------------------------
68/*!
69 * \fn size_t CMesh::nodeIndex (double lon, double lat)
70 * Returns its index if a node exists; otherwise adds the node and returns -1.
71 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
72 * \param [in] lon Node longitude in degrees.
73 * \param [in] lat Node latitude in degress ranged from 0 to 360.
74 * \return node index if a node exists; -1 otherwise
75 */
76  size_t CMesh::nodeIndex (double lon, double lat)
77  {
78    double minBoundLon = 0. ;
79    double maxBoundLon = 360. ;
80    double minBoundLat = -90 ;
81    double maxBoundLat = 90 ;
82    double prec=1e-11 ;
83    double precLon=prec ;
84    double precLat=prec ;
85
86    size_t maxsize_t=numeric_limits<size_t>::max() ;
87    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
88    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
89
90    size_t iMinLon=0 ;
91    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
92    size_t iMinLat=0 ;
93    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
94
95    size_t hash0,hash1,hash2,hash3 ;
96    size_t lon0,lon1,lat0,lat1 ;
97
98    lon0=(lon-minBoundLon)/precLon ;
99    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
100    {
101      if (lon0==iMinLon) lon1=iMaxLon ;
102      else lon1=lon0-1 ;
103    }
104    else
105    {
106      if (lon0==iMaxLon) lon1=iMinLon ;
107      else lon1=lon0+1 ;
108    }
109
110    lat0=(lat-minBoundLat)/precLat ;
111    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
112    {
113      if (lat0==iMinLat) lat1=lat0 ;
114      else lat1=lat0-1 ;
115    }
116    else
117    {
118      if (lat0==iMaxLat) lat1=lat0 ;
119      else lat1=lat0+1 ;
120    }
121
122    hash0=hashPair(lon0,lat0) ;
123    hash1=hashPair(lon0,lat1) ;
124    hash2=hashPair(lon1,lat0) ;
125    hash3=hashPair(lon1,lat1) ;
126
127    boost::unordered_map<size_t, size_t>::iterator end = hashed_map_nodes.end() ;
128    size_t mapSize = hashed_map_nodes.size();
129    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)
130    {
131      hashed_map_nodes[hash0] = mapSize ;
132      hashed_map_nodes[hash1] = mapSize + 1;
133      hashed_map_nodes[hash2] = mapSize + 2;
134      hashed_map_nodes[hash3] = mapSize + 3;
135      return -1;
136    }
137    else
138      return ( (hashed_map_nodes[hash0]+1) / 4 );
139
140  } // nodeIndex()
141
142///----------------------------------------------------------------
143/*!
144 * \fn void CMesh::hashDblDbl (double lon, double lat)
145 * Creates two hash values for each dimension, longitude and latitude.
146 * \param [in] lon Node longitude in degrees.
147 * \param [in] lat Node latitude in degress ranged from 0 to 360.
148 */
149  void hashDblDbl (double lon, double lat)
150  {
151    double minBoundLon = 0. ;
152    double maxBoundLon = 360. ;
153    double minBoundLat = -90 ;
154    double maxBoundLat = 90 ;
155    double prec=1e-11 ;
156    double precLon=prec ;
157    double precLat=prec ;
158
159    size_t maxsize_t=numeric_limits<size_t>::max() ;
160    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
161    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
162
163    size_t iMinLon=0 ;
164    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
165    size_t iMinLat=0 ;
166    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
167
168    size_t hash0,hash1,hash2,hash3 ;
169    size_t lon0,lon1,lat0,lat1 ;
170
171    lon0=(lon-minBoundLon)/precLon ;
172    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
173    {
174      if (lon0==iMinLon) lon1=iMaxLon ;
175      else lon1=lon0-1 ;
176    }
177    else
178    {
179      if (lon0==iMaxLon) lon1=iMinLon ;
180      else lon1=lon0+1 ;
181    }
182
183    lat0=(lat-minBoundLat)/precLat ;
184    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
185    {
186      if (lat0==iMinLat) lat1=lat0 ;
187      else lat1=lat0-1 ;
188    }
189    else
190    {
191      if (lat0==iMaxLat) lat1=lat0 ;
192      else lat1=lat0+1 ;
193    }
194
195    hash0=hashPair(lon0,lat0) ;
196    hash1=hashPair(lon0,lat1) ;
197    hash2=hashPair(lon1,lat0) ;
198    hash3=hashPair(lon1,lat1) ;
199
200  } // hashDblDbl
201
202///----------------------------------------------------------------
203  std::pair<int,int> make_ordered_pair(int a, int b)
204  {
205    if ( a < b )
206      return std::pair<int,int>(a,b);
207    else
208      return std::pair<int,int>(b,a);
209  }
210
211///----------------------------------------------------------------
212/*!
213 * \fn void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
214            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
215 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
216 * \param [in] lonvalue  Array of longitudes.
217 * \param [in] latvalue  Array of latitudes.
218 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
219 * \param [in] bounds_lat Array of boundarry latitudes. Its size depends on the element type.
220 */
221  void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
222            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
223  {
224    nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
225
226    if (nvertex == 1)
227    {
228      nbNodes = lonvalue.numElements();
229      node_lon.resizeAndPreserve(nbNodes);
230      node_lat.resizeAndPreserve(nbNodes);
231      for (int nn = 0; nn < nbNodes; ++nn)
232      {
233        if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end())
234        {
235          map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ;
236          node_lon(nn) = lonvalue(nn);
237          node_lat(nn) = latvalue(nn);
238        }
239      }
240    }
241    else if (nvertex == 2)
242    {
243      nbEdges = bounds_lon.shape()[1];
244
245      // Create nodes and edge_node connectivity
246      node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes
247      node_lat.resizeAndPreserve(nbEdges*nvertex);
248      edge_nodes.resizeAndPreserve(nvertex, nbEdges);
249
250      for (int ne = 0; ne < nbEdges; ++ne)
251      {
252        for (int nv = 0; nv < nvertex; ++nv)
253        {
254          if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end())
255          {
256            map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes ;
257            edge_nodes(nv,ne) = nbNodes ;
258            node_lon(nbNodes) = bounds_lon(nv, ne);
259            node_lat(nbNodes) = bounds_lat(nv, ne);
260            ++nbNodes ;
261          }
262          else
263            edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))];
264        }
265      }
266      node_lon.resizeAndPreserve(nbNodes);
267      node_lat.resizeAndPreserve(nbNodes);
268
269      // Create edges
270      edge_lon.resizeAndPreserve(nbEdges);
271      edge_lat.resizeAndPreserve(nbEdges);
272
273      for (int ne = 0; ne < nbEdges; ++ne)
274      {
275        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
276        {
277          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
278          edge_lon(ne) = lonvalue(ne);
279          edge_lat(ne) = latvalue(ne);
280        }
281
282      }
283      edgesAreWritten = true;
284    }
285    else
286    {
287      nbFaces = bounds_lon.shape()[1];
288 
289      // Create nodes and face_node connectivity
290      node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes
291      node_lat.resizeAndPreserve(nbFaces*nvertex);
292      face_nodes.resize(nvertex, nbFaces);
293 
294      for (int nf = 0; nf < nbFaces; ++nf)
295      {
296        for (int nv = 0; nv < nvertex; ++nv)
297        {
298          if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end())
299          {
300            map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes ;
301            face_nodes(nv,nf) = nbNodes ;
302            node_lon(nbNodes) = bounds_lon(nv, nf);
303            node_lat(nbNodes) = bounds_lat(nv ,nf);
304            ++nbNodes ;
305          }
306          else
307          {
308            face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
309          }
310        }
311      }
312      node_lon.resizeAndPreserve(nbNodes);
313      node_lat.resizeAndPreserve(nbNodes);
314 
315      // Create edges and edge_nodes connectivity
316      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges
317      edge_lat.resizeAndPreserve(nbFaces*nvertex);
318      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex);
319      for (int nf = 0; nf < nbFaces; ++nf)
320      {
321        for (int nv1 = 0; nv1 < nvertex; ++nv1)
322        {
323          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
324          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
325          {
326            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ;
327            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf);
328            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
329                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
330                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);;
331            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
332            ++nbEdges;
333          }
334        }
335      }
336      edge_nodes.resizeAndPreserve(2, nbEdges);
337      edge_lon.resizeAndPreserve(nbEdges);
338      edge_lat.resizeAndPreserve(nbEdges);
339 
340      // Create faces
341      face_lon.resize(nbFaces);
342      face_lat.resize(nbFaces);
343      face_lon = lonvalue;
344      face_lat = latvalue;
345      facesAreWritten = true;
346    } // nvertex > 2
347   
348  } // createMesh()
349
350///----------------------------------------------------------------
351/*!
352 * \fn void CMesh::createMeshEpsilon(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
353            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
354 * Creates or updates a mesh for the three types of mesh elements: nodes, edges, and faces.
355 * Precision check is implemented with two hash values for each dimension, longitude and latitude.
356 * \param [in] lonvalue  Array of longitudes.
357 * \param [in] latvalue  Array of latitudes.
358 * \param [in] bounds_lon Array of boundary longitudes. Its size depends on the element type.
359 * \param [in] bounds_lat Array of boundarry latitudes. Its size depends on the element type.
360 */
361  void CMesh::createMeshEpsilon(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
362            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
363  {
364    nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
365
366    if (nvertex == 1)
367    {
368      nbNodes = lonvalue.numElements();
369      node_lon.resizeAndPreserve(nbNodes);
370      node_lat.resizeAndPreserve(nbNodes);
371      for (int nn = 0; nn < nbNodes; ++nn)
372      {
373        if ( nodeIndex(lonvalue(nn), latvalue(nn)) == -1 )
374        {
375          node_lon(nn) = lonvalue(nn);
376          node_lat(nn) = latvalue(nn);
377        }
378      }
379    nodesAreWritten = true;
380    }
381    else if (nvertex == 2)
382    {
383      nbEdges = bounds_lon.shape()[1];
384
385      // Create nodes and edge_node connectivity
386      node_lon.resizeAndPreserve(nbEdges*nvertex); // Max possible number of nodes
387      node_lat.resizeAndPreserve(nbEdges*nvertex);
388      edge_nodes.resizeAndPreserve(nvertex, nbEdges);
389
390      for (int ne = 0; ne < nbEdges; ++ne)
391      {
392        for (int nv = 0; nv < nvertex; ++nv)
393        {
394          if ( nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne)) == -1)
395          {
396            edge_nodes(nv,ne) = nbNodes ;
397            node_lon(nbNodes) = bounds_lon(nv, ne);
398            node_lat(nbNodes) = bounds_lat(nv, ne);
399            ++nbNodes ;
400          }
401          else
402            edge_nodes(nv,ne) = nodeIndex(bounds_lon(nv, ne), bounds_lat(nv ,ne));
403        }
404      }
405      node_lon.resizeAndPreserve(nbNodes);
406      node_lat.resizeAndPreserve(nbNodes);
407
408      // Create edges
409      edge_lon.resizeAndPreserve(nbEdges);
410      edge_lat.resizeAndPreserve(nbEdges);
411
412      for (int ne = 0; ne < nbEdges; ++ne)
413      {
414        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
415        {
416          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
417          edge_lon(ne) = lonvalue(ne);
418          edge_lat(ne) = latvalue(ne);
419        }
420      }
421      edgesAreWritten = true;
422    } // nvertex = 2
423    else
424    {
425      nbFaces = bounds_lon.shape()[1];
426
427      // Create nodes and face_node connectivity
428      node_lon.resizeAndPreserve(nbFaces*nvertex);  // Max possible number of nodes
429      node_lat.resizeAndPreserve(nbFaces*nvertex);
430      face_nodes.resize(nvertex, nbFaces);
431
432      for (int nf = 0; nf < nbFaces; ++nf)
433      {
434        for (int nv = 0; nv < nvertex; ++nv)
435        {
436          if ( nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf)) == -1)
437          {
438            face_nodes(nv,nf) = nbNodes ;
439            node_lon(nbNodes) = bounds_lon(nv, nf);
440            node_lat(nbNodes) = bounds_lat(nv ,nf);
441            ++nbNodes ;
442          }
443          else
444          {
445            face_nodes(nv,nf) = nodeIndex(bounds_lon(nv, nf), bounds_lat(nv ,nf));
446          }
447        }
448      }
449      node_lon.resizeAndPreserve(nbNodes);
450      node_lat.resizeAndPreserve(nbNodes);
451
452      // Create edges and edge_nodes connectivity
453      edge_lon.resizeAndPreserve(nbFaces*nvertex); // Max possible number of edges
454      edge_lat.resizeAndPreserve(nbFaces*nvertex);
455      edge_nodes.resizeAndPreserve(2, nbFaces*nvertex);
456      for (int nf = 0; nf < nbFaces; ++nf)
457      {
458        for (int nv1 = 0; nv1 < nvertex; ++nv1)
459        {
460          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
461          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
462          {
463            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges ;
464            edge_nodes(Range::all(),nbEdges) = face_nodes(nv1,nf), face_nodes(nv2,nf);
465            edge_lon(nbEdges) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
466                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
467                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);;
468            edge_lat(nbEdges) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
469            ++nbEdges;
470          }
471        }
472      }
473      edge_nodes.resizeAndPreserve(2, nbEdges);
474      edge_lon.resizeAndPreserve(nbEdges);
475      edge_lat.resizeAndPreserve(nbEdges);
476
477      // Create faces
478      face_lon.resize(nbFaces);
479      face_lat.resize(nbFaces);
480      face_lon = lonvalue;
481      face_lat = latvalue;
482      facesAreWritten = true;
483    } // nvertex >= 3
484
485  } // createMeshEpsilon
486
487} // namespace xios
Note: See TracBrowser for help on using the repository browser.