Changeset 931


Ignore:
Timestamp:
09/20/16 16:48:07 (8 years ago)
Author:
oabramkina
Message:

Mesh connectivity:

Functions for determining local and global face-to-face connectivity added. Both types of neighboring cells (those that share a node or those that share an edge) are considered.

Location:
XIOS/trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/src/io/nc4_data_output.cpp

    r929 r931  
    450450      StdString domainName = domain->name; 
    451451      domain->assignMesh(domainName, domain->nvertex); 
    452       //domain->mesh->createMesh(domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv); 
    453452      domain->mesh->createMeshEpsilon(server->intraComm, domain->lonvalue_srv, domain->latvalue_srv, domain->bounds_lon_srv, domain->bounds_lat_srv); 
    454453 

  • XIOS/trunk/src/node/mesh.cpp

    r929 r931  
    15231523  ///---------------------------------------------------------------- 
    15241524  /*! 
    1525    * \fn  void CMesh::getNghbFacesNodeType(const MPI_Comm& comm, 
    1526                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
    1527                                CArray<int, 2>& nghbFaces) 
     1525   * \fn  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx, 
     1526                                              const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1527                                              CArray<int, 2>& nghbFaces) 
    15281528   * Finds neighboring cells of a local domain for node-type of neighbors. 
    15291529   * \param [in] comm 
     1530   * \param [in] face_idx Array with global indexes. 
    15301531   * \param [in] bounds_lon Array of boundary longitudes. 
    15311532   * \param [in] bounds_lat Array of boundary latitudes. 
    … …  
    15331534   */ 
    15341535 
    1535   void CMesh::getNghbFacesNodeType(const MPI_Comm& comm, 
     1536  void CMesh::getGloNghbFacesNodeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx, 
    15361537                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
    15371538                               CArray<int, 2>& nghbFaces) 
    15381539  { 
    1539     int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     1540    int nvertex = bounds_lon.rows(); 
    15401541    int nbFaces = bounds_lon.shape()[1]; 
    15411542    nghbFaces.resize(2, nbFaces*10);    // some estimate on max number of neighbouring cells 
    … …  
    15441545    MPI_Comm_rank(comm, &mpiRank); 
    15451546    MPI_Comm_size(comm, &mpiSize); 
    1546  
    1547     // For determining the global face index 
    1548     size_t nbFacesOnProc = nbFaces; 
    1549     size_t nbFacesAccum; 
    1550     MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    1551     nbFacesAccum -= nbFaces; 
    15521547 
    15531548    // (1) Generating unique node indexes 
    … …  
    16301625        it = nodeIdx2IdxMin.find(myNodeIdx); 
    16311626        size_t nodeIdxMin = (it->second)[0]; 
    1632         size_t faceIdx = nbFacesAccum + nf; 
     1627        size_t faceIdx = face_idx(nf); 
    16331628        if (nodeIdxMin2Face.count(nodeIdxMin) == 0) 
    16341629        { 
    … …  
    16801675    } 
    16811676    nghbFaces.resizeAndPreserve(2, nbNghb); 
    1682   } // getNghbFacesNodeType 
     1677  } // getGloNghbFacesNodeType 
    16831678 
    16841679  ///---------------------------------------------------------------- 
    16851680  /*! 
    1686    * \fn  void CMesh::getNghbFacesEdgeType(const MPI_Comm& comm, 
     1681   * \fn  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx, 
    16871682                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
    16881683                               CArray<int, 2>& nghbFaces) 
    16891684   * Finds neighboring cells of a local domain for edge-type of neighbors. 
    16901685   * \param [in] comm 
     1686   * \param [in] face_idx Array with global indexes. 
    16911687   * \param [in] bounds_lon Array of boundary longitudes. 
    16921688   * \param [in] bounds_lat Array of boundary latitudes. 
    … …  
    16941690   */ 
    16951691 
    1696   void CMesh::getNghbFacesEdgeType(const MPI_Comm& comm, 
     1692  void CMesh::getGloNghbFacesEdgeType(const MPI_Comm& comm, const CArray<int, 1>& face_idx, 
    16971693                               const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
    16981694                               CArray<int, 2>& nghbFaces) 
    16991695  { 
    1700     int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows(); 
     1696    int nvertex = bounds_lon.rows(); 
    17011697    int nbFaces = bounds_lon.shape()[1]; 
    17021698    nghbFaces.resize(2, nbFaces*10);    // estimate of max number of neighbouring cells 
    … …  
    17051701    MPI_Comm_rank(comm, &mpiRank); 
    17061702    MPI_Comm_size(comm, &mpiSize); 
    1707  
    1708     // For determining the global face index 
    1709     size_t nbFacesOnProc = nbFaces; 
    1710     size_t nbFacesAccum; 
    1711     MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm); 
    1712     nbFacesAccum -= nbFaces; 
    17131703 
    17141704    // (1) Generating unique node indexes 
    … …  
    17981788        size_t nodeIdxMin1 = (it1->second)[0]; 
    17991789        size_t nodeIdxMin2 = (it2->second)[0]; 
    1800         size_t faceIdx = nbFacesAccum + nf; 
     1790        size_t faceIdx = face_idx(nf); 
    18011791 
    18021792        if (nodeIdxMin1 != nodeIdxMin2) 
    … …  
    18281818    for (int nf = 0; nf < nbFaces; ++nf) 
    18291819    { 
    1830       nbNghb = 0; 
    18311820      for (int nv1 = 0; nv1 < nvertex; ++nv1) 
    18321821      { 
    … …  
    18661855    } 
    18671856    nghbFaces.resizeAndPreserve(2, nbNghb); 
    1868   } // getNghbFacesEdgeType 
     1857  } // getGloNghbFacesEdgeType 
    18691858 
    18701859  ///---------------------------------------------------------------- 
    18711860  /*! 
    1872    * \fn void getDomainExtended(const int nghbType, const MPI_Comm& comm, 
    1873                                 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
    1874                                 CArray<size_t, 1>& nghbFaces) 
     1861   * \fn void getGlobalNghbFaces (const int nghbType, const MPI_Comm& comm, const CArray<int, 1>& face_idx, 
     1862                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1863                                  CArray<size_t, 1>& nghbFaces) 
     1864   * Finds neighboring faces owned by other procs. 
    18751865   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges. 
    18761866   * \param [in] comm 
     1867   * \param [in] face_idx Array with global indexes. 
    18771868   * \param [in] bounds_lon Array of boundary longitudes. 
    18781869   * \param [in] bounds_lat Array of boundary latitudes. 
    1879    * \param [out] nghbFaces Array containing neighboring faces and their ranks. The ordering is face1, rank1,.... 
     1870   * \param [out] nghbFaces 2D array containing neighboring faces and owner ranks. 
    18801871   */ 
    18811872 
    1882   void CMesh::getDomainExtended(const int nghbType, const MPI_Comm& comm, 
    1883       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
    1884       CArray<int, 2>& nghbFaces) 
     1873  void CMesh::getGlobalNghbFaces(const int nghbType, const MPI_Comm& comm, 
     1874                                 const CArray<int, 1>& face_idx, 
     1875                                 const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1876                                 CArray<int, 2>& nghbFaces) 
    18851877  { 
    18861878    if (nghbType == 0) 
    1887       getNghbFacesNodeType(comm, bounds_lon, bounds_lat, nghbFaces); 
     1879      getGloNghbFacesNodeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces); 
    18881880    else 
    1889       getNghbFacesEdgeType(comm, bounds_lon, bounds_lat, nghbFaces); 
    1890   } // getDomainExtended 
     1881      getGloNghbFacesEdgeType(comm, face_idx, bounds_lon, bounds_lat, nghbFaces); 
     1882  } // getGlobalNghbFaces 
     1883 
     1884  ///---------------------------------------------------------------- 
     1885  /*! 
     1886   * \fn void getLocalNghbFaces (const int nghbType, 
     1887                                  const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1888                                  CArray<size_t, 1>& nghbFaces) 
     1889   * \param [in] nghbType 0 for faces sharing nodes, otherwise for faces sharing edges. 
     1890   * \param [in] bounds_lon Array of boundary longitudes. 
     1891   * \param [in] bounds_lat Array of boundary latitudes. 
     1892   * \param [out] nghbFaces 1D array containing neighboring faces. 
     1893   */ 
     1894 
     1895  void CMesh::getLocalNghbFaces(const int nghbType, const CArray<int, 1>& face_idx, 
     1896                                const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1897                                CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces) 
     1898  { 
     1899    if (nghbType == 0) 
     1900      getLocNghbFacesNodeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces); 
     1901    else 
     1902      getLocNghbFacesEdgeType(face_idx, bounds_lon, bounds_lat, nghbFaces, nbNghbFaces); 
     1903  } // getLocalNghbFaces 
     1904 
     1905  ///---------------------------------------------------------------- 
     1906  /*! 
     1907   * \fn void getLocNghbFacesNodeType (const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1908                                       CArray<int, 2>& nghbFaces) 
     1909   * \param [in] face_idx Array with global face indexing. 
     1910   * \param [in] bounds_lon Array of boundary longitudes. 
     1911   * \param [in] bounds_lat Array of boundary latitudes. 
     1912   * \param [out] nghbFaces 2D array containing neighboring faces. 
     1913   * \param [out] nbNghbFaces Array containing number of neighboring faces. 
     1914   */ 
     1915 
     1916  void CMesh::getLocNghbFacesNodeType (const CArray<int, 1>& face_idx, 
     1917                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1918                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces) 
     1919  { 
     1920    int nvertex = bounds_lon.rows(); 
     1921    int nbFaces = bounds_lon.shape()[1]; 
     1922    int nbNodes = 0; 
     1923    nbNghbFaces.resize(nbFaces); 
     1924    nbNghbFaces = 0; 
     1925 
     1926    // nodeToFaces connectivity 
     1927    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces; 
     1928    for (int nf = 0; nf < nbFaces; ++nf) 
     1929      for (int nv = 0; nv < nvertex; ++nv) 
     1930      { 
     1931        size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0]; 
     1932        nodeToFaces[nodeHash].push_back(face_idx(nf)); 
     1933      } 
     1934 
     1935    // faceToFaces connectivity 
     1936    boost::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant) 
     1937    int maxNb = 20;                            // some assumption on the max possible number of neighboring cells 
     1938    faceToFaces.resize(maxNb, nbFaces); 
     1939    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it; 
     1940    for (it = nodeToFaces.begin(); it != nodeToFaces.end(); ++it) 
     1941    { 
     1942      int size = it->second.size(); 
     1943      for (int i = 0; i < (size-1); ++i) 
     1944      { 
     1945        int face1 = it->second[i]; 
     1946        for (int j = i+1; j < size; ++j) 
     1947        { 
     1948          int face2 = it->second[j]; 
     1949          if (face2 != face1) 
     1950          { 
     1951            int hashFace = hashPairOrdered(face1, face2); 
     1952            if (mapFaces.count(hashFace) == 0) 
     1953            { 
     1954              faceToFaces(nbNghbFaces(face1), face1) = face2; 
     1955              faceToFaces(nbNghbFaces(face2), face2) = face1; 
     1956              ++nbNghbFaces(face1); 
     1957              ++nbNghbFaces(face2); 
     1958              mapFaces[hashFace] = hashFace; 
     1959            } 
     1960          } 
     1961        } 
     1962      } 
     1963    } 
     1964  } //getLocNghbFacesEdgeType 
     1965 
     1966 
     1967  ///---------------------------------------------------------------- 
     1968  /*! 
     1969   * \fn void getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx, 
     1970   *                                   const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1971   *                                   CArray<int, 2>& nghbFaces, CArray<int, 1>& nbNghbFaces) 
     1972   * \param [in] face_idx Array with global face indexing. 
     1973   * \param [in] bounds_lon Array of boundary longitudes. 
     1974   * \param [in] bounds_lat Array of boundary latitudes. 
     1975   * \param [out] nghbFaces 2D array containing neighboring faces. 
     1976   * \param [out] nbNghbFaces Array containing number of neighboring faces. 
     1977   */ 
     1978 
     1979  void CMesh::getLocNghbFacesEdgeType (const CArray<int, 1>& face_idx, 
     1980                                       const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat, 
     1981                                       CArray<int, 2>& faceToFaces, CArray<int, 1>& nbNghbFaces) 
     1982  { 
     1983    int nvertex = bounds_lon.rows(); 
     1984    int nbFaces = bounds_lon.shape()[1]; 
     1985    int nbNodes = 0; 
     1986    int nbEdges = 0; 
     1987    nbNghbFaces.resize(nbFaces); 
     1988    nbNghbFaces = 0; 
     1989 
     1990    // faceToNodes connectivity 
     1991    CArray<double, 2> faceToNodes (nvertex, nbFaces); 
     1992 
     1993    boost::unordered_map <pair<double,double>, int> mapNodes; 
     1994 
     1995    for (int nf = 0; nf < nbFaces; ++nf) 
     1996      for (int nv = 0; nv < nvertex; ++nv) 
     1997      { 
     1998        if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end()) 
     1999        { 
     2000          mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes; 
     2001          faceToNodes(nv,nf) = nbNodes ; 
     2002          ++nbNodes ; 
     2003        } 
     2004        else 
     2005          faceToNodes(nv,nf) = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))]; 
     2006      } 
     2007 
     2008    // faceToFaces connectivity 
     2009    boost::unordered_map <pair<int,int>, int> mapEdges; 
     2010    faceToFaces.resize(nvertex, nbFaces); 
     2011    CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible 
     2012 
     2013    for (int nf = 0; nf < nbFaces; ++nf) 
     2014    { 
     2015      for (int nv1 = 0; nv1 < nvertex; ++nv1) 
     2016      { 
     2017        int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation 
     2018        int face = face_idx(nf); 
     2019        int node1 = faceToNodes(nv1,face); 
     2020        int node2 = faceToNodes(nv2,face); 
     2021        if (node1 != node2) 
     2022        { 
     2023          if (mapEdges.find(make_ordered_pair (node1, node2)) == mapEdges.end()) 
     2024          { 
     2025            mapEdges[make_ordered_pair (node1, node2)] = nbEdges; 
     2026            edgeToFaces(0,nbEdges) = face; 
     2027            ++nbEdges; 
     2028          } 
     2029          else 
     2030          { 
     2031            int edge = mapEdges[make_ordered_pair (node1, node2)]; 
     2032            edgeToFaces(1, edge) = face; 
     2033            int face1 = face; 
     2034            int face2 = edgeToFaces(0,edge); 
     2035            faceToFaces(nbNghbFaces(face1), face1) =  face2; 
     2036            faceToFaces(nbNghbFaces(face2), face2) =  face1; 
     2037            ++nbNghbFaces(face1); 
     2038            ++nbNghbFaces(face2); 
     2039          } 
     2040        } // node1 != node2 
     2041      } // nv 
     2042    } // nf 
     2043 
     2044  } //getLocNghbFacesEdgeType 
     2045 
    18912046 
    18922047} // namespace xios 

  • XIOS/trunk/src/node/mesh.hpp

    r929 r931  
    6464                             const CArray<double, 2>&, const CArray<double, 2>& ); 
    6565 
    66       void getDomainExtended(const int, const MPI_Comm&, 
     66      void getGlobalNghbFaces(const int, const MPI_Comm&, const CArray<int, 1>&, 
    6767                              const CArray<double, 2>&, const CArray<double, 2>&, 
    6868                              CArray<int, 2>&); 
     69 
     70      void getLocalNghbFaces(const int, const CArray<int, 1>&, 
     71                             const CArray<double, 2>&, const CArray<double, 2>&, 
     72                             CArray<int, 2>&, CArray<int, 1>&); 
    6973             
    7074      static CMesh* getMesh(StdString, int); 
    … …  
    8084      CClientClientDHTSizet* pNodeGlobalIndex;                    // pointer to a map <nodeHash, nodeIdxGlo> 
    8185      CClientClientDHTSizet* pEdgeGlobalIndex;                    // pointer to a map <edgeHash, edgeIdxGlo> 
    82       void getNghbFacesNodeType(const MPI_Comm&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 
    83       void getNghbFacesEdgeType(const MPI_Comm&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 
     86      void getGloNghbFacesNodeType(const MPI_Comm&, const CArray<int, 1>&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 
     87      void getGloNghbFacesEdgeType(const MPI_Comm&, const CArray<int, 1>&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&); 
     88      void getLocNghbFacesNodeType(const CArray<int, 1>&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&, CArray<int, 1>&); 
     89      void getLocNghbFacesEdgeType(const CArray<int, 1>&, const CArray<double, 2>&, const CArray<double, 2>&, CArray<int, 2>&, CArray<int, 1>&); 
    8490 
    8591      vector<size_t> createHashes (const double, const double); 
Note: See TracChangeset for help on using the changeset viewer.