moab
|
00001 #include "NCHelperHOMME.hpp" 00002 #include "moab/ReadUtilIface.hpp" 00003 #include "moab/FileOptions.hpp" 00004 #include "moab/SpectralMeshTool.hpp" 00005 00006 #include <cmath> 00007 00008 #define ERRORR(rval, str) \ 00009 if (MB_SUCCESS != rval) {_readNC->readMeshIface->report_error("%s", str); return rval;} 00010 00011 #define ERRORS(err, str) \ 00012 if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_FAILURE;} 00013 00014 namespace moab { 00015 00016 NCHelperHOMME::NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet) 00017 : UcdNCHelper(readNC, fileId, opts, fileSet), 00018 _spectralOrder(-1), connectId(-1), isConnFile(false) 00019 { 00020 // Calculate spectral order 00021 std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("np"); 00022 if (attIt != readNC->globalAtts.end()) { 00023 int success = NCFUNC(get_att_int)(readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &_spectralOrder); 00024 if (success != 0) 00025 readNC->readMeshIface->report_error("%s", "Failed to read np global attribute int data."); 00026 else 00027 _spectralOrder--; // Spectral order is one less than np 00028 } 00029 else { 00030 // As can_read_file() returns true and there is no global attribute "np", it should be a connectivity file 00031 isConnFile = true; 00032 _spectralOrder = 3; // Assume np is 4 00033 } 00034 } 00035 00036 bool NCHelperHOMME::can_read_file(ReadNC* readNC, int fileId) 00037 { 00038 // If global attribute "np" exists then it should be the HOMME grid 00039 if (readNC->globalAtts.find("np") != readNC->globalAtts.end()) { 00040 // Make sure it is CAM grid 00041 std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source"); 00042 if (attIt == readNC->globalAtts.end()) { 00043 readNC->readMeshIface->report_error("%s", "File does not have source global attribute."); 00044 return false; 00045 } 00046 unsigned int sz = attIt->second.attLen; 00047 std::string att_data; 00048 att_data.resize(sz + 1); 00049 att_data[sz] = '\000'; 00050 int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]); 00051 if (success != 0) { 00052 readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data."); 00053 return false; 00054 } 00055 if (att_data.find("CAM") == std::string::npos) 00056 return false; 00057 00058 return true; 00059 } 00060 else { 00061 // If dimension names "ncol" AND "ncorners" AND "ncells" exist, then it should be the HOMME connectivity file 00062 // In this case, the mesh can still be created 00063 std::vector<std::string>& dimNames = readNC->dimNames; 00064 if ((std::find(dimNames.begin(), dimNames.end(), std::string("ncol")) != dimNames.end()) && 00065 (std::find(dimNames.begin(), dimNames.end(), std::string("ncorners")) != dimNames.end()) && 00066 (std::find(dimNames.begin(), dimNames.end(), std::string("ncells")) != dimNames.end())) 00067 return true; 00068 } 00069 00070 return false; 00071 } 00072 00073 ErrorCode NCHelperHOMME::init_mesh_vals() 00074 { 00075 std::vector<std::string>& dimNames = _readNC->dimNames; 00076 std::vector<int>& dimLens = _readNC->dimLens; 00077 std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo; 00078 00079 ErrorCode rval; 00080 unsigned int idx; 00081 std::vector<std::string>::iterator vit; 00082 00083 // Look for time dimension 00084 if (isConnFile) { 00085 // Connectivity file might not have time dimension 00086 } 00087 else { 00088 if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end()) 00089 idx = vit - dimNames.begin(); 00090 else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end()) 00091 idx = vit - dimNames.begin(); 00092 else { 00093 ERRORR(MB_FAILURE, "Couldn't find 'time' or 't' dimension."); 00094 } 00095 tDim = idx; 00096 nTimeSteps = dimLens[idx]; 00097 } 00098 00099 // Get number of vertices (labeled as number of columns) 00100 if ((vit = std::find(dimNames.begin(), dimNames.end(), "ncol")) != dimNames.end()) 00101 idx = vit - dimNames.begin(); 00102 else { 00103 ERRORR(MB_FAILURE, "Couldn't find 'ncol' dimension."); 00104 } 00105 vDim = idx; 00106 nVertices = dimLens[idx]; 00107 00108 // Set number of cells 00109 nCells = nVertices - 2; 00110 00111 // Get number of levels 00112 if (isConnFile) { 00113 // Connectivity file might not have level dimension 00114 } 00115 else { 00116 if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end()) 00117 idx = vit - dimNames.begin(); 00118 else if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end()) 00119 idx = vit - dimNames.begin(); 00120 else { 00121 ERRORR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension."); 00122 } 00123 levDim = idx; 00124 nLevels = dimLens[idx]; 00125 } 00126 00127 // Store lon values in xVertVals 00128 std::map<std::string, ReadNC::VarData>::iterator vmit; 00129 if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { 00130 rval = read_coordinate("lon", 0, nVertices - 1, xVertVals); 00131 ERRORR(rval, "Trouble reading 'lon' variable."); 00132 } 00133 else { 00134 ERRORR(MB_FAILURE, "Couldn't find 'lon' variable."); 00135 } 00136 00137 // Store lat values in yVertVals 00138 if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { 00139 rval = read_coordinate("lat", 0, nVertices - 1, yVertVals); 00140 ERRORR(rval, "Trouble reading 'lat' variable."); 00141 } 00142 else { 00143 ERRORR(MB_FAILURE, "Couldn't find 'lat' variable."); 00144 } 00145 00146 // Store lev values in levVals 00147 if (isConnFile) { 00148 // Connectivity file might not have level variable 00149 } 00150 else { 00151 if ((vmit = varInfo.find("lev")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { 00152 rval = read_coordinate("lev", 0, nLevels - 1, levVals); 00153 ERRORR(rval, "Trouble reading 'lev' variable."); 00154 00155 // Decide whether down is positive 00156 char posval[10]; 00157 int success = NCFUNC(get_att_text)(_fileId, (*vmit).second.varId, "positive", posval); 00158 if (0 == success && !strcmp(posval, "down")) { 00159 for (std::vector<double>::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit) 00160 (*dvit) *= -1.0; 00161 } 00162 } 00163 else { 00164 ERRORR(MB_FAILURE, "Couldn't find 'lev' variable."); 00165 } 00166 } 00167 00168 // Store time coordinate values in tVals 00169 if (isConnFile) { 00170 // Connectivity file might not have time variable 00171 } 00172 else { 00173 if ((vmit = varInfo.find("time")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { 00174 rval = read_coordinate("time", 0, nTimeSteps - 1, tVals); 00175 ERRORR(rval, "Trouble reading 'time' variable."); 00176 } 00177 else if ((vmit = varInfo.find("t")) != varInfo.end() && (*vmit).second.varDims.size() == 1) { 00178 rval = read_coordinate("t", 0, nTimeSteps - 1, tVals); 00179 ERRORR(rval, "Trouble reading 't' variable."); 00180 } 00181 else { 00182 // If expected time variable is not available, set dummy time coordinate values to tVals 00183 for (int t = 0; t < nTimeSteps; t++) 00184 tVals.push_back((double)t); 00185 } 00186 } 00187 00188 // For each variable, determine the entity location type and number of levels 00189 std::map<std::string, ReadNC::VarData>::iterator mit; 00190 for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) { 00191 ReadNC::VarData& vd = (*mit).second; 00192 00193 vd.entLoc = ReadNC::ENTLOCSET; 00194 if (std::find(vd.varDims.begin(), vd.varDims.end(), tDim) != vd.varDims.end()) { 00195 if (std::find(vd.varDims.begin(), vd.varDims.end(), vDim) != vd.varDims.end()) 00196 vd.entLoc = ReadNC::ENTLOCVERT; 00197 } 00198 00199 vd.numLev = 1; 00200 if (std::find(vd.varDims.begin(), vd.varDims.end(), levDim) != vd.varDims.end()) 00201 vd.numLev = nLevels; 00202 } 00203 00204 // Hack: create dummy variables for dimensions (like ncol) with no corresponding coordinate variables 00205 rval = create_dummy_variables(); 00206 ERRORR(rval, "Failed to create dummy variables."); 00207 00208 return MB_SUCCESS; 00209 } 00210 00211 // When noMesh option is used on this read, the old ReadNC class instance for last read can get out 00212 // of scope (and deleted). The old instance initialized localGidVerts properly when the mesh was 00213 // created, but it is now lost. The new instance (will not create the mesh with noMesh option) has 00214 // to restore it based on the existing mesh from last read 00215 ErrorCode NCHelperHOMME::check_existing_mesh() 00216 { 00217 Interface*& mbImpl = _readNC->mbImpl; 00218 Tag& mGlobalIdTag = _readNC->mGlobalIdTag; 00219 bool& noMesh = _readNC->noMesh; 00220 00221 if (noMesh && localGidVerts.empty()) { 00222 // We need to populate localGidVerts range with the gids of vertices from the tmp_set 00223 // localGidVerts is important in reading the variable data into the nodes 00224 // also, for our purposes, localGidVerts is truly the GLOBAL_ID tag data, not other 00225 // file_id tags that could get passed around in other scenarios for parallel reading 00226 // for nodal_partition, this local gid is easier, should be initialized with only 00227 // the owned nodes 00228 00229 // We need to get all vertices from tmp_set (it is the input set in no_mesh scenario) 00230 Range local_verts; 00231 ErrorCode rval = mbImpl->get_entities_by_dimension(_fileSet, 0, local_verts); 00232 if (MB_FAILURE == rval) 00233 return rval; 00234 00235 if (!local_verts.empty()) { 00236 std::vector<int> gids(local_verts.size()); 00237 00238 // !IMPORTANT : this has to be the GLOBAL_ID tag 00239 rval = mbImpl->tag_get_data(mGlobalIdTag, local_verts, &gids[0]); 00240 if (MB_FAILURE == rval) 00241 return rval; 00242 00243 // Restore localGidVerts 00244 std::copy(gids.rbegin(), gids.rend(), range_inserter(localGidVerts)); 00245 nLocalVertices = localGidVerts.size(); 00246 } 00247 } 00248 00249 return MB_SUCCESS; 00250 } 00251 00252 ErrorCode NCHelperHOMME::create_mesh(Range& faces) 00253 { 00254 Interface*& mbImpl = _readNC->mbImpl; 00255 std::string& fileName = _readNC->fileName; 00256 Tag& mGlobalIdTag = _readNC->mGlobalIdTag; 00257 const Tag*& mpFileIdTag = _readNC->mpFileIdTag; 00258 DebugOutput& dbgOut = _readNC->dbgOut; 00259 bool& spectralMesh = _readNC->spectralMesh; 00260 int& gatherSetRank = _readNC->gatherSetRank; 00261 00262 int rank = 0; 00263 int procs = 1; 00264 #ifdef USE_MPI 00265 bool& isParallel = _readNC->isParallel; 00266 if (isParallel) { 00267 ParallelComm*& myPcomm = _readNC->myPcomm; 00268 rank = myPcomm->proc_config().proc_rank(); 00269 procs = myPcomm->proc_config().proc_size(); 00270 } 00271 #endif 00272 00273 ErrorCode rval; 00274 int success = 0; 00275 00276 // Need to get/read connectivity data before creating elements 00277 std::string conn_fname; 00278 00279 if (isConnFile) { 00280 // Connectivity file has already been read 00281 connectId = _readNC->fileId; 00282 } 00283 else { 00284 // Try to open the connectivity file through CONN option, if used 00285 rval = _opts.get_str_option("CONN", conn_fname); 00286 if (MB_SUCCESS != rval) { 00287 // Default convention for reading HOMME is a file HommeMapping.nc in same dir as data file 00288 conn_fname = std::string(fileName); 00289 size_t idx = conn_fname.find_last_of("/"); 00290 if (idx != std::string::npos) 00291 conn_fname = conn_fname.substr(0, idx).append("/HommeMapping.nc"); 00292 else 00293 conn_fname = "HommeMapping.nc"; 00294 } 00295 #ifdef PNETCDF_FILE 00296 #ifdef USE_MPI 00297 if (isParallel) { 00298 ParallelComm*& myPcomm = _readNC->myPcomm; 00299 success = NCFUNC(open)(myPcomm->proc_config().proc_comm(), conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId); 00300 } 00301 else 00302 success = NCFUNC(open)(MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId); 00303 #endif 00304 #else 00305 success = NCFUNC(open)(conn_fname.c_str(), 0, &connectId); 00306 #endif 00307 ERRORS(success, "Failed on open."); 00308 } 00309 00310 std::vector<std::string> conn_names; 00311 std::vector<int> conn_vals; 00312 rval = _readNC->get_dimensions(connectId, conn_names, conn_vals); 00313 ERRORR(rval, "Failed to get dimensions for connectivity."); 00314 00315 // Read connectivity into temporary variable 00316 int num_fine_quads = 0; 00317 int num_coarse_quads = 0; 00318 int start_idx = 0; 00319 std::vector<std::string>::iterator vit; 00320 int idx = 0; 00321 if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncells")) != conn_names.end()) 00322 idx = vit - conn_names.begin(); 00323 else if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncenters")) != conn_names.end()) 00324 idx = vit - conn_names.begin(); 00325 else { 00326 ERRORR(MB_FAILURE, "Failed to get number of quads."); 00327 } 00328 int num_quads = conn_vals[idx]; 00329 if (!isConnFile && num_quads != nCells) { 00330 dbgOut.tprintf(1, "Warning: number of quads from %s and cells from %s are inconsistent; num_quads = %d, nCells = %d.\n", 00331 conn_fname.c_str(), fileName.c_str(), num_quads, nCells); 00332 } 00333 00334 // Get the connectivity into tmp_conn2 and permute into tmp_conn 00335 int cornerVarId; 00336 success = NCFUNC(inq_varid)(connectId, "element_corners", &cornerVarId); 00337 ERRORS(success, "Failed to get variable id."); 00338 NCDF_SIZE tmp_starts[2] = {0, 0}; 00339 NCDF_SIZE tmp_counts[2] = {4, static_cast<NCDF_SIZE>(num_quads)}; 00340 std::vector<int> tmp_conn(4 * num_quads), tmp_conn2(4 * num_quads); 00341 success = NCFUNCAG(_vara_int)(connectId, cornerVarId, tmp_starts, tmp_counts, &tmp_conn2[0]); 00342 ERRORS(success, "Failed to get temporary connectivity."); 00343 if (isConnFile) { 00344 // This data/connectivity file will be closed later in ReadNC::load_file() 00345 } 00346 else { 00347 success = NCFUNC(close)(connectId); 00348 ERRORS(success, "Failed on close."); 00349 } 00350 // Permute the connectivity 00351 for (int i = 0; i < num_quads; i++) { 00352 tmp_conn[4 * i] = tmp_conn2[i]; 00353 tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads]; 00354 tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads]; 00355 tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads]; 00356 } 00357 00358 // Need to know whether we'll be creating gather mesh later, to make sure we allocate enough space 00359 // in one shot 00360 bool create_gathers = false; 00361 if (rank == gatherSetRank) 00362 create_gathers = true; 00363 00364 // Compute the number of local quads, accounting for coarse or fine representation 00365 // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2 00366 int spectral_unit = (spectralMesh ? _spectralOrder * _spectralOrder : 1); 00367 // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh, num_coarse_quads = num_fine_quads 00368 num_coarse_quads = int(std::floor(1.0 * num_quads / (spectral_unit * procs))); 00369 // start_idx is the starting index in the HommeMapping connectivity list for this proc, before converting to coarse quad representation 00370 start_idx = 4 * rank * num_coarse_quads * spectral_unit; 00371 // iextra = # coarse quads extra after equal split over procs 00372 int iextra = num_quads % (procs * spectral_unit); 00373 if (rank < iextra) 00374 num_coarse_quads++; 00375 start_idx += 4 * spectral_unit * std::min(rank, iextra); 00376 // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned to this proc 00377 num_fine_quads = spectral_unit * num_coarse_quads; 00378 00379 // Now create num_coarse_quads 00380 EntityHandle* conn_arr; 00381 EntityHandle start_vertex; 00382 Range tmp_range; 00383 00384 // Read connectivity into that space 00385 EntityHandle* sv_ptr = NULL; 00386 EntityHandle start_quad; 00387 SpectralMeshTool smt(mbImpl, _spectralOrder); 00388 if (!spectralMesh) { 00389 rval = _readNC->readMeshIface->get_element_connect(num_coarse_quads, 4, 00390 MBQUAD, 0, start_quad, conn_arr, 00391 // Might have to create gather mesh later 00392 (create_gathers ? num_coarse_quads + num_quads : num_coarse_quads)); 00393 ERRORR(rval, "Failed to create quads."); 00394 tmp_range.insert(start_quad, start_quad + num_coarse_quads - 1); 00395 std::copy(&tmp_conn[start_idx], &tmp_conn[start_idx + 4 * num_fine_quads], conn_arr); 00396 std::copy(conn_arr, conn_arr + 4 * num_fine_quads, range_inserter(localGidVerts)); 00397 } 00398 else { 00399 rval = smt.create_spectral_elems(&tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx, &localGidVerts); 00400 ERRORR(rval, "Failed to create spectral elements."); 00401 int count, v_per_e; 00402 rval = mbImpl->connect_iterate(tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count); 00403 ERRORR(rval, "Failed to get connectivity of spectral elements."); 00404 rval = mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_range.begin(), tmp_range.end(), 00405 count, (void*&)sv_ptr); 00406 ERRORR(rval, "Failed to get fine connectivity of spectral elements."); 00407 } 00408 00409 // Create vertices 00410 nLocalVertices = localGidVerts.size(); 00411 std::vector<double*> arrays; 00412 rval = _readNC->readMeshIface->get_node_coords(3, nLocalVertices, 0, start_vertex, arrays, 00413 // Might have to create gather mesh later 00414 (create_gathers ? nLocalVertices + nVertices : nLocalVertices)); 00415 ERRORR(rval, "Couldn't create vertices in HOMME mesh."); 00416 00417 // Set vertex coordinates 00418 Range::iterator rit; 00419 double* xptr = arrays[0]; 00420 double* yptr = arrays[1]; 00421 double* zptr = arrays[2]; 00422 int i; 00423 for (i = 0, rit = localGidVerts.begin(); i < nLocalVertices; i++, ++rit) { 00424 assert(*rit < xVertVals.size() + 1); 00425 xptr[i] = xVertVals[(*rit) - 1]; // lon 00426 yptr[i] = yVertVals[(*rit) - 1]; // lat 00427 } 00428 00429 // Convert lon/lat/rad to x/y/z 00430 const double pideg = acos(-1.0) / 180.0; 00431 double rad = (isConnFile) ? 8000.0 : 8000.0 + levVals[0]; 00432 for (i = 0; i < nLocalVertices; i++) { 00433 double cosphi = cos(pideg * yptr[i]); 00434 double zmult = sin(pideg * yptr[i]); 00435 double xmult = cosphi * cos(xptr[i] * pideg); 00436 double ymult = cosphi * sin(xptr[i] * pideg); 00437 xptr[i] = rad * xmult; 00438 yptr[i] = rad * ymult; 00439 zptr[i] = rad * zmult; 00440 } 00441 00442 // Get ptr to gid memory for vertices 00443 Range vert_range(start_vertex, start_vertex + nLocalVertices - 1); 00444 void* data; 00445 int count; 00446 rval = mbImpl->tag_iterate(mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data); 00447 ERRORR(rval, "Failed to get tag iterator."); 00448 assert(count == nLocalVertices); 00449 int* gid_data = (int*) data; 00450 std::copy(localGidVerts.begin(), localGidVerts.end(), gid_data); 00451 00452 // Duplicate global id data, which will be used to resolve sharing 00453 if (mpFileIdTag) { 00454 rval = mbImpl->tag_iterate(*mpFileIdTag, vert_range.begin(), vert_range.end(), count, data); 00455 ERRORR(rval, "Failed to get tag iterator on file id tag."); 00456 assert(count == nLocalVertices); 00457 gid_data = (int*) data; 00458 std::copy(localGidVerts.begin(), localGidVerts.end(), gid_data); 00459 } 00460 00461 // Create map from file ids to vertex handles, used later to set connectivity 00462 std::map<EntityHandle, EntityHandle> vert_handles; 00463 for (rit = localGidVerts.begin(), i = 0; rit != localGidVerts.end(); ++rit, i++) 00464 vert_handles[*rit] = start_vertex + i; 00465 00466 // Compute proper handles in connectivity using offset 00467 for (int q = 0; q < 4 * num_coarse_quads; q++) { 00468 conn_arr[q] = vert_handles[conn_arr[q]]; 00469 assert(conn_arr[q]); 00470 } 00471 if (spectralMesh) { 00472 int verts_per_quad = (_spectralOrder + 1) * (_spectralOrder + 1); 00473 for (int q = 0; q < verts_per_quad * num_coarse_quads; q++) { 00474 sv_ptr[q] = vert_handles[sv_ptr[q]]; 00475 assert(sv_ptr[q]); 00476 } 00477 } 00478 00479 // Add new vertices and elements to the set 00480 faces.merge(tmp_range); 00481 tmp_range.insert(start_vertex, start_vertex + nLocalVertices - 1); 00482 rval = mbImpl->add_entities(_fileSet, tmp_range); 00483 ERRORR(rval, "Couldn't add new vertices and quads to file set."); 00484 00485 // Mark the set with the spectral order 00486 Tag sporder; 00487 rval = mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_CREAT | MB_TAG_SPARSE); 00488 ERRORR(rval, "Couldn't create spectral order tag."); 00489 rval = mbImpl->tag_set_data(sporder, &_fileSet, 1, &_spectralOrder); 00490 ERRORR(rval, "Couldn't set value for spectral order tag."); 00491 00492 if (create_gathers) { 00493 EntityHandle gather_set; 00494 rval = _readNC->readMeshIface->create_gather_set(gather_set); 00495 ERRORR(rval, "Trouble creating gather set."); 00496 00497 // Create vertices 00498 arrays.clear(); 00499 // Don't need to specify allocation number here, because we know enough verts were created before 00500 rval = _readNC->readMeshIface->get_node_coords(3, nVertices, 0, start_vertex, arrays); 00501 ERRORR(rval, "Couldn't create vertices in HOMME mesh for gather set."); 00502 00503 xptr = arrays[0]; 00504 yptr = arrays[1]; 00505 zptr = arrays[2]; 00506 for (i = 0; i < nVertices; i++) { 00507 double cosphi = cos(pideg * yVertVals[i]); 00508 double zmult = sin(pideg * yVertVals[i]); 00509 double xmult = cosphi * cos(xVertVals[i] * pideg); 00510 double ymult = cosphi * sin(xVertVals[i] * pideg); 00511 xptr[i] = rad * xmult; 00512 yptr[i] = rad * ymult; 00513 zptr[i] = rad * zmult; 00514 } 00515 00516 // Get ptr to gid memory for vertices 00517 Range gather_verts(start_vertex, start_vertex + nVertices - 1); 00518 rval = mbImpl->tag_iterate(mGlobalIdTag, gather_verts.begin(), gather_verts.end(), count, data); 00519 ERRORR(rval, "Failed to get tag iterator."); 00520 assert(count == nVertices); 00521 gid_data = (int*) data; 00522 for (int j = 1; j <= nVertices; j++) 00523 gid_data[j - 1] = j; 00524 // Set the file id tag too, it should be bigger something not interfering with global id 00525 if (mpFileIdTag) { 00526 rval = mbImpl->tag_iterate(*mpFileIdTag, gather_verts.begin(), gather_verts.end(), count, data); 00527 ERRORR(rval, "Failed to get tag iterator in file id tag."); 00528 assert(count == nVertices); 00529 gid_data = (int*) data; 00530 for (int j = 1; j <= nVertices; j++) 00531 gid_data[j - 1] = nVertices + j; // bigger than global id tag 00532 } 00533 00534 rval = mbImpl->add_entities(gather_set, gather_verts); 00535 ERRORR(rval, "Couldn't add vertices to gather set."); 00536 00537 // Create quads 00538 Range gather_quads; 00539 // Don't need to specify allocation number here, because we know enough quads were created before 00540 rval = _readNC->readMeshIface->get_element_connect(num_quads, 4, 00541 MBQUAD, 0, start_quad, conn_arr); 00542 ERRORR(rval, "Failed to create quads."); 00543 gather_quads.insert(start_quad, start_quad + num_quads - 1); 00544 std::copy(&tmp_conn[0], &tmp_conn[4 * num_quads], conn_arr); 00545 for (i = 0; i != 4 * num_quads; i++) 00546 conn_arr[i] += start_vertex - 1; // Connectivity array is shifted by where the gather verts start 00547 rval = mbImpl->add_entities(gather_set, gather_quads); 00548 ERRORR(rval, "Couldn't add quads to gather set."); 00549 } 00550 00551 return MB_SUCCESS; 00552 } 00553 00554 ErrorCode NCHelperHOMME::read_ucd_variable_to_nonset_allocate(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums) 00555 { 00556 Interface*& mbImpl = _readNC->mbImpl; 00557 std::vector<int>& dimLens = _readNC->dimLens; 00558 DebugOutput& dbgOut = _readNC->dbgOut; 00559 00560 ErrorCode rval = MB_SUCCESS; 00561 00562 Range* range = NULL; 00563 00564 // Get vertices in set 00565 Range verts; 00566 rval = mbImpl->get_entities_by_dimension(_fileSet, 0, verts); 00567 ERRORR(rval, "Trouble getting vertices in set."); 00568 assert("Should only have a single vertex subrange, since they were read in one shot" && 00569 verts.psize() == 1); 00570 00571 for (unsigned int i = 0; i < vdatas.size(); i++) { 00572 // Support non-set variables with 3 dimensions like (time, lev, ncol) 00573 assert(3 == vdatas[i].varDims.size()); 00574 00575 // For a non-set variable, time should be the first dimension 00576 assert(tDim == vdatas[i].varDims[0]); 00577 00578 // Set up readStarts and readCounts 00579 vdatas[i].readStarts.resize(3); 00580 vdatas[i].readCounts.resize(3); 00581 00582 // First: time 00583 vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later 00584 vdatas[i].readCounts[0] = 1; 00585 00586 // Next: lev 00587 vdatas[i].readStarts[1] = 0; 00588 vdatas[i].readCounts[1] = vdatas[i].numLev; 00589 00590 // Finally: ncol 00591 switch (vdatas[i].entLoc) { 00592 case ReadNC::ENTLOCVERT: 00593 // Vertices 00594 // Start from the first localGidVerts 00595 // Actually, this will be reset later on in a loop 00596 vdatas[i].readStarts[2] = localGidVerts[0] - 1; 00597 vdatas[i].readCounts[2] = nLocalVertices; 00598 range = &verts; 00599 break; 00600 default: 00601 ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable."); 00602 break; 00603 } 00604 00605 // Get variable size 00606 vdatas[i].sz = 1; 00607 for (std::size_t idx = 0; idx != 3; idx++) 00608 vdatas[i].sz *= vdatas[i].readCounts[idx]; 00609 00610 for (unsigned int t = 0; t < tstep_nums.size(); t++) { 00611 dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]); 00612 00613 if (tstep_nums[t] >= dimLens[tDim]) { 00614 ERRORR(MB_INDEX_OUT_OF_RANGE, "Wrong value for a timestep number."); 00615 } 00616 00617 // Get the tag to read into 00618 if (!vdatas[i].varTags[t]) { 00619 rval = get_tag_to_nonset(vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev); 00620 ERRORR(rval, "Trouble getting tag."); 00621 } 00622 00623 // Get ptr to tag space 00624 void* data; 00625 int count; 00626 rval = mbImpl->tag_iterate(vdatas[i].varTags[t], range->begin(), range->end(), count, data); 00627 ERRORR(rval, "Failed to get tag iterator."); 00628 assert((unsigned)count == range->size()); 00629 vdatas[i].varDatas[t] = data; 00630 } 00631 } 00632 00633 return rval; 00634 } 00635 00636 #ifdef PNETCDF_FILE 00637 ErrorCode NCHelperHOMME::read_ucd_variable_to_nonset_async(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums) 00638 { 00639 DebugOutput& dbgOut = _readNC->dbgOut; 00640 00641 ErrorCode rval = read_ucd_variable_to_nonset_allocate(vdatas, tstep_nums); 00642 ERRORR(rval, "Trouble allocating read variables."); 00643 00644 // Finally, read into that space 00645 int success; 00646 00647 for (unsigned int i = 0; i < vdatas.size(); i++) { 00648 std::size_t sz = vdatas[i].sz; 00649 00650 // A typical supported variable: float T(time, lev, ncol) 00651 // For tag values, need transpose (lev, ncol) to (ncol, lev) 00652 size_t ni = vdatas[i].readCounts[2]; // ncol 00653 size_t nj = 1; // Here we should just set nj to 1 00654 size_t nk = vdatas[i].readCounts[1]; // lev 00655 00656 for (unsigned int t = 0; t < tstep_nums.size(); t++) { 00657 // We will synchronize all these reads with the other processors, 00658 // so the wait will be inside this double loop; is it too much? 00659 size_t nb_reads = localGidVerts.psize(); 00660 std::vector<int> requests(nb_reads), statuss(nb_reads); 00661 size_t idxReq = 0; 00662 00663 // Tag data for this timestep 00664 void* data = vdatas[i].varDatas[t]; 00665 00666 // Set readStart for each timestep along time dimension 00667 vdatas[i].readStarts[0] = tstep_nums[t]; 00668 00669 switch (vdatas[i].varDataType) { 00670 case NC_BYTE: 00671 case NC_CHAR: { 00672 ERRORR(MB_FAILURE, "not implemented"); 00673 break; 00674 } 00675 case NC_DOUBLE: { 00676 // Copied from float case 00677 std::vector<double> tmpdoubledata(sz); 00678 00679 // In the case of ucd mesh, and on multiple proc, 00680 // we need to read as many times as subranges we have in the 00681 // localGidVerts range; 00682 // basically, we have to give a different point 00683 // for data to start, for every subrange :( 00684 size_t indexInDoubleArray = 0; 00685 size_t ic = 0; 00686 for (Range::pair_iterator pair_iter = localGidVerts.pair_begin(); 00687 pair_iter != localGidVerts.pair_end(); 00688 pair_iter++, ic++) { 00689 EntityHandle starth = pair_iter->first; 00690 EntityHandle endh = pair_iter->second; // inclusive 00691 vdatas[i].readStarts[2] = (NCDF_SIZE) (starth - 1); 00692 vdatas[i].readCounts[2] = (NCDF_SIZE) (endh - starth + 1); 00693 00694 // Do a partial read, in each subrange 00695 // wait outside this loop 00696 success = NCFUNCREQG(_vara_double)(_fileId, vdatas[i].varId, 00697 &(vdatas[i].readStarts[0]), &(vdatas[i].readCounts[0]), 00698 &(tmpdoubledata[indexInDoubleArray]), &requests[idxReq++]); 00699 ERRORS(success, "Failed to read double data in loop"); 00700 // We need to increment the index in double array for the 00701 // next subrange 00702 indexInDoubleArray += (endh - starth + 1) * 1 * vdatas[i].numLev; 00703 } 00704 assert(ic == localGidVerts.psize()); 00705 00706 success = ncmpi_wait_all(_fileId, requests.size(), &requests[0], &statuss[0]); 00707 ERRORS(success, "Failed on wait_all."); 00708 00709 if (vdatas[i].numLev != 1) 00710 // Transpose (lev, ncol) to (ncol, lev) 00711 success = kji_to_jik_stride(ni, nj, nk, data, &tmpdoubledata[0], localGidVerts); 00712 else { 00713 for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++) 00714 ((double*) data)[idx] = tmpdoubledata[idx]; 00715 } 00716 ERRORS(success, "Failed to read double data."); 00717 break; 00718 } 00719 case NC_FLOAT: { 00720 std::vector<float> tmpfloatdata(sz); 00721 00722 // In the case of ucd mesh, and on multiple proc, 00723 // we need to read as many times as subranges we have in the 00724 // localGidVerts range; 00725 // basically, we have to give a different point 00726 // for data to start, for every subrange :( 00727 size_t indexInFloatArray = 0; 00728 size_t ic = 0; 00729 for (Range::pair_iterator pair_iter = localGidVerts.pair_begin(); 00730 pair_iter != localGidVerts.pair_end(); 00731 pair_iter++, ic++) { 00732 EntityHandle starth = pair_iter->first; 00733 EntityHandle endh = pair_iter->second; // inclusive 00734 vdatas[i].readStarts[2] = (NCDF_SIZE) (starth - 1); 00735 vdatas[i].readCounts[2] = (NCDF_SIZE) (endh - starth + 1); 00736 00737 // Do a partial read, in each subrange 00738 // wait outside this loop 00739 success = NCFUNCREQG(_vara_float)(_fileId, vdatas[i].varId, 00740 &(vdatas[i].readStarts[0]), &(vdatas[i].readCounts[0]), 00741 &(tmpfloatdata[indexInFloatArray]), &requests[idxReq++]); 00742 ERRORS(success, "Failed to read float data in loop"); 00743 // We need to increment the index in float array for the 00744 // next subrange 00745 indexInFloatArray += (endh - starth + 1) * 1 * vdatas[i].numLev; 00746 } 00747 assert(ic == localGidVerts.psize()); 00748 00749 success = ncmpi_wait_all(_fileId, requests.size(), &requests[0], &statuss[0]); 00750 ERRORS(success, "Failed on wait_all."); 00751 00752 if (vdatas[i].numLev != 1) 00753 // Transpose (lev, ncol) to (ncol, lev) 00754 success = kji_to_jik_stride(ni, nj, nk, data, &tmpfloatdata[0], localGidVerts); 00755 else { 00756 for (std::size_t idx = 0; idx != tmpfloatdata.size(); idx++) 00757 ((float*) data)[idx] = tmpfloatdata[idx]; 00758 } 00759 ERRORS(success, "Failed to read float data."); 00760 break; 00761 } 00762 case NC_INT: { 00763 ERRORR(MB_FAILURE, "not implemented"); 00764 break; 00765 } 00766 case NC_SHORT: { 00767 ERRORR(MB_FAILURE, "not implemented"); 00768 break; 00769 } 00770 default: 00771 success = 1; 00772 } 00773 00774 if (success) 00775 ERRORR(MB_FAILURE, "Trouble reading variable."); 00776 } 00777 } 00778 00779 for (unsigned int i = 0; i < vdatas.size(); i++) { 00780 for (unsigned int t = 0; t < tstep_nums.size(); t++) { 00781 dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]); 00782 ErrorCode tmp_rval = convert_variable(vdatas[i], t); 00783 if (MB_SUCCESS != tmp_rval) 00784 rval = tmp_rval; 00785 } 00786 } 00787 00788 // Debug output, if requested 00789 if (1 == dbgOut.get_verbosity()) { 00790 dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str()); 00791 for (unsigned int i = 1; i < vdatas.size(); i++) 00792 dbgOut.printf(1, ", %s ", vdatas[i].varName.c_str()); 00793 dbgOut.tprintf(1, "\n"); 00794 } 00795 00796 return rval; 00797 } 00798 #else 00799 ErrorCode NCHelperHOMME::read_ucd_variable_to_nonset(std::vector<ReadNC::VarData>& vdatas, std::vector<int>& tstep_nums) 00800 { 00801 DebugOutput& dbgOut = _readNC->dbgOut; 00802 00803 ErrorCode rval = read_ucd_variable_to_nonset_allocate(vdatas, tstep_nums); 00804 ERRORR(rval, "Trouble allocating read variables."); 00805 00806 // Finally, read into that space 00807 int success; 00808 for (unsigned int i = 0; i < vdatas.size(); i++) { 00809 std::size_t sz = vdatas[i].sz; 00810 00811 // A typical supported variable: float T(time, lev, ncol) 00812 // For tag values, need transpose (lev, ncol) to (ncol, lev) 00813 size_t ni = vdatas[i].readCounts[2]; // ncol 00814 size_t nj = 1; // Here we should just set nj to 1 00815 size_t nk = vdatas[i].readCounts[1]; // lev 00816 00817 for (unsigned int t = 0; t < tstep_nums.size(); t++) { 00818 // Tag data for this timestep 00819 void* data = vdatas[i].varDatas[t]; 00820 00821 // Set readStart for each timestep along time dimension 00822 vdatas[i].readStarts[0] = tstep_nums[t]; 00823 00824 switch (vdatas[i].varDataType) { 00825 case NC_BYTE: 00826 case NC_CHAR: { 00827 ERRORR(MB_FAILURE, "not implemented"); 00828 break; 00829 } 00830 case NC_DOUBLE: { 00831 // Copied from float case 00832 std::vector<double> tmpdoubledata(sz); 00833 00834 // In the case of ucd mesh, and on multiple proc, 00835 // we need to read as many times as subranges we have in the 00836 // localGidVerts range; 00837 // basically, we have to give a different point 00838 // for data to start, for every subrange :( 00839 size_t indexInDoubleArray = 0; 00840 size_t ic = 0; 00841 for (Range::pair_iterator pair_iter = localGidVerts.pair_begin(); 00842 pair_iter != localGidVerts.pair_end(); 00843 pair_iter++, ic++) { 00844 EntityHandle starth = pair_iter->first; 00845 EntityHandle endh = pair_iter->second; // Inclusive 00846 vdatas[i].readStarts[2] = (NCDF_SIZE) (starth - 1); 00847 vdatas[i].readCounts[2] = (NCDF_SIZE) (endh - starth + 1); 00848 00849 success = NCFUNCAG(_vara_double)(_fileId, vdatas[i].varId, 00850 &(vdatas[i].readStarts[0]), &(vdatas[i].readCounts[0]), 00851 &(tmpdoubledata[indexInDoubleArray])); 00852 ERRORS(success, "Failed to read float data in loop"); 00853 // We need to increment the index in double array for the 00854 // next subrange 00855 indexInDoubleArray += (endh - starth + 1) * 1 * vdatas[i].numLev; 00856 } 00857 assert(ic == localGidVerts.psize()); 00858 00859 if (vdatas[i].numLev != 1) 00860 // Transpose (lev, ncol) to (ncol, lev) 00861 success = kji_to_jik_stride(ni, nj, nk, data, &tmpdoubledata[0], localGidVerts); 00862 else { 00863 for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++) 00864 ((double*) data)[idx] = tmpdoubledata[idx]; 00865 } 00866 ERRORS(success, "Failed to read double data."); 00867 break; 00868 } 00869 case NC_FLOAT: { 00870 std::vector<float> tmpfloatdata(sz); 00871 00872 // In the case of ucd mesh, and on multiple proc, 00873 // we need to read as many times as subranges we have in the 00874 // localGidVerts range; 00875 // basically, we have to give a different point 00876 // for data to start, for every subrange :( 00877 size_t indexInFloatArray = 0; 00878 size_t ic = 0; 00879 for (Range::pair_iterator pair_iter = localGidVerts.pair_begin(); 00880 pair_iter != localGidVerts.pair_end(); 00881 pair_iter++, ic++) { 00882 EntityHandle starth = pair_iter->first; 00883 EntityHandle endh = pair_iter->second; // Inclusive 00884 vdatas[i].readStarts[2] = (NCDF_SIZE) (starth - 1); 00885 vdatas[i].readCounts[2] = (NCDF_SIZE) (endh - starth + 1); 00886 00887 success = NCFUNCAG(_vara_float)(_fileId, vdatas[i].varId, 00888 &(vdatas[i].readStarts[0]), &(vdatas[i].readCounts[0]), 00889 &(tmpfloatdata[indexInFloatArray])); 00890 ERRORS(success, "Failed to read float data in loop"); 00891 // We need to increment the index in float array for the 00892 // next subrange 00893 indexInFloatArray += (endh - starth + 1) * 1 * vdatas[i].numLev; 00894 } 00895 assert(ic == localGidVerts.psize()); 00896 00897 if (vdatas[i].numLev != 1) 00898 // Transpose (lev, ncol) to (ncol, lev) 00899 success = kji_to_jik_stride(ni, nj, nk, data, &tmpfloatdata[0], localGidVerts); 00900 else { 00901 for (std::size_t idx = 0; idx != tmpfloatdata.size(); idx++) 00902 ((float*) data)[idx] = tmpfloatdata[idx]; 00903 } 00904 ERRORS(success, "Failed to read float data."); 00905 break; 00906 } 00907 case NC_INT: { 00908 ERRORR(MB_FAILURE, "not implemented"); 00909 break; 00910 } 00911 case NC_SHORT: { 00912 ERRORR(MB_FAILURE, "not implemented"); 00913 break; 00914 } 00915 default: 00916 success = 1; 00917 } 00918 00919 if (success) 00920 ERRORR(MB_FAILURE, "Trouble reading variable."); 00921 } 00922 } 00923 00924 for (unsigned int i = 0; i < vdatas.size(); i++) { 00925 for (unsigned int t = 0; t < tstep_nums.size(); t++) { 00926 dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]); 00927 ErrorCode tmp_rval = convert_variable(vdatas[i], t); 00928 if (MB_SUCCESS != tmp_rval) 00929 rval = tmp_rval; 00930 } 00931 } 00932 // Debug output, if requested 00933 if (1 == dbgOut.get_verbosity()) { 00934 dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str()); 00935 for (unsigned int i = 1; i < vdatas.size(); i++) 00936 dbgOut.printf(1, ", %s ", vdatas[i].varName.c_str()); 00937 dbgOut.tprintf(1, "\n"); 00938 } 00939 00940 return rval; 00941 } 00942 #endif 00943 00944 } // namespace moab