Changeset 941
- Timestamp:
- 09/22/16 10:59:17 (8 years ago)
- Location:
- XIOS/trunk/src
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/trunk/src/client_client_dht_decl.cpp
r924 r941 14 14 template class CClientClientDHTTemplate<int>; 15 15 template class CClientClientDHTTemplate<size_t>; 16 template class CClientClientDHTTemplate<double>; 16 17 template class CClientClientDHTTemplate<PairIntInt>; 17 18 -
XIOS/trunk/src/client_client_dht_template.hpp
r924 r941 126 126 typedef CClientClientDHTTemplate<int> CClientClientDHTInt; 127 127 typedef CClientClientDHTTemplate<size_t> CClientClientDHTSizet; 128 typedef CClientClientDHTTemplate<double> CClientClientDHTDouble; 128 129 typedef CClientClientDHTTemplate<PairIntInt> CClientClientDHTPairIntInt; 129 130 -
XIOS/trunk/src/config/expand_domain_attribute.conf
r935 r941 1 1 DECLARE_ATTRIBUTE(int, order) 2 3 DECLARE_ENUM2(type,node,edge) -
XIOS/trunk/src/declare_ref_func.hpp
r922 r941 18 18 bool hasDirect##type##Reference(void) const; \ 19 19 C##type* getDirect##type##Reference(void) const; \ 20 const StdString& get##type##OutputName(void) const; 20 const StdString& get##type##OutputName(void) const; \ 21 21 void setAttributesReference(bool apply = true); \ 22 bool hasRefTo(C##type* ref) const; \ 22 23 \ 23 24 private: \ 24 std::vector<C##type*> refObjects; \ 25 StdString outputName_; 25 std::vector<C##type*> refObjects; 26 26 27 27 // Definitions … … 32 32 std::set<C##type*> tmpRefObjects; \ 33 33 C##type* refer_ptr = this; \ 34 std::vector<C##type*>().swap(refObjects); \ 34 35 refObjects.push_back(this); \ 35 36 \ … … 107 108 return getId(); \ 108 109 } \ 110 \ 111 bool C##type::hasRefTo(C##type* ref) const \ 112 { \ 113 bool found = false; \ 114 for (int idx = 0; idx < refObjects.size(); ++idx) \ 115 if (ref == refObjects[idx]) { found = true; break; } \ 116 return found; \ 117 } 109 118 110 119 #endif // __XIOS_DECLARE_REF_FUNC_HPP__ -
XIOS/trunk/src/node/expand_domain.cpp
r939 r941 38 38 void CExpandDomain::checkValid(CDomain* domainDst) 39 39 { 40 if (CDomain::type_attr::unstructured != domainDst->type) 41 { 42 ERROR("CExpandDomain::checkValid(CDomain* domainDst)", 43 << "Domain extension is only supported for unstructured" << std::endl 44 << "Check type of domain destination, id = " << domainDst->getId()); 45 } 46 47 if (this->type.isEmpty()) this->type.setValue(CExpandDomain::type_attr::edge); 40 48 } 41 49 -
XIOS/trunk/src/node/mesh.cpp
- Property svn:executable set to *
r931 r941 1923 1923 nbNghbFaces.resize(nbFaces); 1924 1924 nbNghbFaces = 0; 1925 1926 // nodeToFaces connectivity 1927 CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces; 1925 int maxNb = 10; // assumption on max number of edges sharing the same node 1926 1927 // faceToNodes connectivity 1928 CArray<int, 2> nodeToFaces(maxNb, nbFaces*nvertex); 1929 nodeToFaces = -1; 1930 boost::unordered_map <pair<double,double>, int> mapNodes; 1931 vector<int> countNodes(nbFaces*nvertex); 1932 countNodes.assign(nbFaces*nvertex, 0); 1933 1928 1934 for (int nf = 0; nf < nbFaces; ++nf) 1929 1935 for (int nv = 0; nv < nvertex; ++nv) 1930 1936 { 1931 size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0]; 1932 nodeToFaces[nodeHash].push_back(face_idx(nf)); 1937 if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end()) 1938 { 1939 mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes; 1940 nodeToFaces(countNodes[nbNodes], nbNodes) = face_idx(nf); 1941 ++(countNodes[nbNodes]); 1942 ++nbNodes; 1943 } 1944 else 1945 { 1946 int node = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))]; 1947 nodeToFaces(countNodes[node], node) = face_idx(nf); 1948 ++(countNodes[node]); 1949 } 1933 1950 } 1951 nodeToFaces.resizeAndPreserve(maxNb, nbNodes); 1952 1953 vector<int> tmpVec (nbNodes); 1954 for (int i = 0; i < nbNodes; ++i) 1955 { 1956 for (int j = 0; j < maxNb; ++j) 1957 tmpVec[j] = nodeToFaces(j,i); 1958 } 1959 1934 1960 1935 1961 // faceToFaces connectivity 1936 1962 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 } 1963 faceToFaces.resize(nvertex, nbFaces); 1964 for (int nn = 0; nn < nbNodes; ++nn) 1965 { 1966 int nf1 = 0; 1967 while (nodeToFaces(nf1+1, nn) != -1) 1968 { 1969 int nf2 = nf1+1; 1970 int face1 = nodeToFaces(nf1, nn); 1971 while (nodeToFaces(nf2, nn) != -1) 1972 { 1973 int face2 = nodeToFaces(nf2, nn); 1974 int hash = hashPairOrdered(face1, face2); 1975 if (mapFaces.count(hash) == 0) 1976 { 1977 faceToFaces(nbNghbFaces(face1), face1) = face2; 1978 faceToFaces(nbNghbFaces(face2), face2) = face1; 1979 ++nbNghbFaces(face1); 1980 ++nbNghbFaces(face2); 1981 } 1982 else 1983 mapFaces[hash] = hash; 1984 ++nf2; 1985 } 1986 ++nf1; 1962 1987 } 1963 1988 } 1989 1964 1990 } //getLocNghbFacesEdgeType 1965 1991 -
XIOS/trunk/src/transformation/algo_types.hpp
r934 r941 25 25 #include "domain_algorithm_compute_connectivity.hpp" 26 26 #include "compute_connectivity_domain.hpp" 27 #include "expand_domain.hpp" 28 #include "domain_algorithm_expand.hpp" 27 29 28 30 #endif // __XIOS_ALGORITHM_TRANSFORMATION_TYPES_HPP__ -
XIOS/trunk/src/transformation/domain_algorithm_expand.cpp
r935 r941 3 3 \author Ha NGUYEN 4 4 \since 08 Aug 2016 5 \date 08 Aug20166 7 \brief Algorithm for expand onan domain.5 \date 19 Sep 2016 6 7 \brief Algorithm for expanding an domain. 8 8 */ 9 9 #include "domain_algorithm_expand.hpp" … … 13 13 #include "grid.hpp" 14 14 #include "grid_transformation_factory_impl.hpp" 15 #include "context.hpp" 16 #include "context_client.hpp" 15 17 16 18 namespace xios { … … 45 47 : CDomainAlgorithmTransformation(domainDestination, domainSource) 46 48 { 49 if (domainDestination == domainSource) 50 { 51 ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)", 52 << "Domain source and domain destination are the same. Please make sure domain destination refers to domain source" << std::endl 53 << "Domain source " <<domainSource->getId() << std::endl 54 << "Domain destination " <<domainDestination->getId() << std::endl); 55 } 56 57 if (!domainDestination->hasRefTo(domainSource)) 58 { 59 ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)", 60 << "Domain domain destination must refer to domain source (directly or indirectly) by domain_ref " << std::endl 61 << "Domain source " <<domainSource->getId() << std::endl 62 << "Domain destination " <<domainDestination->getId() << std::endl); 63 } 64 47 65 this->type_ = (ELEMENT_MODIFICATION_WITH_DATA); 48 66 expandDomain->checkValid(domainDestination); 49 } 67 68 switch (expandDomain->type) 69 { 70 case CExpandDomain::type_attr::node : 71 expandDomainNodeConnectivity(domainDestination, 72 domainSource); 73 break; 74 case CExpandDomain::type_attr::edge : 75 expandDomainEdgeConnectivity(domainDestination, 76 domainSource); 77 break; 78 default: 79 break; 80 } 81 } 82 83 /*! 84 * Expand domain with edge-type neighbor 85 * \param[in/out] domainDestination domain destination and will be modified 86 * \param[in] domainSource domain source 87 */ 88 void CDomainAlgorithmExpand::expandDomainEdgeConnectivity(CDomain* domainDestination, 89 CDomain* domainSource) 90 { 91 CContext* context = CContext::getCurrent(); 92 CContextClient* client=context->client; 93 94 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 95 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 96 CArray<int,2> neighborsSrc; 97 98 int type = 1; // For edge 99 CMesh mesh; 100 mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc); 101 updateDomainAttributes(domainDestination, domainSource, neighborsSrc); 102 } 103 104 /*! 105 * Expand domain with node-type neighbor 106 * \param[in/out] domainDestination domain destination and will be modified 107 * \param[in] domainSource domain source 108 */ 109 void CDomainAlgorithmExpand::expandDomainNodeConnectivity(CDomain* domainDestination, 110 CDomain* domainSource) 111 { 112 CContext* context = CContext::getCurrent(); 113 CContextClient* client=context->client; 114 115 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 116 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 117 CArray<int,2> neighborsSrc; 118 119 int type = 0; // For node 120 CMesh mesh; 121 mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc); 122 updateDomainAttributes(domainDestination, domainSource, neighborsSrc); 123 } 124 125 /*! 126 * Extend domain destination and update its attributes 127 * Suppose that domain destination and domain source have the same values for all attributes (by inheritance) 128 * \param [in/out] domainDestination domain destination 129 * \param [in] domainSource domain source 130 * \param [in] neighborsDomainSrc domain extended part 131 */ 132 void CDomainAlgorithmExpand::updateDomainAttributes(CDomain* domainDestination, 133 CDomain* domainSource, 134 CArray<int,2>& neighborsDomainSrc) 135 { 136 CContext* context = CContext::getCurrent(); 137 CContextClient* client=context->client; 138 139 // First of all, "copy" all attributes of domain source to domain destination 140 StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue() 141 : ""; 142 domainDestination->domain_ref.setValue(domainSource->getId()); 143 domainDestination->solveRefInheritance(true); 144 if (domainDstRef.empty()) domainDestination->domain_ref.reset(); 145 else domainDestination->domain_ref.setValue(domainDstRef); 146 147 // Now extend domain destination 148 int niGlob = domainSource->ni_glo; 149 CArray<bool,1>& mask_1d_src = domainSource->mask_1d; 150 CArray<int,1>& i_index_src = domainSource->i_index; 151 CArray<double,1>& lon_src = domainSource->lonvalue_1d; 152 CArray<double,1>& lat_src = domainSource->latvalue_1d; 153 CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d; 154 CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d; 155 CArray<int,1>& data_i_index_src = domainSource->data_i_index; 156 157 int oldNbLocal = i_index_src.numElements(), index, globalIndex; 158 // Uncompress data_i_index 159 CArray<int,1> data_i_index_src_full(oldNbLocal); 160 int nbUnMaskedPointOnLocalDomain = 0; 161 data_i_index_src_full = -1; // Suppose all values are masked 162 for (int idx = 0; idx < data_i_index_src.numElements(); ++idx) 163 { 164 int dataIdx = data_i_index_src(idx); 165 if ((0 <= dataIdx) && (dataIdx < oldNbLocal)) 166 { 167 data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIdx; 168 ++nbUnMaskedPointOnLocalDomain; 169 } 170 } 171 172 CArray<bool,1>& mask_1d_dst = domainDestination->mask_1d; 173 CArray<int,1>& i_index_dst = domainDestination->i_index; 174 CArray<int,1>& j_index_dst = domainDestination->j_index; 175 CArray<double,1>& lon_dst = domainDestination->lonvalue_1d; 176 CArray<double,1>& lat_dst = domainDestination->latvalue_1d; 177 CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d; 178 CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d; 179 CArray<int,1>& data_i_index_dst = domainDestination->data_i_index; 180 CArray<int,1>& data_j_index_dst = domainDestination->data_j_index; 181 182 // Resize all array-like attributes of domain destination 183 int nbNeighbor = neighborsDomainSrc.shape()[1]; 184 int newNbLocalDst = nbNeighbor + oldNbLocal; 185 int nVertex = bounds_lon_dst.shape()[0]; 186 187 mask_1d_dst.resizeAndPreserve(newNbLocalDst); 188 i_index_dst.resizeAndPreserve(newNbLocalDst); 189 j_index_dst.resizeAndPreserve(newNbLocalDst); 190 lon_dst.resizeAndPreserve(newNbLocalDst); 191 lat_dst.resizeAndPreserve(newNbLocalDst); 192 bounds_lon_dst.resizeAndPreserve(nVertex, newNbLocalDst); 193 bounds_lat_dst.resizeAndPreserve(nVertex, newNbLocalDst); 194 CArray<int,1> data_i_index_dst_full(newNbLocalDst); 195 data_i_index_dst_full(Range(0,oldNbLocal-1)) = data_i_index_src_full; 196 data_i_index_dst_full(Range(oldNbLocal,newNbLocalDst-1)) = -1; 197 198 // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...) 199 // Global index mapping between destination and source 200 this->transformationMapping_.resize(1); 201 this->transformationWeight_.resize(1); 202 TransformationIndexMap& transMap = this->transformationMapping_[0]; 203 TransformationWeightMap& transWeight = this->transformationWeight_[0]; 204 205 transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor())); 206 transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor())); 207 // First, index mapping for local domain 208 for (int idx = 0; idx < oldNbLocal; ++idx) 209 { 210 index = i_index_dst(idx); 211 transMap[index].push_back(index); 212 transWeight[index].push_back(1.0); 213 } 214 // Then, index mapping for extended part 215 for (int idx = 0; idx < nbNeighbor; ++idx) 216 { 217 index = idx + oldNbLocal; 218 globalIndex = neighborsDomainSrc(0,idx); 219 i_index_dst(index) = globalIndex; 220 j_index_dst(index) = 0; 221 transMap[globalIndex].push_back(globalIndex); 222 transWeight[globalIndex].push_back(1.0); 223 } 224 225 // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...) 226 CClientClientDHTDouble::Index2VectorInfoTypeMap localData; 227 localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor())); 228 // Information exchanged among domains (attention to their order), number in parentheses presents size of data 229 // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1) 230 int dataPackageSize = 1 + 1 + // lon + lat 231 nVertex + nVertex + //bounds_lon + bounds_lat 232 1 + // mask_1d_dst; 233 1; // data_i_index 234 // Initialize database 235 for (int idx = 0; idx < oldNbLocal; ++idx) 236 { 237 index = i_index_src(idx); 238 localData[index].resize(dataPackageSize); 239 std::vector<double>& data = localData[index]; 240 241 //Pack data 242 int dataIdx = 0; 243 data[dataIdx] = lon_src(idx);++dataIdx; 244 data[dataIdx] = lat_src(idx);++dataIdx; 245 for (int i = 0; i < nVertex; ++i) 246 { 247 data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx; 248 } 249 for (int i = 0; i < nVertex; ++i) 250 { 251 data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx; 252 } 253 data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1; ++dataIdx; 254 data[dataIdx] = data_i_index_src_full(idx); 255 } 256 257 CClientClientDHTDouble dhtData(localData,client->intraComm); 258 CArray<size_t,1> neighborInd(nbNeighbor); 259 for (int idx = 0; idx < nbNeighbor; ++idx) 260 neighborInd(idx) = neighborsDomainSrc(0,idx); 261 262 // Compute local data on other domains 263 dhtData.computeIndexInfoMapping(neighborInd); 264 CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap(); 265 CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it; 266 // Ok get neighbor data 267 size_t nIdx; 268 int nbUnMaskedPointOnExtendedPart = 0; 269 for (int idx = 0; idx < nbNeighbor; ++idx) 270 { 271 nIdx = neighborInd(idx); 272 it = neighborData.find(nIdx); 273 if (ite != it) 274 { 275 index = idx + oldNbLocal; 276 std::vector<double>& data = it->second; 277 // Unpack data 278 int dataIdx = 0; 279 lon_dst(index) = data[dataIdx]; ++dataIdx; 280 lat_dst(index) = data[dataIdx]; ++dataIdx; 281 for (int i = 0; i < nVertex; ++i) 282 { 283 bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx; 284 } 285 for (int i = 0; i < nVertex; ++i) 286 { 287 bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx; 288 } 289 mask_1d_dst(index) = (1.0 == data[dataIdx]) ? true : false; ++dataIdx; 290 data_i_index_dst_full(index) = (int)(data[dataIdx]); 291 if (0 <= data_i_index_dst_full(index)) 292 { 293 data_i_index_dst_full(index) = index; 294 ++nbUnMaskedPointOnExtendedPart; 295 } 296 } 297 } 298 299 300 // Finally, update data_i_index 301 int nbUnMaskedPointOnNewDstDomain = (nbUnMaskedPointOnExtendedPart + nbUnMaskedPointOnLocalDomain); 302 int count = 0, dataIdx; 303 for (int idx = 0; idx < newNbLocalDst; ++idx) 304 { 305 dataIdx = data_i_index_dst_full(idx); 306 if ((0 <= dataIdx) && (dataIdx < newNbLocalDst)) 307 { 308 ++count; 309 } 310 } 311 312 data_i_index_dst.resize(count); 313 data_j_index_dst.resize(count); 314 data_j_index_dst = 0; 315 316 count = 0; 317 for (int idx = 0; idx < newNbLocalDst; ++idx) 318 { 319 dataIdx = data_i_index_dst_full(idx); 320 if ((0 <= dataIdx) && (dataIdx < newNbLocalDst)) 321 { 322 data_i_index_dst(count) = dataIdx; 323 ++count; 324 } 325 } 326 327 // Update ni 328 domainDestination->ni.setValue(newNbLocalDst); 329 330 } 331 50 332 51 333 /*! -
XIOS/trunk/src/transformation/domain_algorithm_expand.hpp
r935 r941 28 28 29 29 static bool registerTrans(); 30 31 protected: 32 void expandDomainEdgeConnectivity(CDomain* domainDestination, CDomain* domainSource); 33 void expandDomainNodeConnectivity(CDomain* domainDestination, CDomain* domainSource); 34 void updateDomainAttributes(CDomain* domainDestination, 35 CDomain* domainSource, 36 CArray<int,2>& neighborsDomainSrc); 37 30 38 protected: 31 39 void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); -
XIOS/trunk/src/transformation/grid_transformation.cpp
r933 r941 340 340 algo = algoTransformation_[std::distance(itb, it)]; 341 341 342 if ((0 != algo) && (CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type())) // Only registered transformation can be executed 342 if ((0 != algo) && 343 ((CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type()) || 344 (CGenericAlgorithmTransformation::ELEMENT_MODIFICATION_WITH_DATA == algo->type()))) // Only registered transformation can be executed 343 345 { 344 346 algo->computeIndexSourceMapping(dataAuxInputs); -
XIOS/trunk/src/transformation/grid_transformation_selector.cpp
r934 r941 325 325 CDomainAlgorithmInterpolate::registerTrans(); 326 326 CDomainAlgorithmZoom::registerTrans(); 327 328 } 329 330 } 327 CDomainAlgorithmExpand::registerTrans(); 328 329 } 330 331 }
Note: See TracChangeset
for help on using the changeset viewer.