Ignore:
Timestamp:
06/18/18 20:32:55 (6 years ago)
Author:
yushan
Message:

branch_openmp merged with trunk r1544

File:
1 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/branch_openmp/src/node/mesh.cpp

    r1328 r1545  
    77#include "mesh.hpp" 
    88using namespace ep_lib; 
     9#include <boost/functional/hash.hpp> 
    910 
    1011namespace xios { 
    … …  
    3233  } 
    3334 
    34   //std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>(); 
    35   //std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >(); 
    36  
    3735  std::map <StdString, CMesh> *CMesh::meshList_ptr = 0; 
    3836  std::map <StdString, vector<int> > *CMesh::domainList_ptr = 0; 
    … …  
    5048    if(CMesh::meshList_ptr == NULL)   CMesh::meshList_ptr   = new std::map <StdString, CMesh>(); 
    5149 
    52     //CMesh::domainList[meshName].push_back(nvertex); 
    5350    CMesh::domainList_ptr->at(meshName).push_back(nvertex); 
    5451 
    55     //if ( CMesh::meshList.begin() != CMesh::meshList.end() ) 
    5652    if ( CMesh::meshList_ptr->begin() != CMesh::meshList_ptr->end() ) 
    5753    { 
    58       //for (std::map<StdString, CMesh>::iterator it=CMesh::meshList.begin(); it!=CMesh::meshList.end(); ++it) 
    5954      for (std::map<StdString, CMesh>::iterator it=CMesh::meshList_ptr->begin(); it!=CMesh::meshList_ptr->end(); ++it) 
    6055      { 
    6156        if (it->first == meshName) 
    62           //return &meshList[meshName]; 
    6357          return &meshList_ptr->at(meshName); 
    6458        else 
    6559        { 
    6660          CMesh newMesh; 
    67           //CMesh::meshList.insert( make_pair(meshName, newMesh) ); 
    6861          CMesh::meshList_ptr->insert( make_pair(meshName, newMesh) ); 
    69           //return &meshList[meshName]; 
    7062          return &meshList_ptr->at(meshName); 
    7163        } 
    … …  
    7567    { 
    7668      CMesh newMesh; 
    77       //CMesh::meshList.insert( make_pair(meshName, newMesh) ); 
    7869      CMesh::meshList_ptr->insert( make_pair(meshName, newMesh) ); 
    79       //return &meshList[meshName]; 
    8070      return &meshList_ptr->at(meshName); 
    8171    } 
    … …  
    151141    } 
    152142 
    153  
    154 ///---------------------------------------------------------------- 
    155 /*! 
    156  * \fn size_t CMesh::nodeIndex (double lon, double lat) 
    157  * Returns its index if a node exists; otherwise adds the node and returns -1. 
    158  * Precision check is implemented with two hash values for each dimension, longitude and latitude. 
    159  * \param [in] lon Node longitude in degrees. 
    160  * \param [in] lat Node latitude in degrees ranged from 0 to 360. 
    161  * \return node index if a node exists; -1 otherwise 
    162  */ 
    163   size_t CMesh::nodeIndex (double lon, double lat) 
    164   { 
    165     double minBoundLon = 0. ; 
    166     double maxBoundLon = 360. ; 
    167     double minBoundLat = -90 ; 
    168     double maxBoundLat = 90 ; 
    169     double prec=1e-11 ; 
    170     double precLon=prec ; 
    171     double precLat=prec ; 
    172  
    173     size_t maxsize_t=numeric_limits<size_t>::max() ; 
    174     if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ; 
    175     if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ; 
    176  
    177     size_t iMinLon=0 ; 
    178     size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ; 
    179     size_t iMinLat=0 ; 
    180     size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ; 
    181  
    182     size_t hash0,hash1,hash2,hash3 ; 
    183     size_t lon0,lon1,lat0,lat1 ; 
    184  
    185     lon0=(lon-minBoundLon)/precLon ; 
    186     if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon) 
    187     { 
    188       if (lon0==iMinLon) lon1=iMaxLon ; 
    189       else lon1=lon0-1 ; 
    190     } 
    191     else 
    192     { 
    193       if (lon0==iMaxLon) lon1=iMinLon ; 
    194       else lon1=lon0+1 ; 
    195     } 
    196  
    197     lat0=(lat-minBoundLat)/precLat ; 
    198     if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat) 
    199     { 
    200       if (lat0==iMinLat) lat1=lat0 ; 
    201       else lat1=lat0-1 ; 
    202     } 
    203     else 
    204     { 
    205       if (lat0==iMaxLat) lat1=lat0 ; 
    206       else lat1=lat0+1 ; 
    207     } 
    208  
    209     hash0=hashPair(lon0,lat0) ; 
    210     hash1=hashPair(lon0,lat1) ; 
    211     hash2=hashPair(lon1,lat0) ; 
    212     hash3=hashPair(lon1,lat1) ; 
    213  
    214     boost::unordered_map<size_t, size_t>::iterator end = hashed_map_nodes.end() ; 
    215     size_t mapSize = hashed_map_nodes.size(); 
    216     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) 
    217     { 
    218       hashed_map_nodes[hash0] = mapSize ; 
    219       hashed_map_nodes[hash1] = mapSize + 1; 
    220       hashed_map_nodes[hash2] = mapSize + 2; 
    221       hashed_map_nodes[hash3] = mapSize + 3; 
    222       return -1; 
    223     } 
    224     else 
    225       return ( (hashed_map_nodes[hash0]+1) / 4 ); 
    226  
    227   } // nodeIndex() 
    228  
    229143///---------------------------------------------------------------- 
    230144/*! 
    … …  
    313227 * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type. 
    314228 */ 
    315   void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue, 
    316             const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat) 
    317   { 
    318     int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
    319  
    320     if (nvertex == 1) 
    321     { 
    322       nbNodes_ = lonvalue.numElements(); 
    323       node_lon.resizeAndPreserve(nbNodes_); 
    324       node_lat.resizeAndPreserve(nbNodes_); 
    325       for (int nn = 0; nn < nbNodes_; ++nn) 
    326       { 
    327         if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end()) 
    328         { 
    329           map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ; 
    330           node_lon(nn) = lonvalue(nn); 
    331           node_lat(nn) = latvalue(nn); 
    332         } 
    333       } 
    334     } 
    335     else if (nvertex == 2) 
    336     { 
    337       nbEdges_ = bounds_lon.shape()[1]; 
    338  
    339       // Create nodes and edge_node connectivity 
    340       node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes 
    341       node_lat.resizeAndPreserve(nbEdges_*nvertex); 
    342       edge_nodes.resizeAndPreserve(nvertex, nbEdges_); 
    343  
    344       for (int ne = 0; ne < nbEdges_; ++ne) 
    345       { 
    346         for (int nv = 0; nv < nvertex; ++nv) 
    347         { 
    348           if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end()) 
    349           { 
    350             map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ; 
    351             edge_nodes(nv,ne) = nbNodes_ ; 
    352             node_lon(nbNodes_) = bounds_lon(nv, ne); 
    353             node_lat(nbNodes_) = bounds_lat(nv, ne); 
    354             ++nbNodes_ ; 
    355           } 
    356           else 
    357             edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))]; 
    358         } 
    359       } 
    360       node_lon.resizeAndPreserve(nbNodes_); 
    361       node_lat.resizeAndPreserve(nbNodes_); 
    362  
    363       // Create edges 
    364       edge_lon.resizeAndPreserve(nbEdges_); 
    365       edge_lat.resizeAndPreserve(nbEdges_); 
    366  
    367       for (int ne = 0; ne < nbEdges_; ++ne) 
    368       { 
    369         if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end()) 
    370         { 
    371           map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ; 
    372           edge_lon(ne) = lonvalue(ne); 
    373           edge_lat(ne) = latvalue(ne); 
    374         } 
    375  
    376       } 
    377       edgesAreWritten = true; 
    378     } 
    379     else 
    380     { 
    381       nbFaces_ = bounds_lon.shape()[1]; 
    382    
    383       // Create nodes and face_node connectivity 
    384       node_lon.resizeAndPreserve(nbFaces_*nvertex);  // Max possible number of nodes 
    385       node_lat.resizeAndPreserve(nbFaces_*nvertex); 
    386       face_nodes.resize(nvertex, nbFaces_); 
    387    
    388       for (int nf = 0; nf < nbFaces_; ++nf) 
    389       { 
    390         for (int nv = 0; nv < nvertex; ++nv) 
    391         { 
    392           if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end()) 
    393           { 
    394             map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ; 
    395             face_nodes(nv,nf) = nbNodes_ ; 
    396             node_lon(nbNodes_) = bounds_lon(nv, nf); 
    397             node_lat(nbNodes_) = bounds_lat(nv ,nf); 
    398             ++nbNodes_ ; 
    399           } 
    400           else 
    401           { 
    402             face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))]; 
    403           } 
    404         } 
    405       } 
    406       node_lon.resizeAndPreserve(nbNodes_); 
    407       node_lat.resizeAndPreserve(nbNodes_); 
    408    
    409       // Create edges and edge_nodes connectivity 
    410       edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges 
    411       edge_lat.resizeAndPreserve(nbFaces_*nvertex); 
    412       edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex); 
    413       edge_faces.resize(2, nbFaces_*nvertex); 
    414       face_edges.resize(nvertex, nbFaces_); 
    415       face_faces.resize(nvertex, nbFaces_); 
    416  
    417       vector<int> countEdges(nbFaces_*nvertex);   // needed in case if edges have been already generated 
    418       vector<int> countFaces(nbFaces_); 
    419       countEdges.assign(nbFaces_*nvertex, 0); 
    420       countFaces.assign(nbFaces_, 0); 
    421       int edge; 
    422       for (int nf = 0; nf < nbFaces_; ++nf) 
    423       { 
    424         for (int nv1 = 0; nv1 < nvertex; ++nv1) 
    425         { 
    426           int nv = 0; 
    427           int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
    428           if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end()) 
    429           { 
    430             map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ; 
    431             face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))]; 
    432             edge_faces(0,nbEdges_) = nf; 
    433             edge_faces(1,nbEdges_) = -999; 
    434             face_faces(nv1,nf) = 999999; 
    435             edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf); 
    436             edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ? 
    437                         (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) : 
    438                         (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.); 
    439             edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5; 
    440             ++nbEdges_; 
    441           } 
    442           else 
    443           { 
    444             edge = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))]; 
    445             face_edges(nv1,nf) = edge; 
    446             if (edgesAreWritten) 
    447             { 
    448               edge_faces(countEdges[edge], edge) = nf; 
    449               if (countEdges[edge]==0) 
    450               { 
    451                 face_faces(nv1,nf) = 999999; 
    452               } 
    453               else 
    454               { 
    455                 int face1 = nf; // = edge_faces(1,edge) 
    456                 int face2 = edge_faces(0,edge); 
    457                 face_faces(countFaces[face1], face1) =  face2; 
    458                 face_faces(countFaces[face2], face2) =  face1; 
    459                 ++(countFaces[face1]); 
    460                 ++(countFaces[face2]); 
    461               } 
    462             } 
    463             else 
    464             { 
    465               edge_faces(1,edge) = nf; 
    466               int face1 = nf; // = edge_faces(1,edge) 
    467               int face2 = edge_faces(0,edge); 
    468               face_faces(countFaces[face1], face1) =  face2; 
    469               face_faces(countFaces[face2], face2) =  face1; 
    470               ++(countFaces[face1]); 
    471               ++(countFaces[face2]); 
    472             } 
    473             ++(countEdges[edge]); 
    474           } 
    475         } 
    476       } 
    477       edge_nodes.resizeAndPreserve(2, nbEdges_); 
    478       edge_faces.resizeAndPreserve(2, nbEdges_); 
    479       edge_lon.resizeAndPreserve(nbEdges_); 
    480       edge_lat.resizeAndPreserve(nbEdges_); 
    481  
    482       // Create faces 
    483       face_lon.resize(nbFaces_); 
    484       face_lat.resize(nbFaces_); 
    485       face_lon = lonvalue; 
    486       face_lat = latvalue; 
    487       facesAreWritten = true; 
    488  
    489     } // nvertex > 2 
    490      
    491   } // createMesh() 
    492229 
    493230///---------------------------------------------------------------- 
    … …  
    547284 
    548285      // For determining the global edge index 
    549       size_t nbEdgesOnProc = nbEdges_; 
    550       size_t nbEdgesAccum; 
     286      unsigned long nbEdgesOnProc = nbEdges_; 
     287      unsigned long nbEdgesAccum; 
    551288      MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    552289      nbEdgesAccum -= nbEdges_; 
    … …  
    677414         // Maps modified in this step: 
    678415         // nodeIdx2Idx = <idx, idxGlo> 
    679          int nodeCount = nodeIdx2Idx.size(); 
    680          int nodeStart, nbNodes; 
     416         unsigned long nodeCount = nodeIdx2Idx.size(); 
     417         unsigned long nodeStart, nbNodes; 
    681418         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    682419         int nNodes = nodeStart; 
    … …  
    770507 
    771508      // For determining the global face index 
    772       size_t nbFacesOnProc = nbFaces_; 
    773       size_t nbFacesAccum; 
     509      unsigned long nbFacesOnProc = nbFaces_; 
     510      unsigned long nbFacesAccum; 
    774511      MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    775512      nbFacesAccum -= nbFaces_; 
    … …  
    895632        } 
    896633 
    897         int edgeStart, nbEdges; 
     634        unsigned long edgeStart, nbEdges; 
    898635        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    899636        int nEdges = edgeStart; 
    … …  
    1115852        } 
    1116853 
    1117         int edgeCount = edgeIdx2Idx.size(); 
    1118         int edgeStart, nbEdges; 
     854        unsigned long edgeCount = edgeIdx2Idx.size(); 
     855        unsigned long edgeStart, nbEdges; 
    1119856        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    1120857        int nEdges = edgeStart; 
    … …  
    13851122        // Maps modified in this step: 
    13861123        // nodeIdx2Idx = <idx, idxGlo> 
    1387         int nodeCount = nodeIdx2Idx.size(); 
    1388         int nodeStart, nbNodes; 
     1124        unsigned long nodeCount = nodeIdx2Idx.size(); 
     1125        unsigned long nodeStart, nbNodes; 
    13891126        MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    13901127        int nNodes = nodeStart; 
    … …  
    15051242        } 
    15061243 
    1507         int edgeCount = edgeIdx2Idx.size(); 
    1508         int edgeStart, nbEdges; 
     1244        unsigned long edgeCount = edgeIdx2Idx.size(); 
     1245        unsigned long edgeStart, nbEdges; 
    15091246        MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    15101247        int nEdges = edgeStart; 
    … …  
    21031840 
    21041841    // faceToFaces connectivity 
    2105     boost::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant) 
     1842    std::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant) 
    21061843    int maxNb = 20;                            // some assumption on the max possible number of neighboring cells 
    21071844    faceToFaces.resize(maxNb, nbFaces); 
    … …  
    21601897    CArray<double, 2> faceToNodes (nvertex, nbFaces); 
    21611898 
    2162     boost::unordered_map <pair<double,double>, int> mapNodes; 
     1899    std::unordered_map <pairDouble, int, boost::hash<pairDouble> > mapNodes; 
    21631900 
    21641901    for (int nf = 0; nf < nbFaces; ++nf) 
    … …  
    21761913 
    21771914    // faceToFaces connectivity 
    2178     boost::unordered_map <pair<int,int>, int> mapEdges; 
     1915    std::unordered_map <pairInt, int, boost::hash<pairInt> > mapEdges; 
    21791916    faceToFaces.resize(nvertex, nbFaces); 
    21801917    CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible 
Note: See TracChangeset for help on using the changeset viewer.