Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
WriteHDF5Parallel.cpp
Go to the documentation of this file.
1 #undef DEBUG
2 #undef TIME_DEBUG
3 
4 #include <cstdarg>
5 #include <ctime>
6 #include <cstdlib>
7 
8 #include <cstring>
9 #include <cassert>
10 
11 #include <vector>
12 #include <set>
13 #include <map>
14 #include <utility>
15 #include <iostream>
16 #include <sstream>
17 #include <string>
18 
19 #include "moab/Interface.hpp"
20 #include "Internals.hpp"
21 #include "MBTagConventions.hpp"
22 #include "MBParallelConventions.h"
23 #include "moab/ParallelComm.hpp"
24 #include "moab/CN.hpp"
25 #include "moab/Range.hpp"
26 #include "moab/CpuTimer.hpp"
27 
28 #include "WriteHDF5Parallel.hpp"
29 
30 #ifndef MOAB_HAVE_HDF5
31 #error Attempt to compile WriteHDF5Parallel with HDF5 support disabled
32 #endif
33 
34 #include <H5Tpublic.h>
35 #include <H5Ppublic.h>
36 #include <H5FDmpi.h>
37 #include <H5FDmpio.h>
38 
39 #include "mhdf.h"
40 
41 #include "IODebugTrack.hpp"
42 #include "moab/FileOptions.hpp"
43 
44 namespace
45 {
46 template < bool Condition >
47 struct STATIC_ASSERTION;
48 template <>
49 struct STATIC_ASSERTION< true >
50 {
51 };
52 } // namespace
53 
54 #define PP_CAT_( a, b ) a##b
55 #define PP_CAT( a, b ) PP_CAT_( a, b )
56 #define STATIC_ASSERT( Condition ) \
57  enum \
58  { \
59  PP_CAT( dummy, __LINE__ ) = sizeof( ::STATIC_ASSERTION< (bool)( Condition ) > ) \
60  }
61 
62 namespace moab
63 {
64 
65 #ifndef _WIN32 // problematic for windows
66 // Need an MPI type that we can put handles in
67 STATIC_ASSERT( sizeof( unsigned long ) >= sizeof( EntityHandle ) );
68 
69 // Need an MPI type that we can put file IDs in
70 STATIC_ASSERT( sizeof( unsigned long ) >= sizeof( WriteHDF5::wid_t ) );
71 #endif
72 
73 // This function doesn't do anything useful. It's just a nice
74 // place to set a break point to determine why the reader fails.
75 static inline ErrorCode error( ErrorCode rval )
76 {
77  return rval;
78 }
79 
80 const char* mpi_err_str( int errorcode )
81 {
82  static char buffer[2048];
83  int len = sizeof( buffer );
84  MPI_Error_string( errorcode, buffer, &len );
85  buffer[std::min( (size_t)len, sizeof( buffer ) - 1 )] = '\0';
86  return buffer;
87 }
88 
89 #define MPI_FAILURE_MSG( A ) \
90  "MPI Failure at " __FILE__ ":%d : (Code %d) %s\n", __LINE__, (int)( A ), mpi_err_str( ( A ) )
91 
92 #define CHECK_MPI( A ) \
93  do \
94  { \
95  if( MPI_SUCCESS != ( A ) ) \
96  { \
97  MB_SET_ERR_CONT( "MPI Failure : (Code " << (int)( A ) << ") " << mpi_err_str( ( A ) ) ); \
98  dbgOut.printf( 1, MPI_FAILURE_MSG( ( A ) ) ); \
99  return error( MB_FAILURE ); \
100  } \
101  } while( false )
102 
103 #define MB_FAILURE_MSG( A ) "MOAB_Failure at " __FILE__ ":%d : %s (%d)\n", __LINE__, ErrorCodeStr[( A )], (int)( A )
104 
105 #define CHECK_MB( A ) \
106  do \
107  { \
108  if( MB_SUCCESS != ( A ) ) \
109  { \
110  MB_SET_ERR_CONT( "MOAB Failure : " << ErrorCodeStr[( A )] ); \
111  dbgOut.printf( 1, MB_FAILURE_MSG( ( A ) ) ); \
112  return error( A ); \
113  } \
114  } while( false )
115 
116 #define HDF_FAILURE_MSG( A ) "MHDF Failure at " __FILE__ ":%d : %s\n", __LINE__, mhdf_message( &( A ) )
117 
118 #define CHECK_HDF( A ) \
119  do \
120  { \
121  if( mhdf_isError( &( A ) ) ) \
122  { \
123  MB_SET_ERR_CONT( "MHDF Failure : " << mhdf_message( &( A ) ) ); \
124  dbgOut.printf( 1, HDF_FAILURE_MSG( ( A ) ) ); \
125  return error( MB_FAILURE ); \
126  } \
127  } while( false )
128 
129 #define CHECK_HDFN( A ) \
130  do \
131  { \
132  if( mhdf_isError( &( A ) ) ) \
133  { \
134  MB_SET_ERR_CONT( "MHDF Failure : " << mhdf_message( &( A ) ) ); \
135  return error( MB_FAILURE ); \
136  } \
137  } while( false )
138 
139 #ifdef VALGRIND
140 #include <valgrind/memcheck.h>
141 
142 template < typename T >
143 inline void VALGRIND_MAKE_VEC_UNDEFINED( std::vector< T >& v )
144 {
145  if( v.size() )
146  {
147  }
148  (void)VALGRIND_MAKE_MEM_UNDEFINED( v.data(), v.size() * sizeof( T ) );
149 }
150 
151 #else
152 #ifndef VALGRIND_CHECK_MEM_IS_DEFINED
153 #define VALGRIND_CHECK_MEM_IS_DEFINED( a, b ) ( (void)0 )
154 #endif
155 #ifndef VALGRIND_CHECK_MEM_IS_ADDRESSABLE
156 #define VALGRIND_CHECK_MEM_IS_ADDRESSABLE( a, b ) ( (void)0 )
157 #endif
158 #ifndef VALGRIND_MAKE_MEM_UNDEFINED
159 #define VALGRIND_MAKE_MEM_UNDEFINED( a, b ) ( (void)0 )
160 #endif
161 
162 template < typename T >
163 inline void VALGRIND_MAKE_VEC_UNDEFINED( std::vector< T >& )
164 {
165  /* Nothing to do */
166 }
167 
168 #endif
169 
170 #ifndef NDEBUG
171 #define START_SERIAL \
172  for( unsigned _x = 0; _x < myPcomm->proc_config().proc_size(); ++_x ) \
173  { \
174  MPI_Barrier( myPcomm->proc_config().proc_comm() ); \
175  if( _x != myPcomm->proc_config().proc_rank() ) continue
176 #define END_SERIAL \
177  } \
178  MPI_Barrier( myPcomm->proc_config().proc_comm() )
179 #else
180 #define START_SERIAL
181 #define END_SERIAL
182 #endif
183 
184 static int my_Gatherv( void* sendbuf,
185  int sendcount,
186  MPI_Datatype sendtype,
187  std::vector< unsigned char >& recvbuf,
188  std::vector< int >& recvcounts,
189  int root,
190  MPI_Comm comm )
191 {
192  int nproc, rank, bytes, err;
193  MPI_Comm_size( comm, &nproc );
194  MPI_Comm_rank( comm, &rank );
195  MPI_Type_size( sendtype, &bytes );
196 
197  recvcounts.resize( rank == root ? nproc : 0 );
198  err = MPI_Gather( &sendcount, 1, MPI_INT, recvcounts.data(), 1, MPI_INT, root, comm );
199  if( MPI_SUCCESS != err ) return err;
200 
201  std::vector< int > disp( recvcounts.size() );
202  if( root == rank )
203  {
204  disp[0] = 0;
205  for( int i = 1; i < nproc; ++i )
206  disp[i] = disp[i - 1] + recvcounts[i - 1];
207  recvbuf.resize( bytes * ( disp.back() + recvcounts.back() ) );
208  }
209 
210  return MPI_Gatherv( sendbuf, sendcount, sendtype, recvbuf.data(), recvcounts.data(), disp.data(), sendtype, root,
211  comm );
212 }
213 
214 static void print_type_sets( Interface* iFace, DebugOutput* str, Range& sets )
215 {
216  const unsigned VB = 2;
217  if( str->get_verbosity() < VB ) return;
218 
219  Tag gid, did, bid, sid, nid;
220  gid = iFace->globalId_tag();
221  iFace->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, did );
222  iFace->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, bid );
223  iFace->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, nid );
224  iFace->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, sid );
225  Range typesets[10];
226  const char* typenames[] = { "Block ", "Sideset ", "NodeSet", "Vertex", "Curve",
227  "Surface", "Volume", "Body", "Other" };
228  for( Range::iterator riter = sets.begin(); riter != sets.end(); ++riter )
229  {
230  unsigned dim, id; //, oldsize;
231  if( MB_SUCCESS == iFace->tag_get_data( bid, &*riter, 1, &id ) )
232  dim = 0;
233  else if( MB_SUCCESS == iFace->tag_get_data( sid, &*riter, 1, &id ) )
234  dim = 1;
235  else if( MB_SUCCESS == iFace->tag_get_data( nid, &*riter, 1, &id ) )
236  dim = 2;
237  else if( MB_SUCCESS == iFace->tag_get_data( did, &*riter, 1, &dim ) )
238  {
239  id = 0;
240  iFace->tag_get_data( gid, &*riter, 1, &id );
241  dim += 3;
242  }
243  else
244  {
245  id = *riter;
246  dim = 9;
247  }
248 
249  // oldsize = typesets[dim].size();
250  typesets[dim].insert( id );
251  // assert(typesets[dim].size() - oldsize == 1);
252  }
253  for( int ii = 0; ii < 9; ++ii )
254  {
255  char tmp[64];
256  snprintf( tmp, 64, "%s (%lu) ", typenames[ii], (unsigned long)typesets[ii].size() );
257  str->print( VB, tmp, typesets[ii] );
258  }
259  str->printf( VB, "Total: %lu\n", (unsigned long)sets.size() );
260 }
261 
262 #define debug_barrier() debug_barrier_line( __LINE__ )
263 
264 void WriteHDF5Parallel::debug_barrier_line( int lineno )
265 {
266  const unsigned threshold = 2;
267  static unsigned long count = 0;
268  if( dbgOut.get_verbosity() >= threshold && myPcomm )
269  {
270  dbgOut.printf( threshold, "*********** Debug Barrier %lu (@%d)***********\n", ++count, lineno );
271  MPI_Barrier( myPcomm->proc_config().proc_comm() );
272  }
273 }
274 
275 WriterIface* WriteHDF5Parallel::factory( Interface* iface )
276 {
277  return new WriteHDF5Parallel( iface );
278 }
279 
280 WriteHDF5Parallel::WriteHDF5Parallel( Interface* iface )
281  : WriteHDF5( iface ), myPcomm( NULL ), pcommAllocated( false ), hslabOp( H5S_SELECT_OR )
282 {
283 }
284 
285 WriteHDF5Parallel::~WriteHDF5Parallel()
286 {
287  if( pcommAllocated && myPcomm ) delete myPcomm;
288 }
289 
290 // The parent WriteHDF5 class has ExportSet structs that are
291 // populated with the entities to be written, grouped by type
292 // (and for elements, connectivity length). This function:
293 // o determines which entities are to be written by a remote processor
294 // o removes those entities from the ExportSet structs in WriteMesh
295 // o passes them back in a Range
296 ErrorCode WriteHDF5Parallel::gather_interface_meshes( Range& nonowned )
297 {
298  ErrorCode result;
299 
300  // START_SERIAL;
301  dbgOut.print( 3, "Pre-interface mesh:\n" );
302  dbgOut.print( 3, nodeSet.range );
303  for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
304  dbgOut.print( 3, eiter->range );
305  dbgOut.print( 3, setSet.range );
306 
307  // Move handles of non-owned entities from lists of entities
308  // that this processor will write to the 'nonowned' list.
309 
310  nonowned.clear();
311  result = myPcomm->filter_pstatus( nodeSet.range, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonowned );
312  if( MB_SUCCESS != result ) return error( result );
313  nodeSet.range = subtract( nodeSet.range, nonowned );
314 
315  for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
316  {
317  Range tmpset;
318  result = myPcomm->filter_pstatus( eiter->range, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &tmpset );
319  if( MB_SUCCESS != result ) return error( result );
320  eiter->range = subtract( eiter->range, tmpset );
321  nonowned.merge( tmpset );
322  }
323 
324  dbgOut.print( 3, "Post-interface mesh:\n" );
325  dbgOut.print( 3, nodeSet.range );
326  for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
327  dbgOut.print( 3, eiter->range );
328  dbgOut.print( 3, setSet.range );
329 
330  // END_SERIAL;
331 
332  return MB_SUCCESS;
333 }
334 
335 ErrorCode WriteHDF5Parallel::parallel_create_file( const char* filename,
336  bool overwrite,
337  const std::vector< std::string >& qa_records,
338  const FileOptions& opts,
339  const Tag* user_tag_list,
340  int user_tag_count,
341  int dimension,
342  double* times )
343 {
344  ErrorCode rval;
345  mhdf_Status status;
346 
347  int pcomm_no = 0;
348  opts.get_int_option( "PARALLEL_COMM", pcomm_no );
349 
350  myPcomm = ParallelComm::get_pcomm( iFace, pcomm_no );
351  if( 0 == myPcomm )
352  {
353  myPcomm = new ParallelComm( iFace, MPI_COMM_WORLD );
354  pcommAllocated = true;
355  }
356 
357  MPI_Info info = MPI_INFO_NULL;
358  std::string cb_size;
359  rval = opts.get_str_option( "CB_BUFFER_SIZE", cb_size );
360  if( MB_SUCCESS == rval )
361  {
362  MPI_Info_create( &info );
363  MPI_Info_set( info, const_cast< char* >( "cb_buffer_size" ), const_cast< char* >( cb_size.c_str() ) );
364  }
365 
366  dbgOut.set_rank( myPcomm->proc_config().proc_rank() );
367  dbgOut.limit_output_to_first_N_procs( 32 );
368  Range nonlocal;
369  debug_barrier();
370  dbgOut.tprint( 1, "Gathering interface meshes\n" );
371  rval = gather_interface_meshes( nonlocal );
372  if( MB_SUCCESS != rval ) return error( rval );
373 
374  /**************** Get tag names for sets likely to be shared ***********/
375  // debug_barrier();
376  // dbgOut.tprint(1, "Getting shared entity sets\n");
377  // rval = get_sharedset_tags();
378  // if (MB_SUCCESS != rval) return error(rval);
379 
380  /**************** Create actual file and write meta info ***************/
381 
382  debug_barrier();
383  if( myPcomm->proc_config().proc_rank() == 0 )
384  {
385  dbgOut.tprintf( 1, "Creating file: %s\n", filename );
386 
387  // Create the file
388  const char* type_names[MBMAXTYPE];
389  memset( type_names, 0, MBMAXTYPE * sizeof( char* ) );
390  for( EntityType i = MBEDGE; i < MBENTITYSET; ++i )
391  type_names[i] = CN::EntityTypeName( i );
392 
393  dbgOut.tprint( 1, "call mhdf_createFile\n" );
394  filePtr = mhdf_createFile( filename, overwrite, type_names, MBMAXTYPE, id_type, &status );
395  if( !filePtr )
396  {
397  MB_SET_ERR( MB_FAILURE, mhdf_message( &status ) );
398  }
399 
400  dbgOut.tprint( 1, "call write_qa\n" );
401  rval = write_qa( qa_records );
402  if( MB_SUCCESS != rval ) return error( rval );
403  }
404 
405  /**************** Create node coordinate table ***************/
406  CpuTimer timer;
407  debug_barrier();
408  dbgOut.tprint( 1, "creating node table\n" );
409  topState.start( "creating node table" );
410  rval = create_node_table( dimension );
411  topState.end( rval );
412  if( MB_SUCCESS != rval ) return error( rval );
413  if( times ) times[CREATE_NODE_TIME] = timer.time_elapsed();
414 
415  /**************** Create element tables ***************/
416 
417  debug_barrier();
418  dbgOut.tprint( 1, "negotiating element types\n" );
419  topState.start( "negotiating element types" );
420  rval = negotiate_type_list();
421  topState.end( rval );
422  if( MB_SUCCESS != rval ) return error( rval );
423  if( times ) times[NEGOTIATE_TYPES_TIME] = timer.time_elapsed();
424  dbgOut.tprint( 1, "creating element table\n" );
425  topState.start( "creating element tables" );
426  rval = create_element_tables();
427  topState.end( rval );
428  if( MB_SUCCESS != rval ) return error( rval );
429  if( times ) times[CREATE_ELEM_TIME] = timer.time_elapsed();
430 
431  /*************** Exchange file IDs *****************/
432 
433  debug_barrier();
434  dbgOut.tprint( 1, "communicating file ids\n" );
435  topState.start( "communicating file ids" );
436  rval = exchange_file_ids( nonlocal );
437  topState.end( rval );
438  if( MB_SUCCESS != rval ) return error( rval );
439  if( times ) times[FILEID_EXCHANGE_TIME] = timer.time_elapsed();
440 
441  /**************** Create meshset tables *********************/
442 
443  debug_barrier();
444  dbgOut.tprint( 1, "creating meshset table\n" );
445  topState.start( "creating meshset tables" );
446  rval = create_meshset_tables( times );
447  topState.end( rval );
448  if( MB_SUCCESS != rval ) return error( rval );
449  if( times ) times[CREATE_SET_TIME] = timer.time_elapsed();
450 
451  /**************** Create adjacency tables *********************/
452 
453  debug_barrier();
454  dbgOut.tprint( 1, "creating adjacency table\n" );
455  topState.start( "creating adjacency tables" );
456  rval = create_adjacency_tables();
457  topState.end( rval );
458  if( MB_SUCCESS != rval ) return error( rval );
459  if( times ) times[CREATE_ADJ_TIME] = timer.time_elapsed();
460 
461  /**************** Create tag data *********************/
462 
463  debug_barrier();
464  dbgOut.tprint( 1, "creating tag tables\n" );
465  topState.start( "creating tag tables" );
466  rval = gather_tags( user_tag_list, user_tag_count );
467  if( MB_SUCCESS != rval ) return error( rval );
468  rval = create_tag_tables();
469  topState.end( rval );
470  if( MB_SUCCESS != rval ) return error( rval );
471  if( times ) times[CREATE_TAG_TIME] = timer.time_elapsed();
472 
473  /************** Close serial file and reopen parallel *****************/
474 
475  if( 0 == myPcomm->proc_config().proc_rank() ) mhdf_closeFile( filePtr, &status );
476 
477  MPI_Barrier( myPcomm->proc_config().proc_comm() );
478  dbgOut.tprint( 1, "(re)opening file in parallel mode\n" );
479  unsigned long junk;
480  hid_t hdf_opt = H5Pcreate( H5P_FILE_ACCESS );
481  H5Pset_fapl_mpio( hdf_opt, myPcomm->proc_config().proc_comm(), info );
482  filePtr = mhdf_openFileWithOpt( filename, 1, &junk, id_type, hdf_opt, &status );
483  H5Pclose( hdf_opt );
484  if( !filePtr )
485  {
486  MB_SET_ERR( MB_FAILURE, mhdf_message( &status ) );
487  }
488 
489  if( collectiveIO )
490  {
491  dbgOut.print( 1, "USING COLLECTIVE IO\n" );
492  writeProp = H5Pcreate( H5P_DATASET_XFER );
493  H5Pset_dxpl_mpio( writeProp, H5FD_MPIO_COLLECTIVE );
494  }
495 
496  /* Test if we can use H5S_APPEND when selecting hyperslabs */
497  if( MB_SUCCESS != opts.get_null_option( "HYPERSLAB_OR" ) &&
498  ( MB_SUCCESS == opts.get_null_option( "HYPERSLAB_APPEND" ) || HDF5_can_append_hyperslabs() ) )
499  {
500  dbgOut.print( 1, "HDF5 library supports H5Sselect_hyperlsab with H5S_SELECT_APPEND\n" );
501  hslabOp = H5S_SELECT_APPEND;
502  }
503 
504  dbgOut.tprint( 1, "Exiting parallel_create_file\n" );
505  return MB_SUCCESS;
506 }
507 
508 class TagNameCompare
509 {
510  Interface* iFace;
511  std::string name1, name2;
512 
513  public:
514  TagNameCompare( Interface* iface ) : iFace( iface ) {}
515  bool operator()( const WriteHDF5::TagDesc& t1, const WriteHDF5::TagDesc& t2 );
516 };
517 
518 bool TagNameCompare::operator()( const WriteHDF5::TagDesc& t1, const WriteHDF5::TagDesc& t2 )
519 {
520  iFace->tag_get_name( t1.tag_id, name1 );
521  iFace->tag_get_name( t2.tag_id, name2 );
522  return name1 < name2;
523 }
524 
525 struct serial_tag_data
526 {
527  TagType storage;
528  DataType type;
529  int size;
530  int name_len;
531  int def_val_len;
532  char name[sizeof( unsigned long )];
533 
534  static size_t pad( size_t len )
535  {
536  if( len % sizeof( unsigned long ) )
537  return len + sizeof( unsigned long ) - len % sizeof( unsigned long );
538  else
539  return len;
540  }
541 
542  static size_t def_val_bytes( int def_val_len, DataType type )
543  {
544  switch( type )
545  {
546  case MB_TYPE_BIT:
547  return def_val_len ? 1 : 0;
548  case MB_TYPE_OPAQUE:
549  return def_val_len;
550  case MB_TYPE_INTEGER:
551  return def_val_len * sizeof( int );
552  case MB_TYPE_DOUBLE:
553  return def_val_len * sizeof( double );
554  case MB_TYPE_HANDLE:
555  return def_val_len * sizeof( EntityHandle );
556  }
557  return 0;
558  }
559 
560  static size_t len( int name_len, int def_val_len, DataType type )
561  {
562  return sizeof( serial_tag_data ) + pad( name_len + def_val_bytes( def_val_len, type ) ) -
563  sizeof( unsigned long );
564  }
565  size_t len() const
566  {
567  return len( name_len, def_val_len, type );
568  }
569  void* default_value()
570  {
571  return def_val_len ? name + name_len : 0;
572  }
573  const void* default_value() const
574  {
575  return const_cast< serial_tag_data* >( this )->default_value();
576  }
577  void set_default_value( const void* val )
578  {
579  memcpy( default_value(), val, def_val_bytes( def_val_len, type ) );
580  }
581 };
582 
583 ErrorCode WriteHDF5Parallel::append_serial_tag_data( std::vector< unsigned char >& buffer,
584  const WriteHDF5::TagDesc& tag )
585 {
586  ErrorCode rval;
587 
588  std::string name;
589  rval = iFace->tag_get_name( tag.tag_id, name );
590  if( MB_SUCCESS != rval ) return error( rval );
591 
592  // Get name length, including space for null char
593  size_t name_len = name.size() + 1;
594  if( name_len == 1 ) return MB_SUCCESS; // Skip tags with no name
595 
596  DataType data_type;
597  rval = iFace->tag_get_data_type( tag.tag_id, data_type );
598  if( MB_SUCCESS != rval ) return error( rval );
599 
600  // Get default value
601  int def_val_len;
602  const void* def_val;
603  if( MB_SUCCESS != iFace->tag_get_default_value( tag.tag_id, def_val, def_val_len ) )
604  {
605  def_val_len = 0;
606  def_val = 0;
607  }
608 
609  // Allocate struct within buffer
610  size_t init_size = buffer.size();
611  buffer.resize( init_size + serial_tag_data::len( name_len, def_val_len, data_type ) );
612  serial_tag_data* ptr = reinterpret_cast< serial_tag_data* >( &buffer[init_size] );
613 
614  // Populate struct
615  rval = iFace->tag_get_type( tag.tag_id, ptr->storage );
616  if( MB_SUCCESS != rval ) return error( rval );
617  ptr->type = data_type;
618  rval = iFace->tag_get_length( tag.tag_id, ptr->size );
619  if( MB_VARIABLE_DATA_LENGTH == rval )
620  ptr->size = MB_VARIABLE_LENGTH;
621  else if( MB_SUCCESS != rval )
622  return error( rval );
623  ptr->name_len = name_len;
624  Range range;
625  memset( ptr->name, 0, ptr->name_len );
626  memcpy( ptr->name, name.data(), name.size() );
627  ptr->def_val_len = def_val_len;
628  ptr->set_default_value( def_val );
629 
630  return MB_SUCCESS;
631 }
632 
633 ErrorCode WriteHDF5Parallel::check_serial_tag_data( const std::vector< unsigned char >& buffer,
634  std::vector< TagDesc* >* missing,
635  std::vector< TagDesc* >* newlist )
636 {
637  ErrorCode rval;
638 
639  // Use 'write_sparse' field as a 'visited' mark
640  std::list< TagDesc >::iterator tag_iter;
641  if( missing )
642  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
643  tag_iter->write_sparse = true;
644 
645  // Use a set as a temporary for what will ultimately go in
646  // newlist because we need to pass back newlist in the order
647  // of the tagList member.
648  std::set< TagDesc* > newset;
649 
650  // Iterate over data from, updating the local list of tags.
651  // Be careful to keep tagList sorted such that in the end all
652  // procs have the same list in the same order.
653  std::vector< unsigned char >::const_iterator diter = buffer.begin();
654  tag_iter = tagList.begin();
655  while( diter < buffer.end() )
656  {
657  // Get struct from buffer
658  const serial_tag_data* ptr = reinterpret_cast< const serial_tag_data* >( &*diter );
659 
660  // Find local struct for tag
661  std::string name( ptr->name );
662  std::string n;
663  iFace->tag_get_name( tag_iter->tag_id, n ); // Second time we've called, so shouldn't fail
664  if( n > name )
665  {
666  tag_iter = tagList.begin(); // New proc, start search from beginning
667  }
668  iFace->tag_get_name( tag_iter->tag_id, n );
669  while( n < name )
670  {
671  ++tag_iter;
672  if( tag_iter == tagList.end() ) break;
673  iFace->tag_get_name( tag_iter->tag_id, n );
674  }
675  if( tag_iter == tagList.end() || n != name )
676  { // New tag
677  TagDesc newtag;
678 
679  if( ptr->size == MB_VARIABLE_LENGTH )
680  rval = iFace->tag_get_handle( name.c_str(), ptr->def_val_len, ptr->type, newtag.tag_id,
681  MB_TAG_VARLEN | MB_TAG_CREAT | ptr->storage, ptr->default_value() );
682  else
683  rval = iFace->tag_get_handle( name.c_str(), ptr->size, ptr->type, newtag.tag_id,
684  MB_TAG_CREAT | ptr->storage, ptr->default_value() );
685  if( MB_SUCCESS != rval ) return error( rval );
686 
687  newtag.sparse_offset = 0;
688  newtag.var_data_offset = 0;
689  newtag.write_sparse = false;
690  newtag.max_num_ents = 0;
691  newtag.max_num_vals = 0;
692 
693  tag_iter = tagList.insert( tag_iter, newtag );
694  if( newlist ) newset.insert( &*tag_iter );
695  }
696  else
697  { // Check that tag is as expected
698  DataType type;
699  iFace->tag_get_data_type( tag_iter->tag_id, type );
700  if( type != ptr->type )
701  {
702  MB_SET_ERR( MB_FAILURE, "Processes have inconsistent data type for tag \"" << name << "\"" );
703  }
704  int size;
705  iFace->tag_get_length( tag_iter->tag_id, size );
706  if( size != ptr->size )
707  {
708  MB_SET_ERR( MB_FAILURE, "Processes have inconsistent size for tag \"" << name << "\"" );
709  }
710  tag_iter->write_sparse = false;
711  }
712 
713  // Step to next variable-length struct.
714  diter += ptr->len();
715  }
716 
717  // Now pass back any local tags that weren't in the buffer
718  if( missing )
719  {
720  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
721  {
722  if( tag_iter->write_sparse )
723  {
724  tag_iter->write_sparse = false;
725  missing->push_back( &*tag_iter );
726  }
727  }
728  }
729 
730  // Be careful to populate newlist in the same, sorted, order as tagList
731  if( newlist )
732  {
733  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
734  if( newset.find( &*tag_iter ) != newset.end() ) newlist->push_back( &*tag_iter );
735  }
736 
737  return MB_SUCCESS;
738 }
739 
740 static void set_bit( int position, unsigned char* bytes )
741 {
742  int byte = position / 8;
743  int bit = position % 8;
744  bytes[byte] |= ( ( (unsigned char)1 ) << bit );
745 }
746 
747 static bool get_bit( int position, const unsigned char* bytes )
748 {
749  int byte = position / 8;
750  int bit = position % 8;
751  return 0 != ( bytes[byte] & ( ( (unsigned char)1 ) << bit ) );
752 }
753 
754 ErrorCode WriteHDF5Parallel::create_tag_tables()
755 {
756  std::list< TagDesc >::iterator tag_iter;
757  ErrorCode rval;
758  int err;
759  const int rank = myPcomm->proc_config().proc_rank();
760  const MPI_Comm comm = myPcomm->proc_config().proc_comm();
761 
762  subState.start( "negotiating tag list" );
763 
764  dbgOut.tprint( 1, "communicating tag metadata\n" );
765 
766  dbgOut.printf( 2, "Exchanging tag data for %d tags.\n", (int)tagList.size() );
767 
768  // Sort tagList contents in alphabetical order by tag name
769  tagList.sort( TagNameCompare( iFace ) );
770 
771  // Negotiate total list of tags to write
772 
773  // Build concatenated list of all tag data
774  std::vector< unsigned char > tag_buffer;
775  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
776  {
777  rval = append_serial_tag_data( tag_buffer, *tag_iter );
778  CHECK_MB( rval );
779  }
780 
781  // Broadcast list from root to all other procs
782  unsigned long size = tag_buffer.size();
783  err = MPI_Bcast( &size, 1, MPI_UNSIGNED_LONG, 0, comm );
784  CHECK_MPI( err );
785  tag_buffer.resize( size );
786  err = MPI_Bcast( tag_buffer.data(), size, MPI_UNSIGNED_CHAR, 0, comm );
787  CHECK_MPI( err );
788 
789  // Update local tag list
790  std::vector< TagDesc* > missing;
791  rval = check_serial_tag_data( tag_buffer, &missing, 0 );
792  CHECK_MB( rval );
793 
794  // Check if we're done (0->done, 1->more, 2+->error)
795  int code, lcode = ( MB_SUCCESS != rval ) ? rval + 2 : missing.empty() ? 0 : 1;
796  err = MPI_Allreduce( &lcode, &code, 1, MPI_INT, MPI_MAX, comm );
797  CHECK_MPI( err );
798  if( code > 1 )
799  {
800  MB_SET_ERR_CONT( "Inconsistent tag definitions between procs" );
801  return error( (ErrorCode)( code - 2 ) );
802  }
803 
804  // If not done...
805  if( code )
806  {
807  dbgOut.print( 1, "Not all procs had same tag definitions, negotiating...\n" );
808 
809  // Get tags defined on this proc but not on root proc
810  tag_buffer.clear();
811  for( size_t i = 0; i < missing.size(); ++i )
812  {
813  rval = append_serial_tag_data( tag_buffer, *missing[i] );
814  CHECK_MB( rval );
815  }
816 
817  // Gather extra tag definitions on root processor
818  std::vector< int > junk; // don't care how many from each proc
819  assert( rank || tag_buffer.empty() ); // must be empty on root
820  err = my_Gatherv( tag_buffer.data(), tag_buffer.size(), MPI_UNSIGNED_CHAR, tag_buffer, junk, 0, comm );
821  CHECK_MPI( err );
822 
823  // Process serialized tag descriptions on root, and
824  rval = MB_SUCCESS;
825  if( 0 == rank )
826  {
827  // Process serialized tag descriptions on root, and
828  std::vector< TagDesc* > newlist;
829  rval = check_serial_tag_data( tag_buffer, 0, &newlist );
830  tag_buffer.clear();
831  // re-serialize a unique list of new tag definitions
832  for( size_t i = 0; MB_SUCCESS == rval && i != newlist.size(); ++i )
833  {
834  rval = append_serial_tag_data( tag_buffer, *newlist[i] );
835  CHECK_MB( rval );
836  }
837  }
838 
839  // Broadcast any new tag definitions from root to other procs
840  long this_size = tag_buffer.size();
841  if( MB_SUCCESS != rval ) this_size = -rval;
842  err = MPI_Bcast( &this_size, 1, MPI_LONG, 0, comm );
843  CHECK_MPI( err );
844  if( this_size < 0 )
845  {
846  MB_SET_ERR_CONT( "Inconsistent tag definitions between procs" );
847  return error( (ErrorCode)-this_size );
848  }
849  tag_buffer.resize( this_size );
850  err = MPI_Bcast( tag_buffer.data(), this_size, MPI_UNSIGNED_CHAR, 0, comm );
851  CHECK_MPI( err );
852 
853  // Process new tag definitions
854  rval = check_serial_tag_data( tag_buffer, 0, 0 );
855  CHECK_MB( rval );
856  }
857 
858  subState.end();
859  subState.start( "negotiate which element/tag combinations are dense" );
860 
861  // Figure out for which tag/element combinations we can
862  // write dense tag data.
863 
864  // Construct a table of bits,
865  // where each row of the table corresponds to a tag
866  // and each column to an element group.
867 
868  // Two extra, because first is nodes and last is sets.
869  // (n+7)/8 is ceil(n/8)
870  const int bytes_per_tag = ( exportList.size() + 9 ) / 8;
871  std::vector< unsigned char > data( bytes_per_tag * tagList.size(), 0 );
872  std::vector< unsigned char > recv( data.size(), 0 );
873  unsigned char* iter = data.data();
874  if( writeTagDense && !data.empty() )
875  {
876  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, iter += bytes_per_tag )
877  {
878 
879  Range tagged;
880  rval = get_sparse_tagged_entities( *tag_iter, tagged );
881  CHECK_MB( rval );
882 
883  int s;
884  if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) ) continue;
885 
886  std::string n;
887  iFace->tag_get_name( tag_iter->tag_id,
888  n ); // Second time we've called, so shouldn't fail
889 
890  // Check if we want to write this tag in dense format even if not
891  // all of the entities have a tag value. The criterion of this
892  // is that the tag be dense, have a default value, and have at
893  // least 2/3 of the entities tagged.
894  bool prefer_dense = false;
895  TagType type;
896  rval = iFace->tag_get_type( tag_iter->tag_id, type );
897  CHECK_MB( rval );
898  if( MB_TAG_DENSE == type )
899  {
900  const void* defval = 0;
901  rval = iFace->tag_get_default_value( tag_iter->tag_id, defval, s );
902  if( MB_SUCCESS == rval ) prefer_dense = true;
903  }
904 
905  int i = 0;
906  if( check_dense_format_tag( nodeSet, tagged, prefer_dense ) )
907  {
908  set_bit( i, iter );
909  dbgOut.printf( 2, "Can write dense data for \"%s\"/Nodes\n", n.c_str() );
910  }
911  std::list< ExportSet >::const_iterator ex_iter = exportList.begin();
912  for( ++i; ex_iter != exportList.end(); ++i, ++ex_iter )
913  {
914  // when writing in parallel, on some partitions, some of these element ranges might
915  // be empty so do not turn this tag as sparse, just because of that, leave it dense,
916  // if we prefer dense
917  if( ( prefer_dense && ex_iter->range.empty() ) ||
918  check_dense_format_tag( *ex_iter, tagged, prefer_dense ) )
919  {
920  set_bit( i, iter );
921  dbgOut.printf( 2, "Can write dense data for \"%s\"/%s\n", n.c_str(), ex_iter->name() );
922  }
923  }
924  if( check_dense_format_tag( setSet, tagged, prefer_dense ) )
925  {
926  set_bit( i, iter );
927  dbgOut.printf( 2, "Can write dense data for \"%s\"/Sets\n", n.c_str() );
928  }
929  }
930 
931  // Do bit-wise AND of list over all processors (only write dense format
932  // if all proccesses want dense format for this group of entities).
933  err = MPI_Allreduce( data.data(), recv.data(), data.size(), MPI_UNSIGNED_CHAR, MPI_BAND,
934  myPcomm->proc_config().proc_comm() );
935  CHECK_MPI( err );
936  } // if (writeTagDense)
937 
938  // Store initial counts for sparse-formatted tag data.
939  // The total number of values to send and receive will be the number of
940  // tags plus the number of var-len tags because we need to negotiate
941  // offsets into two different tables for the var-len tags.
942  std::vector< long > counts;
943 
944  // Record dense tag/element combinations
945  iter = recv.data();
946  const unsigned char* iter2 = data.data();
947  for( tag_iter = tagList.begin(); tag_iter != tagList.end();
948  ++tag_iter, iter += bytes_per_tag, iter2 += bytes_per_tag )
949  {
950 
951  Range tagged;
952  rval = get_sparse_tagged_entities( *tag_iter, tagged );
953  CHECK_MB( rval );
954 
955  std::string n;
956  iFace->tag_get_name( tag_iter->tag_id, n ); // Second time we've called, so shouldn't fail
957 
958  int i = 0;
959  if( get_bit( i, iter ) )
960  {
961  assert( get_bit( i, iter2 ) );
962  tag_iter->dense_list.push_back( nodeSet );
963  tagged -= nodeSet.range;
964  dbgOut.printf( 2, "Will write dense data for \"%s\"/Nodes\n", n.c_str() );
965  }
966  std::list< ExportSet >::const_iterator ex_iter = exportList.begin();
967  for( ++i; ex_iter != exportList.end(); ++i, ++ex_iter )
968  {
969  if( get_bit( i, iter ) )
970  {
971  assert( get_bit( i, iter2 ) );
972  tag_iter->dense_list.push_back( *ex_iter );
973  dbgOut.printf( 2, "WIll write dense data for \"%s\"/%s\n", n.c_str(), ex_iter->name() );
974  tagged -= ex_iter->range;
975  }
976  }
977  if( get_bit( i, iter ) )
978  {
979  assert( get_bit( i, iter2 ) );
980  tag_iter->dense_list.push_back( setSet );
981  dbgOut.printf( 2, "Will write dense data for \"%s\"/Sets\n", n.c_str() );
982  tagged -= setSet.range;
983  }
984 
985  counts.push_back( tagged.size() );
986 
987  int s;
988  if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
989  {
990  unsigned long data_len;
991  rval = get_tag_data_length( *tag_iter, tagged, data_len );
992  CHECK_MB( rval );
993  counts.push_back( data_len );
994  }
995  }
996 
997  subState.end();
998  subState.start( "Negotiate offsets for sparse tag info" );
999 
1000  std::vector< long > offsets( counts.size() ), maxima( counts.size() ), totals( counts.size() );
1001  rval = create_dataset( counts.size(), counts.data(), offsets.data(), maxima.data(), totals.data() );
1002  CHECK_MB( rval );
1003 
1004  // Copy values into local structs and if root then create tables
1005  size_t idx = 0;
1006  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, ++idx )
1007  {
1008  assert( idx < counts.size() );
1009  tag_iter->sparse_offset = offsets[idx];
1010  tag_iter->max_num_ents = maxima[idx];
1011  tag_iter->write_sparse = ( 0 != totals[idx] );
1012  int s;
1013  if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
1014  {
1015  ++idx;
1016  assert( idx < counts.size() );
1017  tag_iter->var_data_offset = offsets[idx];
1018  tag_iter->max_num_vals = maxima[idx];
1019  }
1020  else
1021  {
1022  tag_iter->var_data_offset = 0;
1023  tag_iter->max_num_vals = 0;
1024  }
1025  }
1026 
1027  subState.end();
1028 
1029  // Create tag tables on root process
1030  if( 0 == myPcomm->proc_config().proc_rank() )
1031  {
1032  size_t iidx = 0;
1033  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, ++iidx )
1034  {
1035  assert( iidx < totals.size() );
1036  unsigned long num_ents = totals[iidx];
1037  unsigned long num_val = 0;
1038  int s;
1039  if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
1040  {
1041  ++iidx;
1042  assert( iidx < totals.size() );
1043  num_val = totals[iidx];
1044  }
1045  dbgOut.printf( 2, "Writing tag description for tag 0x%lx with %lu values\n",
1046  (unsigned long)tag_iter->tag_id, num_val ? num_val : num_ents );
1047 
1048  rval = create_tag( *tag_iter, num_ents, num_val );
1049  if( MB_SUCCESS != rval ) return error( rval );
1050  }
1051  }
1052 
1053  if( dbgOut.get_verbosity() > 1 )
1054  {
1055  dbgOut.printf( 2, "Tags: %12s %8s %8s %8s %8s %8s\n", "Name", "Count", "Offset", "Var Off", "Max Ent",
1056  "Handle" );
1057 
1058  for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
1059  {
1060  std::string name;
1061  iFace->tag_get_name( tag_iter->tag_id, name );
1062  size_t this_size;
1063  get_num_sparse_tagged_entities( *tag_iter, this_size );
1064  dbgOut.printf( 2, "%18s %8lu %8lu %8lu %8lu 0x%7lx\n", name.c_str(), (unsigned long)this_size,
1065  (unsigned long)tag_iter->sparse_offset, (unsigned long)tag_iter->var_data_offset,
1066  (unsigned long)tag_iter->max_num_ents, (unsigned long)tag_iter->tag_id );
1067  }
1068  }
1069 
1070  return MB_SUCCESS;
1071 }
1072 
1073 struct DatasetVals
1074 {
1075  long start_id;
1076  long max_count;
1077  long total;
1078 };
1079 STATIC_ASSERT( ( sizeof( DatasetVals ) == 3 * sizeof( long ) ) );
1080 
1081 ErrorCode WriteHDF5Parallel::create_dataset( int num_datasets,
1082  const long* num_owned,
1083  long* offsets_out,
1084  long* max_proc_entities,
1085  long* total_entities,
1086  const DataSetCreator& creator,
1087  ExportSet* groups[],
1088  wid_t* first_ids_out )
1089 {
1090  int result;
1091  ErrorCode rval;
1092  const unsigned rank = myPcomm->proc_config().proc_rank();
1093  const unsigned nproc = myPcomm->proc_config().proc_size();
1094  const MPI_Comm comm = myPcomm->proc_config().proc_comm();
1095 
1096  // Gather entity counts for each processor on root
1097  std::vector< long > counts( rank ? 0 : nproc * num_datasets );
1098  (void)VALGRIND_CHECK_MEM_IS_DEFINED( &num_owned, sizeof( long ) );
1099  result = MPI_Gather( const_cast< long* >( num_owned ), num_datasets, MPI_LONG, counts.data(), num_datasets,
1100  MPI_LONG, 0, comm );
1101  CHECK_MPI( result );
1102 
1103  // Create node data in file
1104  DatasetVals zero_val = { 0, 0, 0 };
1105  std::vector< DatasetVals > cumulative( num_datasets, zero_val );
1106  if( rank == 0 )
1107  {
1108  for( unsigned i = 0; i < nproc; i++ )
1109  {
1110  const long* proc_data = &counts[i * num_datasets];
1111  for( int index = 0; index < num_datasets; ++index )
1112  {
1113  cumulative[index].total += proc_data[index];
1114  if( proc_data[index] > cumulative[index].max_count ) cumulative[index].max_count = proc_data[index];
1115  }
1116  }
1117 
1118  for( int index = 0; index < num_datasets; ++index )
1119  {
1120  if( cumulative[index].total )
1121  {
1122  rval = creator( this, cumulative[index].total, groups ? groups[index] : 0, cumulative[index].start_id );
1123  CHECK_MB( rval );
1124  }
1125  else
1126  {
1127  cumulative[index].start_id = -1;
1128  }
1129  }
1130  }
1131 
1132  // Send id offset to every proc
1133  result = MPI_Bcast( (void*)cumulative.data(), 3 * num_datasets, MPI_LONG, 0, comm );
1134  CHECK_MPI( result );
1135  for( int index = 0; index < num_datasets; ++index )
1136  {
1137  if( first_ids_out ) first_ids_out[index] = (wid_t)cumulative[index].start_id;
1138  max_proc_entities[index] = cumulative[index].max_count;
1139  total_entities[index] = cumulative[index].total;
1140  }
1141 
1142  // Convert array of per-process counts to per-process offsets
1143  if( rank == 0 )
1144  {
1145  // Initialize prev_size with data sizes for root process
1146  std::vector< long > prev_size( counts.begin(), counts.begin() + num_datasets );
1147  // Root process gets offset zero
1148  std::fill( counts.begin(), counts.begin() + num_datasets, 0L );
1149  // For each proc other than this one (root)
1150  for( unsigned i = 1; i < nproc; ++i )
1151  {
1152  // Get pointer to offsets for previous process in list
1153  long* prev_data = &counts[( i - 1 ) * num_datasets];
1154  // Get pointer to offsets for this process in list
1155  long* proc_data = &counts[i * num_datasets];
1156  // For each data set
1157  for( int j = 0; j < num_datasets; ++j )
1158  {
1159  // Get size of data in dataset from process i
1160  long mysize = proc_data[j];
1161  // Offset for process i is offset of previous process plus
1162  // number of values previous process will write
1163  proc_data[j] = prev_data[j] + prev_size[j];
1164  // Store my size, as it is no longer available in 'counts'
1165  prev_size[j] = mysize;
1166  }
1167  }
1168  }
1169 
1170  // Send each proc it's offset in the table
1171  if( rank == 0 )
1172  {
1173  (void)VALGRIND_CHECK_MEM_IS_DEFINED( counts.data(), num_datasets * nproc * sizeof( long ) );
1174  }
1175  result = MPI_Scatter( counts.data(), num_datasets, MPI_LONG, offsets_out, num_datasets, MPI_LONG, 0, comm );
1176  CHECK_MPI( result );
1177 
1178  return MB_SUCCESS;
1179 }
1180 
1181 ErrorCode WriteHDF5Parallel::create_node_table( int dimension )
1182 {
1183  nodeSet.num_nodes = dimension; // Put it here so NodeSetCreator can access it
1184  struct NodeSetCreator : public DataSetCreator
1185  {
1186  ErrorCode operator()( WriteHDF5* file, long count, const ExportSet* group, long& start_id ) const
1187  {
1188  mhdf_Status status;
1189  hid_t handle = mhdf_createNodeCoords( file->file_ptr(), group->num_nodes, count, &start_id, &status );
1190  CHECK_HDFN( status );
1191  mhdf_closeData( file->file_ptr(), handle, &status );
1192  CHECK_HDFN( status );
1193  return MB_SUCCESS;
1194  }
1195  };
1196 
1197  const long count = nodeSet.range.size();
1198  ExportSet* array[] = { &nodeSet };
1199  ErrorCode rval = create_dataset( 1, &count, &nodeSet.offset, &nodeSet.max_num_ents, &nodeSet.total_num_ents,
1200  NodeSetCreator(), array, &nodeSet.first_id );
1201  CHECK_MB( rval );
1202  return assign_ids( nodeSet.range, nodeSet.first_id + nodeSet.offset );
1203 }
1204 
1205 struct elemtype
1206 {
1207  int mbtype;
1208  int numnode;
1209 
1210  elemtype( int vals[2] ) : mbtype( vals[0] ), numnode( vals[1] ) {}
1211  elemtype( int t, int n ) : mbtype( t ), numnode( n ) {}
1212 
1213  bool operator==( const elemtype& other ) const
1214  {
1215  return mbtype == other.mbtype && ( mbtype == MBENTITYSET || numnode == other.numnode );
1216  }
1217  bool operator<( const elemtype& other ) const
1218  {
1219  if( mbtype > other.mbtype ) return false;
1220 
1221  return mbtype < other.mbtype || ( mbtype != MBENTITYSET && numnode < other.numnode );
1222  }
1223  bool operator!=( const elemtype& other ) const
1224  {
1225  return !this->operator==( other );
1226  }
1227 };
1228 
1229 ErrorCode WriteHDF5Parallel::negotiate_type_list()
1230 {
1231  int result;
1232  const MPI_Comm comm = myPcomm->proc_config().proc_comm();
1233 
1234  exportList.sort();
1235  int num_types = exportList.size();
1236 
1237  // Get list of types on this processor
1238  typedef std::vector< std::pair< int, int > > typelist;
1239  typelist my_types( num_types );
1240  (void)VALGRIND_MAKE_VEC_UNDEFINED( my_types );
1241  typelist::iterator viter = my_types.begin();
1242  for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
1243  {
1244  viter->first = eiter->type;
1245  viter->second = eiter->num_nodes;
1246  ++viter;
1247  }
1248 
1249  dbgOut.print( 2, "Local Element Types:\n" );
1250  for( viter = my_types.begin(); viter != my_types.end(); ++viter )
1251  {
1252  int type = viter->first;
1253  int count = viter->second;
1254  dbgOut.printf( 2, " %s : %d\n", CN::EntityTypeName( (EntityType)type ), count );
1255  }
1256 
1257  // Broadcast number of types from root to all nodes
1258  int num_types0 = num_types;
1259  result = MPI_Bcast( &num_types0, 1, MPI_INT, 0, comm );
1260  CHECK_MPI( result );
1261  // Broadcast type list from root to all nodes
1262  typelist root_types( num_types0 );
1263  if( 0 == myPcomm->proc_config().proc_rank() ) root_types = my_types;
1264  result = MPI_Bcast( (void*)root_types.data(), 2 * num_types0, MPI_INT, 0, comm );
1265  CHECK_MPI( result );
1266 
1267  // Build local list of any types that root did not know about
1268  typelist non_root_types;
1269  viter = root_types.begin();
1270  for( typelist::iterator iter = my_types.begin(); iter != my_types.end(); ++iter )
1271  {
1272  if( viter == root_types.end() || *viter != *iter )
1273  non_root_types.push_back( *iter );
1274  else
1275  ++viter;
1276  }
1277 
1278  // Determine if any process had types not defined on the root
1279  int non_root_count = non_root_types.size();
1280  int not_done;
1281  result = MPI_Allreduce( &non_root_count, &not_done, 1, MPI_INT, MPI_LOR, comm );
1282  CHECK_MPI( result );
1283  if( not_done )
1284  {
1285  // Get number of types each processor has that root does not
1286  std::vector< int > counts( myPcomm->proc_config().proc_size() );
1287  int two_count = 2 * non_root_count;
1288  result = MPI_Gather( &two_count, 1, MPI_INT, counts.data(), 1, MPI_INT, 0, comm );
1289  CHECK_MPI( result );
1290 
1291  // Get list of types from each processor
1292  std::vector< int > displs( myPcomm->proc_config().proc_size() + 1 );
1293  (void)VALGRIND_MAKE_VEC_UNDEFINED( displs );
1294  displs[0] = 0;
1295  for( unsigned long i = 1; i <= myPcomm->proc_config().proc_size(); ++i )
1296  displs[i] = displs[i - 1] + counts[i - 1];
1297  int total = displs[myPcomm->proc_config().proc_size()];
1298  typelist alltypes( total / 2 );
1299  (void)VALGRIND_MAKE_VEC_UNDEFINED( alltypes );
1300  (void)VALGRIND_CHECK_MEM_IS_DEFINED( non_root_types.data(), non_root_types.size() * sizeof( int ) );
1301  result = MPI_Gatherv( (void*)non_root_types.data(), 2 * non_root_count, MPI_INT, (int*)alltypes.data(),
1302  counts.data(), displs.data(), MPI_INT, 0, comm );
1303  CHECK_MPI( result );
1304 
1305  // Merge type lists.
1306  // Prefer O(n) insertions with O(ln n) search time because
1307  // we expect data from a potentially large number of processes,
1308  // but with a small total number of element types.
1309  if( 0 == myPcomm->proc_config().proc_rank() )
1310  {
1311  for( viter = alltypes.begin(); viter != alltypes.end(); ++viter )
1312  {
1313  typelist::iterator titer = std::lower_bound( my_types.begin(), my_types.end(), *viter );
1314  if( titer == my_types.end() || *titer != *viter ) my_types.insert( titer, *viter );
1315  }
1316 
1317  dbgOut.print( 2, "Global Element Types:\n" );
1318  for( viter = my_types.begin(); viter != my_types.end(); ++viter )
1319  dbgOut.printf( 2, " %s : %d\n", CN::EntityTypeName( (EntityType)viter->first ), viter->second );
1320  }
1321 
1322  // Send total number of types to each processor
1323  total = my_types.size();
1324  result = MPI_Bcast( &total, 1, MPI_INT, 0, comm );
1325  CHECK_MPI( result );
1326 
1327  // Send list of types to each processor
1328  my_types.resize( total );
1329  result = MPI_Bcast( (void*)my_types.data(), 2 * total, MPI_INT, 0, comm );
1330  CHECK_MPI( result );
1331  }
1332  else
1333  {
1334  // Special case: if root had types but some subset of procs did not
1335  // have those types, but there are no types that the root doesn't
1336  // know about then we still need to update processes that are missing
1337  // types.
1338  my_types.swap( root_types );
1339  }
1340 
1341  // Insert missing types into exportList, with an empty
1342  // range of entities to export.
1343  std::list< ExportSet >::iterator ex_iter = exportList.begin();
1344  for( viter = my_types.begin(); viter != my_types.end(); ++viter )
1345  {
1346  while( ex_iter != exportList.end() && *ex_iter < *viter )
1347  ++ex_iter;
1348 
1349  if( ex_iter == exportList.end() || !( *ex_iter == *viter ) )
1350  {
1351  ExportSet insert;
1352  insert.type = (EntityType)viter->first;
1353  insert.num_nodes = viter->second;
1354  insert.first_id = 0;
1355  insert.offset = 0;
1356  insert.adj_offset = 0;
1357  ex_iter = exportList.insert( ex_iter, insert );
1358  }
1359  }
1360 
1361  return MB_SUCCESS;
1362 }
1363 
1364 ErrorCode WriteHDF5Parallel::create_element_tables()
1365 {
1366  struct ElemSetCreator : public DataSetCreator
1367  {
1368  ErrorCode operator()( WriteHDF5* file, long size, const ExportSet* ex, long& start_id ) const
1369  {
1370  return file->create_elem_table( *ex, size, start_id );
1371  }
1372  };
1373 
1374  const int numtypes = exportList.size();
1375  if( numtypes > 0 )
1376  {
1377  std::vector< ExportSet* > groups( numtypes );
1378  std::vector< long > counts( numtypes ), offsets( numtypes ), max_ents( numtypes ), total_ents( numtypes );
1379  std::vector< wid_t > start_ids( numtypes );
1380 
1381  size_t idx = 0;
1382  std::list< ExportSet >::iterator ex_iter;
1383  for( ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter, ++idx )
1384  {
1385  groups[idx] = &*ex_iter;
1386  counts[idx] = ex_iter->range.size();
1387  }
1388  ErrorCode rval = create_dataset( numtypes, counts.data(), offsets.data(), max_ents.data(), total_ents.data(),
1389  ElemSetCreator(), groups.data(), start_ids.data() );
1390  CHECK_MB( rval );
1391 
1392  for( idx = 0, ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter, ++idx )
1393  {
1394  ex_iter->first_id = start_ids[idx];
1395  ex_iter->offset = offsets[idx];
1396  ex_iter->max_num_ents = max_ents[idx];
1397  ex_iter->total_num_ents = total_ents[idx];
1398  rval = assign_ids( ex_iter->range, ex_iter->first_id + ex_iter->offset );
1399  CHECK_MB( rval );
1400  }
1401  }
1402 
1403  return MB_SUCCESS;
1404 }
1405 
1406 ErrorCode WriteHDF5Parallel::create_adjacency_tables()
1407 {
1408  struct AdjSetCreator : public DataSetCreator
1409  {
1410  ErrorCode operator()( WriteHDF5* file, long size, const ExportSet* ex, long& start_id ) const
1411  {
1412  mhdf_Status status;
1413  hid_t handle = mhdf_createAdjacency( file->file_ptr(), ex->name(), size, &status );
1414  CHECK_HDFN( status );
1415  mhdf_closeData( file->file_ptr(), handle, &status );
1416  CHECK_HDFN( status );
1417  start_id = -1;
1418  return MB_SUCCESS;
1419  }
1420  };
1421 
1422  std::vector< ExportSet* > groups;
1423 #ifdef WRITE_NODE_ADJACENCIES
1424  groups.push_back( &nodeSet );
1425 #endif
1426  for( std::list< ExportSet >::iterator ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter )
1427  groups.push_back( &*ex_iter );
1428 
1429  ErrorCode rval;
1430  const int numtypes = groups.size();
1431  if( numtypes > 0 )
1432  {
1433  std::vector< long > counts( numtypes );
1434  std::vector< long > offsets( numtypes );
1435  std::vector< long > max_ents( numtypes );
1436  std::vector< long > totals( numtypes );
1437  for( int i = 0; i < numtypes; ++i )
1438  {
1439  wid_t count;
1440  rval = count_adjacencies( groups[i]->range, count );
1441  CHECK_MB( rval );
1442  counts[i] = count;
1443  }
1444 
1445  rval = create_dataset( numtypes, counts.data(), offsets.data(), max_ents.data(), totals.data(), AdjSetCreator(),
1446  groups.data() );
1447  CHECK_MB( rval );
1448 
1449  // Cppcheck warning (false positive): variable groups is assigned a value that is never used
1450  for( int i = 0; i < numtypes; ++i )
1451  {
1452  groups[i]->max_num_adjs = max_ents[i];
1453  groups[i]->adj_offset = offsets[i];
1454  }
1455  }
1456  return MB_SUCCESS;
1457 }
1458 
1459 const unsigned SSVB = 3;
1460 
1461 void WriteHDF5Parallel::print_set_sharing_data( const Range& range, const char* label, Tag idt )
1462 {
1463  dbgOut.printf( SSVB, "set\tid\towner\t%-*s\tfid\tshared\n", (int)( sizeof( EntityHandle ) * 2 ), "handle" );
1464  for( Range::iterator it = range.begin(); it != range.end(); ++it )
1465  {
1466  int id;
1467  iFace->tag_get_data( idt, &*it, 1, &id );
1468  EntityHandle handle = 0;
1469  unsigned owner = 0;
1470  wid_t file_id = 0;
1471  myPcomm->get_entityset_owner( *it, owner, &handle );
1472  if( !idMap.find( *it, file_id ) ) file_id = 0;
1473  dbgOut.printf( SSVB, "%s\t%d\t%u\t%lx\t%lu\t", label, id, owner, (unsigned long)handle,
1474  (unsigned long)file_id );
1475  std::vector< unsigned > procs;
1476  myPcomm->get_entityset_procs( *it, procs );
1477  if( procs.empty() )
1478  dbgOut.print( SSVB, "<none>\n" );
1479  else
1480  {
1481  for( unsigned i = 0; i < procs.size() - 1; ++i )
1482  dbgOut.printf( SSVB, "%u,", procs[i] );
1483  dbgOut.printf( SSVB, "%u\n", procs.back() );
1484  }
1485  }
1486 }
1487 
1488 void WriteHDF5Parallel::print_shared_sets()
1489 {
1490  const char* tag_names[][2] = { { MATERIAL_SET_TAG_NAME, "block" },
1491  { DIRICHLET_SET_TAG_NAME, "nodeset" },
1492  { NEUMANN_SET_TAG_NAME, "sideset" },
1493  { 0, 0 } };
1494 
1495  for( int i = 0; tag_names[i][0]; ++i )
1496  {
1497  Tag tag;
1498  if( MB_SUCCESS != iFace->tag_get_handle( tag_names[i][0], 1, MB_TYPE_INTEGER, tag ) ) continue;
1499 
1500  Range tagged;
1501  iFace->get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, 0, 1, tagged );
1502  print_set_sharing_data( tagged, tag_names[i][1], tag );
1503  }
1504 
1505  Tag geom, id;
1506  if( MB_SUCCESS != iFace->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom ) ) return;
1507  id = iFace->globalId_tag();
1508 
1509  const char* geom_names[] = { "vertex", "curve", "surface", "volume" };
1510  for( int d = 0; d <= 3; ++d )
1511  {
1512  Range tagged;
1513  const void* vals[] = { &d };
1514  iFace->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom, vals, 1, tagged );
1515  print_set_sharing_data( tagged, geom_names[d], id );
1516  }
1517 }
1518 
1519 ErrorCode WriteHDF5Parallel::communicate_shared_set_ids( const Range& owned, const Range& remote )
1520 {
1521  ErrorCode rval;
1522  int mperr;
1523  const int TAG = 0xD0E;
1524  // const unsigned rank = myPcomm->proc_config().proc_rank();
1525  const MPI_Comm comm = myPcomm->proc_config().proc_comm();
1526 
1527  dbgOut.tprint( 1, "COMMUNICATING SHARED SET IDS\n" );
1528  dbgOut.print( 6, "Owned, shared sets: ", owned );
1529 
1530  // Post receive buffers for all procs for which we share sets
1531 
1532  std::vector< unsigned > procs;
1533  rval = myPcomm->get_entityset_owners( procs );
1534  CHECK_MB( rval );
1535  std::vector< unsigned >::iterator it = std::find( procs.begin(), procs.end(), myPcomm->proc_config().proc_rank() );
1536  if( it != procs.end() ) procs.erase( it );
1537 
1538  std::vector< MPI_Request > recv_req( procs.size(), MPI_REQUEST_NULL );
1539  std::vector< std::vector< unsigned long > > recv_buf( procs.size() );
1540 
1541  size_t recv_count = 0;
1542  for( size_t i = 0; i < procs.size(); ++i )
1543  {
1544  Range tmp;
1545  rval = myPcomm->get_owned_sets( procs[i], tmp );
1546  CHECK_MB( rval );
1547  size_t count =
1548  intersect( tmp, remote ).size(); // Necessary because we might not be writing all of the database
1549  if( count )
1550  {
1551  dbgOut.printf( 6, "Sets owned by proc %u (remote handles): ", procs[i] );
1552  if( dbgOut.get_verbosity() >= 6 )
1553  {
1554  Range remote_handles;
1555  tmp = intersect( tmp, remote );
1556  for( Range::iterator j = tmp.begin(); j != tmp.end(); ++j )
1557  {
1558  unsigned r;
1559  EntityHandle h;
1560  myPcomm->get_entityset_owner( *j, r, &h );
1561  assert( r == procs[i] );
1562  remote_handles.insert( h );
1563  }
1564  dbgOut.print( 6, remote_handles );
1565  }
1566  recv_count++;
1567  recv_buf[i].resize( 2 * count + 1 );
1568  dbgOut.printf( 5, "Posting receive buffer of size %lu for proc %u (%lu of %lu owned sets)\n",
1569  (unsigned long)recv_buf[i].size(), procs[i], count, tmp.size() );
1570  mperr = MPI_Irecv( recv_buf[i].data(), recv_buf[i].size(), MPI_UNSIGNED_LONG, procs[i], TAG, comm,
1571  &recv_req[i] );
1572  CHECK_MPI( mperr );
1573  }
1574  }
1575 
1576  // Send set ids to all procs with which we share them
1577 
1578  // First build per-process lists of sets for which we need to send data
1579  std::map< unsigned, Range > send_sets;
1580  std::vector< unsigned > set_procs;
1581  for( Range::reverse_iterator i = owned.rbegin(); i != owned.rend(); ++i )
1582  {
1583  set_procs.clear();
1584  rval = myPcomm->get_entityset_procs( *i, set_procs );
1585  CHECK_MB( rval );
1586  for( size_t j = 0; j < set_procs.size(); ++j )
1587  if( set_procs[j] != myPcomm->proc_config().proc_rank() ) send_sets[set_procs[j]].insert( *i );
1588  }
1589  assert( send_sets.find( myPcomm->proc_config().proc_rank() ) == send_sets.end() );
1590 
1591  // Now send the data
1592  std::vector< std::vector< unsigned long > > send_buf( send_sets.size() );
1593  std::vector< MPI_Request > send_req( send_sets.size() );
1594  std::map< unsigned, Range >::iterator si = send_sets.begin();
1595  for( size_t i = 0; si != send_sets.end(); ++si, ++i )
1596  {
1597  dbgOut.printf( 6, "Sending data for shared sets to proc %u: ", si->first );
1598  dbgOut.print( 6, si->second );
1599 
1600  send_buf[i].reserve( 2 * si->second.size() + 1 );
1601  send_buf[i].push_back( si->second.size() );
1602  for( Range::iterator j = si->second.begin(); j != si->second.end(); ++j )
1603  {
1604  send_buf[i].push_back( *j );
1605  send_buf[i].push_back( idMap.find( *j ) );
1606  }
1607  dbgOut.printf( 5, "Sending buffer of size %lu to proc %u (%lu of %lu owned sets)\n",
1608  (unsigned long)send_buf[i].size(), si->first, si->second.size(), owned.size() );
1609  mperr =
1610  MPI_Isend( send_buf[i].data(), send_buf[i].size(), MPI_UNSIGNED_LONG, si->first, TAG, comm, &send_req[i] );
1611  }
1612 
1613  // Process received data
1614  MPI_Status status;
1615  int idx;
1616  while( recv_count-- )
1617  {
1618  mperr = MPI_Waitany( recv_req.size(), recv_req.data(), &idx, &status );
1619  CHECK_MPI( mperr );
1620 
1621  assert( (unsigned)status.MPI_SOURCE == procs[idx] );
1622  assert( 2 * recv_buf[idx].front() + 1 == recv_buf[idx].size() );
1623  const size_t n = std::min< size_t >( recv_buf[idx].front(), ( recv_buf[idx].size() - 1 ) / 2 );
1624  dbgOut.printf( 5, "Received buffer of size %lu from proc %d\n", (unsigned long)( 2 * n + 1 ),
1625  (int)status.MPI_SOURCE );
1626 
1627  for( size_t i = 0; i < n; ++i )
1628  {
1629  EntityHandle handle = 0;
1630  rval = myPcomm->get_entityset_local_handle( procs[idx], recv_buf[idx][2 * i + 1], handle );
1631  CHECK_MB( rval );
1632  assert( handle != 0 );
1633  if( !idMap.insert( handle, recv_buf[idx][2 * i + 2], 1 ).second )
1634  error( MB_FAILURE ); // Conflicting IDs??????
1635  }
1636 
1637  recv_req[idx] = MPI_REQUEST_NULL;
1638  }
1639  assert( MPI_SUCCESS == MPI_Waitany( recv_req.size(), recv_req.data(), &idx, &status ) &&
1640  MPI_UNDEFINED == idx ); // Check that we got them all
1641 
1642  // Wait for all sends to complete before we release send
1643  // buffers (implicitly releases when we return from this function)
1644 
1645  std::vector< MPI_Status > stats( send_req.size() );
1646  mperr = MPI_Waitall( send_req.size(), send_req.data(), stats.data() );
1647  CHECK_MPI( mperr );
1648 
1649  if( dbgOut.get_verbosity() >= SSVB ) print_shared_sets();
1650 
1651  return MB_SUCCESS;
1652 }
1653 
1654 // void get_global_ids(Interface* iFace, const unsigned long* ptr,
1655 // size_t len, unsigned flags,
1656 // std::vector<int>& ids)
1657 //{
1658 // Tag idtag;
1659 // iFace->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, idtag);
1660 // for (size_t i = 0; i < len; ++i) {
1661 // if (flags & MESHSET_ORDERED) {
1662 // int tmp;
1663 // iFace->tag_get_data(idtag, ptr + i, 1, &tmp);
1664 // ids.push_back(tmp);
1665 // continue;
1666 // }
1667 //
1668 // EntityHandle s = ptr[i];
1669 // EntityHandle e = ptr[++i];
1670 // for (; s <= e; ++s) {
1671 // int tmp;
1672 // iFace->tag_get_data(idtag, &s, 1, &tmp);
1673 // ids.push_back(tmp);
1674 // }
1675 // }
1676 //}
1677 
1678 ErrorCode WriteHDF5Parallel::pack_set( Range::const_iterator it, unsigned long* buffer, size_t buffer_size )
1679 {
1680  ErrorCode rval;
1681  const EntityHandle* ptr;
1682  int len;
1683  unsigned char flags;
1684  std::vector< wid_t > tmp;
1685  size_t newlen;
1686 
1687  // Buffer must always contain at least flags and desired sizes
1688  assert( buffer_size >= 4 );
1689  buffer_size -= 4;
1690 
1691  Range::const_iterator nd = it;
1692  ++nd;
1693  rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::CONTENTS, &len, &flags );
1694  CHECK_MB( rval );
1695 
1696  // Tag mattag;
1697  // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
1698  // int block;
1699  // if (MB_SUCCESS != iFace->tag_get_data(mattag, &*it, 1, &block))
1700  // block = 0;
1701  //
1702  // if (block) {
1703  // std::vector<int> ids;
1704  // get_global_ids(iFace, ptr, len, flags, ids);
1705  //}
1706 
1707  if( len && !( flags & MESHSET_ORDERED ) )
1708  {
1709  tmp.clear();
1710  bool blocked = false;
1711  assert( ( 0 == len % 2 ) );
1712  rval = range_to_blocked_list( ptr, len / 2, tmp, blocked );
1713  CHECK_MB( rval );
1714  if( blocked ) flags |= mhdf_SET_RANGE_BIT;
1715  }
1716  else
1717  {
1718  tmp.resize( len );
1719  rval = vector_to_id_list( ptr, len, tmp.data(), newlen, true );
1720  CHECK_MB( rval );
1721  tmp.resize( newlen );
1722  }
1723 
1724  buffer[0] = flags;
1725  buffer[1] = tmp.size();
1726  if( tmp.size() <= buffer_size ) std::copy( tmp.begin(), tmp.end(), buffer + 4 );
1727 
1728  rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::CHILDREN, &len );
1729  CHECK_MB( rval );
1730  tmp.resize( len );
1731  rval = vector_to_id_list( ptr, len, tmp.data(), newlen, true );
1732  tmp.resize( newlen );
1733  buffer[2] = tmp.size();
1734  if( tmp.size() <= buffer_size - buffer[1] ) std::copy( tmp.begin(), tmp.end(), buffer + 4 + buffer[1] );
1735 
1736  rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::PARENTS, &len );
1737  CHECK_MB( rval );
1738  tmp.resize( len );
1739  rval = vector_to_id_list( ptr, len, tmp.data(), newlen, true );
1740  tmp.resize( newlen );
1741  buffer[3] = tmp.size();
1742  if( tmp.size() <= buffer_size - buffer[1] - buffer[2] )
1743  std::copy( tmp.begin(), tmp.end(), buffer + 4 + buffer[1] + buffer[2] );
1744 
1745  return MB_SUCCESS;
1746 }
1747 
1748 template < typename TYPE >
1749 static void convert_to_ranged_ids( const TYPE* buffer, size_t len, std::vector< WriteHDF5::wid_t >& result )
1750 {
1751  if( !len )
1752  {
1753  result.clear();
1754  return;
1755  }
1756 
1757  result.resize( len * 2 );
1758  Range tmp;
1759  for( size_t i = 0; i < len; i++ )
1760  tmp.insert( (EntityHandle)buffer[i] );
1761  result.resize( tmp.psize() * 2 );
1762  int j = 0;
1763  for( Range::const_pair_iterator pit = tmp.const_pair_begin(); pit != tmp.const_pair_end(); ++pit, j++ )
1764  {
1765  result[2 * j] = pit->first;
1766  result[2 * j + 1] = pit->second - pit->first + 1;
1767  }
1768 }
1769 
1770 static void merge_ranged_ids( const unsigned long* range_list, size_t len, std::vector< WriteHDF5::wid_t >& result )
1771 {
1772  typedef WriteHDF5::wid_t wid_t;
1773  assert( 0 == len % 2 );
1774  assert( 0 == result.size() % 2 );
1775  STATIC_ASSERT( sizeof( std::pair< wid_t, wid_t > ) == 2 * sizeof( wid_t ) );
1776 
1777  result.insert( result.end(), range_list, range_list + len );
1778  size_t plen = result.size() / 2;
1779  Range tmp;
1780  for( size_t i = 0; i < plen; i++ )
1781  {
1782  EntityHandle starth = (EntityHandle)result[2 * i];
1783  EntityHandle endh = (EntityHandle)result[2 * i] + (wid_t)result[2 * i + 1] - 1; // id + count - 1
1784  tmp.insert( starth, endh );
1785  }
1786  // Now convert back to std::vector<WriteHDF5::wid_t>, compressed range format
1787  result.resize( tmp.psize() * 2 );
1788  int j = 0;
1789  for( Range::const_pair_iterator pit = tmp.const_pair_begin(); pit != tmp.const_pair_end(); ++pit, j++ )
1790  {
1791  result[2 * j] = pit->first;
1792  result[2 * j + 1] = pit->second - pit->first + 1;
1793  }
1794 }
1795 
1796 static void merge_vector_ids( const unsigned long* list, size_t len, std::vector< WriteHDF5::wid_t >& result )
1797 {
1798  result.insert( result.end(), list, list + len );
1799 }
1800 
1801 ErrorCode WriteHDF5Parallel::unpack_set( EntityHandle set, const unsigned long* buffer, size_t buffer_size )
1802 {
1803  // Use local variables for readability
1804  if( buffer_size < 4 ) return MB_FAILURE;
1805  assert( buffer[1] + buffer[2] + buffer[3] <= buffer_size );
1806 #ifdef NDEBUG
1807  UNUSED( buffer_size );
1808 #endif
1809  const unsigned long flags = buffer[0];
1810  unsigned long num_content = buffer[1];
1811  const unsigned long num_child = buffer[2];
1812  const unsigned long num_parent = buffer[3];
1813  const unsigned long* contents = buffer + 4;
1814  const unsigned long* children = contents + num_content;
1815  const unsigned long* parents = children + num_child;
1816 
1817  SpecialSetData* data = find_set_data( set );
1818  assert( NULL != data );
1819  if( NULL == data ) return MB_FAILURE;
1820 
1821  // Tag mattag;
1822  // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
1823  // int block;
1824  // if (MB_SUCCESS != iFace->tag_get_data(mattag, &set, 1, &block))
1825  // block = 0;
1826 
1827  // If either the current data or the new data is in ranged format,
1828  // then change the other to ranged format if it isn't already
1829  // in both cases when they differ, the data will end up "compressed range"
1830  std::vector< wid_t > tmp;
1831  if( ( flags & mhdf_SET_RANGE_BIT ) != ( data->setFlags & mhdf_SET_RANGE_BIT ) )
1832  {
1833  if( flags & mhdf_SET_RANGE_BIT )
1834  {
1835  tmp = data->contentIds;
1836  convert_to_ranged_ids( tmp.data(), tmp.size(), data->contentIds );
1837  data->setFlags |= mhdf_SET_RANGE_BIT;
1838  }
1839  else
1840  {
1841  tmp.clear();
1842  convert_to_ranged_ids( contents, num_content, tmp );
1843  num_content = tmp.size();
1844  if( sizeof( wid_t ) < sizeof( long ) )
1845  {
1846  size_t old_size = tmp.size();
1847  tmp.resize( sizeof( long ) * old_size / sizeof( wid_t ) );
1848  unsigned long* array = reinterpret_cast< unsigned long* >( tmp.data() );
1849  for( long i = ( (long)old_size ) - 1; i >= 0; --i )
1850  array[i] = tmp[i];
1851  contents = array;
1852  }
1853  else if( sizeof( wid_t ) > sizeof( long ) )
1854  {
1855  unsigned long* array = reinterpret_cast< unsigned long* >( tmp.data() );
1856  std::copy( tmp.begin(), tmp.end(), array );
1857  }
1858  contents = reinterpret_cast< unsigned long* >( tmp.data() );
1859  }
1860  }
1861 
1862  if( data->setFlags & mhdf_SET_RANGE_BIT )
1863  merge_ranged_ids( contents, num_content, data->contentIds );
1864  else
1865  merge_vector_ids( contents, num_content, data->contentIds );
1866 
1867  merge_vector_ids( children, num_child, data->childIds );
1868  merge_vector_ids( parents, num_parent, data->parentIds );
1869  return MB_SUCCESS;
1870 }
1871 
1872 ErrorCode WriteHDF5Parallel::communicate_shared_set_data( const Range& owned, const Range& remote )
1873 {
1874  ErrorCode rval;
1875  int mperr;
1876  const unsigned rank = myPcomm->proc_config().proc_rank();
1877  const MPI_Comm comm = myPcomm->proc_config().proc_comm();
1878 
1879  dbgOut.tprintf( 1, "COMMUNICATING SHARED SET DATA (%lu owned & %lu remote)\n", (unsigned long)owned.size(),
1880  (unsigned long)remote.size() );
1881 
1882  // Calculate the total number of messages to be in transit (send and receive)
1883  size_t nummess = 0;
1884  std::vector< unsigned > procs;
1885  ;
1886  Range shared( owned );
1887  shared.merge( remote );
1888  for( Range::iterator i = shared.begin(); i != shared.end(); ++i )
1889  {
1890  procs.clear();
1891  rval = myPcomm->get_entityset_procs( *i, procs );
1892  CHECK_MB( rval );
1893  nummess += procs.size();
1894  }
1895 
1896  // Choose a receive buffer size. We need 4*sizeof(long) minimum,
1897  // but that is almost useless so use 16*sizeof(long) as the minimum
1898  // instead. Choose an upper limit such that we don't exceed 32 MB
1899  // of allocated memory (unless we absolutely must to meet the minimum.)
1900  // Also, don't initially choose buffers larger than 128*sizeof(long).
1901  const size_t MAX_BUFFER_MEM = 32 * 1024 * 1024 / sizeof( long );
1902  // const size_t INIT_BUFFER_SIZE = 128;
1903  const size_t INIT_BUFFER_SIZE = 1024;
1904  const size_t MIN_BUFFER_SIZE = 16;
1905  size_t init_buff_size = INIT_BUFFER_SIZE;
1906  if( init_buff_size * nummess > MAX_BUFFER_MEM ) init_buff_size = MAX_BUFFER_MEM / nummess;
1907  if( init_buff_size < MIN_BUFFER_SIZE ) init_buff_size = MIN_BUFFER_SIZE;
1908 
1909  dbgOut.printf( 2, "Using buffer size of %lu for an expected message count of %lu\n", (unsigned long)init_buff_size,
1910  (unsigned long)nummess );
1911 
1912  // Count number of recvs
1913  size_t numrecv = 0;
1914  for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
1915  {
1916  procs.clear();
1917  rval = myPcomm->get_entityset_procs( *i, procs );
1918  CHECK_MB( rval );
1919  numrecv += procs.size();
1920  if( std::find( procs.begin(), procs.end(), rank ) != procs.end() ) --numrecv;
1921  }
1922 
1923  // Post receive buffers for all owned sets for all sharing procs
1924  std::vector< MPI_Request > recv_req( numrecv, MPI_REQUEST_NULL );
1925  std::vector< MPI_Request > lrecv_req( numrecv, MPI_REQUEST_NULL );
1926 
1927  std::vector< std::vector< unsigned long > > recv_buf( numrecv, std::vector< unsigned long >( init_buff_size ) );
1928  int idx = 0;
1929  for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
1930  {
1931  procs.clear();
1932  rval = myPcomm->get_entityset_procs( *i, procs );
1933  CHECK_MB( rval );
1934  for( size_t j = 0; j < procs.size(); ++j )
1935  {
1936  if( procs[j] == rank ) continue;
1937  int tag = ID_FROM_HANDLE( *i );
1938  if( *i != CREATE_HANDLE( MBENTITYSET, tag ) )
1939  {
1940 #ifndef NDEBUG
1941  abort();
1942 #endif
1943  CHECK_MB( MB_FAILURE );
1944  }
1945  dbgOut.printf( 5, "Posting buffer to receive set %d from proc %u\n", tag, procs[j] );
1946  mperr = MPI_Irecv( recv_buf[idx].data(), init_buff_size, MPI_UNSIGNED_LONG, procs[j], tag, comm,
1947  &recv_req[idx] );
1948  CHECK_MPI( mperr );
1949  ++idx;
1950  }
1951  }
1952  assert( (size_t)idx == numrecv );
1953 
1954  // Now send set data for all remote sets that I know about
1955  std::vector< MPI_Request > send_req( remote.size() );
1956  std::vector< std::vector< unsigned long > > send_buf( remote.size() );
1957  idx = 0;
1958  for( Range::iterator i = remote.begin(); i != remote.end(); ++i, ++idx )
1959  {
1960  send_buf[idx].resize( init_buff_size );
1961  rval = pack_set( i, send_buf[idx].data(), init_buff_size );
1962  CHECK_MB( rval );
1963  EntityHandle remote_handle;
1964  unsigned owner;
1965  rval = myPcomm->get_entityset_owner( *i, owner, &remote_handle );
1966  CHECK_MB( rval );
1967 
1968  int tag = ID_FROM_HANDLE( remote_handle );
1969  assert( remote_handle == CREATE_HANDLE( MBENTITYSET, tag ) );
1970  dbgOut.printf( 5, "Sending %lu values for set %d to proc %u\n",
1971  send_buf[idx][1] + send_buf[idx][2] + send_buf[idx][3] + 4, tag, owner );
1972  mperr = MPI_Isend( send_buf[idx].data(), init_buff_size, MPI_UNSIGNED_LONG, owner, tag, comm, &send_req[idx] );
1973  CHECK_MPI( mperr );
1974  }
1975 
1976  // Tag mattag;
1977  // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
1978 
1979  // Now initialize local data for managing contents of owned, shared sets
1980  assert( specialSets.empty() );
1981  specialSets.clear();
1982  specialSets.reserve( owned.size() );
1983  for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
1984  {
1985  // int block;
1986  // if (MB_SUCCESS != iFace->tag_get_data(mattag, &*i, 1, &block))
1987  // block = 0;
1988  // std::vector<int> ids;
1989 
1990  SpecialSetData data;
1991  data.setHandle = *i;
1992  rval = iFace->get_meshset_options( *i, data.setFlags );
1993  CHECK_MB( rval );
1994  specialSets.push_back( data );
1995  std::vector< EntityHandle > list;
1996  if( data.setFlags & MESHSET_ORDERED )
1997  {
1998  list.clear();
1999  rval = iFace->get_entities_by_handle( *i, list );
2000  CHECK_MB( rval );
2001  rval = vector_to_id_list( list, specialSets.back().contentIds, true );
2002  CHECK_MB( rval );
2003  // if (block)
2004  // get_global_ids(iFace, &list[0], list.size(), MESHSET_ORDERED, ids);
2005  }
2006  else
2007  {
2008  Range range;
2009  rval = iFace->get_entities_by_handle( *i, range );
2010  CHECK_MB( rval );
2011  bool ranged;
2012  rval = range_to_blocked_list( range, specialSets.back().contentIds, ranged );
2013  if( ranged ) specialSets.back().setFlags |= mhdf_SET_RANGE_BIT;
2014  // if (block) {
2015  // std::vector<EntityHandle> tmp;
2016  // for (Range::const_pair_iterator pi = range.const_pair_begin(); pi !=
2017  // range.const_pair_end(); ++pi) {
2018  // tmp.push_back(pi->first);
2019  // tmp.push_back(pi->second);
2020  // }
2021  // get_global_ids(iFace, &tmp[0], tmp.size(), ranged ? 0 : MESHSET_ORDERED, ids);
2022  //}
2023  }
2024 
2025  list.clear();
2026  rval = iFace->get_parent_meshsets( *i, list );
2027  CHECK_MB( rval );
2028  rval = vector_to_id_list( list, specialSets.back().parentIds, true );
2029  CHECK_MB( rval );
2030  rval = iFace->get_child_meshsets( *i, list );
2031  CHECK_MB( rval );
2032  rval = vector_to_id_list( list, specialSets.back().childIds, true );
2033  CHECK_MB( rval );
2034  }
2035 
2036  // Process received buffers, repost larger buffers where necessary
2037  size_t remaining = numrecv;
2038  numrecv = 0;
2039  while( remaining-- )
2040  {
2041  std::vector< unsigned long > dead;
2042  MPI_Status status;
2043  mperr = MPI_Waitany( recv_req.size(), recv_req.data(), &idx, &status );
2044  CHECK_MPI( mperr );
2045  EntityHandle handle = CREATE_HANDLE( MBENTITYSET, status.MPI_TAG );
2046  std::vector< unsigned long >& buff = recv_buf[idx];
2047  size_t size = buff[1] + buff[2] + buff[3] + 4;
2048  dbgOut.printf( 5, "Received %lu values for set %d from proc %d\n", (unsigned long)size, status.MPI_TAG,
2049  status.MPI_SOURCE );
2050  if( size <= init_buff_size )
2051  {
2052  rval = unpack_set( handle, buff.data(), init_buff_size );
2053  CHECK_MB( rval );
2054  dead.swap( buff ); // Release memory
2055  }
2056  else
2057  {
2058  // Data was too big for init_buff_size
2059  // repost with larger buffer
2060  buff.resize( size );
2061  dbgOut.printf( 5, "Re-Posting buffer to receive set %d from proc %d with size %lu\n", status.MPI_TAG,
2062  status.MPI_SOURCE, (unsigned long)size );
2063  mperr = MPI_Irecv( buff.data(), size, MPI_UNSIGNED_LONG, status.MPI_SOURCE, status.MPI_TAG, comm,
2064  &lrecv_req[idx] );
2065  CHECK_MPI( mperr );
2066  ++numrecv;
2067  }
2068  recv_req[idx] = MPI_REQUEST_NULL;
2069  }
2070 
2071  // Wait for sends to complete
2072  MPI_Waitall( send_req.size(), send_req.data(), MPI_STATUSES_IGNORE );
2073 
2074  // Re-send sets that didn't fit initial buffer size
2075  idx = 0;
2076  for( Range::iterator i = remote.begin(); i != remote.end(); ++i, ++idx )
2077  {
2078  std::vector< unsigned long >& buff = send_buf[idx];
2079  size_t size = buff[1] + buff[2] + buff[3] + 4;
2080  if( size <= init_buff_size ) continue;
2081 
2082  buff.resize( size );
2083  rval = pack_set( i, buff.data(), size );
2084  CHECK_MB( rval );
2085  EntityHandle remote_handle;
2086  unsigned owner;
2087  rval = myPcomm->get_entityset_owner( *i, owner, &remote_handle );
2088  CHECK_MB( rval );
2089 
2090  int tag = ID_FROM_HANDLE( remote_handle );
2091  assert( remote_handle == CREATE_HANDLE( MBENTITYSET, tag ) );
2092  dbgOut.printf( 5, "Sending %lu values for set %d to proc %u\n", (unsigned long)size, tag, owner );
2093  mperr = MPI_Isend( buff.data(), size, MPI_UNSIGNED_LONG, owner, tag, comm, &send_req[idx] );
2094  CHECK_MPI( mperr );
2095  }
2096 
2097  // Process received buffers
2098  remaining = numrecv;
2099  while( remaining-- )
2100  {
2101  std::vector< unsigned long > dead;
2102  MPI_Status status;
2103  mperr = MPI_Waitany( lrecv_req.size(), lrecv_req.data(), &idx, &status );
2104  CHECK_MPI( mperr );
2105  EntityHandle handle = CREATE_HANDLE( MBENTITYSET, status.MPI_TAG );
2106  std::vector< unsigned long >& buff = recv_buf[idx];
2107  dbgOut.printf( 5, "Received %lu values for set %d from proc %d\n", 4 + buff[1] + buff[2] + buff[3],
2108  status.MPI_TAG, status.MPI_SOURCE );
2109  rval = unpack_set( handle, buff.data(), buff.size() );
2110  CHECK_MB( rval );
2111  dead.swap( buff ); // Release memory
2112 
2113  lrecv_req[idx] = MPI_REQUEST_NULL;
2114  }
2115 
2116  // Wait for sends to complete
2117  MPI_Waitall( send_req.size(), send_req.data(), MPI_STATUSES_IGNORE );
2118 
2119  return MB_SUCCESS;
2120 }
2121 
2122 ErrorCode WriteHDF5Parallel::create_meshset_tables( double* times )
2123 {
2124  Range::const_iterator riter;
2125  const unsigned rank = myPcomm->proc_config().proc_rank();
2126 
2127  START_SERIAL;
2128  print_type_sets( iFace, &dbgOut, setSet.range );
2129  END_SERIAL;
2130  CpuTimer timer;
2131 
2132  // Remove remote sets from setSets
2133  Range shared, owned, remote;
2134  ErrorCode rval = myPcomm->get_shared_sets( shared );
2135  CHECK_MB( rval );
2136  shared = intersect( shared, setSet.range );
2137  rval = myPcomm->get_owned_sets( rank, owned );
2138  CHECK_MB( rval );
2139  owned = intersect( owned, setSet.range );
2140  remote = subtract( shared, owned );
2141  setSet.range = subtract( setSet.range, remote );
2142 
2143  // Create set meta table
2144  struct SetDescCreator : public DataSetCreator
2145  {
2146  ErrorCode operator()( WriteHDF5* writer, long size, const ExportSet*, long& start_id ) const
2147  {
2148  return writer->create_set_meta( size, start_id );
2149  }
2150  };
2151  long count = setSet.range.size();
2152  rval = create_dataset( 1, &count, &setSet.offset, &setSet.max_num_ents, &setSet.total_num_ents, SetDescCreator(),
2153  NULL, &setSet.first_id );
2154  CHECK_MB( rval );
2155  writeSets = setSet.max_num_ents > 0;
2156 
2157  rval = assign_ids( setSet.range, setSet.first_id + setSet.offset );
2158  CHECK_MB( rval );
2159  if( times ) times[SET_OFFSET_TIME] = timer.time_elapsed();
2160 
2161  // Exchange file IDS for sets between all procs
2162  rval = communicate_shared_set_ids( owned, remote );
2163  CHECK_MB( rval );
2164  if( times ) times[SHARED_SET_IDS] = timer.time_elapsed();
2165 
2166  // Communicate remote set contents, children, etc.
2167  rval = communicate_shared_set_data( owned, remote );
2168  CHECK_MB( rval );
2169  if( times ) times[SHARED_SET_CONTENTS] = timer.time_elapsed();
2170 
2171  // Communicate counts for owned sets
2172  long data_counts[3]; // { #contents, #children, #parents }
2173  rval = count_set_size( setSet.range, data_counts[0], data_counts[1], data_counts[2] );
2174  CHECK_MB( rval );
2175  if( times ) times[SET_OFFSET_TIME] += timer.time_elapsed();
2176 
2177  long offsets[3], max_counts[3], totals[3];
2178  rval = create_dataset( 3, data_counts, offsets, max_counts, totals );
2179  CHECK_MB( rval );
2180 
2181  // Create the datasets
2182  if( 0 == myPcomm->proc_config().proc_rank() )
2183  {
2184  rval = create_set_tables( totals[0], totals[1], totals[2] );
2185  CHECK_MB( rval );
2186  }
2187 
2188  // Store communicated global data
2189  setContentsOffset = offsets[0];
2190  setChildrenOffset = offsets[1];
2191  setParentsOffset = offsets[2];
2192  maxNumSetContents = max_counts[0];
2193  maxNumSetChildren = max_counts[1];
2194  maxNumSetParents = max_counts[2];
2195  writeSetContents = totals[0] > 0;
2196  writeSetChildren = totals[1] > 0;
2197  writeSetParents = totals[2] > 0;
2198 
2199  dbgOut.printf( 2, "set contents: %ld local, %ld global, offset = %ld\n", data_counts[0], totals[0], offsets[0] );
2200  dbgOut.printf( 2, "set children: %ld local, %ld global, offset = %ld\n", data_counts[1], totals[1], offsets[1] );
2201  dbgOut.printf( 2, "set parents: %ld local, %ld global, offset = %ld\n", data_counts[2], totals[2], offsets[2] );
2202 
2203  return MB_SUCCESS;
2204 }
2205 
2206 void WriteHDF5Parallel::remove_remote_entities( EntityHandle relative, Range& range )
2207 {
2208  Range result;
2209  result.merge( intersect( range, nodeSet.range ) );
2210  result.merge( intersect( range, setSet.range ) );
2211  for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
2212  result.merge( intersect( range, eiter->range ) );
2213 
2214  // result.merge(intersect(range, myParallelSets));
2215  Range sets;
2216  int junk;
2217  sets.merge( Range::lower_bound( range.begin(), range.end(), CREATE_HANDLE( MBENTITYSET, 0, junk ) ), range.end() );
2218  remove_remote_sets( relative, sets );
2219  result.merge( sets );
2220  range.swap( result );
2221 }
2222 
2223 void WriteHDF5Parallel::remove_remote_sets( EntityHandle /* relative */, Range& range )
2224 {
2225  Range result( intersect( range, setSet.range ) );
2226  // Store the non-intersecting entities separately if needed
2227  // Range remaining(subtract(range, result));
2228  range.swap( result );
2229 }
2230 
2231 void WriteHDF5Parallel::remove_remote_entities( EntityHandle relative, std::vector< EntityHandle >& vect )
2232 {
2233  Range intrsct;
2234  for( std::vector< EntityHandle >::const_iterator iter = vect.begin(); iter != vect.end(); ++iter )
2235  intrsct.insert( *iter );
2236  remove_remote_entities( relative, intrsct );
2237 
2238  unsigned int read, write;
2239  for( read = write = 0; read < vect.size(); ++read )
2240  {
2241  if( intrsct.find( vect[read] ) != intrsct.end() )
2242  {
2243  if( read != write ) vect[write] = vect[read];
2244  ++write;
2245  }
2246  }
2247  if( write != vect.size() ) vect.resize( write );
2248 }
2249 
2250 void WriteHDF5Parallel::remove_remote_sets( EntityHandle relative, std::vector< EntityHandle >& vect )
2251 {
2252  Range intrsct;
2253  for( std::vector< EntityHandle >::const_iterator iter = vect.begin(); iter != vect.end(); ++iter )
2254  intrsct.insert( *iter );
2255  remove_remote_sets( relative, intrsct );
2256 
2257  unsigned int read, write;
2258  for( read = write = 0; read < vect.size(); ++read )
2259  {
2260  if( intrsct.find( vect[read] ) != intrsct.end() )
2261  {
2262  if( read != write ) vect[write] = vect[read];
2263  ++write;
2264  }
2265  }
2266  if( write != vect.size() ) vect.resize( write );
2267 }
2268 
2269 ErrorCode WriteHDF5Parallel::exchange_file_ids( const Range& nonlocal )
2270 {
2271  ErrorCode rval;
2272 
2273  // For each entity owned on the interface, write its file id to
2274  // a tag. The sets of entities to be written should already contain
2275  // only owned entities so by intersecting with them we not only
2276  // filter by entities to be written, but also restrict to entities
2277  // owned by the proc
2278 
2279  // Get list of interface entities
2280  Range imesh, tmp;
2281  for( std::list< ExportSet >::reverse_iterator i = exportList.rbegin(); i != exportList.rend(); ++i )
2282  {
2283  tmp.clear();
2284  rval = myPcomm->filter_pstatus( i->range, PSTATUS_SHARED, PSTATUS_AND, -1, &tmp );
2285  if( MB_SUCCESS != rval ) return error( rval );
2286  imesh.merge( tmp );
2287  }
2288  tmp.clear();
2289  rval = myPcomm->filter_pstatus( nodeSet.range, PSTATUS_SHARED, PSTATUS_AND, -1, &tmp );
2290  if( MB_SUCCESS != rval ) return error( rval );
2291  imesh.merge( tmp );
2292 
2293  // Create tag to store file IDs
2294  EntityHandle default_val = 0;
2295  Tag file_id_tag = 0;
2296  rval = iFace->tag_get_handle( "__hdf5_ll_fileid", 1, MB_TYPE_HANDLE, file_id_tag, MB_TAG_DENSE | MB_TAG_CREAT,
2297  &default_val );
2298  if( MB_SUCCESS != rval ) return error( rval );
2299 
2300  // Copy file IDs into tag
2301  std::vector< EntityHandle > file_id_vect( imesh.size() );
2302  Range::const_iterator i;
2303  std::vector< EntityHandle >::iterator j = file_id_vect.begin();
2304  for( i = imesh.begin(); i != imesh.end(); ++i, ++j )
2305  {
2306  *j = idMap.find( *i );
2307  if( !*j )
2308  {
2309  iFace->tag_delete( file_id_tag );
2310  return error( MB_FAILURE );
2311  }
2312  }
2313  rval = iFace->tag_set_data( file_id_tag, imesh, file_id_vect.data() );
2314  if( MB_SUCCESS != rval )
2315  {
2316  iFace->tag_delete( file_id_tag );
2317  return error( rval );
2318  }
2319 
2320  // Do communication
2321  rval = myPcomm->exchange_tags( file_id_tag, imesh );
2322  if( MB_SUCCESS != rval )
2323  {
2324  iFace->tag_delete( file_id_tag );
2325  return error( rval );
2326  }
2327 
2328  // Copy file IDs from tag into idMap for remote entities
2329  file_id_vect.resize( nonlocal.size() );
2330  rval = iFace->tag_get_data( file_id_tag, nonlocal, file_id_vect.data() );
2331  if( MB_SUCCESS != rval )
2332  {
2333  iFace->tag_delete( file_id_tag );
2334  return error( rval );
2335  }
2336 
2337  j = file_id_vect.begin();
2338  for( i = nonlocal.begin(); i != nonlocal.end(); ++i, ++j )
2339  {
2340  if( *j == 0 )
2341  {
2342  int owner = -1;
2343  myPcomm->get_owner( *i, owner );
2344  const char* name = CN::EntityTypeName( TYPE_FROM_HANDLE( *i ) );
2345  int id = ID_FROM_HANDLE( *i );
2346  MB_SET_ERR_CONT( "Process " << myPcomm->proc_config().proc_rank()
2347  << " did not receive valid id handle for shared " << name << " " << id
2348  << " owned by process " << owner );
2349  dbgOut.printf( 1,
2350  "Did not receive valid remote id for "
2351  "shared %s %d owned by process %d",
2352  name, id, owner );
2353  iFace->tag_delete( file_id_tag );
2354  return error( MB_FAILURE );
2355  }
2356  else
2357  {
2358  if( !idMap.insert( *i, *j, 1 ).second )
2359  {
2360  iFace->tag_delete( file_id_tag );
2361  return error( MB_FAILURE );
2362  }
2363  }
2364  }
2365 
2366 #ifndef NDEBUG
2367  // Check that writer is correct with regards to which entities
2368  // that it owns by verifying that the file ids that we thought
2369  // we were sending where not received instead
2370  file_id_vect.resize( imesh.size() );
2371  rval = iFace->tag_get_data( file_id_tag, imesh, file_id_vect.data() );
2372  if( MB_SUCCESS != rval )
2373  {
2374  iFace->tag_delete( file_id_tag );
2375  return error( rval );
2376  }
2377  int invalid_count = 0;
2378  j = file_id_vect.begin();
2379  for( i = imesh.begin(); i != imesh.end(); ++i, ++j )
2380  {
2381  EntityHandle h = idMap.find( *i );
2382  if( *j != h )
2383  {
2384  ++invalid_count;
2385  dbgOut.printf( 1, "Conflicting ownership for %s %ld\n", CN::EntityTypeName( TYPE_FROM_HANDLE( *i ) ),
2386  (long)ID_FROM_HANDLE( *i ) );
2387  }
2388  }
2389  if( invalid_count )
2390  {
2391  iFace->tag_delete( file_id_tag );
2392  MB_SET_ERR( MB_FAILURE, invalid_count << " entities with conflicting ownership found by process "
2393  << myPcomm->proc_config().proc_rank()
2394  << ". This will result in duplicate entities written to file" );
2395  }
2396 #endif
2397 
2398  return iFace->tag_delete( file_id_tag );
2399 }
2400 
2401 void WriteHDF5Parallel::print_times( const double* times ) const
2402 {
2403  if( !myPcomm )
2404  {
2405  WriteHDF5::print_times( times );
2406  }
2407  else
2408  {
2409  double recv[NUM_TIMES];
2410  MPI_Reduce( (void*)times, recv, NUM_TIMES, MPI_DOUBLE, MPI_MAX, 0, myPcomm->proc_config().proc_comm() );
2411  if( 0 == myPcomm->proc_config().proc_rank() ) WriteHDF5::print_times( recv );
2412  }
2413 }
2414 
2415 } // namespace moab