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