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

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

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