moab
|
00001 00016 //------------------------------------------------------------------------- 00017 // Filename : ReadHDF5.cpp 00018 // 00019 // Purpose : HDF5 Writer 00020 // 00021 // Creator : Jason Kraftcheck 00022 // 00023 // Creation Date : 04/18/04 00024 //------------------------------------------------------------------------- 00025 00026 #include <assert.h> 00027 /* include our MPI header before any HDF5 because otherwise 00028 it will get included indirectly by HDF5 */ 00029 #ifdef USE_MPI 00030 # include "moab_mpi.h" 00031 # include "moab/ParallelComm.hpp" 00032 #endif 00033 #include <H5Tpublic.h> 00034 #include <H5Ppublic.h> 00035 #include <H5Epublic.h> 00036 #include "moab/Interface.hpp" 00037 #include "Internals.hpp" 00038 #include "MBTagConventions.hpp" 00039 #include "ReadHDF5.hpp" 00040 #include "moab/CN.hpp" 00041 #include "moab/FileOptions.hpp" 00042 #ifdef HDF5_PARALLEL 00043 # include <H5FDmpi.h> 00044 # include <H5FDmpio.h> 00045 #endif 00046 //#include "WriteHDF5.hpp" 00047 00048 #include <stdlib.h> 00049 #include <string.h> 00050 #include <limits> 00051 #include <functional> 00052 #include <iostream> 00053 00054 #ifdef BLUEGENE 00055 # include <sys/vfs.h> 00056 // Magic number for GPFS file system (ASCII for 'G' 'P' 'F' 'S') 00057 # define BG_LOCKLESS_GPFS 0x47504653 00058 #endif 00059 00060 #include "IODebugTrack.hpp" 00061 #include "ReadHDF5Dataset.hpp" 00062 #include "ReadHDF5VarLen.hpp" 00063 #include "moab_mpe.h" 00064 00065 namespace moab { 00066 00067 /* If true, coordinates are read in blocked format (all X values before 00068 * Y values before Z values.) If undefined, then all coordinates for a 00069 * given vertex are read at the same time. 00070 */ 00071 const bool DEFAULT_BLOCKED_COORDINATE_IO = false; 00072 00073 /* If true, file is opened first by root node only to read summary, 00074 * file is the closed and the summary is broadcast to all nodes, after 00075 * which all nodes open file in parallel to read data. If undefined, 00076 * file is opened once in parallel and all nodes read summary data. 00077 */ 00078 const bool DEFAULT_BCAST_SUMMARY = true; 00079 00080 /* If true and all processors are to read the same block of data, 00081 * read it on one and broadcast to others rather than using collective 00082 * io 00083 */ 00084 const bool DEFAULT_BCAST_DUPLICATE_READS = true; 00085 00086 #define READ_HDF5_BUFFER_SIZE (128*1024*1024) 00087 00088 #define assert_range( PTR, CNT ) \ 00089 assert( (PTR) >= (void*)dataBuffer ); assert( ((PTR)+(CNT)) <= (void*)(dataBuffer + bufferSize) ); 00090 00091 #define error(A) process_error(A,&dbgOut,__FILE__,__LINE__) 00092 00093 // This function doesn't do anything useful. It's just a nice 00094 // place to set a break point to determine why the reader fails. 00095 //static inline ErrorCode process_error( ErrorCode rval ) 00096 // { return rval; } 00097 static inline ErrorCode process_error(ErrorCode code, DebugOutput* dbgOut, const char* file, int line) 00098 { 00099 if (MB_SUCCESS != code) { 00100 if (dbgOut) 00101 dbgOut->printf(1,"Failure with error code %s at %s:%d\n", ErrorCodeStr[code], file, line); 00102 #if defined(WITH_MPI) && !defined(NDEBUG) 00103 MPI_Abort( MPI_COMM_WORLD ): 00104 #endif 00105 } 00106 return code; 00107 } 00108 00109 // Call \c error function during HDF5 library errors to make 00110 // it easier to trap such errors in the debugger. This function 00111 // gets registered with the HDF5 library as a callback. It 00112 // works the same as the default (H5Eprint), except that it 00113 // also calls the \c error fuction as a no-op. 00114 #if defined(H5E_auto_t_vers) && H5E_auto_t_vers > 1 00115 static herr_t handle_hdf5_error( hid_t stack, void* data ) 00116 { 00117 ReadHDF5::HDF5ErrorHandler* h = reinterpret_cast<ReadHDF5::HDF5ErrorHandler*>(data); 00118 herr_t result = 0; 00119 if (h->func) 00120 result = (*h->func)(stack,h->data); 00121 process_error(MB_FAILURE,0,__FILE__,__LINE__); 00122 return result; 00123 } 00124 #else 00125 static herr_t handle_hdf5_error( void* data ) 00126 { 00127 ReadHDF5::HDF5ErrorHandler* h = reinterpret_cast<ReadHDF5::HDF5ErrorHandler*>(data); 00128 herr_t result = 0; 00129 if (h->func) 00130 result = (*h->func)(h->data); 00131 process_error(MB_FAILURE,0,__FILE__,__LINE__); 00132 return result; 00133 } 00134 #endif 00135 00136 static void copy_sorted_file_ids( const EntityHandle* sorted_ids, 00137 long num_ids, 00138 Range& results ) 00139 { 00140 Range::iterator hint = results.begin(); 00141 long i = 0; 00142 while (i < num_ids) { 00143 EntityHandle start = sorted_ids[i]; 00144 for (++i; i < num_ids && sorted_ids[i] == 1+sorted_ids[i-1]; ++i); 00145 hint = results.insert( hint, start, sorted_ids[i-1] ); 00146 } 00147 } 00148 00149 static void intersect( const mhdf_EntDesc& group, const Range& range, Range& result ) 00150 { 00151 Range::const_iterator s, e; 00152 s = Range::lower_bound( range.begin(), range.end(), group.start_id ); 00153 e = Range::lower_bound( s, range.end(), group.start_id + group.count ); 00154 result.merge( s, e ); 00155 } 00156 00157 #define debug_barrier() debug_barrier_line(__LINE__) 00158 void ReadHDF5::debug_barrier_line(int lineno) 00159 { 00160 #ifdef USE_MPI 00161 if (mpiComm) { 00162 const unsigned threshold = 2; 00163 static unsigned long count = 0; 00164 if (dbgOut.get_verbosity() >= threshold) { 00165 dbgOut.printf( threshold, "*********** Debug Barrier %lu (@%d)***********\n", ++count, lineno); 00166 MPI_Barrier( *mpiComm ); 00167 } 00168 } 00169 #else 00170 if (lineno) {} 00171 #endif 00172 } 00173 00174 class CheckOpenReadHDF5Handles 00175 { 00176 int fileline; 00177 mhdf_FileHandle handle; 00178 int enter_count; 00179 public: 00180 CheckOpenReadHDF5Handles(mhdf_FileHandle file, int line) 00181 : fileline(line), handle(file), 00182 enter_count(mhdf_countOpenHandles(file)) 00183 {} 00184 ~CheckOpenReadHDF5Handles() 00185 { 00186 int new_count = mhdf_countOpenHandles(handle); 00187 if (new_count != enter_count) { 00188 std::cout << "Leaked HDF5 object handle in function at " 00189 << __FILE__ << ":" << fileline << std::endl 00190 << "Open at entrance: " << enter_count << std::endl 00191 << "Open at exit: " << new_count << std::endl; 00192 } 00193 } 00194 }; 00195 00196 00197 #ifdef NDEBUG 00198 #define CHECK_OPEN_HANDLES 00199 #else 00200 #define CHECK_OPEN_HANDLES \ 00201 CheckOpenReadHDF5Handles check_open_handles_(filePtr,__LINE__) 00202 #endif 00203 00204 ReaderIface* ReadHDF5::factory( Interface* iface ) 00205 { return new ReadHDF5( iface ); } 00206 00207 ReadHDF5::ReadHDF5( Interface* iface ) 00208 : bufferSize( READ_HDF5_BUFFER_SIZE ), 00209 dataBuffer( 0 ), 00210 iFace( iface ), 00211 filePtr( 0 ), 00212 fileInfo( 0 ), 00213 readUtil( 0 ), 00214 handleType( 0 ), 00215 indepIO( H5P_DEFAULT ), 00216 collIO( H5P_DEFAULT ), 00217 debugTrack( false ), 00218 dbgOut(stderr), 00219 blockedCoordinateIO(DEFAULT_BLOCKED_COORDINATE_IO), 00220 bcastSummary(DEFAULT_BCAST_SUMMARY), 00221 bcastDuplicateReads(DEFAULT_BCAST_DUPLICATE_READS), 00222 setMeta(0) 00223 { 00224 } 00225 00226 ErrorCode ReadHDF5::init() 00227 { 00228 ErrorCode rval; 00229 00230 if (readUtil) 00231 return MB_SUCCESS; 00232 00233 indepIO = collIO = H5P_DEFAULT; 00234 //WriteHDF5::register_known_tag_types( iFace ); 00235 00236 handleType = H5Tcopy( H5T_NATIVE_ULONG ); 00237 if (handleType < 0) 00238 return error(MB_FAILURE); 00239 00240 if (H5Tset_size( handleType, sizeof(EntityHandle)) < 0) 00241 { 00242 H5Tclose( handleType ); 00243 return error(MB_FAILURE); 00244 } 00245 00246 rval = iFace->query_interface( readUtil ); 00247 if (MB_SUCCESS != rval) 00248 { 00249 H5Tclose( handleType ); 00250 return error(rval); 00251 } 00252 00253 idMap.clear(); 00254 fileInfo = 0; 00255 debugTrack = false; 00256 myPcomm = 0; 00257 00258 return MB_SUCCESS; 00259 } 00260 00261 00262 ReadHDF5::~ReadHDF5() 00263 { 00264 if (!readUtil) // init() failed. 00265 return; 00266 00267 delete [] setMeta; 00268 setMeta = 0; 00269 iFace->release_interface( readUtil ); 00270 H5Tclose( handleType ); 00271 } 00272 00273 ErrorCode ReadHDF5::set_up_read( const char* filename, 00274 const FileOptions& opts ) 00275 { 00276 ErrorCode rval; 00277 mhdf_Status status; 00278 indepIO = collIO = H5P_DEFAULT; 00279 mpiComm = 0; 00280 00281 if (MB_SUCCESS != init()) 00282 return error(MB_FAILURE); 00283 00284 #if defined(H5Eget_auto_vers) && H5Eget_auto_vers > 1 00285 herr_t err = H5Eget_auto( H5E_DEFAULT, &errorHandler.func, &errorHandler.data ); 00286 #else 00287 herr_t err = H5Eget_auto( &errorHandler.func, &errorHandler.data ); 00288 #endif 00289 if (err < 0) { 00290 errorHandler.func = 0; 00291 errorHandler.data = 0; 00292 } 00293 else { 00294 #if defined(H5Eset_auto_vers) && H5Eset_auto_vers > 1 00295 err = H5Eset_auto( H5E_DEFAULT, &handle_hdf5_error, &errorHandler ); 00296 #else 00297 err = H5Eset_auto( &handle_hdf5_error, &errorHandler ); 00298 #endif 00299 if (err < 0) { 00300 errorHandler.func = 0; 00301 errorHandler.data = 0; 00302 } 00303 } 00304 00305 00306 // Set up debug output 00307 int tmpval; 00308 if (MB_SUCCESS == opts.get_int_option("DEBUG_IO", 1, tmpval)) { 00309 dbgOut.set_verbosity(tmpval); 00310 dbgOut.set_prefix("H5M "); 00311 } 00312 dbgOut.limit_output_to_first_N_procs( 32 ); 00313 00314 // Enable some extra checks for reads. Note: amongst other things this 00315 // will print errors if the entire file is not read, so if doing a 00316 // partial read that is not a parallel read, this should be disabled. 00317 debugTrack = (MB_SUCCESS == opts.get_null_option("DEBUG_BINIO")); 00318 00319 opts.get_toggle_option("BLOCKED_COORDINATE_IO",DEFAULT_BLOCKED_COORDINATE_IO,blockedCoordinateIO); 00320 opts.get_toggle_option("BCAST_SUMMARY", DEFAULT_BCAST_SUMMARY, bcastSummary); 00321 opts.get_toggle_option("BCAST_DUPLICATE_READS",DEFAULT_BCAST_DUPLICATE_READS,bcastDuplicateReads); 00322 bool bglockless = (MB_SUCCESS == opts.get_null_option("BGLOCKLESS")); 00323 00324 // Handle parallel options 00325 std::string junk; 00326 bool use_mpio = (MB_SUCCESS == opts.get_null_option("USE_MPIO")); 00327 rval = opts.match_option("PARALLEL", "READ_PART"); 00328 bool parallel = (rval != MB_ENTITY_NOT_FOUND); 00329 nativeParallel = (rval == MB_SUCCESS); 00330 if (use_mpio && !parallel) { 00331 readUtil->report_error( "'USE_MPIO' option specified w/out 'PARALLEL' option" ); 00332 return MB_NOT_IMPLEMENTED; 00333 } 00334 00335 // This option is intended for testing purposes only, and thus 00336 // is not documented anywhere. Decreasing the buffer size can 00337 // expose bugs that would otherwise only be seen when reading 00338 // very large files. 00339 rval = opts.get_int_option( "BUFFER_SIZE", bufferSize ); 00340 if (MB_SUCCESS != rval) { 00341 bufferSize = READ_HDF5_BUFFER_SIZE; 00342 } 00343 else if (bufferSize < (int)std::max( sizeof(EntityHandle), sizeof(void*) )) { 00344 return error(MB_INVALID_SIZE); 00345 } 00346 00347 dataBuffer = (char*)malloc( bufferSize ); 00348 if (!dataBuffer) 00349 return error(MB_MEMORY_ALLOCATION_FAILED); 00350 00351 if (use_mpio || nativeParallel) { 00352 00353 // lockless file IO on IBM BlueGene 00354 std::string pfilename(filename); 00355 #ifdef BLUEGENE 00356 if (!bglockless && 0 != pfilename.find("bglockless:")) { 00357 // check for GPFS file system 00358 struct statfs fsdata; 00359 statfs( filename, &fsdata ); 00360 if (fsdata.f_type == BG_LOCKLESS_GPFS) { 00361 bglockless = true; 00362 } 00363 } 00364 #endif 00365 if (bglockless) { 00366 pfilename = std::string("bglockless:") + pfilename; 00367 } 00368 00369 #ifndef HDF5_PARALLEL 00370 readUtil->report_error("MOAB not configured with parallel HDF5 support"); 00371 free(dataBuffer); 00372 dataBuffer = NULL; 00373 return MB_NOT_IMPLEMENTED; 00374 #else 00375 00376 MPI_Info info = MPI_INFO_NULL; 00377 std::string cb_size; 00378 rval = opts.get_str_option("CB_BUFFER_SIZE", cb_size); 00379 if (MB_SUCCESS == rval) { 00380 MPI_Info_create (&info); 00381 MPI_Info_set (info, const_cast<char*>("cb_buffer_size"), const_cast<char*>(cb_size.c_str())); 00382 } 00383 00384 int pcomm_no = 0; 00385 rval = opts.get_int_option("PARALLEL_COMM", pcomm_no); 00386 if (rval == MB_TYPE_OUT_OF_RANGE) { 00387 readUtil->report_error("Invalid value for PARALLEL_COMM option"); 00388 return rval; 00389 } 00390 myPcomm = ParallelComm::get_pcomm(iFace, pcomm_no); 00391 if (0 == myPcomm) { 00392 myPcomm = new ParallelComm(iFace, MPI_COMM_WORLD); 00393 } 00394 const int rank = myPcomm->proc_config().proc_rank(); 00395 dbgOut.set_rank(rank); 00396 dbgOut.limit_output_to_first_N_procs( 32 ); 00397 mpiComm = new MPI_Comm(myPcomm->proc_config().proc_comm()); 00398 if (bglockless) { 00399 dbgOut.printf( 1, "Enabling lockless IO for BlueGene (filename: \"%s\")\n", pfilename.c_str() ); 00400 } 00401 00402 #ifndef H5_MPI_COMPLEX_DERIVED_DATATYPE_WORKS 00403 dbgOut.print(1,"H5_MPI_COMPLEX_DERIVED_DATATYPE_WORKS is not defined\n"); 00404 #endif 00405 00406 // Open the file in serial on root to read summary 00407 dbgOut.tprint( 1, "Getting file summary\n" ); 00408 fileInfo = 0; 00409 00410 00411 hid_t file_prop; 00412 if (bcastSummary) { 00413 unsigned long size = 0; 00414 if (rank == 0) { 00415 00416 file_prop = H5Pcreate(H5P_FILE_ACCESS); 00417 err = H5Pset_fapl_mpio(file_prop, MPI_COMM_SELF, MPI_INFO_NULL); 00418 assert(file_prop >= 0); 00419 assert(err >= 0); 00420 filePtr = mhdf_openFileWithOpt( pfilename.c_str(), 0, NULL, handleType, file_prop, &status ); 00421 H5Pclose( file_prop ); 00422 00423 if (filePtr) { 00424 fileInfo = mhdf_getFileSummary( filePtr, handleType, &status ); 00425 if (!is_error(status)) { 00426 size = fileInfo->total_size; 00427 fileInfo->offset = (unsigned char*)fileInfo; 00428 } 00429 } 00430 mhdf_closeFile( filePtr, &status ); 00431 if (fileInfo && mhdf_isError(&status)) { 00432 free(fileInfo); 00433 fileInfo = NULL; 00434 } 00435 } 00436 00437 dbgOut.tprint( 1, "Communicating file summary\n" ); 00438 int mpi_err = MPI_Bcast( &size, 1, MPI_UNSIGNED_LONG, 0, myPcomm->proc_config().proc_comm() ); 00439 if (mpi_err || !size) 00440 return MB_FAILURE; 00441 00442 00443 if (rank != 0) 00444 fileInfo = reinterpret_cast<mhdf_FileDesc*>( malloc( size ) ); 00445 00446 MPI_Bcast( fileInfo, size, MPI_BYTE, 0, myPcomm->proc_config().proc_comm() ); 00447 00448 if (rank != 0) 00449 mhdf_fixFileDesc( fileInfo, reinterpret_cast<mhdf_FileDesc*>(fileInfo->offset) ); 00450 } 00451 00452 file_prop = H5Pcreate(H5P_FILE_ACCESS); 00453 err = H5Pset_fapl_mpio(file_prop, myPcomm->proc_config().proc_comm(), info); 00454 assert(file_prop >= 0); 00455 assert(err >= 0); 00456 00457 collIO = H5Pcreate(H5P_DATASET_XFER); 00458 assert(collIO > 0); 00459 err = H5Pset_dxpl_mpio(collIO, H5FD_MPIO_COLLECTIVE); 00460 assert(err >= 0); 00461 indepIO = nativeParallel ? H5P_DEFAULT : collIO; 00462 00463 // re-open file in parallel 00464 dbgOut.tprintf( 1, "Opening \"%s\" for parallel IO\n", pfilename.c_str() ); 00465 filePtr = mhdf_openFileWithOpt( pfilename.c_str(), 0, NULL, handleType, file_prop, &status ); 00466 00467 H5Pclose( file_prop ); 00468 if (!filePtr) 00469 { 00470 readUtil->report_error("%s", mhdf_message( &status )); 00471 free( dataBuffer ); 00472 dataBuffer = NULL; 00473 H5Pclose( indepIO ); 00474 if (collIO != indepIO) 00475 H5Pclose( collIO ); 00476 collIO = indepIO = H5P_DEFAULT; 00477 return error(MB_FAILURE); 00478 } 00479 00480 if (!bcastSummary) { 00481 fileInfo = mhdf_getFileSummary( filePtr, handleType, &status ); 00482 if (is_error(status)) { 00483 readUtil->report_error( "%s", mhdf_message( &status ) ); 00484 free( dataBuffer ); 00485 dataBuffer = NULL; 00486 mhdf_closeFile( filePtr, &status ); 00487 return error(MB_FAILURE); 00488 } 00489 } 00490 00491 #endif // HDF5_PARALLEL 00492 } 00493 else { 00494 00495 // Open the file 00496 filePtr = mhdf_openFile( filename, 0, NULL, handleType, &status ); 00497 if (!filePtr) 00498 { 00499 readUtil->report_error( "%s", mhdf_message( &status )); 00500 free( dataBuffer ); 00501 dataBuffer = NULL; 00502 return error(MB_FAILURE); 00503 } 00504 00505 // get file info 00506 fileInfo = mhdf_getFileSummary( filePtr, handleType, &status ); 00507 if (is_error(status)) { 00508 free( dataBuffer ); 00509 dataBuffer = NULL; 00510 mhdf_closeFile( filePtr, &status ); 00511 return error(MB_FAILURE); 00512 } 00513 } 00514 00515 ReadHDF5Dataset::default_hyperslab_selection_limit(); 00516 int hslimit; 00517 rval = opts.get_int_option( "HYPERSLAB_SELECT_LIMIT", hslimit ); 00518 if (MB_SUCCESS == rval && hslimit > 0) 00519 ReadHDF5Dataset::set_hyperslab_selection_limit( hslimit ); 00520 else 00521 ReadHDF5Dataset::default_hyperslab_selection_limit(); 00522 if (MB_SUCCESS != opts.get_null_option("HYPERSLAB_OR") && 00523 (MB_SUCCESS == opts.get_null_option( "HYPERSLAB_APPEND" ) 00524 || HDF5_can_append_hyperslabs())) { 00525 ReadHDF5Dataset::append_hyperslabs(); 00526 if (MB_SUCCESS != opts.get_int_option( "HYPERSLAB_SELECT_LIMIT", hslimit )) 00527 ReadHDF5Dataset::set_hyperslab_selection_limit( std::numeric_limits<int>::max() ); 00528 dbgOut.print(1,"Using H5S_APPEND for hyperslab selection\n"); 00529 } 00530 00531 return MB_SUCCESS; 00532 } 00533 00534 ErrorCode ReadHDF5::clean_up_read( const FileOptions& ) 00535 { 00536 HDF5ErrorHandler handler; 00537 #if defined(H5Eget_auto_vers) && H5Eget_auto_vers > 1 00538 herr_t err = H5Eget_auto( H5E_DEFAULT, &handler.func, &handler.data ); 00539 #else 00540 herr_t err = H5Eget_auto( &handler.func, &handler.data ); 00541 #endif 00542 if (err >= 0 && handler.func == &handle_hdf5_error) { 00543 assert(handler.data = &errorHandler); 00544 #if defined(H5Eget_auto_vers) && H5Eget_auto_vers > 1 00545 H5Eset_auto( H5E_DEFAULT, errorHandler.func, errorHandler.data ); 00546 #else 00547 H5Eset_auto( errorHandler.func, errorHandler.data ); 00548 #endif 00549 } 00550 00551 free( dataBuffer ); 00552 dataBuffer = NULL; 00553 free( fileInfo ); 00554 fileInfo = NULL; 00555 delete mpiComm; 00556 mpiComm = 0; 00557 00558 if (indepIO != H5P_DEFAULT) 00559 H5Pclose( indepIO ); 00560 if (collIO != indepIO) 00561 H5Pclose( collIO ); 00562 collIO = indepIO = H5P_DEFAULT; 00563 00564 delete [] setMeta; 00565 setMeta = 0; 00566 00567 mhdf_Status status; 00568 mhdf_closeFile( filePtr, &status ); 00569 filePtr = 0; 00570 return is_error(status) ? MB_FAILURE : MB_SUCCESS; 00571 } 00572 00573 ErrorCode ReadHDF5::load_file( const char* filename, 00574 const EntityHandle* file_set, 00575 const FileOptions& opts, 00576 const ReaderIface::SubsetList* subset_list, 00577 const Tag* file_id_tag ) 00578 { 00579 ErrorCode rval; 00580 00581 rval = set_up_read( filename, opts ); 00582 if (MB_SUCCESS != rval) { 00583 clean_up_read(opts); 00584 return rval; 00585 } 00586 00587 // We read the entire set description table regarless of partial 00588 // or complete reads or serial vs parallel reads 00589 rval = read_all_set_meta(); 00590 00591 if (subset_list && MB_SUCCESS == rval) 00592 rval = load_file_partial( subset_list->tag_list, 00593 subset_list->tag_list_length, 00594 subset_list->num_parts, 00595 subset_list->part_number, 00596 opts ); 00597 else 00598 rval = load_file_impl( opts ); 00599 00600 if (MB_SUCCESS == rval && file_id_tag) { 00601 dbgOut.tprint( 1, "Storing file IDs in tag\n" ); 00602 rval = store_file_ids( *file_id_tag ); 00603 } 00604 00605 if (MB_SUCCESS == rval && 0 != file_set) { 00606 dbgOut.tprint( 1, "Reading QA records\n" ); 00607 rval = read_qa( *file_set ); 00608 } 00609 00610 00611 dbgOut.tprint( 1, "Cleaning up\n" ); 00612 ErrorCode rval2 = clean_up_read( opts ); 00613 if (rval == MB_SUCCESS && rval2 != MB_SUCCESS) 00614 rval = rval2; 00615 00616 if (MB_SUCCESS == rval) 00617 dbgOut.tprint(1, "Read finished.\n"); 00618 else { 00619 std::string msg; 00620 iFace->get_last_error(msg); 00621 dbgOut.tprintf(1,"READ FAILED (ERROR CODE %s): %s\n", ErrorCodeStr[rval], msg.c_str()); 00622 } 00623 00624 if (H5P_DEFAULT != collIO) 00625 H5Pclose( collIO ); 00626 if (H5P_DEFAULT != indepIO) 00627 H5Pclose( indepIO ); 00628 collIO = indepIO = H5P_DEFAULT; 00629 00630 return rval; 00631 } 00632 00633 00634 00635 ErrorCode ReadHDF5::load_file_impl( const FileOptions& ) 00636 { 00637 ErrorCode rval; 00638 mhdf_Status status; 00639 std::string tagname; 00640 int i; 00641 00642 CHECK_OPEN_HANDLES; 00643 00644 dbgOut.tprint(1, "Reading all nodes...\n"); 00645 Range ids; 00646 if (fileInfo->nodes.count) { 00647 ids.insert( fileInfo->nodes.start_id, 00648 fileInfo->nodes.start_id + fileInfo->nodes.count - 1); 00649 rval = read_nodes( ids ); 00650 if (MB_SUCCESS != rval) 00651 return error(rval); 00652 } 00653 00654 00655 dbgOut.tprint(1, "Reading all element connectivity...\n"); 00656 std::vector<int> polyhedra; // need to do these last so that faces are loaded 00657 for (i = 0; i < fileInfo->num_elem_desc; ++i) { 00658 if (CN::EntityTypeFromName(fileInfo->elems[i].type) == MBPOLYHEDRON) { 00659 polyhedra.push_back(i); 00660 continue; 00661 } 00662 00663 rval = read_elems( i ); 00664 if (MB_SUCCESS != rval) 00665 return error(rval); 00666 } 00667 for (std::vector<int>::iterator it = polyhedra.begin(); 00668 it != polyhedra.end(); ++it) { 00669 rval = read_elems( *it ); 00670 if (MB_SUCCESS != rval) 00671 return error(rval); 00672 } 00673 00674 dbgOut.tprint(1, "Reading all sets...\n"); 00675 ids.clear(); 00676 if (fileInfo->sets.count) { 00677 ids.insert( fileInfo->sets.start_id, 00678 fileInfo->sets.start_id + fileInfo->sets.count - 1); 00679 rval = read_sets( ids ); 00680 if (rval != MB_SUCCESS) { 00681 return error(rval); 00682 } 00683 } 00684 00685 dbgOut.tprint(1, "Reading all adjacencies...\n"); 00686 for (i = 0; i < fileInfo->num_elem_desc; ++i) { 00687 if (!fileInfo->elems[i].have_adj) 00688 continue; 00689 00690 long table_len; 00691 hid_t table = mhdf_openAdjacency( filePtr, 00692 fileInfo->elems[i].handle, 00693 &table_len, 00694 &status ); 00695 if (is_error(status)) 00696 return error(MB_FAILURE); 00697 00698 rval = read_adjacencies( table, table_len ); 00699 mhdf_closeData( filePtr, table, &status ); 00700 if (MB_SUCCESS != rval) 00701 return error(rval); 00702 if (is_error(status)) 00703 return error(MB_FAILURE); 00704 } 00705 00706 dbgOut.tprint(1, "Reading all tags...\n"); 00707 for (i = 0; i < fileInfo->num_tag_desc; ++i) { 00708 rval = read_tag( i ); 00709 if (MB_SUCCESS != rval) 00710 return error(rval); 00711 } 00712 00713 dbgOut.tprint(1, "Core read finished. Cleaning up...\n"); 00714 return MB_SUCCESS; 00715 } 00716 00717 ErrorCode ReadHDF5::find_int_tag( const char* name, int& index ) 00718 { 00719 for (index = 0; index < fileInfo->num_tag_desc; ++index) 00720 if (!strcmp( name, fileInfo->tags[index].name)) 00721 break; 00722 00723 if (index == fileInfo->num_tag_desc) { 00724 readUtil->report_error( "File does not contain subset tag '%s'", name ); 00725 return error(MB_TAG_NOT_FOUND); 00726 } 00727 00728 if (fileInfo->tags[index].type != mhdf_INTEGER || 00729 fileInfo->tags[index].size != 1) { 00730 readUtil->report_error( "Tag '%s' does not containa single integer value", name ); 00731 return error(MB_TYPE_OUT_OF_RANGE); 00732 } 00733 00734 return MB_SUCCESS; 00735 } 00736 00737 ErrorCode ReadHDF5::get_subset_ids( const ReaderIface::IDTag* subset_list, 00738 int subset_list_length, 00739 Range& file_ids ) 00740 { 00741 ErrorCode rval; 00742 00743 for (int i = 0; i < subset_list_length; ++i) { 00744 00745 int tag_index; 00746 rval = find_int_tag( subset_list[i].tag_name, tag_index ); 00747 if (MB_SUCCESS != rval) 00748 return error(rval); 00749 00750 Range tmp_file_ids; 00751 if (!subset_list[i].num_tag_values) { 00752 rval = get_tagged_entities( tag_index, tmp_file_ids ); 00753 } 00754 else { 00755 std::vector<int> ids( subset_list[i].tag_values, 00756 subset_list[i].tag_values + subset_list[i].num_tag_values ); 00757 std::sort( ids.begin(), ids.end() ); 00758 rval = search_tag_values( tag_index, ids, tmp_file_ids ); 00759 if (MB_SUCCESS != rval) 00760 return error(rval); 00761 } 00762 00763 if (tmp_file_ids.empty()) 00764 return error(MB_ENTITY_NOT_FOUND); 00765 00766 if (i == 0) 00767 file_ids.swap( tmp_file_ids ); 00768 else 00769 file_ids = intersect( tmp_file_ids, file_ids ); 00770 } 00771 00772 return MB_SUCCESS; 00773 } 00774 00775 ErrorCode ReadHDF5::get_partition( Range& tmp_file_ids, int num_parts, int part_number ) 00776 { 00777 00778 CHECK_OPEN_HANDLES; 00779 00780 // check that the tag only identified sets 00781 if ((unsigned long)fileInfo->sets.start_id > tmp_file_ids.front()) { 00782 dbgOut.print(2,"Ignoreing non-set entities with partition set tag\n"); 00783 tmp_file_ids.erase( tmp_file_ids.begin(), 00784 tmp_file_ids.lower_bound( 00785 (EntityHandle)fileInfo->sets.start_id ) ); 00786 } 00787 unsigned long set_end = (unsigned long)fileInfo->sets.start_id + fileInfo->sets.count; 00788 if (tmp_file_ids.back() >= set_end) { 00789 dbgOut.print(2,"Ignoreing non-set entities with partition set tag\n"); 00790 tmp_file_ids.erase( tmp_file_ids.upper_bound( (EntityHandle)set_end ), 00791 tmp_file_ids.end() ); 00792 } 00793 00794 Range::iterator s = tmp_file_ids.begin(); 00795 size_t num_per_proc = tmp_file_ids.size() / num_parts; 00796 size_t num_extra = tmp_file_ids.size() % num_parts; 00797 Range::iterator e; 00798 if (part_number < (long)num_extra) { 00799 s += (num_per_proc+1) * part_number; 00800 e = s; 00801 e += (num_per_proc+1); 00802 } 00803 else { 00804 s += num_per_proc * part_number + num_extra; 00805 e = s; 00806 e += num_per_proc; 00807 } 00808 tmp_file_ids.erase(e, tmp_file_ids.end()); 00809 tmp_file_ids.erase(tmp_file_ids.begin(), s); 00810 00811 return MB_SUCCESS; 00812 } 00813 00814 00815 ErrorCode ReadHDF5::load_file_partial( const ReaderIface::IDTag* subset_list, 00816 int subset_list_length, 00817 int num_parts, 00818 int part_number, 00819 const FileOptions& opts ) 00820 { 00821 mhdf_Status status; 00822 00823 static MPEState mpe_event( "ReadHDF5", "yellow" ); 00824 00825 mpe_event.start( "gather parts" ); 00826 00827 CHECK_OPEN_HANDLES; 00828 00829 for (int i = 0; i < subset_list_length; ++i) { 00830 dbgOut.printf( 2, "Select by \"%s\" with num_tag_values = %d\n", 00831 subset_list[i].tag_name, subset_list[i].num_tag_values ); 00832 if (subset_list[i].num_tag_values) { 00833 assert(0 != subset_list[i].tag_values); 00834 dbgOut.printf( 2, " \"%s\" values = { %d", 00835 subset_list[i].tag_name, subset_list[i].tag_values[0] ); 00836 for (int j = 1; j < subset_list[i].num_tag_values; ++j) 00837 dbgOut.printf( 2, ", %d", subset_list[i].tag_values[j] ); 00838 dbgOut.printf(2," }\n"); 00839 } 00840 } 00841 if (num_parts) 00842 dbgOut.printf( 2, "Partition with num_parts = %d and part_number = %d\n", 00843 num_parts, part_number ); 00844 00845 dbgOut.tprint( 1, "RETRIEVING TAGGED ENTITIES\n" ); 00846 00847 Range file_ids; 00848 ErrorCode rval = get_subset_ids( subset_list, subset_list_length, file_ids ); 00849 if (MB_SUCCESS != rval) 00850 return error(rval); 00851 00852 if (num_parts) { 00853 rval = get_partition( file_ids, num_parts, part_number ); 00854 if (MB_SUCCESS != rval) 00855 return error(rval); 00856 } 00857 00858 dbgOut.print_ints( 4, "Set file IDs for partial read: ", file_ids ); 00859 mpe_event.end(); 00860 mpe_event.start( "gather related sets" ); 00861 dbgOut.tprint( 1, "GATHERING ADDITIONAL ENTITIES\n" ); 00862 00863 enum RecusiveSetMode { RSM_NONE, RSM_SETS, RSM_CONTENTS }; 00864 const char* const set_opts[] = { "NONE", "SETS", "CONTENTS", NULL }; 00865 int child_mode; 00866 rval = opts.match_option( "CHILDREN", set_opts, child_mode ); 00867 if (MB_ENTITY_NOT_FOUND == rval) 00868 child_mode = RSM_CONTENTS; 00869 else if (MB_SUCCESS != rval) { 00870 readUtil->report_error( "Invalid value for 'CHILDREN' option" ); 00871 return error(rval); 00872 } 00873 int content_mode; 00874 rval = opts.match_option( "SETS", set_opts, content_mode ); 00875 if (MB_ENTITY_NOT_FOUND == rval) 00876 content_mode = RSM_CONTENTS; 00877 else if (MB_SUCCESS != rval) { 00878 readUtil->report_error( "Invalid value for 'SETS' option" ); 00879 return error(rval); 00880 } 00881 00882 // If we want the contents of contained/child sets, 00883 // search for them now (before gathering the non-set contents 00884 // of the sets.) 00885 Range sets; 00886 intersect( fileInfo->sets, file_ids, sets ); 00887 if (content_mode == RSM_CONTENTS || child_mode == RSM_CONTENTS) { 00888 dbgOut.tprint( 1, " doing read_set_ids_recursive\n" ); 00889 rval = read_set_ids_recursive( sets, content_mode == RSM_CONTENTS, child_mode == RSM_CONTENTS ); 00890 if (MB_SUCCESS != rval) 00891 return error(rval); 00892 } 00893 00894 debug_barrier(); 00895 00896 // get elements and vertices contained in sets 00897 dbgOut.tprint( 1, " doing get_set_contents\n" ); 00898 rval = get_set_contents( sets, file_ids ); 00899 if (MB_SUCCESS != rval) 00900 return error(rval); 00901 00902 dbgOut.print_ints( 5, "File IDs for partial read: ", file_ids ); 00903 debug_barrier(); 00904 mpe_event.end(); 00905 dbgOut.tprint( 1, "GATHERING NODE IDS\n" ); 00906 00907 // Figure out the maximum dimension of entity to be read 00908 int max_dim = 0; 00909 for (int i = 0; i < fileInfo->num_elem_desc; ++i) { 00910 EntityType type = CN::EntityTypeFromName( fileInfo->elems[i].type ); 00911 if (type <= MBVERTEX || type >= MBENTITYSET) { 00912 assert( false ); // for debug code die for unknown element tyoes 00913 continue; // for release code, skip unknown element types 00914 } 00915 int dim = CN::Dimension(type); 00916 if (dim > max_dim) { 00917 Range subset; 00918 intersect( fileInfo->elems[i].desc, file_ids, subset ); 00919 if (!subset.empty()) 00920 max_dim = dim; 00921 } 00922 } 00923 #ifdef USE_MPI 00924 if (nativeParallel) { 00925 int send = max_dim; 00926 MPI_Allreduce( &send, &max_dim, 1, MPI_INT, MPI_MAX, *mpiComm ); 00927 } 00928 #endif 00929 00930 // if input contained any polyhedra, then need to get faces 00931 // of the polyhedra before the next loop because we need to 00932 // read said faces in that loop. 00933 for (int i = 0; i < fileInfo->num_elem_desc; ++i) { 00934 EntityType type = CN::EntityTypeFromName( fileInfo->elems[i].type ); 00935 if (type != MBPOLYHEDRON) 00936 continue; 00937 00938 debug_barrier(); 00939 dbgOut.print( 2, " Getting polyhedra faces\n" ); 00940 mpe_event.start( "reading connectivity for ", fileInfo->elems[i].handle ); 00941 00942 Range polyhedra; 00943 intersect( fileInfo->elems[i].desc, file_ids, polyhedra ); 00944 rval = read_elems( i, polyhedra, &file_ids ); 00945 mpe_event.end(rval); 00946 if (MB_SUCCESS != rval) 00947 return error(rval); 00948 } 00949 00950 // get node file ids for all elements 00951 Range nodes; 00952 intersect( fileInfo->nodes, file_ids, nodes ); 00953 for (int i = 0; i < fileInfo->num_elem_desc; ++i) { 00954 EntityType type = CN::EntityTypeFromName( fileInfo->elems[i].type ); 00955 if (type <= MBVERTEX || type >= MBENTITYSET) { 00956 assert( false ); // for debug code die for unknown element tyoes 00957 continue; // for release code, skip unknown element types 00958 } 00959 if (MBPOLYHEDRON == type) 00960 continue; 00961 00962 debug_barrier(); 00963 dbgOut.printf( 2, " Getting element node IDs for: %s\n", fileInfo->elems[i].handle ); 00964 00965 Range subset; 00966 intersect( fileInfo->elems[i].desc, file_ids, subset ); 00967 mpe_event.start( "reading connectivity for ", fileInfo->elems[i].handle ); 00968 00969 // If dimension is max_dim, then we can create the elements now 00970 // so we don't have to read the table again later (connectivity 00971 // will be fixed up after nodes are created when update_connectivity()) 00972 // is called. For elements of a smaller dimension, we just build 00973 // the node ID range now because a) we'll have to read the whole 00974 // connectivity table again later, and b) we don't want to worry 00975 // about accidentally creating multiple copies of the same element. 00976 if (CN::Dimension(type) == max_dim) 00977 rval = read_elems( i, subset, &nodes ); 00978 else 00979 rval = read_elems( i, subset, nodes ); 00980 mpe_event.end(rval); 00981 if (MB_SUCCESS != rval) 00982 return error(rval); 00983 } 00984 00985 debug_barrier(); 00986 mpe_event.start( "read coords" ); 00987 dbgOut.tprintf( 1, "READING NODE COORDINATES (%lu nodes in %lu selects)\n", 00988 (unsigned long)nodes.size(), (unsigned long)nodes.psize() ); 00989 00990 // Read node coordinates and create vertices in MOAB 00991 // NOTE: This populates the RangeMap with node file ids, 00992 // which is expected by read_node_adj_elems. 00993 rval = read_nodes( nodes ); 00994 mpe_event.end(rval); 00995 if (MB_SUCCESS != rval) 00996 return error(rval); 00997 00998 debug_barrier(); 00999 dbgOut.tprint( 1, "READING ELEMENTS\n" ); 01000 01001 // decide if we need to read additional elements 01002 enum SideMode { SM_EXPLICIT, SM_NODES, SM_SIDES }; 01003 int side_mode; 01004 const char* const options[] = { "EXPLICIT", "NODES", "SIDES", 0 }; 01005 rval = opts.match_option( "ELEMENTS", options, side_mode ); 01006 if (MB_ENTITY_NOT_FOUND == rval) { 01007 // If only nodes were specified, then default to "NODES", otherwise 01008 // default to "SIDES". 01009 if (0 == max_dim) 01010 side_mode = SM_NODES; 01011 else 01012 side_mode = SM_SIDES; 01013 } 01014 else if (MB_SUCCESS != rval) { 01015 readUtil->report_error( "Invalid value for 'ELEMENTS' option" ); 01016 return error(rval); 01017 } 01018 01019 if (side_mode == SM_SIDES /*ELEMENTS=SIDES*/ && max_dim == 0 /*node-based*/) { 01020 // Read elements until we find something. Once we find someting, 01021 // read only elements of the same dimension. NOTE: loop termination 01022 // criterion changes on both sides (max_dim can be changed in loop 01023 // body). 01024 for (int dim = 3; dim >= max_dim; --dim) { 01025 for (int i = 0; i < fileInfo->num_elem_desc; ++i) { 01026 EntityType type = CN::EntityTypeFromName( fileInfo->elems[i].type ); 01027 if (CN::Dimension(type) == dim) { 01028 debug_barrier(); 01029 dbgOut.tprintf( 2, " Reading node-adjacent elements for: %s\n", fileInfo->elems[i].handle ); 01030 mpe_event.start( "reading connectivity for ", fileInfo->elems[i].handle ); 01031 Range ents; 01032 rval = read_node_adj_elems( fileInfo->elems[i] ); 01033 mpe_event.end(rval); 01034 if (MB_SUCCESS != rval) 01035 return error(rval); 01036 if (!ents.empty()) 01037 max_dim = 3; 01038 } 01039 } 01040 } 01041 } 01042 01043 Range side_entities; 01044 if (side_mode != SM_EXPLICIT /*ELEMENTS=NODES || ELEMENTS=SIDES*/) { 01045 if (0 == max_dim) 01046 max_dim = 4; 01047 // now read any additional elements for which we've already read all 01048 // of the nodes. 01049 for (int dim = max_dim - 1; dim > 0; --dim) { 01050 for (int i = 0; i < fileInfo->num_elem_desc; ++i) { 01051 EntityType type = CN::EntityTypeFromName( fileInfo->elems[i].type ); 01052 if (CN::Dimension(type) == dim) { 01053 debug_barrier(); 01054 dbgOut.tprintf( 2, " Reading node-adjacent elements for: %s\n", fileInfo->elems[i].handle ); 01055 mpe_event.start( "reading connectivity for ", fileInfo->elems[i].handle ); 01056 rval = read_node_adj_elems( fileInfo->elems[i], &side_entities ); 01057 mpe_event.end(rval); 01058 if (MB_SUCCESS != rval) 01059 return error(rval); 01060 } 01061 } 01062 } 01063 } 01064 01065 // We need to do this here for polyhedra to be handled correctly. 01066 // We have to wait until the faces are read in the above code block, 01067 // but need to create the connectivity before doing update_connectivity, 01068 // which might otherwise delete polyhedra faces. 01069 debug_barrier(); 01070 dbgOut.tprint( 1, "UPDATING CONNECTIVITY ARRAYS FOR READ ELEMENTS\n" ); 01071 mpe_event.start( "updating connectivity for elements read before vertices"); 01072 rval = update_connectivity(); 01073 mpe_event.end(); 01074 if (MB_SUCCESS != rval) 01075 return error(rval); 01076 01077 01078 dbgOut.tprint( 1, "READING ADJACENCIES\n" ); 01079 for (int i = 0; i < fileInfo->num_elem_desc; ++i) { 01080 if (fileInfo->elems[i].have_adj && 01081 idMap.intersects( fileInfo->elems[i].desc.start_id, fileInfo->elems[i].desc.count )) { 01082 mpe_event.start( "reading adjacencies for ", fileInfo->elems[i].handle ); 01083 long len; 01084 hid_t th = mhdf_openAdjacency( filePtr, fileInfo->elems[i].handle, &len, &status ); 01085 if (is_error(status)) 01086 return error(MB_FAILURE); 01087 01088 rval = read_adjacencies( th, len ); 01089 mhdf_closeData( filePtr, th, &status ); 01090 mpe_event.end(rval); 01091 if (MB_SUCCESS != rval) 01092 return error(rval); 01093 } 01094 } 01095 01096 // If doing ELEMENTS=SIDES then we need to delete any entities 01097 // that we read that aren't actually sides (e.g. an interior face 01098 // that connects two disjoint portions of the part). Both 01099 // update_connectivity and reading of any explicit adjacencies must 01100 // happen before this. 01101 if (side_mode == SM_SIDES) { 01102 debug_barrier(); 01103 mpe_event.start( "cleaning up non-side lower-dim elements" ); 01104 dbgOut.tprint( 1, "CHECKING FOR AND DELETING NON-SIDE ELEMENTS\n" ); 01105 rval = delete_non_side_elements( side_entities ); 01106 mpe_event.end(rval); 01107 if (MB_SUCCESS != rval) 01108 return error(rval); 01109 } 01110 01111 debug_barrier(); 01112 dbgOut.tprint( 1, "READING SETS\n" ); 01113 01114 // If reading contained/child sets but not their contents then find 01115 // them now. If we were also reading their contents we would 01116 // have found them already. 01117 if (content_mode == RSM_SETS || child_mode == RSM_SETS) { 01118 dbgOut.tprint( 1, " doing read_set_ids_recursive\n" ); 01119 mpe_event.start( "finding recursively contained sets" ); 01120 rval = read_set_ids_recursive( sets, content_mode == RSM_SETS, child_mode == RSM_SETS ); 01121 mpe_event.end(rval); 01122 if (MB_SUCCESS != rval) 01123 return error(rval); 01124 } 01125 01126 dbgOut.tprint( 1, " doing find_sets_containing\n" ); 01127 mpe_event.start( "finding sets containing any read entities" ); 01128 01129 // decide whether to read set-containing parents 01130 bool read_set_containing_parents = true; 01131 std::string tmp_opt; 01132 rval = opts.get_option( "NO_SET_CONTAINING_PARENTS", tmp_opt ); 01133 if (MB_SUCCESS == rval) 01134 read_set_containing_parents = false; 01135 01136 // Append file IDs of sets containing any of the nodes or elements 01137 // we've read up to this point. 01138 rval = find_sets_containing( sets, read_set_containing_parents ); 01139 mpe_event.end(rval); 01140 if (MB_SUCCESS != rval) 01141 return error(rval); 01142 // Now actually read all set data and instantiate sets in MOAB. 01143 // Get any contained sets out of file_ids. 01144 mpe_event.start( "reading set contents/parents/children" ); 01145 EntityHandle first_set = fileInfo->sets.start_id; 01146 sets.merge( file_ids.lower_bound( first_set ), 01147 file_ids.lower_bound( first_set + fileInfo->sets.count ) ); 01148 dbgOut.tprint( 1, " doing read_sets\n" ); 01149 rval = read_sets( sets ); 01150 mpe_event.end(rval); 01151 if (MB_SUCCESS != rval) 01152 return error(rval); 01153 01154 dbgOut.tprint( 1, "READING TAGS\n" ); 01155 01156 for (int i = 0; i < fileInfo->num_tag_desc; ++i) { 01157 mpe_event.start( "reading tag: ", fileInfo->tags[i].name ); 01158 rval = read_tag( i ); 01159 if (MB_SUCCESS != rval) 01160 return error(rval); 01161 } 01162 01163 dbgOut.tprint( 1, "PARTIAL READ COMPLETE.\n" ); 01164 01165 return MB_SUCCESS; 01166 } 01167 01168 ErrorCode ReadHDF5::search_tag_values( int tag_index, 01169 const std::vector<int>& sorted_values, 01170 Range& file_ids, 01171 bool sets_only ) 01172 { 01173 ErrorCode rval; 01174 mhdf_Status status; 01175 std::vector<EntityHandle>::iterator iter; 01176 const mhdf_TagDesc& tag = fileInfo->tags[tag_index]; 01177 long size; 01178 long start_id; 01179 01180 CHECK_OPEN_HANDLES; 01181 01182 debug_barrier(); 01183 01184 // do dense data 01185 01186 hid_t table; 01187 const char* name; 01188 std::vector<EntityHandle> indices; 01189 // These are probably in order of dimension, so iterate 01190 // in reverse order to make Range insertions more efficient. 01191 std::vector<int> grp_indices( tag.dense_elem_indices, tag.dense_elem_indices+tag.num_dense_indices ); 01192 for (std::vector<int>::reverse_iterator i = grp_indices.rbegin(); i != grp_indices.rend(); ++i) 01193 { 01194 int idx = *i; 01195 if (idx == -2) { 01196 name = mhdf_set_type_handle(); 01197 start_id = fileInfo->sets.start_id; 01198 } 01199 else if (sets_only) { 01200 continue; 01201 } 01202 else if (idx == -1) { 01203 name = mhdf_node_type_handle(); 01204 start_id = fileInfo->nodes.start_id; 01205 } 01206 else { 01207 if (idx < 0 || idx >= fileInfo->num_elem_desc) 01208 return error(MB_FAILURE); 01209 name = fileInfo->elems[idx].handle; 01210 start_id = fileInfo->elems[idx].desc.start_id; 01211 } 01212 table = mhdf_openDenseTagData( filePtr, tag.name, name, &size, &status ); 01213 if (is_error(status)) 01214 return error(MB_FAILURE); 01215 rval = search_tag_values( table, size, sorted_values, indices ); 01216 mhdf_closeData( filePtr, table, &status ); 01217 if (MB_SUCCESS != rval || is_error(status)) 01218 return error(MB_FAILURE); 01219 // Convert from table indices to file IDs and add to result list 01220 std::sort( indices.begin(), indices.end(), std::greater<EntityHandle>() ); 01221 std::transform( indices.begin(), indices.end(), range_inserter(file_ids), 01222 std::bind1st( std::plus<long>(), start_id ) ); 01223 indices.clear(); 01224 } 01225 01226 if (!tag.have_sparse) 01227 return MB_SUCCESS; 01228 01229 // do sparse data 01230 01231 hid_t tables[2]; 01232 long junk; // redundant value for non-variable-length tags 01233 mhdf_openSparseTagData( filePtr, tag.name, &size, &junk, tables, &status ); 01234 if (is_error(status)) 01235 return error(MB_FAILURE); 01236 rval = search_tag_values( tables[1], size, sorted_values, indices ); 01237 mhdf_closeData( filePtr, tables[1], &status ); 01238 if (MB_SUCCESS != rval || is_error(status)) { 01239 mhdf_closeData( filePtr, tables[0], &status ); 01240 return error(MB_FAILURE); 01241 } 01242 // convert to ranges 01243 std::sort( indices.begin(), indices.end() ); 01244 std::vector<EntityHandle> ranges; 01245 iter = indices.begin(); 01246 while (iter != indices.end()) { 01247 ranges.push_back( *iter ); 01248 EntityHandle last = *iter; 01249 for (++iter; iter != indices.end() && (last + 1) == *iter; ++iter, ++last); 01250 ranges.push_back( last ); 01251 } 01252 // read file ids 01253 iter = ranges.begin(); 01254 unsigned long offset = 0; 01255 while (iter != ranges.end()) { 01256 long begin = *iter; ++iter; 01257 long end = *iter; ++iter; 01258 mhdf_readSparseTagEntitiesWithOpt( tables[0], begin, end - begin + 1, 01259 handleType, &indices[offset], indepIO, &status ); 01260 if (is_error(status)) { 01261 mhdf_closeData( filePtr, tables[0], &status ); 01262 return error(MB_FAILURE); 01263 } 01264 offset += end - begin + 1; 01265 } 01266 mhdf_closeData( filePtr, tables[0], &status ); 01267 if (is_error(status)) 01268 return error(MB_FAILURE); 01269 assert( offset == indices.size() ); 01270 std::sort( indices.begin(), indices.end() ); 01271 01272 if (sets_only) { 01273 iter = std::lower_bound( indices.begin(), indices.end(), 01274 (EntityHandle)(fileInfo->sets.start_id + fileInfo->sets.count) ); 01275 indices.erase( iter, indices.end() ); 01276 iter = std::lower_bound( indices.begin(), indices.end(), 01277 fileInfo->sets.start_id ); 01278 indices.erase( indices.begin(), iter ); 01279 } 01280 copy_sorted_file_ids( &indices[0], indices.size(), file_ids ); 01281 01282 return MB_SUCCESS; 01283 } 01284 01285 ErrorCode ReadHDF5::get_tagged_entities( int tag_index, Range& file_ids ) 01286 { 01287 const mhdf_TagDesc& tag = fileInfo->tags[tag_index]; 01288 01289 CHECK_OPEN_HANDLES; 01290 01291 // do dense data 01292 Range::iterator hint = file_ids.begin(); 01293 for (int i = 0; i < tag.num_dense_indices; ++i) 01294 { 01295 int idx = tag.dense_elem_indices[i]; 01296 mhdf_EntDesc* ents; 01297 if (idx == -2) 01298 ents = &fileInfo->sets; 01299 else if (idx == -1) 01300 ents = &fileInfo->nodes; 01301 else { 01302 if (idx < 0 || idx >= fileInfo->num_elem_desc) 01303 return error(MB_FAILURE); 01304 ents = &(fileInfo->elems[idx].desc); 01305 } 01306 01307 EntityHandle h = (EntityHandle)ents->start_id; 01308 hint = file_ids.insert( hint, h, h + ents->count - 1 ); 01309 } 01310 01311 if (!tag.have_sparse) 01312 return MB_SUCCESS; 01313 01314 // do sparse data 01315 01316 mhdf_Status status; 01317 hid_t tables[2]; 01318 long size, junk; 01319 mhdf_openSparseTagData( filePtr, tag.name, &size, &junk, tables, &status ); 01320 if (is_error(status)) 01321 return error(MB_FAILURE); 01322 mhdf_closeData( filePtr, tables[1], &status ); 01323 if (is_error(status)) { 01324 mhdf_closeData( filePtr, tables[0], &status ); 01325 return error(MB_FAILURE); 01326 } 01327 01328 hid_t file_type = H5Dget_type( tables[0] ); 01329 if (file_type < 0) 01330 return error(MB_FAILURE); 01331 01332 hint = file_ids.begin(); 01333 EntityHandle* buffer = reinterpret_cast<EntityHandle*>(dataBuffer); 01334 const long buffer_size = bufferSize / std::max(sizeof(EntityHandle),H5Tget_size(file_type)); 01335 long remaining = size, offset = 0; 01336 while (remaining) { 01337 long count = std::min( buffer_size, remaining ); 01338 assert_range( buffer, count ); 01339 mhdf_readSparseTagEntitiesWithOpt( *tables, offset, count, 01340 file_type, buffer, collIO, &status ); 01341 if (is_error(status)) { 01342 H5Tclose(file_type); 01343 mhdf_closeData( filePtr, *tables, &status ); 01344 return error(MB_FAILURE); 01345 } 01346 H5Tconvert( file_type, handleType, count, buffer, NULL, H5P_DEFAULT ); 01347 01348 std::sort( buffer, buffer + count ); 01349 for (long i = 0; i < count; ++i) 01350 hint = file_ids.insert( hint, buffer[i], buffer[i] ); 01351 01352 remaining -= count; 01353 offset += count; 01354 } 01355 01356 H5Tclose(file_type); 01357 mhdf_closeData( filePtr, *tables, &status ); 01358 if (is_error(status)) 01359 return error(MB_FAILURE); 01360 01361 return MB_SUCCESS; 01362 } 01363 01364 ErrorCode ReadHDF5::search_tag_values( hid_t tag_table, 01365 unsigned long table_size, 01366 const std::vector<int>& sorted_values, 01367 std::vector<EntityHandle>& value_indices ) 01368 { 01369 01370 debug_barrier(); 01371 01372 CHECK_OPEN_HANDLES; 01373 01374 mhdf_Status status; 01375 size_t chunk_size = bufferSize / sizeof(int); 01376 int * buffer = reinterpret_cast<int*>(dataBuffer); 01377 size_t remaining = table_size, offset = 0; 01378 while (remaining) { 01379 // Get a block of tag values 01380 size_t count = std::min( chunk_size, remaining ); 01381 assert_range( buffer, count ); 01382 mhdf_readTagValuesWithOpt( tag_table, offset, count, H5T_NATIVE_INT, buffer, collIO, &status ); 01383 if (is_error(status)) 01384 return error(MB_FAILURE); 01385 01386 // search tag values 01387 for (size_t i = 0; i < count; ++i) 01388 if (std::binary_search( sorted_values.begin(), sorted_values.end(), (int)buffer[i] )) 01389 value_indices.push_back( i + offset ); 01390 01391 offset += count; 01392 remaining -= count; 01393 } 01394 01395 return MB_SUCCESS; 01396 } 01397 01398 ErrorCode ReadHDF5::read_nodes( const Range& node_file_ids ) 01399 { 01400 ErrorCode rval; 01401 mhdf_Status status; 01402 const int dim = fileInfo->nodes.vals_per_ent; 01403 Range range; 01404 01405 CHECK_OPEN_HANDLES; 01406 01407 if (node_file_ids.empty()) 01408 return MB_SUCCESS; 01409 01410 int cdim; 01411 rval = iFace->get_dimension( cdim ); 01412 if (MB_SUCCESS != rval) 01413 return error(rval); 01414 01415 if (cdim < dim) 01416 { 01417 rval = iFace->set_dimension( dim ); 01418 if (MB_SUCCESS != rval) 01419 return error(rval); 01420 } 01421 01422 hid_t data_id = mhdf_openNodeCoordsSimple( filePtr, &status ); 01423 if (is_error(status)) 01424 return error(MB_FAILURE); 01425 01426 EntityHandle handle; 01427 std::vector<double*> arrays(dim); 01428 const size_t num_nodes = node_file_ids.size(); 01429 rval = readUtil->get_node_coords( dim, (int)num_nodes, 0, handle, arrays ); 01430 if (MB_SUCCESS != rval) 01431 { 01432 mhdf_closeData( filePtr, data_id, &status ); 01433 return error(rval); 01434 } 01435 01436 if (blockedCoordinateIO) { 01437 try { 01438 for (int d = 0; d < dim; ++d) { 01439 ReadHDF5Dataset reader( "blocked coords", data_id, nativeParallel, mpiComm, false ); 01440 reader.set_column( d ); 01441 reader.set_file_ids( node_file_ids, fileInfo->nodes.start_id, num_nodes, H5T_NATIVE_DOUBLE ); 01442 dbgOut.printf( 3, "Reading %lu chunks for coordinate dimension %d\n", reader.get_read_count(), d ); 01443 // should normally only have one read call, unless sparse nature 01444 // of file_ids caused reader to do something strange 01445 size_t count, offset = 0; 01446 int nn = 0; 01447 while (!reader.done()) { 01448 dbgOut.printf(3,"Reading chunk %d for dimension %d\n", ++nn, d ); 01449 reader.read( arrays[d]+offset, count ); 01450 offset += count; 01451 } 01452 if (offset != num_nodes) { 01453 mhdf_closeData( filePtr, data_id, &status ); 01454 assert(false); 01455 return MB_FAILURE; 01456 } 01457 } 01458 } 01459 catch (ReadHDF5Dataset::Exception) { 01460 mhdf_closeData( filePtr, data_id, &status ); 01461 return error(MB_FAILURE); 01462 } 01463 } 01464 else { // !blockedCoordinateIO 01465 double* buffer = (double*)dataBuffer; 01466 long chunk_size = bufferSize / (3*sizeof(double)); 01467 long coffset = 0; 01468 int nn = 0; 01469 try { 01470 ReadHDF5Dataset reader( "interleaved coords", data_id, nativeParallel, mpiComm, false ); 01471 reader.set_file_ids( node_file_ids, fileInfo->nodes.start_id, chunk_size, H5T_NATIVE_DOUBLE ); 01472 dbgOut.printf( 3, "Reading %lu chunks for coordinate coordinates\n", reader.get_read_count() ); 01473 while (!reader.done()) { 01474 dbgOut.tprintf(3,"Reading chunk %d of node coords\n", ++nn); 01475 01476 size_t count; 01477 reader.read( buffer, count ); 01478 01479 for (size_t i = 0; i < count; ++i) 01480 for (int d = 0; d < dim; ++d) 01481 arrays[d][coffset+i] = buffer[dim*i+d]; 01482 coffset += count; 01483 } 01484 } 01485 catch (ReadHDF5Dataset::Exception) { 01486 mhdf_closeData( filePtr, data_id, &status ); 01487 return error(MB_FAILURE); 01488 } 01489 } 01490 01491 dbgOut.print(3,"Closing node coordinate table\n"); 01492 mhdf_closeData( filePtr, data_id, &status ); 01493 for (int d = dim; d < cdim; ++d) 01494 memset( arrays[d], 0, num_nodes*sizeof(double) ); 01495 01496 dbgOut.printf(3,"Updating ID to handle map for %lu nodes\n", (unsigned long)node_file_ids.size()); 01497 return insert_in_id_map( node_file_ids, handle ); 01498 } 01499 01500 ErrorCode ReadHDF5::read_elems( int i ) 01501 { 01502 Range ids; 01503 ids.insert( fileInfo->elems[i].desc.start_id, 01504 fileInfo->elems[i].desc.start_id + fileInfo->elems[i].desc.count - 1); 01505 return read_elems( i, ids ); 01506 } 01507 01508 ErrorCode ReadHDF5::read_elems( int i, const Range& file_ids, Range* node_ids ) 01509 { 01510 if (fileInfo->elems[i].desc.vals_per_ent < 0) { 01511 if (node_ids != 0) // not implemented for version 3 format of poly data 01512 return error(MB_TYPE_OUT_OF_RANGE); 01513 return read_poly( fileInfo->elems[i], file_ids ); 01514 } 01515 else 01516 return read_elems( fileInfo->elems[i], file_ids, node_ids ); 01517 } 01518 01519 ErrorCode ReadHDF5::read_elems( const mhdf_ElemDesc& elems, const Range& file_ids, Range* node_ids ) 01520 { 01521 01522 CHECK_OPEN_HANDLES; 01523 01524 debug_barrier(); 01525 dbgOut.tprintf( 1, "READING %s CONNECTIVITY (%lu elems in %lu selects)\n", 01526 elems.handle, (unsigned long)file_ids.size(), (unsigned long)file_ids.psize() ); 01527 01528 ErrorCode rval = MB_SUCCESS; 01529 mhdf_Status status; 01530 01531 EntityType type = CN::EntityTypeFromName( elems.type ); 01532 if (type == MBMAXTYPE) 01533 { 01534 readUtil->report_error( "Unknown element type: \"%s\".\n", elems.type ); 01535 return error(MB_FAILURE); 01536 } 01537 01538 const int nodes_per_elem = elems.desc.vals_per_ent; 01539 const size_t count = file_ids.size(); 01540 hid_t data_id = mhdf_openConnectivitySimple( filePtr, elems.handle, &status ); 01541 if (is_error(status)) 01542 return error(MB_FAILURE); 01543 01544 EntityHandle handle; 01545 EntityHandle* array = 0; 01546 if (count>0) 01547 rval = readUtil->get_element_connect( count, nodes_per_elem, type, 01548 0, handle, array ); 01549 if (MB_SUCCESS != rval) 01550 return error(rval); 01551 01552 try { 01553 EntityHandle* buffer = reinterpret_cast<EntityHandle*>(dataBuffer); 01554 const size_t buffer_size = bufferSize/(sizeof(EntityHandle)*nodes_per_elem); 01555 ReadHDF5Dataset reader( elems.handle, data_id, nativeParallel, mpiComm ); 01556 reader.set_file_ids( file_ids, elems.desc.start_id, buffer_size, handleType ); 01557 dbgOut.printf( 3, "Reading connectivity in %lu chunks for element group \"%s\"\n", 01558 reader.get_read_count(), elems.handle ); 01559 EntityHandle* iter = array; 01560 int nn = 0; 01561 while (!reader.done()) { 01562 dbgOut.printf( 3, "Reading chunk %d for \"%s\"\n", ++nn, elems.handle ); 01563 01564 size_t num_read; 01565 reader.read( buffer, num_read ); 01566 iter = std::copy( buffer, buffer+num_read*nodes_per_elem, iter ); 01567 01568 if (node_ids) { 01569 std::sort( buffer, buffer + num_read*nodes_per_elem ); 01570 num_read = std::unique( buffer, buffer + num_read*nodes_per_elem ) - buffer; 01571 copy_sorted_file_ids( buffer, num_read, *node_ids ); 01572 } 01573 } 01574 assert(iter - array == (ptrdiff_t)count * nodes_per_elem); 01575 } 01576 catch (ReadHDF5Dataset::Exception) { 01577 return error(MB_FAILURE); 01578 } 01579 01580 if (!node_ids) { 01581 rval = convert_id_to_handle( array, count*nodes_per_elem ); 01582 if (MB_SUCCESS != rval) 01583 return error(rval); 01584 01585 rval = readUtil->update_adjacencies( handle, count, nodes_per_elem, array ); 01586 if (MB_SUCCESS != rval) 01587 return error(rval); 01588 } 01589 else { 01590 IDConnectivity t; 01591 t.handle = handle; 01592 t.count = count; 01593 t.nodes_per_elem = nodes_per_elem; 01594 t.array = array; 01595 idConnectivityList.push_back(t); 01596 } 01597 01598 return insert_in_id_map( file_ids, handle ); 01599 } 01600 01601 ErrorCode ReadHDF5::update_connectivity() 01602 { 01603 ErrorCode rval; 01604 std::vector<IDConnectivity>::iterator i; 01605 for (i = idConnectivityList.begin(); i != idConnectivityList.end(); ++i) { 01606 rval = convert_id_to_handle( i->array, i->count * i->nodes_per_elem ); 01607 if (MB_SUCCESS != rval) 01608 return error(rval); 01609 01610 rval = readUtil->update_adjacencies( i->handle, i->count, i->nodes_per_elem, i->array ); 01611 if (MB_SUCCESS != rval) 01612 return error(rval); 01613 } 01614 idConnectivityList.clear(); 01615 return MB_SUCCESS; 01616 } 01617 01618 ErrorCode ReadHDF5::read_node_adj_elems( const mhdf_ElemDesc& group, Range* handles_out ) 01619 { 01620 mhdf_Status status; 01621 ErrorCode rval; 01622 01623 CHECK_OPEN_HANDLES; 01624 01625 hid_t table = mhdf_openConnectivitySimple( filePtr, group.handle, &status ); 01626 if (is_error(status)) 01627 return error(MB_FAILURE); 01628 01629 rval = read_node_adj_elems( group, table, handles_out ); 01630 01631 mhdf_closeData( filePtr, table, &status ); 01632 if (MB_SUCCESS == rval && is_error(status)) 01633 return error(rval = MB_FAILURE); 01634 return rval; 01635 } 01636 01637 ErrorCode ReadHDF5::read_node_adj_elems( const mhdf_ElemDesc& group, 01638 hid_t table_handle, 01639 Range* handles_out ) 01640 { 01641 01642 CHECK_OPEN_HANDLES; 01643 01644 debug_barrier(); 01645 01646 mhdf_Status status; 01647 ErrorCode rval; 01648 IODebugTrack debug_track( debugTrack, std::string(group.handle) ); 01649 01650 // copy data to local variables (makes other code clearer) 01651 const int node_per_elem = group.desc.vals_per_ent; 01652 long start_id = group.desc.start_id; 01653 long remaining = group.desc.count; 01654 const EntityType type = CN::EntityTypeFromName( group.type ); 01655 01656 // figure out how many elements we can read in each pass 01657 long* const buffer = reinterpret_cast<long*>( dataBuffer ); 01658 const long buffer_size = bufferSize / (node_per_elem * sizeof(buffer[0])); 01659 // read all element connectivity in buffer_size blocks 01660 long offset = 0; 01661 dbgOut.printf( 3, "Reading node-adjacent elements from \"%s\" in %ld chunks\n", 01662 group.handle, (remaining + buffer_size - 1) / buffer_size ); 01663 int nn = 0; 01664 Range::iterator hint; 01665 if (handles_out) 01666 hint = handles_out->begin(); 01667 while (remaining) { 01668 dbgOut.printf( 3, "Reading chunk %d of connectivity data for \"%s\"\n", ++nn, group.handle ); 01669 01670 // read a block of connectivity data 01671 const long count = std::min( remaining, buffer_size ); 01672 debug_track.record_io( offset, count ); 01673 assert_range( buffer, count*node_per_elem ); 01674 mhdf_readConnectivityWithOpt( table_handle, offset, count, H5T_NATIVE_LONG, buffer, collIO, &status ); 01675 if (is_error(status)) 01676 return error(MB_FAILURE); 01677 offset += count; 01678 remaining -= count; 01679 01680 // count the number of elements in the block that we want, 01681 // zero connectivity for other elements 01682 long num_elem = 0; 01683 long* iter = buffer; 01684 for (long i = 0; i < count; ++i) { 01685 for (int j = 0; j < node_per_elem; ++j) { 01686 iter[j] = (long)idMap.find( iter[j] ); 01687 if (!iter[j]) { 01688 iter[0] = 0; 01689 break; 01690 } 01691 } 01692 if (iter[0]) 01693 ++num_elem; 01694 iter += node_per_elem; 01695 } 01696 01697 if (!num_elem) { 01698 start_id += count; 01699 continue; 01700 } 01701 01702 // create elements 01703 EntityHandle handle; 01704 EntityHandle* array; 01705 rval = readUtil->get_element_connect( (int)num_elem, 01706 node_per_elem, 01707 type, 01708 0, 01709 handle, 01710 array ); 01711 if (MB_SUCCESS != rval) 01712 return error(rval); 01713 01714 // copy all non-zero connectivity values 01715 iter = buffer; 01716 EntityHandle* iter2 = array; 01717 EntityHandle h = handle; 01718 for (long i = 0; i < count; ++i) { 01719 if (!*iter) { 01720 iter += node_per_elem; 01721 continue; 01722 } 01723 if (!idMap.insert( start_id + i, h++, 1 ).second) 01724 return error(MB_FAILURE); 01725 01726 long* const end = iter + node_per_elem; 01727 for (; iter != end; ++iter, ++iter2) 01728 *iter2 = (EntityHandle)*iter; 01729 } 01730 assert( iter2 - array == num_elem * node_per_elem ); 01731 start_id += count; 01732 01733 rval = readUtil->update_adjacencies( handle, num_elem, node_per_elem, array ); 01734 if (MB_SUCCESS != rval) return error(rval); 01735 if (handles_out) 01736 hint = handles_out->insert( hint, handle, handle + num_elem - 1 ); 01737 } 01738 01739 debug_track.all_reduce(); 01740 return MB_SUCCESS; 01741 } 01742 01743 01744 ErrorCode ReadHDF5::read_elems( int i, const Range& elems_in, Range& nodes ) 01745 { 01746 01747 CHECK_OPEN_HANDLES; 01748 01749 debug_barrier(); 01750 dbgOut.tprintf( 1, "READING %s CONNECTIVITY (%lu elems in %lu selects)\n", fileInfo->elems[i].handle, 01751 (unsigned long)elems_in.size(), (unsigned long)elems_in.psize() ); 01752 01753 EntityHandle* const buffer = reinterpret_cast<EntityHandle*>(dataBuffer); 01754 const int node_per_elem = fileInfo->elems[i].desc.vals_per_ent; 01755 const size_t buffer_size = bufferSize / (node_per_elem*sizeof(EntityHandle)); 01756 01757 if (elems_in.empty()) 01758 return MB_SUCCESS; 01759 01760 assert( (long)elems_in.front() >= fileInfo->elems[i].desc.start_id ); 01761 assert( (long)elems_in.back() - fileInfo->elems[i].desc.start_id < fileInfo->elems[i].desc.count ); 01762 01763 // we don't support version 3 style poly element data 01764 if (fileInfo->elems[i].desc.vals_per_ent <= 0) 01765 return error(MB_TYPE_OUT_OF_RANGE); 01766 01767 mhdf_Status status; 01768 hid_t table = mhdf_openConnectivitySimple( filePtr, fileInfo->elems[i].handle, &status ); 01769 if (is_error(status)) 01770 return error(MB_FAILURE); 01771 01772 try { 01773 ReadHDF5Dataset reader( fileInfo->elems[i].handle, table, nativeParallel, mpiComm ); 01774 reader.set_file_ids( elems_in, fileInfo->elems[i].desc.start_id, 01775 buffer_size, handleType ); 01776 dbgOut.printf( 3, "Reading node list in %lu chunks for \"%s\"\n", reader.get_read_count(), fileInfo->elems[i].handle ); 01777 int nn = 0; 01778 while (!reader.done()) { 01779 dbgOut.printf( 3, "Reading chunk %d of \"%s\" connectivity\n", ++nn, fileInfo->elems[i].handle ); 01780 size_t num_read; 01781 reader.read( buffer, num_read ); 01782 std::sort( buffer, buffer + num_read*node_per_elem ); 01783 num_read = std::unique( buffer, buffer + num_read*node_per_elem ) - buffer; 01784 copy_sorted_file_ids( buffer, num_read, nodes ); 01785 } 01786 } 01787 catch (ReadHDF5Dataset::Exception) { 01788 return error(MB_FAILURE); 01789 } 01790 01791 return MB_SUCCESS; 01792 } 01793 01794 01795 ErrorCode ReadHDF5::read_poly( const mhdf_ElemDesc& elems, const Range& file_ids ) 01796 { 01797 class PolyReader : public ReadHDF5VarLen { 01798 private: 01799 const EntityType type; 01800 ReadHDF5* readHDF5; 01801 public: 01802 PolyReader( EntityType elem_type, void* buffer, size_t buffer_size, 01803 ReadHDF5* owner, DebugOutput& dbg ) 01804 : ReadHDF5VarLen( dbg, buffer, buffer_size ), 01805 type(elem_type), readHDF5(owner) 01806 {} 01807 virtual ~PolyReader() {} 01808 ErrorCode store_data( EntityHandle file_id, void* data, long len, bool ) 01809 { 01810 size_t valid; 01811 EntityHandle* conn = reinterpret_cast<EntityHandle*>(data); 01812 readHDF5->convert_id_to_handle( conn, len, valid ); 01813 if (valid != (size_t)len) 01814 return error(MB_ENTITY_NOT_FOUND); 01815 EntityHandle handle; 01816 ErrorCode rval = readHDF5->moab()->create_element( type, conn, len, handle ); 01817 if (MB_SUCCESS != rval) 01818 return error(rval); 01819 01820 rval = readHDF5->insert_in_id_map( file_id, handle ); 01821 return rval; 01822 } 01823 }; 01824 01825 CHECK_OPEN_HANDLES; 01826 01827 debug_barrier(); 01828 01829 EntityType type = CN::EntityTypeFromName( elems.type ); 01830 if (type == MBMAXTYPE) 01831 { 01832 readUtil->report_error( "Unknown element type: \"%s\".\n", elems.type ); 01833 return error(MB_FAILURE); 01834 } 01835 01836 hid_t handles[2]; 01837 mhdf_Status status; 01838 long num_poly, num_conn, first_id; 01839 mhdf_openPolyConnectivity( filePtr, elems.handle, &num_poly, &num_conn, &first_id, 01840 handles, &status ); 01841 if (is_error(status)) 01842 return error(MB_FAILURE); 01843 01844 std::string nm(elems.handle); 01845 ReadHDF5Dataset offset_reader( (nm + " offsets").c_str(), handles[0], nativeParallel, mpiComm, true ); 01846 ReadHDF5Dataset connect_reader( (nm + " data").c_str(), handles[1], nativeParallel, mpiComm, true ); 01847 01848 PolyReader tool( type, dataBuffer, bufferSize, this, dbgOut ); 01849 return tool.read( offset_reader, connect_reader, file_ids, first_id, handleType ); 01850 } 01851 01852 01853 ErrorCode ReadHDF5::delete_non_side_elements( const Range& side_ents ) 01854 { 01855 ErrorCode rval; 01856 01857 // build list of entities that we need to find the sides of 01858 Range explicit_ents; 01859 Range::iterator hint = explicit_ents.begin(); 01860 for (IDMap::iterator i = idMap.begin(); i != idMap.end(); ++i) { 01861 EntityHandle start = i->value; 01862 EntityHandle end = i->value + i->count - 1; 01863 EntityType type = TYPE_FROM_HANDLE(start); 01864 assert( type == TYPE_FROM_HANDLE(end) ); // otherwise handle space entirely full!! 01865 if (type != MBVERTEX && type != MBENTITYSET) 01866 hint = explicit_ents.insert( hint, start, end ); 01867 } 01868 explicit_ents = subtract( explicit_ents, side_ents ); 01869 01870 // figure out which entities we want to delete 01871 Range dead_ents( side_ents ); 01872 Range::iterator ds, de, es; 01873 ds = dead_ents.lower_bound( CN::TypeDimensionMap[1].first ); 01874 de = dead_ents.lower_bound( CN::TypeDimensionMap[2].first, ds ); 01875 if (ds != de) { 01876 // get subset of explicit ents of dimension greater than 1 01877 es = explicit_ents.lower_bound( CN::TypeDimensionMap[2].first ); 01878 Range subset, adj; 01879 subset.insert( es, explicit_ents.end() ); 01880 rval = iFace->get_adjacencies( subset, 1, false, adj, Interface::UNION ); 01881 if (MB_SUCCESS != rval) 01882 return rval; 01883 dead_ents = subtract( dead_ents, adj ); 01884 } 01885 ds = dead_ents.lower_bound( CN::TypeDimensionMap[2].first ); 01886 de = dead_ents.lower_bound( CN::TypeDimensionMap[3].first, ds ); 01887 assert(de == dead_ents.end()); 01888 if (ds != de) { 01889 // get subset of explicit ents of dimension 3 01890 es = explicit_ents.lower_bound( CN::TypeDimensionMap[3].first ); 01891 Range subset, adj; 01892 subset.insert( es, explicit_ents.end() ); 01893 rval = iFace->get_adjacencies( subset, 2, false, adj, Interface::UNION ); 01894 if (MB_SUCCESS != rval) 01895 return rval; 01896 dead_ents = subtract( dead_ents, adj ); 01897 } 01898 01899 // now delete anything remaining in dead_ents 01900 dbgOut.printf( 2, "Deleting %lu elements\n", (unsigned long)dead_ents.size() ); 01901 dbgOut.print( 4, "\tDead entities: ", dead_ents ); 01902 rval = iFace->delete_entities( dead_ents ); 01903 if (MB_SUCCESS != rval) 01904 return error(rval); 01905 01906 // remove dead entities from ID map 01907 while (!dead_ents.empty()) { 01908 EntityHandle start = dead_ents.front(); 01909 EntityID count = dead_ents.const_pair_begin()->second - start + 1; 01910 IDMap::iterator rit; 01911 for (rit = idMap.begin(); rit != idMap.end(); ++rit) 01912 if (rit->value <= start && (long)(start - rit->value) < rit->count) 01913 break; 01914 if (rit == idMap.end()) 01915 return error(MB_FAILURE); 01916 01917 EntityID offset = start - rit->value; 01918 EntityID avail = rit->count - offset; 01919 if (avail < count) 01920 count = avail; 01921 01922 dead_ents.erase( dead_ents.begin(), dead_ents.begin() + count ); 01923 idMap.erase( rit->begin + offset, count ); 01924 } 01925 01926 return MB_SUCCESS; 01927 } 01928 01929 ErrorCode ReadHDF5::read_sets( const Range& file_ids ) 01930 { 01931 01932 CHECK_OPEN_HANDLES; 01933 01934 debug_barrier(); 01935 01936 mhdf_Status status; 01937 ErrorCode rval; 01938 01939 const size_t num_sets = fileInfo->sets.count; 01940 if (!num_sets) // If no sets at all! 01941 return MB_SUCCESS; 01942 01943 // create sets 01944 std::vector<unsigned> flags(file_ids.size()); 01945 Range::iterator si = file_ids.begin(); 01946 for (size_t i = 0; i < flags.size(); ++i, ++si) 01947 flags[i] = setMeta[*si - fileInfo->sets.start_id][3] & ~(long)mhdf_SET_RANGE_BIT; 01948 EntityHandle start_handle; 01949 rval = readUtil->create_entity_sets( flags.size(), &flags[0], 0, start_handle ); 01950 if (MB_SUCCESS != rval) 01951 return error(rval); 01952 rval = insert_in_id_map( file_ids, start_handle ); 01953 if (MB_SUCCESS != rval) 01954 return error(rval); 01955 01956 // read contents 01957 if (fileInfo->have_set_contents) { 01958 long len = 0; 01959 hid_t handle = mhdf_openSetData( filePtr, &len, &status ); 01960 if (is_error(status)) 01961 return error(MB_FAILURE); 01962 01963 ReadHDF5Dataset dat( "set contents", handle, nativeParallel, mpiComm, true ); 01964 rval = read_set_data( file_ids, start_handle, dat, CONTENT ); 01965 if (MB_SUCCESS != rval) 01966 return error(rval); 01967 } 01968 01969 // read set child lists 01970 if (fileInfo->have_set_children) { 01971 long len = 0; 01972 hid_t handle = mhdf_openSetChildren( filePtr, &len, &status ); 01973 if (is_error(status)) 01974 return error(MB_FAILURE); 01975 01976 ReadHDF5Dataset dat( "set children", handle, nativeParallel, mpiComm, true ); 01977 rval = read_set_data( file_ids, start_handle, dat, CHILD ); 01978 if (MB_SUCCESS != rval) 01979 return error(rval); 01980 } 01981 01982 // read set parent lists 01983 if (fileInfo->have_set_parents) { 01984 long len = 0; 01985 hid_t handle = mhdf_openSetParents( filePtr, &len, &status ); 01986 if (is_error(status)) 01987 return error(MB_FAILURE); 01988 01989 ReadHDF5Dataset dat( "set parents", handle, nativeParallel, mpiComm, true ); 01990 rval = read_set_data( file_ids, start_handle, dat, PARENT ); 01991 if (MB_SUCCESS != rval) 01992 return error(rval); 01993 } 01994 01995 return MB_SUCCESS; 01996 } 01997 01998 ErrorCode ReadHDF5::read_all_set_meta() 01999 { 02000 CHECK_OPEN_HANDLES; 02001 02002 assert(!setMeta); 02003 const long num_sets = fileInfo->sets.count; 02004 if (!num_sets) 02005 return MB_SUCCESS; 02006 02007 mhdf_Status status; 02008 hid_t handle = mhdf_openSetMetaSimple( filePtr, &status ); 02009 if (is_error(status)) { 02010 return error(MB_FAILURE); 02011 } 02012 02013 // Allocate extra space if we need it for data conversion 02014 hid_t meta_type = H5Dget_type( handle ); 02015 size_t size = H5Tget_size( meta_type ); 02016 if (size > sizeof(long)) 02017 setMeta = new long[(num_sets * size + (sizeof(long)-1)) / sizeof(long)][4]; 02018 else 02019 setMeta = new long[num_sets][4]; 02020 02021 // set some parameters based on whether or not each proc reads the 02022 // table or only the root reads it and bcasts it to the others 02023 int rank = 0; 02024 bool bcast = false; 02025 hid_t ioprop = H5P_DEFAULT; 02026 #ifdef USE_MPI 02027 MPI_Comm comm = 0; 02028 if (nativeParallel) { 02029 rank = myPcomm->proc_config().proc_rank(); 02030 comm = myPcomm->proc_config().proc_comm(); 02031 bcast = bcastDuplicateReads; 02032 if (!bcast) 02033 ioprop = collIO; 02034 } 02035 #endif 02036 02037 if (!bcast || 0 == rank) { 02038 mhdf_readSetMetaWithOpt( handle, 0, num_sets, meta_type, setMeta, ioprop, &status ); 02039 if (is_error(status)) { 02040 H5Tclose( meta_type ); 02041 mhdf_closeData( filePtr, handle, &status ); 02042 return error(MB_FAILURE); 02043 } 02044 02045 H5Tconvert( meta_type, H5T_NATIVE_LONG, num_sets*4, setMeta, 0, H5P_DEFAULT ); 02046 } 02047 mhdf_closeData( filePtr, handle, &status ); 02048 if (is_error(status)) 02049 return error(MB_FAILURE); 02050 H5Tclose( meta_type ); 02051 02052 if (bcast) { 02053 #ifdef USE_MPI 02054 int ierr = MPI_Bcast( setMeta, num_sets*4, MPI_LONG, 0, comm ); 02055 if (MPI_SUCCESS != ierr) 02056 return error(MB_FAILURE); 02057 #else 02058 assert(rank == 0); // if not MPI, then only one proc 02059 #endif 02060 } 02061 02062 return MB_SUCCESS; 02063 } 02064 02065 02066 ErrorCode ReadHDF5::read_set_ids_recursive( Range& sets_in_out, 02067 bool contained_sets, 02068 bool child_sets ) 02069 { 02070 02071 CHECK_OPEN_HANDLES; 02072 mhdf_Status status; 02073 02074 if (!fileInfo->have_set_children) 02075 child_sets = false; 02076 if (!fileInfo->have_set_contents) 02077 contained_sets = false; 02078 if (!child_sets && !contained_sets) 02079 return MB_SUCCESS; 02080 02081 // open data tables 02082 if (fileInfo->sets.count == 0) { 02083 assert( sets_in_out.empty() ); 02084 return MB_SUCCESS; 02085 } 02086 02087 if (!contained_sets && !child_sets) 02088 return MB_SUCCESS; 02089 02090 ReadHDF5Dataset cont( "set contents", false, mpiComm ); 02091 ReadHDF5Dataset child( "set children", false, mpiComm ); 02092 02093 if (contained_sets) { 02094 long content_len = 0; 02095 hid_t content_handle = mhdf_openSetData( filePtr, &content_len, &status ); 02096 if (is_error(status)) 02097 return error(MB_FAILURE); 02098 try { 02099 cont.init( content_handle, true ); 02100 } 02101 catch ( ReadHDF5Dataset::Exception ) { 02102 return error(MB_FAILURE); 02103 } 02104 } 02105 02106 if (child_sets) { 02107 long child_len = 0; 02108 hid_t child_handle = mhdf_openSetChildren( filePtr, &child_len, &status ); 02109 if (is_error(status)) 02110 return error(MB_FAILURE); 02111 try { 02112 child.init( child_handle, true ); 02113 } 02114 catch ( ReadHDF5Dataset::Exception ) { 02115 return error(MB_FAILURE); 02116 } 02117 } 02118 02119 ErrorCode rval = MB_SUCCESS; 02120 Range children, new_children(sets_in_out); 02121 int iteration_count = 0; 02122 do { 02123 ++iteration_count; 02124 dbgOut.tprintf(2,"Iteration %d of read_set_ids_recursive\n",iteration_count); 02125 children.clear(); 02126 if (child_sets) { 02127 rval = read_set_data( new_children, 0, child, CHILD, &children ); 02128 if (MB_SUCCESS != rval) 02129 break; 02130 } 02131 if (contained_sets) { 02132 rval = read_set_data( new_children, 0, cont, CONTENT, &children ); 02133 // remove any non-set values 02134 Range::iterator it = children.lower_bound( fileInfo->sets.start_id ); 02135 children.erase( children.begin(), it ); 02136 it = children.lower_bound( fileInfo->sets.start_id + fileInfo->sets.count ); 02137 children.erase( it, children.end() ); 02138 if (MB_SUCCESS != rval) 02139 break; 02140 } 02141 new_children = subtract( children, sets_in_out ); 02142 dbgOut.print_ints( 2, "Adding additional contained/child sets", new_children ); 02143 sets_in_out.merge( new_children ); 02144 } while (!new_children.empty()); 02145 02146 return MB_SUCCESS; 02147 } 02148 02149 ErrorCode ReadHDF5::find_sets_containing( Range& sets_out, 02150 bool read_set_containing_parents) 02151 { 02152 ErrorCode rval; 02153 mhdf_Status status; 02154 02155 CHECK_OPEN_HANDLES; 02156 02157 if (!fileInfo->have_set_contents) 02158 return MB_SUCCESS; 02159 assert( fileInfo->sets.count ); 02160 02161 // open data tables 02162 long content_len = 0; 02163 hid_t content_handle = mhdf_openSetData( filePtr, &content_len, &status ); 02164 if (is_error(status)) 02165 return error(MB_FAILURE); 02166 02167 hid_t data_type = H5Dget_type( content_handle ); 02168 02169 rval = find_sets_containing( content_handle, data_type, content_len, 02170 read_set_containing_parents, sets_out ); 02171 02172 H5Tclose( data_type ); 02173 02174 mhdf_closeData( filePtr, content_handle, &status ); 02175 if(MB_SUCCESS == rval && is_error(status)) 02176 return error(MB_FAILURE); 02177 02178 return rval; 02179 } 02180 02181 static bool set_map_intersect( bool ranged, 02182 const long* contents, 02183 int content_len, 02184 const RangeMap<long,EntityHandle>& id_map ) 02185 { 02186 if (ranged) { 02187 if (!content_len || id_map.empty()) 02188 return false; 02189 02190 const long* j = contents; 02191 const long* const end = contents + content_len; 02192 assert(content_len % 2 == 0); 02193 while (j != end) { 02194 long start = *(j++); 02195 long count = *(j++); 02196 if (id_map.intersects( start, count )) 02197 return true; 02198 } 02199 } 02200 else { 02201 const long* const end = contents + content_len; 02202 for (const long* i = contents; i != end; ++i) 02203 if (id_map.exists( *i )) 02204 return true; 02205 } 02206 return false; 02207 } 02208 02209 struct SetContOffComp { 02210 bool operator()(const long a1[4], const long a2) 02211 { return a1[ReadHDF5::CONTENT] < a2; } 02212 }; 02213 02214 02215 ErrorCode ReadHDF5::find_sets_containing( hid_t contents_handle, 02216 hid_t content_type, 02217 long contents_len, 02218 bool read_set_containing_parents, 02219 Range& file_ids ) 02220 { 02221 02222 CHECK_OPEN_HANDLES; 02223 02224 // Scan all set contents data 02225 02226 const size_t content_size = H5Tget_size( content_type ); 02227 const long num_sets = fileInfo->sets.count; 02228 dbgOut.printf( 2, "Searching contents of %ld\n", num_sets ); 02229 mhdf_Status status; 02230 02231 int rank = 0; 02232 bool bcast = false; 02233 #ifdef USE_MPI 02234 MPI_Comm comm = 0; 02235 if (nativeParallel) { 02236 rank = myPcomm->proc_config().proc_rank(); 02237 comm = myPcomm->proc_config().proc_comm(); 02238 bcast = bcastDuplicateReads; 02239 } 02240 #endif 02241 02242 // check offsets so that we don't read past end of table or 02243 // walk off end of array. 02244 long prev = -1; 02245 for (long i = 0; i < num_sets; ++i) { 02246 if (setMeta[i][CONTENT] < prev) { 02247 std::cerr << "Invalid data in set contents offsets at position " 02248 << i << ": index " << setMeta[i][CONTENT] 02249 << " is less than previous index " << prev << std::endl; 02250 std::cerr.flush(); 02251 return error(MB_FAILURE); 02252 } 02253 prev = setMeta[i][CONTENT]; 02254 } 02255 if (setMeta[num_sets-1][CONTENT] >= contents_len) { 02256 std::cerr << "Maximum set content index " << setMeta[num_sets-1][CONTENT] 02257 << " exceeds contents table length of " << contents_len 02258 << std::endl; 02259 std::cerr.flush(); 02260 return error(MB_FAILURE); 02261 } 02262 02263 // set up buffer for reading set contents 02264 long* const content_buffer = (long*)dataBuffer; 02265 const long content_len = bufferSize / std::max( content_size, sizeof(long) ); 02266 02267 // scan set table 02268 Range::iterator hint = file_ids.begin(); 02269 Range tmp_range; 02270 long prev_idx = -1; 02271 int mm = 0; 02272 long sets_offset = 0; 02273 while (sets_offset < num_sets) { 02274 long sets_count = std::lower_bound( setMeta + sets_offset, 02275 setMeta + num_sets, 02276 content_len + prev_idx, 02277 SetContOffComp() 02278 ) - setMeta - sets_offset; 02279 assert(sets_count >= 0 && sets_offset + sets_count <= num_sets); 02280 if (!sets_count) { // contents of single set don't fit in buffer 02281 long content_remaining = setMeta[sets_offset][CONTENT] - prev_idx; 02282 long content_offset = prev_idx+1; 02283 while (content_remaining) { 02284 long content_count = content_len < content_remaining ? 02285 2*(content_len/2) : content_remaining; 02286 assert_range( content_buffer, content_count ); 02287 dbgOut.printf( 3, "Reading chunk %d (%ld values) from set contents table\n", ++mm, content_count); 02288 if (!bcast || 0 == rank) { 02289 if (!bcast) 02290 mhdf_readSetDataWithOpt( contents_handle, content_offset, 02291 content_count, content_type, 02292 content_buffer, collIO, &status ); 02293 else 02294 mhdf_readSetData( contents_handle, content_offset, 02295 content_count, content_type, 02296 content_buffer, &status ); 02297 if (is_error(status)) 02298 return error(MB_FAILURE); 02299 02300 H5Tconvert( content_type, H5T_NATIVE_LONG, content_count, content_buffer, 0, H5P_DEFAULT ); 02301 } 02302 if (bcast) { 02303 #ifdef USE_MPI 02304 int ierr = MPI_Bcast( content_buffer, content_count, MPI_LONG, 0, comm ); 02305 if (MPI_SUCCESS != ierr) 02306 return error(MB_FAILURE); 02307 #else 02308 assert(rank == 0); // if not MPI, then only one proc 02309 #endif 02310 } 02311 02312 if (read_set_containing_parents) { 02313 tmp_range.clear(); 02314 if (setMeta[sets_offset][3] & mhdf_SET_RANGE_BIT) tmp_range.insert(*content_buffer, *(content_buffer+1)); 02315 else std::copy(content_buffer, content_buffer+content_count, range_inserter(tmp_range)); 02316 tmp_range = intersect(tmp_range, file_ids); 02317 } 02318 02319 if (!tmp_range.empty() || 02320 set_map_intersect( setMeta[sets_offset][3] & mhdf_SET_RANGE_BIT, 02321 content_buffer, content_count, idMap )) { 02322 long id = fileInfo->sets.start_id + sets_offset; 02323 hint = file_ids.insert( hint, id, id ); 02324 if (!nativeParallel) // don't stop if doing READ_PART because we need to read collectively 02325 break; 02326 } 02327 content_remaining -= content_count; 02328 content_offset += content_count; 02329 } 02330 prev_idx = setMeta[sets_offset][CONTENT]; 02331 sets_count = 1; 02332 } 02333 else if (long read_num = setMeta[sets_offset + sets_count - 1][CONTENT] - prev_idx) { 02334 assert(sets_count > 0); 02335 assert_range( content_buffer, read_num ); 02336 dbgOut.printf( 3, "Reading chunk %d (%ld values) from set contents table\n", ++mm, read_num); 02337 if (!bcast || 0 == rank) { 02338 if (!bcast) 02339 mhdf_readSetDataWithOpt( contents_handle, prev_idx+1, read_num, 02340 content_type, content_buffer, collIO, &status ); 02341 else 02342 mhdf_readSetData( contents_handle, prev_idx+1, read_num, 02343 content_type, content_buffer, &status ); 02344 if (is_error(status)) 02345 return error(MB_FAILURE); 02346 02347 H5Tconvert( content_type, H5T_NATIVE_LONG, read_num, content_buffer, 0, H5P_DEFAULT ); 02348 } 02349 if (bcast) { 02350 #ifdef USE_MPI 02351 int ierr = MPI_Bcast( content_buffer, read_num, MPI_LONG, 0, comm ); 02352 if (MPI_SUCCESS != ierr) 02353 return error(MB_FAILURE); 02354 #else 02355 assert(rank == 0); // if not MPI, then only one proc 02356 #endif 02357 } 02358 02359 long* buff_iter = content_buffer; 02360 for (long i = 0; i < sets_count; ++i) { 02361 long set_size = setMeta[i+sets_offset][CONTENT] - prev_idx; 02362 prev_idx += set_size; 02363 02364 // check whether contents include set already being loaded 02365 if (read_set_containing_parents) { 02366 tmp_range.clear(); 02367 if (setMeta[sets_offset+i][3] & mhdf_SET_RANGE_BIT) tmp_range.insert(*buff_iter, *(buff_iter+1)); 02368 else std::copy(buff_iter, buff_iter+set_size, range_inserter(tmp_range)); 02369 tmp_range = intersect(tmp_range, file_ids); 02370 } 02371 02372 if (!tmp_range.empty() || 02373 set_map_intersect( setMeta[sets_offset+i][3] & mhdf_SET_RANGE_BIT, 02374 buff_iter, set_size, idMap )) { 02375 long id = fileInfo->sets.start_id + sets_offset + i; 02376 hint = file_ids.insert( hint, id, id ); 02377 } 02378 buff_iter += set_size; 02379 } 02380 } 02381 02382 sets_offset += sets_count; 02383 } 02384 02385 return MB_SUCCESS; 02386 } 02387 02388 static Range::iterator copy_set_contents( 02389 Range::iterator hint, 02390 int ranged, 02391 EntityHandle* contents, 02392 long length, 02393 Range& results ) 02394 { 02395 if (ranged) { 02396 assert( length%2 == 0 ); 02397 for (long i = 0; i < length; i += 2) 02398 hint = results.insert( hint, contents[i], contents[i] + contents[i+1] - 1 ); 02399 } 02400 else { 02401 std::sort( contents, contents+length ); 02402 for (long i = 0; i < length; ++i) 02403 hint = results.insert( hint, contents[i] ); 02404 } 02405 return hint; 02406 } 02407 02408 ErrorCode ReadHDF5::read_set_data( const Range& set_file_ids, 02409 EntityHandle start_handle, 02410 ReadHDF5Dataset& data, 02411 SetMode mode, 02412 Range* file_ids_out ) 02413 { 02414 ErrorCode rval; 02415 Range::const_pair_iterator pi; 02416 Range::iterator out_hint; 02417 if (file_ids_out) 02418 out_hint = file_ids_out->begin(); 02419 02420 // Construct range of offsets into data table at which to read 02421 // Note: all offsets are incremented by TWEAK because Range cannot 02422 // store zeros. 02423 const long TWEAK = 1; 02424 Range data_offsets; 02425 Range::iterator hint = data_offsets.begin(); 02426 pi = set_file_ids.const_pair_begin(); 02427 if ((long)pi->first == fileInfo->sets.start_id) { 02428 long second = pi->second - fileInfo->sets.start_id; 02429 if (setMeta[second][mode] >= 0) 02430 hint = data_offsets.insert( hint, TWEAK, setMeta[second][mode]+TWEAK ); 02431 ++pi; 02432 } 02433 for ( ; pi != set_file_ids.const_pair_end(); ++pi) { 02434 long first = pi->first - fileInfo->sets.start_id; 02435 long second = pi->second - fileInfo->sets.start_id; 02436 long idx1 = setMeta[first-1][mode]+1; 02437 long idx2 = setMeta[second][mode]; 02438 if (idx2 >= idx1) 02439 hint = data_offsets.insert( hint, idx1+TWEAK, idx2+TWEAK ); 02440 } 02441 try { 02442 data.set_file_ids( data_offsets, TWEAK, bufferSize/sizeof(EntityHandle), handleType ); 02443 } 02444 catch (ReadHDF5Dataset::Exception ) { 02445 return MB_FAILURE; 02446 } 02447 02448 // we need to increment this for each processed set because 02449 // the sets were created in the order of the ids in file_ids. 02450 EntityHandle h = start_handle; 02451 02452 const long ranged_flag = (mode == CONTENT) ? mhdf_SET_RANGE_BIT : 0; 02453 02454 std::vector<EntityHandle> partial; // for when we read only part of the contents of a set/entity 02455 Range::const_iterator fileid_iter = set_file_ids.begin(); 02456 EntityHandle* buffer = reinterpret_cast<EntityHandle*>(dataBuffer); 02457 size_t count, offset; 02458 02459 int nn = 0; 02460 while (!data.done()) { 02461 dbgOut.printf( 3, "Reading chunk %d of %s\n", ++nn, data.get_debug_desc() ); 02462 try { 02463 data.read( buffer, count ); 02464 } 02465 catch (ReadHDF5Dataset::Exception ) { 02466 return MB_FAILURE; 02467 } 02468 02469 // assert not appropriate here - I might have treated all my file ids, but maybe 02470 // another proc hasn't; for me, count will be zero, so I won't do anything, but 02471 // I still need to go through the motions to make the read work 02472 02473 // Handle 'special' case where we read some, but not all 02474 // of the data for an entity during the last iteration. 02475 offset = 0; 02476 if (!partial.empty()) { // didn't read all of previous entity 02477 assert( fileid_iter != set_file_ids.end() ); 02478 size_t num_prev = partial.size(); 02479 size_t idx = *fileid_iter - fileInfo->sets.start_id; 02480 size_t len = idx ? setMeta[idx][mode] - setMeta[idx-1][mode] : setMeta[idx][mode] + 1; 02481 offset = len - num_prev; 02482 if (offset > count) { // still don't have all 02483 partial.insert( partial.end(), buffer, buffer+count ); 02484 continue; 02485 } 02486 02487 partial.insert( partial.end(), buffer, buffer+offset ); 02488 if (file_ids_out) { 02489 out_hint = copy_set_contents( out_hint, setMeta[idx][3] & ranged_flag, 02490 &partial[0], partial.size(), *file_ids_out); 02491 } 02492 else { 02493 switch (mode) { 02494 size_t valid; 02495 case CONTENT: 02496 if (setMeta[idx][3] & ranged_flag) { 02497 if (len % 2) 02498 return error(MB_INDEX_OUT_OF_RANGE); 02499 Range range; 02500 convert_range_to_handle( &partial[0], len/2, range ); 02501 rval = moab()->add_entities( h, range ); 02502 } 02503 else { 02504 convert_id_to_handle( &partial[0], len, valid ); 02505 rval = moab()->add_entities( h, &partial[0], valid ); 02506 } 02507 break; 02508 case CHILD: 02509 convert_id_to_handle( &partial[0], len, valid ); 02510 rval = moab()->add_child_meshsets( h, &partial[0], valid ); 02511 break; 02512 case PARENT: 02513 convert_id_to_handle( &partial[0], len, valid ); 02514 rval = moab()->add_parent_meshsets( h, &partial[0], valid ); 02515 break; 02516 } 02517 if (MB_SUCCESS != rval) 02518 return error(rval); 02519 } 02520 02521 ++fileid_iter; 02522 ++h; 02523 partial.clear(); 02524 } 02525 02526 // Process contents for all entities for which we 02527 // have read the complete list 02528 while (offset < count) { 02529 assert( fileid_iter != set_file_ids.end() ); 02530 size_t idx = *fileid_iter - fileInfo->sets.start_id; 02531 size_t len = idx ? setMeta[idx][mode] - setMeta[idx-1][mode] : setMeta[idx][mode] + 1; 02532 // If we did not read all of the final entity, 02533 // store what we did read to be processed in the 02534 // next iteration 02535 if (offset + len > count) { 02536 partial.insert( partial.end(), buffer + offset, buffer + count ); 02537 break; 02538 } 02539 02540 if (file_ids_out) { 02541 out_hint = copy_set_contents( out_hint, setMeta[idx][3] & ranged_flag, 02542 buffer + offset, len, *file_ids_out); 02543 } 02544 else { 02545 switch (mode) { 02546 size_t valid; 02547 case CONTENT: 02548 if (setMeta[idx][3] & ranged_flag) { 02549 if (len % 2) 02550 return error(MB_INDEX_OUT_OF_RANGE); 02551 Range range; 02552 convert_range_to_handle( buffer+offset, len/2, range ); 02553 rval = moab()->add_entities( h, range ); 02554 } 02555 else { 02556 convert_id_to_handle( buffer+offset, len, valid ); 02557 rval = moab()->add_entities( h, buffer+offset, valid ); 02558 } 02559 break; 02560 case CHILD: 02561 convert_id_to_handle( buffer+offset, len, valid ); 02562 rval = moab()->add_child_meshsets( h, buffer+offset, valid ); 02563 break; 02564 case PARENT: 02565 convert_id_to_handle( buffer+offset, len, valid ); 02566 rval = moab()->add_parent_meshsets( h, buffer+offset, valid ); 02567 break; 02568 } 02569 if (MB_SUCCESS != rval) 02570 return error(rval); 02571 } 02572 02573 ++fileid_iter; 02574 ++h; 02575 offset += len; 02576 } 02577 } 02578 02579 return MB_SUCCESS; 02580 } 02581 02582 ErrorCode ReadHDF5::get_set_contents( const Range& sets, Range& file_ids ) 02583 { 02584 02585 CHECK_OPEN_HANDLES; 02586 02587 if (!fileInfo->have_set_contents) 02588 return MB_SUCCESS; 02589 dbgOut.tprint(2,"Reading set contained file IDs\n"); 02590 try { 02591 mhdf_Status status; 02592 long content_len; 02593 hid_t contents = mhdf_openSetData( filePtr, &content_len, &status ); 02594 if (is_error(status)) 02595 return error(MB_FAILURE); 02596 ReadHDF5Dataset data( "set contents", contents, nativeParallel, mpiComm, true ); 02597 02598 02599 return read_set_data( sets, 0, data, CONTENT, &file_ids ); 02600 } 02601 catch (ReadHDF5Dataset::Exception) { 02602 return error(MB_FAILURE); 02603 } 02604 } 02605 02606 ErrorCode ReadHDF5::read_adjacencies( hid_t table, long table_len ) 02607 { 02608 02609 CHECK_OPEN_HANDLES; 02610 02611 ErrorCode rval; 02612 mhdf_Status status; 02613 02614 debug_barrier(); 02615 02616 hid_t read_type = H5Dget_type( table ); 02617 if (read_type < 0) 02618 return error(MB_FAILURE); 02619 const bool convert = !H5Tequal( read_type, handleType ); 02620 02621 EntityHandle* buffer = (EntityHandle*)dataBuffer; 02622 size_t chunk_size = bufferSize / H5Tget_size(read_type); 02623 size_t remaining = table_len; 02624 size_t left_over = 0; 02625 size_t offset = 0; 02626 dbgOut.printf( 3, "Reading adjacency list in %lu chunks\n", 02627 (unsigned long)(remaining + chunk_size - 1)/chunk_size ); 02628 int nn = 0; 02629 while (remaining) 02630 { 02631 dbgOut.printf( 3, "Reading chunk %d of adjacency list\n", ++nn ); 02632 02633 size_t count = std::min( chunk_size, remaining ); 02634 count -= left_over; 02635 remaining -= count; 02636 02637 assert_range( buffer + left_over, count ); 02638 mhdf_readAdjacencyWithOpt( table, offset, count, read_type, buffer + left_over, 02639 collIO, &status ); 02640 if (is_error(status)) 02641 return error(MB_FAILURE); 02642 02643 if (convert) { 02644 herr_t err = H5Tconvert( read_type, handleType, count, buffer + left_over, 0, H5P_DEFAULT ); 02645 if (err < 0) 02646 return error(MB_FAILURE); 02647 } 02648 02649 EntityHandle* iter = buffer; 02650 EntityHandle* end = buffer + count + left_over; 02651 while (end - iter >= 3) 02652 { 02653 EntityHandle h = idMap.find( *iter++ ); 02654 EntityHandle count2 = *iter++; 02655 if (!h) { 02656 iter += count2; 02657 continue; 02658 } 02659 02660 if (count2 < 1) 02661 return error(MB_FAILURE); 02662 02663 if (end < count2 + iter) 02664 { 02665 iter -= 2; 02666 break; 02667 } 02668 02669 size_t valid; 02670 convert_id_to_handle( iter, count2, valid, idMap ); 02671 rval = iFace->add_adjacencies( h, iter, valid, false ); 02672 if (MB_SUCCESS != rval) 02673 return error(rval); 02674 02675 iter += count2; 02676 } 02677 02678 left_over = end - iter; 02679 assert_range( (char*)buffer, left_over ); 02680 assert_range( (char*)iter, left_over ); 02681 memmove( buffer, iter, left_over ); 02682 } 02683 02684 assert(!left_over); // unexpected truncation of data 02685 02686 return MB_SUCCESS; 02687 } 02688 02689 02690 ErrorCode ReadHDF5::read_tag( int tag_index ) 02691 { 02692 02693 CHECK_OPEN_HANDLES; 02694 02695 dbgOut.tprintf(2, "Reading tag \"%s\"\n", fileInfo->tags[tag_index].name ); 02696 02697 debug_barrier(); 02698 02699 02700 ErrorCode rval; 02701 mhdf_Status status; 02702 Tag tag = 0; 02703 hid_t read_type = -1; 02704 bool table_type; 02705 rval = create_tag( fileInfo->tags[tag_index], tag, read_type ); 02706 if (MB_SUCCESS != rval) 02707 return error(rval); 02708 02709 if (fileInfo->tags[tag_index].have_sparse) { 02710 hid_t handles[3]; 02711 long num_ent, num_val; 02712 mhdf_openSparseTagData( filePtr, 02713 fileInfo->tags[tag_index].name, 02714 &num_ent, &num_val, 02715 handles, &status ); 02716 if (is_error(status)) { 02717 if (read_type) H5Tclose( read_type ); 02718 return error(MB_FAILURE); 02719 } 02720 02721 table_type = false; 02722 if (read_type == 0) { 02723 read_type = H5Dget_type( handles[1] ); 02724 if (read_type == 0) { 02725 mhdf_closeData( filePtr, handles[0], &status ); 02726 mhdf_closeData( filePtr, handles[0], &status ); 02727 if (fileInfo->tags[tag_index].size <= 0) 02728 mhdf_closeData( filePtr, handles[2], &status ); 02729 return error(MB_FAILURE); 02730 } 02731 table_type = true; 02732 } 02733 02734 if (fileInfo->tags[tag_index].size > 0) { 02735 dbgOut.printf(2, "Reading sparse data for tag \"%s\"\n", fileInfo->tags[tag_index].name ); 02736 rval = read_sparse_tag( tag, read_type, handles[0], handles[1], num_ent ); 02737 } 02738 else { 02739 dbgOut.printf(2, "Reading var-len sparse data for tag \"%s\"\n", fileInfo->tags[tag_index].name ); 02740 rval = read_var_len_tag( tag, read_type, handles[0], handles[1], handles[2], num_ent, num_val ); 02741 } 02742 02743 if (table_type) { 02744 H5Tclose(read_type); 02745 read_type = 0; 02746 } 02747 02748 mhdf_closeData( filePtr, handles[0], &status ); 02749 if (MB_SUCCESS == rval && is_error(status)) 02750 rval = MB_FAILURE; 02751 mhdf_closeData( filePtr, handles[1], &status ); 02752 if (MB_SUCCESS == rval && is_error(status)) 02753 rval = MB_FAILURE; 02754 if (fileInfo->tags[tag_index].size <= 0) { 02755 mhdf_closeData( filePtr, handles[2], &status ); 02756 if (MB_SUCCESS == rval && is_error(status)) 02757 rval = MB_FAILURE; 02758 } 02759 if (MB_SUCCESS != rval) { 02760 if (read_type) H5Tclose( read_type ); 02761 return error(rval); 02762 } 02763 } 02764 02765 for (int j = 0; j < fileInfo->tags[tag_index].num_dense_indices; ++j) { 02766 long count; 02767 const char* name = 0; 02768 mhdf_EntDesc* desc; 02769 int elem_idx = fileInfo->tags[tag_index].dense_elem_indices[j]; 02770 if (elem_idx == -2) { 02771 desc = &fileInfo->sets; 02772 name = mhdf_set_type_handle(); 02773 } 02774 else if (elem_idx == -1) { 02775 desc = &fileInfo->nodes; 02776 name = mhdf_node_type_handle(); 02777 } 02778 else if (elem_idx >= 0 && elem_idx < fileInfo->num_elem_desc) { 02779 desc = &fileInfo->elems[elem_idx].desc; 02780 name = fileInfo->elems[elem_idx].handle; 02781 } 02782 else { 02783 return error(MB_FAILURE); 02784 } 02785 02786 dbgOut.printf(2, "Read dense data block for tag \"%s\" on \"%s\"\n", fileInfo->tags[tag_index].name, name ); 02787 02788 hid_t handle = mhdf_openDenseTagData( filePtr, 02789 fileInfo->tags[tag_index].name, 02790 name, 02791 &count, &status ); 02792 if (is_error(status)) { 02793 rval = error(MB_FAILURE); 02794 break; 02795 } 02796 02797 if (count > desc->count) { 02798 readUtil->report_error( "Invalid data length for dense tag data: %s/%s\n", 02799 name, fileInfo->tags[tag_index].name ); 02800 mhdf_closeData( filePtr, handle, &status ); 02801 rval = error(MB_FAILURE); 02802 break; 02803 } 02804 02805 table_type = false; 02806 if (read_type == 0) { 02807 read_type = H5Dget_type( handle ); 02808 if (read_type == 0) { 02809 mhdf_closeData( filePtr, handle, &status ); 02810 return error(MB_FAILURE); 02811 } 02812 table_type = true; 02813 } 02814 02815 rval = read_dense_tag( tag, name, read_type, handle, desc->start_id, count ); 02816 02817 if (table_type) { 02818 H5Tclose( read_type ); 02819 read_type = 0; 02820 } 02821 02822 mhdf_closeData( filePtr, handle, &status ); 02823 if (MB_SUCCESS != rval) 02824 break; 02825 if (is_error(status)) { 02826 rval = error(MB_FAILURE); 02827 break; 02828 } 02829 } 02830 02831 if (read_type) 02832 H5Tclose( read_type ); 02833 return rval; 02834 } 02835 02836 02837 02838 02839 02840 ErrorCode ReadHDF5::create_tag( const mhdf_TagDesc& info, 02841 Tag& handle, 02842 hid_t& hdf_type ) 02843 { 02844 02845 CHECK_OPEN_HANDLES; 02846 02847 ErrorCode rval; 02848 mhdf_Status status; 02849 TagType storage; 02850 DataType mb_type; 02851 bool re_read_default = false; 02852 02853 switch (info.storage) { 02854 case mhdf_DENSE_TYPE : storage = MB_TAG_DENSE ; break; 02855 case mhdf_SPARSE_TYPE: storage = MB_TAG_SPARSE; break; 02856 case mhdf_BIT_TYPE : storage = MB_TAG_BIT; break; 02857 case mhdf_MESH_TYPE : storage = MB_TAG_MESH; break; 02858 default: 02859 readUtil->report_error( "Invalid storage type for tag '%s': %d\n", info.name, info.storage ); 02860 return error(MB_FAILURE); 02861 } 02862 02863 // Type-specific stuff 02864 if (info.type == mhdf_BITFIELD) { 02865 if (info.size < 1 || info.size > 8) 02866 { 02867 readUtil->report_error( "Invalid bit tag: class is MB_TAG_BIT, num bits = %d\n", info.size ); 02868 return error(MB_FAILURE); 02869 } 02870 hdf_type = H5Tcopy(H5T_NATIVE_B8); 02871 mb_type = MB_TYPE_BIT; 02872 if (hdf_type < 0) 02873 return error(MB_FAILURE); 02874 } 02875 else if (info.type == mhdf_OPAQUE) { 02876 mb_type = MB_TYPE_OPAQUE; 02877 02878 // Check for user-provided type 02879 Tag type_handle; 02880 std::string tag_type_name = "__hdf5_tag_type_"; 02881 tag_type_name += info.name; 02882 rval = iFace->tag_get_handle( tag_type_name.c_str(), sizeof(hid_t), MB_TYPE_OPAQUE, type_handle ); 02883 if (MB_SUCCESS == rval) { 02884 EntityHandle root = 0; 02885 rval = iFace->tag_get_data( type_handle, &root, 1, &hdf_type ); 02886 if (MB_SUCCESS != rval) 02887 return error(rval); 02888 hdf_type = H5Tcopy( hdf_type ); 02889 re_read_default = true; 02890 } 02891 else if (MB_TAG_NOT_FOUND == rval) { 02892 hdf_type = 0; 02893 } 02894 else 02895 return error(rval); 02896 02897 if (hdf_type < 0) 02898 return error(MB_FAILURE); 02899 } 02900 else { 02901 switch (info.type) 02902 { 02903 case mhdf_INTEGER: 02904 hdf_type = H5T_NATIVE_INT; 02905 mb_type = MB_TYPE_INTEGER; 02906 break; 02907 02908 case mhdf_FLOAT: 02909 hdf_type = H5T_NATIVE_DOUBLE; 02910 mb_type = MB_TYPE_DOUBLE; 02911 break; 02912 02913 case mhdf_BOOLEAN: 02914 hdf_type = H5T_NATIVE_UINT; 02915 mb_type = MB_TYPE_INTEGER; 02916 break; 02917 02918 case mhdf_ENTITY_ID: 02919 hdf_type = handleType; 02920 mb_type = MB_TYPE_HANDLE; 02921 break; 02922 02923 default: 02924 return error(MB_FAILURE); 02925 } 02926 02927 if (info.size > 1) { // array 02928 hsize_t tmpsize = info.size; 02929 #if defined(H5Tarray_create_vers) && H5Tarray_create_vers > 1 02930 hdf_type = H5Tarray_create2( hdf_type, 1, &tmpsize ); 02931 #else 02932 hdf_type = H5Tarray_create( hdf_type, 1, &tmpsize, NULL ); 02933 #endif 02934 } 02935 else { 02936 hdf_type = H5Tcopy( hdf_type ); 02937 } 02938 if (hdf_type < 0) 02939 return error(MB_FAILURE); 02940 } 02941 02942 02943 // If default or global/mesh value in file, read it. 02944 if (info.default_value || info.global_value) 02945 { 02946 if (re_read_default) { 02947 mhdf_getTagValues( filePtr, info.name, hdf_type, info.default_value, info.global_value, &status ); 02948 if (mhdf_isError( &status )) 02949 { 02950 readUtil->report_error( "%s", mhdf_message( &status ) ); 02951 if (hdf_type) H5Tclose( hdf_type ); 02952 return error(MB_FAILURE); 02953 } 02954 } 02955 02956 if (MB_TYPE_HANDLE == mb_type) { 02957 if (info.default_value) { 02958 rval = convert_id_to_handle( (EntityHandle*)info.default_value, info.default_value_size ); 02959 if (MB_SUCCESS != rval) { 02960 if (hdf_type) H5Tclose( hdf_type ); 02961 return error(rval); 02962 } 02963 } 02964 if (info.global_value) { 02965 rval = convert_id_to_handle( (EntityHandle*)info.global_value, info.global_value_size ); 02966 if (MB_SUCCESS != rval) { 02967 if (hdf_type) H5Tclose( hdf_type ); 02968 return error(rval); 02969 } 02970 } 02971 } 02972 } 02973 02974 // get tag handle, creating if necessary 02975 if (info.size < 0) 02976 rval = iFace->tag_get_handle( info.name, info.default_value_size, 02977 mb_type, handle, storage|MB_TAG_CREAT|MB_TAG_VARLEN|MB_TAG_DFTOK, 02978 info.default_value ); 02979 else 02980 rval = iFace->tag_get_handle( info.name, info.size, mb_type, handle, 02981 storage|MB_TAG_CREAT|MB_TAG_DFTOK, info.default_value ); 02982 if (MB_SUCCESS != rval) 02983 { 02984 readUtil->report_error( "Tag type in file does not match type in " 02985 "database for \"%s\"\n", info.name ); 02986 if (hdf_type) H5Tclose( hdf_type ); 02987 return error(MB_FAILURE); 02988 } 02989 02990 if (info.global_value) { 02991 EntityHandle root = 0; 02992 if (info.size > 0) { // fixed-length tag 02993 rval = iFace->tag_set_data( handle, &root, 1, info.global_value ); 02994 } 02995 else { 02996 int tag_size = info.global_value_size; 02997 rval = iFace->tag_set_by_ptr( handle, &root, 1, &info.global_value, &tag_size ); 02998 } 02999 if (MB_SUCCESS != rval) { 03000 if (hdf_type) H5Tclose( hdf_type ); 03001 return error(rval); 03002 } 03003 } 03004 03005 return MB_SUCCESS; 03006 } 03007 03008 ErrorCode ReadHDF5::read_dense_tag( Tag tag_handle, 03009 const char* ent_name, 03010 hid_t hdf_read_type, 03011 hid_t data, 03012 long start_id, 03013 long num_values ) 03014 { 03015 03016 CHECK_OPEN_HANDLES; 03017 03018 ErrorCode rval; 03019 DataType mb_type; 03020 03021 rval = iFace->tag_get_data_type( tag_handle, mb_type ); 03022 if (MB_SUCCESS != rval) 03023 return error(rval); 03024 03025 03026 int read_size; 03027 rval = iFace->tag_get_bytes( tag_handle, read_size ); 03028 if (MB_SUCCESS != rval) // wrong function for variable-length tags 03029 return error(rval); 03030 //if (MB_TYPE_BIT == mb_type) 03031 // read_size = (read_size + 7)/8; // convert bits to bytes, plus 7 for ceiling 03032 03033 if (hdf_read_type) { // if not opaque 03034 hsize_t hdf_size = H5Tget_size( hdf_read_type ); 03035 if (hdf_size != (hsize_t)read_size) 03036 return error(MB_FAILURE); 03037 } 03038 03039 // get actual entities read from file 03040 Range file_ids, handles; 03041 Range::iterator f_ins = file_ids.begin(), h_ins = handles.begin(); 03042 IDMap::iterator l, u; 03043 l = idMap.lower_bound( start_id ); 03044 u = idMap.lower_bound( start_id + num_values - 1 ); 03045 if (l != idMap.end() && start_id + num_values > l->begin) { 03046 if (l == u) { 03047 size_t beg = std::max(start_id, l->begin); 03048 size_t end = std::min(start_id + num_values, u->begin + u->count) - 1; 03049 f_ins = file_ids.insert( f_ins, beg, end ); 03050 h_ins = handles.insert( h_ins, l->value + (beg - l->begin), 03051 l->value + (end - l->begin) ); 03052 } 03053 else { 03054 size_t beg = std::max(start_id, l->begin); 03055 f_ins = file_ids.insert( f_ins, beg, l->begin + l->count - 1 ); 03056 h_ins = handles.insert( h_ins, l->value + (beg - l->begin), l->value + l->count - 1 ); 03057 for (++l; l != u; ++l) { 03058 f_ins = file_ids.insert( f_ins, l->begin, l->begin + l->count - 1 ); 03059 h_ins = handles.insert( h_ins, l->value, l->value + l->count - 1 ); 03060 } 03061 if (u != idMap.end() && u->begin < start_id + num_values) { 03062 size_t end = std::min( start_id + num_values, u->begin + u->count - 1 ); 03063 f_ins = file_ids.insert( f_ins, u->begin, end ); 03064 h_ins = handles.insert( h_ins, u->value, u->value + end - u->begin ); 03065 } 03066 } 03067 } 03068 03069 // Given that all of the entities for this dense tag data should 03070 // have been created as a single contiguous block, the resulting 03071 // MOAB handle range should be contiguous. 03072 // THE ABOVE IS NOT NECESSARILY TRUE. SOMETIMES LOWER-DIMENSION 03073 // ENTS ARE READ AND THEN DELETED FOR PARTIAL READS. 03074 //assert( handles.empty() || handles.size() == (handles.back() - handles.front() + 1)); 03075 03076 std::string tn("<error>"); 03077 iFace->tag_get_name( tag_handle, tn ); 03078 tn += " data for "; 03079 tn += ent_name; 03080 try { 03081 h_ins = handles.begin(); 03082 ReadHDF5Dataset reader( tn.c_str(), data, nativeParallel, mpiComm, false ); 03083 long buffer_size = bufferSize / read_size; 03084 reader.set_file_ids( file_ids, start_id, buffer_size, hdf_read_type ); 03085 dbgOut.printf( 3, "Reading dense data for tag \"%s\" and group \"%s\" in %lu chunks\n", 03086 tn.c_str(), ent_name, reader.get_read_count() ); 03087 int nn = 0; 03088 while (!reader.done()) { 03089 dbgOut.printf( 3, "Reading chunk %d of \"%s\" data\n", ++nn, tn.c_str() ); 03090 03091 size_t count; 03092 reader.read( dataBuffer, count ); 03093 03094 if (MB_TYPE_HANDLE == mb_type) { 03095 rval = convert_id_to_handle( (EntityHandle*)dataBuffer, count * read_size / sizeof(EntityHandle) ); 03096 if (MB_SUCCESS != rval) 03097 return error(rval); 03098 } 03099 03100 Range ents; 03101 Range::iterator end = h_ins; 03102 end += count; 03103 ents.insert( h_ins, end ); 03104 h_ins = end; 03105 03106 rval = iFace->tag_set_data( tag_handle, ents, dataBuffer ); 03107 if (MB_SUCCESS != rval) { 03108 dbgOut.printf(1,"Internal error setting data for tag \"%s\"\n", tn.c_str()); 03109 return error(rval); 03110 } 03111 } 03112 } 03113 catch (ReadHDF5Dataset::Exception) { 03114 dbgOut.printf(1,"Internal error reading dense data for tag \"%s\"\n",tn.c_str()); 03115 return error(MB_FAILURE); 03116 } 03117 03118 return MB_SUCCESS; 03119 } 03120 03121 03122 // Read entire ID table and for those file IDs corresponding 03123 // to entities that we have read from the file add both the 03124 // offset into the offset range and the handle into the handle 03125 // range. If handles are not ordered, switch to using a vector. 03126 ErrorCode ReadHDF5::read_sparse_tag_indices( const char* name, 03127 hid_t id_table, 03128 EntityHandle start_offset,// can't put zero in a Range 03129 Range& offset_range, 03130 Range& handle_range, 03131 std::vector<EntityHandle>& handle_vect ) 03132 { 03133 03134 CHECK_OPEN_HANDLES; 03135 03136 offset_range.clear(); 03137 handle_range.clear(); 03138 handle_vect.clear(); 03139 03140 ErrorCode rval; 03141 Range::iterator handle_hint = handle_range.begin(); 03142 Range::iterator offset_hint = offset_range.begin(); 03143 03144 EntityHandle* idbuf = (EntityHandle*)dataBuffer; 03145 size_t idbuf_size = bufferSize / sizeof(EntityHandle); 03146 03147 std::string tn(name); 03148 tn += " indices"; 03149 03150 assert(start_offset > 0); // can't put zero in a Range 03151 try { 03152 ReadHDF5Dataset id_reader( tn.c_str(), id_table, nativeParallel, mpiComm, false ); 03153 id_reader.set_all_file_ids( idbuf_size, handleType ); 03154 size_t offset = start_offset; 03155 dbgOut.printf( 3, "Reading file ids for sparse tag \"%s\" in %lu chunks\n", name, id_reader.get_read_count() ); 03156 int nn = 0; 03157 while (!id_reader.done()) {\ 03158 dbgOut.printf( 3, "Reading chunk %d of \"%s\" IDs\n", ++nn, name ); 03159 size_t count; 03160 id_reader.read( idbuf, count ); 03161 03162 rval = convert_id_to_handle( idbuf, count ); 03163 if (MB_SUCCESS != rval) 03164 return error(rval); 03165 03166 // idbuf will now contain zero-valued handles for those 03167 // tag values that correspond to entities we are not reading 03168 // from the file. 03169 for (size_t i = 0; i < count; ++i) { 03170 if (idbuf[i]) { 03171 offset_hint = offset_range.insert( offset_hint, offset+i ); 03172 if (!handle_vect.empty()) { 03173 handle_vect.push_back( idbuf[i] ); 03174 } 03175 else if (handle_range.empty() || idbuf[i] > handle_range.back()) { 03176 handle_hint = handle_range.insert( handle_hint, idbuf[i] ); 03177 } 03178 else { 03179 handle_vect.resize( handle_range.size() ); 03180 std::copy( handle_range.begin(), handle_range.end(), handle_vect.begin() ); 03181 handle_range.clear(); 03182 handle_vect.push_back( idbuf[i] ); 03183 dbgOut.print(2,"Switching to unordered list for tag handle list\n"); 03184 } 03185 } 03186 } 03187 03188 offset += count; 03189 } 03190 } 03191 catch (ReadHDF5Dataset::Exception) { 03192 return error(MB_FAILURE); 03193 } 03194 03195 return MB_SUCCESS; 03196 } 03197 03198 03199 03200 ErrorCode ReadHDF5::read_sparse_tag( Tag tag_handle, 03201 hid_t hdf_read_type, 03202 hid_t id_table, 03203 hid_t value_table, 03204 long /*num_values*/ ) 03205 { 03206 03207 CHECK_OPEN_HANDLES; 03208 03209 // Read entire ID table and for those file IDs corresponding 03210 // to entities that we have read from the file add both the 03211 // offset into the offset range and the handle into the handle 03212 // range. If handles are not ordered, switch to using a vector. 03213 const EntityHandle base_offset = 1; // can't put zero in a Range 03214 std::vector<EntityHandle> handle_vect; 03215 Range handle_range, offset_range; 03216 std::string tn("<error>"); 03217 iFace->tag_get_name( tag_handle, tn ); 03218 ErrorCode rval = read_sparse_tag_indices( tn.c_str(), 03219 id_table, base_offset, 03220 offset_range, handle_range, 03221 handle_vect ); 03222 03223 DataType mbtype; 03224 rval = iFace->tag_get_data_type( tag_handle, mbtype ); 03225 if (MB_SUCCESS != rval) 03226 return error(rval); 03227 03228 int read_size; 03229 rval = iFace->tag_get_bytes( tag_handle, read_size ); 03230 if (MB_SUCCESS != rval) // wrong function for variable-length tags 03231 return error(rval); 03232 //if (MB_TYPE_BIT == mbtype) 03233 // read_size = (read_size + 7)/8; // convert bits to bytes, plus 7 for ceiling 03234 03235 if (hdf_read_type) { // if not opaque 03236 hsize_t hdf_size = H5Tget_size( hdf_read_type ); 03237 if (hdf_size != (hsize_t)read_size) 03238 return error(MB_FAILURE); 03239 } 03240 03241 const int handles_per_tag = read_size/sizeof(EntityHandle); 03242 03243 // Now read data values 03244 size_t chunk_size = bufferSize / read_size; 03245 try { 03246 ReadHDF5Dataset val_reader( (tn + " values").c_str(), value_table, nativeParallel, mpiComm, false ); 03247 val_reader.set_file_ids( offset_range, base_offset, chunk_size, hdf_read_type ); 03248 dbgOut.printf( 3, "Reading sparse values for tag \"%s\" in %lu chunks\n", tn.c_str(), val_reader.get_read_count() ); 03249 int nn = 0; 03250 size_t offset = 0; 03251 while (!val_reader.done()) { 03252 dbgOut.printf( 3, "Reading chunk %d of \"%s\" values\n", ++nn, tn.c_str() ); 03253 size_t count; 03254 val_reader.read( dataBuffer, count ); 03255 if (MB_TYPE_HANDLE == mbtype) { 03256 rval = convert_id_to_handle( (EntityHandle*)dataBuffer, count*handles_per_tag ); 03257 if (MB_SUCCESS != rval) 03258 return error(rval); 03259 } 03260 03261 if (!handle_vect.empty()) { 03262 rval = iFace->tag_set_data( tag_handle, &handle_vect[offset], count, dataBuffer ); 03263 offset += count; 03264 } 03265 else { 03266 Range r; 03267 r.merge( handle_range.begin(), handle_range.begin() + count ); 03268 handle_range.erase( handle_range.begin(), handle_range.begin() + count ); 03269 rval = iFace->tag_set_data( tag_handle, r, dataBuffer ); 03270 } 03271 if (MB_SUCCESS != rval) 03272 return error(rval); 03273 } 03274 } 03275 catch (ReadHDF5Dataset::Exception) { 03276 return error(MB_FAILURE); 03277 } 03278 03279 return MB_SUCCESS; 03280 } 03281 03282 ErrorCode ReadHDF5::read_var_len_tag( Tag tag_handle, 03283 hid_t hdf_read_type, 03284 hid_t ent_table, 03285 hid_t val_table, 03286 hid_t off_table, 03287 long /*num_entities*/, 03288 long /*num_values*/ ) 03289 { 03290 03291 CHECK_OPEN_HANDLES; 03292 03293 ErrorCode rval; 03294 DataType mbtype; 03295 03296 rval = iFace->tag_get_data_type( tag_handle, mbtype ); 03297 if (MB_SUCCESS != rval) 03298 return error(rval); 03299 03300 // can't do variable-length bit tags 03301 if (MB_TYPE_BIT == mbtype) 03302 return error(MB_VARIABLE_DATA_LENGTH); 03303 03304 // if here, MOAB tag must be variable-length 03305 int mbsize; 03306 if (MB_VARIABLE_DATA_LENGTH != iFace->tag_get_bytes( tag_handle, mbsize )) { 03307 assert(false); 03308 return error(MB_VARIABLE_DATA_LENGTH); 03309 } 03310 03311 int read_size; 03312 if (hdf_read_type) { 03313 hsize_t hdf_size = H5Tget_size( hdf_read_type ); 03314 if (hdf_size < 1) 03315 return error(MB_FAILURE); 03316 read_size = hdf_size; 03317 } 03318 else { 03319 // opaque 03320 read_size = 1; 03321 } 03322 03323 std::string tn("<error>"); 03324 iFace->tag_get_name( tag_handle, tn ); 03325 03326 // Read entire ID table and for those file IDs corresponding 03327 // to entities that we have read from the file add both the 03328 // offset into the offset range and the handle into the handle 03329 // range. If handles are not ordered, switch to using a vector. 03330 const EntityHandle base_offset = 1; // can't put zero in a Range 03331 std::vector<EntityHandle> handle_vect; 03332 Range handle_range, offset_range; 03333 rval = read_sparse_tag_indices( tn.c_str(), 03334 ent_table, base_offset, 03335 offset_range, handle_range, 03336 handle_vect ); 03337 03338 // This code only works if the id_table is an ordered list. 03339 // This assumption was also true for the previous iteration 03340 // of this code, but wasn't checked. MOAB's file writer 03341 // always writes an ordered list for id_table. 03342 if (!handle_vect.empty()) { 03343 readUtil->report_error("Unordered file ids for variable length tag not supported.\n"); 03344 return MB_FAILURE; 03345 } 03346 03347 class VTReader : public ReadHDF5VarLen { 03348 Tag tagHandle; 03349 bool isHandle; 03350 size_t readSize; 03351 ReadHDF5* readHDF5; 03352 public: 03353 ErrorCode store_data( EntityHandle file_id, void* data, long count, bool ) 03354 { 03355 ErrorCode rval1; 03356 if (isHandle) { 03357 assert(readSize == sizeof(EntityHandle)); 03358 rval1 = readHDF5->convert_id_to_handle( (EntityHandle*)data, count ); 03359 if (MB_SUCCESS != rval1) 03360 return error(rval1); 03361 } 03362 int n = count; 03363 return readHDF5->moab()->tag_set_by_ptr( tagHandle, &file_id, 1, &data, &n ); 03364 } 03365 VTReader( DebugOutput& debug_output, void* buffer, size_t buffer_size, 03366 Tag tag, bool is_handle_tag, size_t read_size1, ReadHDF5* owner ) 03367 : ReadHDF5VarLen( debug_output, buffer, buffer_size ), 03368 tagHandle(tag), 03369 isHandle(is_handle_tag), 03370 readSize(read_size1), 03371 readHDF5(owner) 03372 {} 03373 }; 03374 03375 VTReader tool( dbgOut, dataBuffer, bufferSize, tag_handle, 03376 MB_TYPE_HANDLE == mbtype, read_size, this ); 03377 try { 03378 // Read offsets into value table. 03379 std::vector<unsigned> counts; 03380 Range offsets; 03381 ReadHDF5Dataset off_reader( (tn + " offsets").c_str(), off_table, nativeParallel, mpiComm, false ); 03382 rval = tool.read_offsets( off_reader, offset_range, base_offset, 03383 base_offset, offsets, counts ); 03384 if (MB_SUCCESS != rval) 03385 return error(rval); 03386 03387 // Read tag values 03388 Range empty; 03389 ReadHDF5Dataset val_reader( (tn + " values").c_str(), val_table, nativeParallel, mpiComm, false ); 03390 rval = tool.read_data( val_reader, offsets, base_offset, hdf_read_type, 03391 handle_range, counts, empty ); 03392 if (MB_SUCCESS != rval) 03393 return error(rval); 03394 } 03395 catch (ReadHDF5Dataset::Exception) { 03396 return error(MB_FAILURE); 03397 } 03398 03399 return MB_SUCCESS; 03400 } 03401 03402 ErrorCode ReadHDF5::convert_id_to_handle( EntityHandle* array, 03403 size_t size ) 03404 { 03405 convert_id_to_handle( array, size, idMap ); 03406 return MB_SUCCESS; 03407 } 03408 03409 void ReadHDF5::convert_id_to_handle( EntityHandle* array, 03410 size_t size, 03411 const RangeMap<long,EntityHandle>& id_map ) 03412 { 03413 for (EntityHandle* const end = array + size; array != end; ++array) 03414 *array = id_map.find( *array ); 03415 } 03416 03417 void ReadHDF5::convert_id_to_handle( EntityHandle* array, 03418 size_t size, size_t& new_size, 03419 const RangeMap<long,EntityHandle>& id_map ) 03420 { 03421 RangeMap<long,EntityHandle>::const_iterator it; 03422 new_size = 0; 03423 for (size_t i = 0; i < size; ++i) { 03424 it = id_map.lower_bound( array[i] ); 03425 if (it != id_map.end() && it->begin <= (long)array[i]) 03426 array[new_size++] = it->value + (array[i] - it->begin); 03427 } 03428 } 03429 03430 void ReadHDF5::convert_range_to_handle( const EntityHandle* ranges, 03431 size_t num_ranges, 03432 const RangeMap<long,EntityHandle>& id_map, 03433 Range& merge ) 03434 { 03435 RangeMap<long,EntityHandle>::iterator it = id_map.begin(); 03436 Range::iterator hint = merge.begin(); 03437 for (size_t i = 0; i < num_ranges; ++i) { 03438 long id = ranges[2*i]; 03439 const long end = id + ranges[2*i+1]; 03440 // we assume that 'ranges' is sorted, but check just in case it isn't. 03441 if (it == id_map.end() || it->begin > id) 03442 it = id_map.begin(); 03443 it = id_map.lower_bound( it, id_map.end(), id ); 03444 if (it == id_map.end()) 03445 continue; 03446 if (id < it->begin) 03447 id = it->begin; 03448 while (id < end) { 03449 if (id < it->begin) id = it->begin; 03450 const long off = id - it->begin; 03451 long count = std::min( it->count - off, end - id ); 03452 // it is possible that this new subrange is starting after the end 03453 // it will result in negative count, which does not make sense 03454 // we are done with this range, go to the next one 03455 if (count <= 0) 03456 break; 03457 hint = merge.insert( hint, it->value + off, it->value + off + count - 1 ); 03458 id += count; 03459 if (id < end) 03460 { 03461 if (++it == id_map.end()) 03462 break; 03463 if (it->begin > end) 03464 break; // 03465 } 03466 } 03467 } 03468 } 03469 03470 ErrorCode ReadHDF5::convert_range_to_handle( const EntityHandle* array, 03471 size_t num_ranges, 03472 Range& range ) 03473 { 03474 convert_range_to_handle( array, num_ranges, idMap, range ); 03475 return MB_SUCCESS; 03476 } 03477 03478 03479 ErrorCode ReadHDF5::insert_in_id_map( const Range& file_ids, 03480 EntityHandle start_id ) 03481 { 03482 IDMap tmp_map; 03483 bool merge = !idMap.empty() && !file_ids.empty() && idMap.back().begin > (long)file_ids.front(); 03484 IDMap& map = merge ? tmp_map : idMap; 03485 Range::const_pair_iterator p; 03486 for (p = file_ids.const_pair_begin(); p != file_ids.const_pair_end(); ++p) { 03487 size_t count = p->second - p->first + 1; 03488 if (!map.insert( p->first, start_id, count ).second) 03489 return error(MB_FAILURE); 03490 start_id += count; 03491 } 03492 if (merge && !idMap.merge( tmp_map )) 03493 return error(MB_FAILURE); 03494 03495 return MB_SUCCESS; 03496 } 03497 03498 ErrorCode ReadHDF5::insert_in_id_map( long file_id, 03499 EntityHandle handle ) 03500 { 03501 if (!idMap.insert( file_id, handle, 1 ).second) 03502 return error(MB_FAILURE); 03503 return MB_SUCCESS; 03504 } 03505 03506 03507 ErrorCode ReadHDF5::read_qa( EntityHandle ) 03508 { 03509 03510 CHECK_OPEN_HANDLES; 03511 03512 mhdf_Status status; 03513 std::vector<std::string> qa_list; 03514 03515 int qa_len; 03516 char** qa = mhdf_readHistory( filePtr, &qa_len, &status ); 03517 if (mhdf_isError( &status )) 03518 { 03519 readUtil->report_error( "%s", mhdf_message( &status ) ); 03520 return error(MB_FAILURE); 03521 } 03522 qa_list.resize(qa_len); 03523 for (int i = 0; i < qa_len; i++) 03524 { 03525 qa_list[i] = qa[i]; 03526 free( qa[i] ); 03527 } 03528 free( qa ); 03529 03532 return MB_SUCCESS; 03533 } 03534 03535 ErrorCode ReadHDF5::store_file_ids( Tag tag ) 03536 { 03537 03538 CHECK_OPEN_HANDLES; 03539 03540 typedef int tag_type; 03541 tag_type* buffer = reinterpret_cast<tag_type*>(dataBuffer); 03542 const long buffer_size = bufferSize / sizeof(tag_type); 03543 for (IDMap::iterator i = idMap.begin(); i != idMap.end(); ++i) { 03544 IDMap::Range range = *i; 03545 03546 // make sure the values will fit in the tag type 03547 IDMap::key_type rv = range.begin + (range.count - 1); 03548 tag_type tv = (tag_type)rv; 03549 if ((IDMap::key_type)tv != rv) { 03550 assert(false); 03551 return MB_INDEX_OUT_OF_RANGE; 03552 } 03553 03554 while (range.count) { 03555 long count = buffer_size < range.count ? buffer_size : range.count; 03556 03557 Range handles; 03558 handles.insert( range.value, range.value + count - 1 ); 03559 range.value += count; 03560 range.count -= count; 03561 for (long j = 0; j < count; ++j) 03562 buffer[j] = (tag_type)range.begin++; 03563 03564 ErrorCode rval = iFace->tag_set_data( tag, handles, buffer ); 03565 if (MB_SUCCESS != rval) 03566 return rval; 03567 } 03568 } 03569 return MB_SUCCESS; 03570 } 03571 03572 ErrorCode ReadHDF5::read_tag_values( const char* file_name, 03573 const char* tag_name, 03574 const FileOptions& opts, 03575 std::vector<int>& tag_values_out, 03576 const SubsetList* subset_list ) 03577 { 03578 ErrorCode rval; 03579 03580 rval = set_up_read( file_name, opts ); 03581 if (MB_SUCCESS != rval) 03582 return error(rval); 03583 03584 int tag_index; 03585 rval = find_int_tag( tag_name, tag_index ); 03586 if (MB_SUCCESS != rval) { 03587 clean_up_read( opts ); 03588 return error(rval); 03589 } 03590 03591 if (subset_list) { 03592 Range file_ids; 03593 rval = get_subset_ids( subset_list->tag_list, subset_list->tag_list_length, file_ids ); 03594 if (MB_SUCCESS != rval) { 03595 clean_up_read( opts ); 03596 return error(rval); 03597 } 03598 03599 rval = read_tag_values_partial( tag_index, file_ids, tag_values_out ); 03600 if (MB_SUCCESS != rval) { 03601 clean_up_read( opts ); 03602 return error(rval); 03603 } 03604 } 03605 else { 03606 rval = read_tag_values_all( tag_index, tag_values_out ); 03607 if (MB_SUCCESS != rval) { 03608 clean_up_read( opts ); 03609 return error(rval); 03610 } 03611 } 03612 03613 return clean_up_read( opts ); 03614 } 03615 03616 ErrorCode ReadHDF5::read_tag_values_partial( int tag_index, 03617 const Range& file_ids, 03618 std::vector<int>& tag_values ) 03619 { 03620 03621 CHECK_OPEN_HANDLES; 03622 03623 mhdf_Status status; 03624 const mhdf_TagDesc& tag = fileInfo->tags[tag_index]; 03625 long num_ent, num_val; 03626 size_t count; 03627 std::string tn(tag.name); 03628 03629 // read sparse values 03630 if (tag.have_sparse) { 03631 hid_t handles[3]; 03632 mhdf_openSparseTagData( filePtr, tag.name, &num_ent, &num_val, handles, &status ); 03633 if (mhdf_isError( &status )) { 03634 readUtil->report_error( "%s", mhdf_message( &status ) ); 03635 return error(MB_FAILURE); 03636 } 03637 03638 try { 03639 // read all entity handles and fill 'offsets' with ranges of 03640 // offsets into the data table for entities that we want. 03641 Range offsets; 03642 long* buffer = reinterpret_cast<long*>(dataBuffer); 03643 const long buffer_size = bufferSize/sizeof(long); 03644 ReadHDF5Dataset ids( (tn + " ids").c_str(), handles[0], nativeParallel, mpiComm ); 03645 ids.set_all_file_ids( buffer_size, H5T_NATIVE_LONG ); 03646 size_t offset = 0; 03647 dbgOut.printf( 3, "Reading sparse IDs for tag \"%s\" in %lu chunks\n", 03648 tag.name, ids.get_read_count() ); 03649 int nn = 0; 03650 while (!ids.done()) { 03651 dbgOut.printf( 3, "Reading chunk %d of IDs for \"%s\"\n", ++nn, tag.name ); 03652 ids.read( buffer, count ); 03653 03654 std::sort( buffer, buffer+count ); 03655 Range::iterator ins = offsets.begin(); 03656 Range::const_iterator i = file_ids.begin(); 03657 for (size_t j = 0; j < count; ++j) { 03658 while (i != file_ids.end() && (long)*i < buffer[j]) 03659 ++i; 03660 if (i == file_ids.end()) 03661 break; 03662 if ((long)*i == buffer[j]) { 03663 ins = offsets.insert( ins, j+offset, j+offset ); 03664 } 03665 } 03666 03667 offset += count; 03668 } 03669 03670 tag_values.clear(); 03671 tag_values.reserve( offsets.size() ); 03672 const size_t data_buffer_size = bufferSize/sizeof(int); 03673 int* data_buffer = reinterpret_cast<int*>(dataBuffer); 03674 ReadHDF5Dataset vals( (tn + " sparse vals").c_str(), handles[1], nativeParallel, mpiComm ); 03675 vals.set_file_ids( offsets, 0, data_buffer_size, H5T_NATIVE_INT ); 03676 dbgOut.printf( 3, "Reading sparse values for tag \"%s\" in %lu chunks\n", 03677 tag.name, vals.get_read_count() ); 03678 nn = 0; 03679 // should normally only have one read call, unless sparse nature 03680 // of file_ids caused reader to do something strange 03681 while (!vals.done()) { 03682 dbgOut.printf( 3, "Reading chunk %d of values for \"%s\"\n", ++nn, tag.name ); 03683 vals.read( data_buffer, count ); 03684 tag_values.insert( tag_values.end(), data_buffer, data_buffer+count ); 03685 } 03686 } 03687 catch (ReadHDF5Dataset::Exception) { 03688 return error(MB_FAILURE); 03689 } 03690 } 03691 03692 std::sort( tag_values.begin(), tag_values.end() ); 03693 tag_values.erase( std::unique(tag_values.begin(), tag_values.end()), tag_values.end() ); 03694 03695 // read dense values 03696 std::vector<int> prev_data, curr_data; 03697 for (int i = 0; i < tag.num_dense_indices; ++i) { 03698 int grp = tag.dense_elem_indices[i]; 03699 const char* gname = 0; 03700 mhdf_EntDesc* desc = 0; 03701 if (grp == -1) { 03702 gname = mhdf_node_type_handle(); 03703 desc = &fileInfo->nodes; 03704 } 03705 else if (grp == -2) { 03706 gname = mhdf_set_type_handle(); 03707 desc = &fileInfo->sets; 03708 } 03709 else { 03710 assert(grp >= 0 && grp < fileInfo->num_elem_desc); 03711 gname = fileInfo->elems[grp].handle; 03712 desc = &fileInfo->elems[grp].desc; 03713 } 03714 03715 Range::iterator s = file_ids.lower_bound( (EntityHandle)(desc->start_id) ); 03716 Range::iterator e = Range::lower_bound( s, file_ids.end(), 03717 (EntityHandle)(desc->start_id) + desc->count ); 03718 Range subset; 03719 subset.merge( s, e ); 03720 03721 hid_t handle = mhdf_openDenseTagData( filePtr, tag.name, gname, &num_val, &status ); 03722 if (mhdf_isError( &status )) { 03723 readUtil->report_error( "%s", mhdf_message( &status ) ); 03724 return error(MB_FAILURE); 03725 } 03726 03727 try { 03728 curr_data.clear(); 03729 tag_values.reserve( subset.size() ); 03730 const size_t data_buffer_size = bufferSize/sizeof(int); 03731 int* data_buffer = reinterpret_cast<int*>(dataBuffer); 03732 03733 ReadHDF5Dataset reader( (tn + " dense vals").c_str(), handle, nativeParallel, mpiComm ); 03734 reader.set_file_ids( subset, desc->start_id, data_buffer_size, H5T_NATIVE_INT ); 03735 dbgOut.printf( 3, "Reading dense data for tag \"%s\" and group \"%s\" in %lu chunks\n", 03736 tag.name, fileInfo->elems[grp].handle, reader.get_read_count() ); 03737 int nn = 0; 03738 // should normally only have one read call, unless sparse nature 03739 // of file_ids caused reader to do something strange 03740 while (!reader.done()) { 03741 dbgOut.printf( 3, "Reading chunk %d of \"%s\"/\"%s\"\n", ++nn, tag.name, fileInfo->elems[grp].handle ); 03742 reader.read( data_buffer, count ); 03743 curr_data.insert( curr_data.end(), data_buffer, data_buffer + count ); 03744 } 03745 } 03746 catch (ReadHDF5Dataset::Exception) { 03747 return error(MB_FAILURE); 03748 } 03749 03750 std::sort( curr_data.begin(), curr_data.end() ); 03751 curr_data.erase( std::unique(curr_data.begin(), curr_data.end()), curr_data.end() ); 03752 prev_data.clear(); 03753 tag_values.swap( prev_data ); 03754 std::set_union( prev_data.begin(), prev_data.end(), 03755 curr_data.begin(), curr_data.end(), 03756 std::back_inserter( tag_values ) ); 03757 } 03758 03759 return MB_SUCCESS; 03760 } 03761 03762 ErrorCode ReadHDF5::read_tag_values_all( int tag_index, 03763 std::vector<int>& tag_values ) 03764 { 03765 03766 CHECK_OPEN_HANDLES; 03767 03768 mhdf_Status status; 03769 const mhdf_TagDesc& tag = fileInfo->tags[tag_index]; 03770 long junk, num_val; 03771 03772 // read sparse values 03773 if (tag.have_sparse) { 03774 hid_t handles[3]; 03775 mhdf_openSparseTagData( filePtr, tag.name, &junk, &num_val, handles, &status ); 03776 if (mhdf_isError( &status )) { 03777 readUtil->report_error( "%s", mhdf_message( &status ) ); 03778 return error(MB_FAILURE); 03779 } 03780 03781 mhdf_closeData( filePtr, handles[0], &status ); 03782 if (mhdf_isError( &status )) { 03783 readUtil->report_error( "%s", mhdf_message( &status ) ); 03784 mhdf_closeData( filePtr, handles[1], &status ); 03785 return error(MB_FAILURE); 03786 } 03787 03788 hid_t file_type = H5Dget_type( handles[1] ); 03789 tag_values.resize( num_val ); 03790 mhdf_readTagValuesWithOpt( handles[1], 0, num_val, file_type, 03791 &tag_values[0], collIO, &status ); 03792 if (mhdf_isError( &status )) { 03793 readUtil->report_error( "%s", mhdf_message( &status ) ); 03794 H5Tclose( file_type ); 03795 mhdf_closeData( filePtr, handles[1], &status ); 03796 return error(MB_FAILURE); 03797 } 03798 H5Tconvert( file_type, H5T_NATIVE_INT, num_val, &tag_values[0], 0, H5P_DEFAULT ); 03799 H5Tclose( file_type ); 03800 03801 mhdf_closeData( filePtr, handles[1], &status ); 03802 if (mhdf_isError( &status )) { 03803 readUtil->report_error( "%s", mhdf_message( &status ) ); 03804 return error(MB_FAILURE); 03805 } 03806 } 03807 03808 std::sort( tag_values.begin(), tag_values.end() ); 03809 tag_values.erase( std::unique(tag_values.begin(), tag_values.end()), tag_values.end() ); 03810 03811 // read dense values 03812 std::vector<int> prev_data, curr_data; 03813 for (int i = 0; i < tag.num_dense_indices; ++i) { 03814 int grp = tag.dense_elem_indices[i]; 03815 const char* gname = 0; 03816 if (grp == -1) 03817 gname = mhdf_node_type_handle(); 03818 else if (grp == -2) 03819 gname = mhdf_set_type_handle(); 03820 else 03821 gname = fileInfo->elems[grp].handle; 03822 hid_t handle = mhdf_openDenseTagData( filePtr, tag.name, gname, &num_val, &status ); 03823 if (mhdf_isError( &status )) { 03824 readUtil->report_error( "%s", mhdf_message( &status ) ); 03825 return error(MB_FAILURE); 03826 } 03827 03828 hid_t file_type = H5Dget_type( handle ); 03829 curr_data.resize( num_val ); 03830 mhdf_readTagValuesWithOpt( handle, 0, num_val, file_type, &curr_data[0], collIO, &status ); 03831 if (mhdf_isError( &status )) { 03832 readUtil->report_error( "%s", mhdf_message( &status ) ); 03833 H5Tclose( file_type ); 03834 mhdf_closeData( filePtr, handle, &status ); 03835 return error(MB_FAILURE); 03836 } 03837 03838 H5Tconvert( file_type, H5T_NATIVE_INT, num_val, &curr_data[0], 0, H5P_DEFAULT ); 03839 H5Tclose( file_type ); 03840 mhdf_closeData( filePtr, handle, &status ); 03841 if (mhdf_isError( &status )) { 03842 readUtil->report_error( "%s", mhdf_message( &status ) ); 03843 return error(MB_FAILURE); 03844 } 03845 03846 std::sort( curr_data.begin(), curr_data.end() ); 03847 curr_data.erase( std::unique(curr_data.begin(), curr_data.end()), curr_data.end() ); 03848 03849 prev_data.clear(); 03850 tag_values.swap( prev_data ); 03851 std::set_union( prev_data.begin(), prev_data.end(), 03852 curr_data.begin(), curr_data.end(), 03853 std::back_inserter( tag_values ) ); 03854 } 03855 03856 return MB_SUCCESS; 03857 } 03858 03859 } // namespace moab