moab
ReadHDF5.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines