[879] | 1 | /*! |
---|
| 2 | \file mesh.cpp |
---|
| 3 | \author Olga Abramkina |
---|
| 4 | \brief Definition of class CMesh. |
---|
| 5 | */ |
---|
| 6 | |
---|
| 7 | #include "mesh.hpp" |
---|
| 8 | |
---|
| 9 | namespace 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 |
---|