Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReorderTool.cpp
Go to the documentation of this file.
1 /** \file ReorderTool.cpp
2  * \author Jason Kraftcheck
3  * \date 2011-05-23
4  */
5 
6 #include "moab/ReorderTool.hpp"
7 #include "moab/Core.hpp"
8 #include "moab/Range.hpp"
9 
10 #include "SequenceManager.hpp"
11 #include "TypeSequenceManager.hpp"
12 #include "EntitySequence.hpp"
13 
14 #include <algorithm>
15 #include <numeric>
16 #include <set>
17 #include <iostream>
18 
19 namespace moab
20 {
21 
22 // no-op function as a convenient spot to set a breakpoint
23 static inline ErrorCode error( ErrorCode val )
24 {
25  return val;
26 }
27 
28 #define CHKERR \
29  if( MB_SUCCESS != rval ) return error( rval )
30 
31 #define UNRECOVERABLE( ERRCODE ) \
32  do \
33  { \
34  if( MB_SUCCESS != ( ERRCODE ) ) \
35  { \
36  error( ( ERRCODE ) ); \
37  std::cerr << "Unreconverable error during mesh reorder." << std::endl \
38  << "Error Code " << ( ERRCODE ) << " at " << __FILE__ << ":" << __LINE__ << std::endl; \
39  std::cerr.flush(); \
40  abort(); \
41  } \
42  } while( false )
43 
44 static ErrorCode check_tag_type( Interface* moab, Tag tag, DataType exp_type, int exp_size )
45 {
46  ErrorCode rval;
47  DataType act_type;
48  int act_size;
49 
50  rval = moab->tag_get_data_type( tag, act_type );CHKERR;
51 
52  rval = moab->tag_get_bytes( tag, act_size );CHKERR;
53 
54  if( act_type != exp_type || act_size != exp_size ) return MB_TYPE_OUT_OF_RANGE;
55 
56  return MB_SUCCESS;
57 }
58 
59 ErrorCode ReorderTool::handle_order_from_int_tag( Tag tag, int skip_value, Tag& new_handles )
60 {
61  ErrorCode rval;
62 
63  // check that input tag handles are what we expect
64  rval = check_tag_type( mMB, tag, MB_TYPE_INTEGER, sizeof( int ) );CHKERR;
65  EntityHandle zero = 0;
66  rval = mMB->tag_get_handle( 0, 1, MB_TYPE_HANDLE, new_handles, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_EXCL, &zero );CHKERR;
67 
68  // Can only reorder within same type/connectivity (or vertex/num_coords)
69  // groupings. Call helper function for each such grouping.
70  for( EntityType t = MBVERTEX; t < MBENTITYSET; ++t )
71  {
72  // Get list of all connectivity lengths (or vertex dimensions)
73  // that exist for type t.
76  std::set< int > values;
77  for( i = seqs.begin(); i != seqs.end(); ++i )
78  {
79  EntitySequence* seq = *i;
80  // 0 values per entity implies structured data, which
81  // we cannot reorder. Skip those.
82  if( t == MBVERTEX || 0 < seq->values_per_entity() ) values.insert( seq->values_per_entity() );
83  }
84 
85  // Call helper function for each (type,size) tuple.
86  std::set< int >::iterator j;
87  for( j = values.begin(); j != values.end(); ++j )
88  {
89  rval = handle_order_from_int_tag( t, *j, tag, skip_value, new_handles );
90  if( MB_SUCCESS != rval )
91  {
92  mMB->tag_delete( new_handles );
93  return error( rval );
94  }
95  }
96  }
97 
98  return MB_SUCCESS;
99 }
100 
101 void ReorderTool::get_entities( EntityType t, int vals_per_ent, Range& entities )
102 {
103  Range::iterator hint = entities.begin();
104  ;
107  for( s = seqs.begin(); s != seqs.end(); ++s )
108  {
109  EntitySequence* seq = *s;
110  if( seq->values_per_entity() == vals_per_ent )
111  hint = entities.insert( hint, seq->start_handle(), seq->end_handle() );
112  }
113 }
114 
116  int vals_per_ent,
117  Tag tag,
118  int skip_value,
119  Tag new_handles )
120 {
121  ErrorCode rval;
122 
123  // check that input tag handles are what we expect
124  rval = check_tag_type( mMB, tag, MB_TYPE_INTEGER, sizeof( int ) );CHKERR;
125  rval = check_tag_type( mMB, new_handles, MB_TYPE_HANDLE, sizeof( EntityHandle ) );CHKERR;
126 
127  // get entities to re-order
128  Range entities;
129  get_entities( t, vals_per_ent, entities );
130 
131  // get entity order values
132  std::vector< int > sortvals( entities.size() );
133  rval = mMB->tag_get_data( tag, entities, &sortvals[0] );CHKERR;
134 
135  // remove any entities for which value is skip_value
136  size_t r = 0, w = 0;
137  for( Range::iterator i = entities.begin(); i != entities.end(); ++r )
138  {
139  if( sortvals[r] == skip_value )
140  i = entities.erase( i );
141  else
142  {
143  sortvals[w++] = sortvals[r];
144  ++i;
145  }
146  }
147  sortvals.resize( w );
148 
149  // sort
150  std::sort( sortvals.begin(), sortvals.end() );
151  // Convert to unique list and offsets for each value
152  // When done, sortvals will contain unique, sortvals and
153  // offsets will contain, for each unique value in sortvals,
154  // the number of values that occured in the non-unique list
155  // before the first instance of that value.
156  std::vector< size_t > offsets;
157  offsets.push_back( 0 );
158  offsets.push_back( 1 );
159  for( w = 0, r = 1; r < sortvals.size(); ++r )
160  {
161  if( sortvals[r] == sortvals[w] )
162  {
163  ++offsets.back();
164  }
165  else
166  {
167  ++w;
168  sortvals[w] = sortvals[r];
169  offsets.push_back( offsets.back() + 1 );
170  }
171  }
172  ++w;
173  assert( w + 1 == offsets.size() );
174  sortvals.resize( w );
175 
176  // Tag each entity with its new handle
177  for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
178  {
179  int val;
180  rval = mMB->tag_get_data( tag, &*i, 1, &val );CHKERR;
181  w = std::lower_bound( sortvals.begin(), sortvals.end(), val ) - sortvals.begin();
182  assert( w < sortvals.size() );
183  size_t offset = offsets[w];
184  ++offsets[w];
185  // should maybe copy range into array to avoid possibly n^2 behavior here
186  EntityHandle h = *( entities.begin() + offset );
187  rval = mMB->tag_set_data( new_handles, &*i, 1, &h );CHKERR;
188  }
189 
190  return MB_SUCCESS;
191 }
192 
194 {
195  ErrorCode rval;
196 
197  if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
198 
199  Tag order_tag;
200  const int negone = -1;
201  rval = mMB->tag_get_handle( 0, 1, MB_TYPE_INTEGER, order_tag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_EXCL, &negone );
202  if( MB_SUCCESS != rval )
203  {
204  mMB->tag_delete( handle_tag );
205  handle_tag = 0;
206  return error( rval );
207  }
208 
209  std::vector< std::vector< EntityHandle >* > data;
210  rval = int_order_from_sets_and_adj( sets, order_tag, negone, data );
211  for( size_t i = 0; i < data.size(); ++i )
212  delete data[i];
213  if( MB_SUCCESS != rval )
214  {
215  mMB->tag_delete( order_tag );
216  return error( rval );
217  }
218 
219  rval = handle_order_from_int_tag( order_tag, negone, handle_tag );
220  if( MB_SUCCESS != rval )
221  {
222  mMB->tag_delete( order_tag );
223  return error( rval );
224  }
225 
226  rval = mMB->tag_delete( order_tag );
227  if( MB_SUCCESS != rval ) return error( rval );
228 
229  return MB_SUCCESS;
230 }
231 
232 // Compare function to use for a map keyed on pointers to sorted vectors
234 {
235  bool operator()( const std::vector< EntityHandle >* v1, const std::vector< EntityHandle >* v2 ) const
236  {
237  std::vector< EntityHandle >::const_iterator i1, i2;
238  for( i1 = v1->begin(), i2 = v2->begin(); i1 != v1->end(); ++i1, ++i2 )
239  {
240  if( i2 == v2->end() || *i1 > *i2 )
241  return false;
242  else if( *i1 < *i2 )
243  return true;
244  }
245  return i2 != v2->end();
246  }
247 };
248 
250  Tag order_tag,
251  int skip_val,
252  std::vector< std::vector< EntityHandle >* >& revMap )
253 {
254  ErrorCode rval;
255 
256  if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
257 
258  rval = check_tag_type( mMB, order_tag, MB_TYPE_INTEGER, sizeof( int ) );CHKERR;
259 
260  // Compare function to use for a map keyed on pointers to sorted vectors
261  CompSortedVect lessthan;
262  // Map from sorted list of handles to index
263  std::map< std::vector< EntityHandle >*, int, CompSortedVect > forMap;
264  std::map< std::vector< EntityHandle >*, int, CompSortedVect >::iterator fiter, fiter2;
265  std::vector< EntityHandle > sharing; // tmp storage for entity
266 
267  // for each set
268  for( Range::iterator s = sets.begin(); s != sets.end(); ++s )
269  {
270 
271  // gather up all entities and adjacencies
272  Range tmp, ents, adj[4]; // indexed by dimension
273  for( int dim = 0; dim < 4; ++dim )
274  {
275  rval = mMB->get_entities_by_dimension( *s, dim, tmp );CHKERR;
276  for( int ldim = 0; ldim < dim; ++ldim )
277  {
278  rval = mMB->get_adjacencies( tmp, ldim, false, adj[ldim], Interface::UNION );CHKERR;
279  }
280  for( int udim = dim + 1; udim <= 3; ++udim )
281  {
282  rval = mMB->get_adjacencies( tmp, udim, false, adj[udim], Interface::UNION );CHKERR;
283  }
284  ents.merge( tmp );
285  tmp.clear();
286  }
287  for( int dim = 0; dim < 4; ++dim )
288  {
289  ents.merge( adj[dim] );
290  adj[dim].clear();
291  }
292 
293  // process each entity
294  for( Range::iterator e = ents.begin(); e != ents.end(); ++e )
295  {
296  int val;
297  rval = mMB->tag_get_data( order_tag, &*e, 1, &val );CHKERR;
298 
299  // If this entity is already in one or more of the sets (either
300  // directly or through adjacency) then get the existing list of
301  // sets and append this set handle (we are iterating over sets
302  // in sorted order, so appending should aways maintain a sorted
303  // list.)
304  sharing.clear();
305  if( val != skip_val )
306  {
307  sharing = *revMap[val];
308  assert( std::lower_bound( sharing.begin(), sharing.end(), *s ) == sharing.end() );
309  }
310  sharing.push_back( *s );
311 
312  // Check if new sharing list already exists in forward map
313  fiter = forMap.lower_bound( &sharing );
314  if( fiter == forMap.end() || lessthan( fiter->first, &sharing ) )
315  {
316  // Add new sharing list to forward and reverse maps.
317  std::vector< EntityHandle >* newvec = new std::vector< EntityHandle >;
318  newvec->swap( sharing );
319  if( (int)revMap.size() == skip_val ) revMap.push_back( 0 );
320  fiter2 =
321  forMap.insert( fiter, std::pair< std::vector< EntityHandle >*, int >( newvec, revMap.size() ) );
322  assert( fiter2 != fiter );
323  fiter = fiter2;
324  revMap.push_back( newvec );
325  }
326 
327  // Update index on entity
328  val = fiter->second;
329  rval = mMB->tag_set_data( order_tag, &*e, 1, &val );CHKERR;
330  }
331  }
332 
333  return MB_SUCCESS;
334 }
335 
337  const Range& old_handles,
338  std::vector< EntityHandle >& new_handles )
339 {
340  new_handles.resize( old_handles.size() );
341  ErrorCode rval = mMB->tag_get_data( tag, old_handles, ( new_handles.empty() ) ? NULL : &new_handles[0] );CHKERR;
342 
343  Range::const_iterator it1 = old_handles.begin();
344  std::vector< EntityHandle >::iterator it2 = new_handles.begin();
345  for( ; it1 != old_handles.end(); ++it1, ++it2 )
346  if( 0 == *it2 ) *it2 = *it1;
347 
348  return MB_SUCCESS;
349 }
350 
352  const std::vector< EntityHandle >& old_handles,
353  std::vector< EntityHandle >& new_handles )
354 {
355  new_handles.resize( old_handles.size() );
356  return get_reordered_handles( tag, &old_handles[0], &new_handles[0], old_handles.size() );
357 }
358 
360  const EntityHandle* old_handles,
361  EntityHandle* new_handles,
362  size_t num_handles )
363 {
364  ErrorCode rval = mMB->tag_get_data( tag, old_handles, num_handles, new_handles );CHKERR;
365 
366  for( size_t i = 0; i < num_handles; ++i )
367  if( 0 == new_handles[i] ) new_handles[i] = old_handles[i];
368 
369  return MB_SUCCESS;
370 }
371 
372 ErrorCode ReorderTool::get_new_handles( Tag tag, Range& old_handles, std::vector< EntityHandle >& newhandles )
373 {
374  // get new handles for tagged entities
375  newhandles.resize( old_handles.size() );
376  ErrorCode rval = mMB->tag_get_data( tag, old_handles, ( newhandles.empty() ) ? NULL : &newhandles[0] );CHKERR;
377 
378  // remove entities that were not reordered
379  Range::iterator i = old_handles.begin();
380  size_t w = 0;
381  for( size_t r = 0; r < newhandles.size(); ++r )
382  {
383  if( 0 != newhandles[r] )
384  {
385  newhandles[w] = newhandles[r];
386  ++w;
387  ++i;
388  }
389  else
390  {
391  i = old_handles.erase( i );
392  }
393  }
394  newhandles.resize( w );
395  assert( newhandles.size() == old_handles.size() );
396  return MB_SUCCESS;
397 }
398 
400 {
401  ErrorCode rval;
402 
403  rval = check_tag_type( mMB, new_handles, MB_TYPE_HANDLE, sizeof( EntityHandle ) );CHKERR;
404  EntityHandle defval;
405  rval = mMB->tag_get_default_value( new_handles, &defval );CHKERR;
406  if( 0 != defval ) return error( MB_INDEX_OUT_OF_RANGE );
407 
408  // Can only reorder within same type/connectivity (or vertex/num_coords)
409  // groupings.
410  for( EntityType t = MBVERTEX; t < MBENTITYSET; ++t )
411  {
412  // Get list of all connectivity lengths (or vertex dimensions)
413  // that exist for type t.
416  std::set< int > values;
417  for( i = seqs.begin(); i != seqs.end(); ++i )
418  {
419  EntitySequence* seq = *i;
420  // 0 values per entity implies structured data, which
421  // we cannot reorder. Skip those.
422  if( t == MBVERTEX || 0 < seq->values_per_entity() ) values.insert( seq->values_per_entity() );
423  }
424 
425  // reorder primary data for each (type,size) tuple.
426  std::set< int >::iterator j;
427  for( j = values.begin(); j != values.end(); ++j )
428  {
429  Range entities;
430  get_entities( t, *j, entities );
431  std::vector< EntityHandle > handles;
432  rval = get_reordered_handles( new_handles, entities, handles );
433  UNRECOVERABLE( rval );
434 
435  if( t == MBVERTEX )
436  {
437  std::vector< double > coords( entities.size() * 3 );
438  rval = mMB->get_coords( entities, &coords[0] );
439  UNRECOVERABLE( rval );
440  rval = mMB->set_coords( &handles[0], handles.size(), &coords[0] );
441  UNRECOVERABLE( rval );
442  }
443  else
444  {
445  std::vector< EntityHandle > conn;
446  conn.reserve( entities.size() * *j );
447  std::vector< EntityHandle > old_handles;
448  old_handles.resize( entities.size() );
449  std::copy( entities.begin(), entities.end(), old_handles.begin() );
450  rval = mMB->get_connectivity( &old_handles[0], old_handles.size(), conn, false );
451  UNRECOVERABLE( rval );
452  old_handles.clear();
453  old_handles = conn;
454  rval = get_reordered_handles( new_handles, old_handles, conn );
455  UNRECOVERABLE( rval );
456  for( unsigned int h = 0; h < handles.size(); ++h )
457  {
458  rval = mMB->set_connectivity( handles[h], &conn[h * *j], *j );
459  UNRECOVERABLE( rval );
460  }
461  }
462  }
463  }
464 
465  // now update tag data
466  std::vector< Tag > tag_handles;
467  mMB->tag_get_tags( tag_handles );
468  for( size_t i = 0; i < tag_handles.size(); ++i )
469  {
470  Tag tag = tag_handles[i];
471  if( tag == new_handles ) // don't mess up mapping from old to new handles
472  continue;
473 
474  for( EntityType t = MBVERTEX; t <= MBENTITYSET; ++t )
475  {
476  rval = reorder_tag_data( t, new_handles, tag );
477  UNRECOVERABLE( rval );
478  }
479  }
480 
481  rval = update_set_contents( new_handles );
482  UNRECOVERABLE( rval );
483 
484  return MB_SUCCESS;
485 }
486 
487 ErrorCode ReorderTool::reorder_tag_data( EntityType etype, Tag new_handles, Tag tag )
488 {
489  ErrorCode rval;
490 
491  int tagsize;
492  DataType tagtype;
493  rval = mMB->tag_get_data_type( tag, tagtype );
494  if( MB_SUCCESS != rval ) return error( rval );
495  if( MB_TYPE_BIT == tagtype )
496  tagsize = 1;
497  else
498  {
499  rval = mMB->tag_get_bytes( tag, tagsize );
500  if( MB_VARIABLE_DATA_LENGTH == rval )
501  tagsize = -1;
502  else if( MB_SUCCESS != rval )
503  return error( rval );
504  }
505 
506  // we don't re-order set handles, so we don't care about sets
507  // unless the tag contains handles that need to be updated
508  if( MBENTITYSET == etype && MB_TYPE_HANDLE != tagtype ) return MB_SUCCESS;
509 
510  // get tagged entities
511  Range old_tagged;
512  rval = mMB->get_entities_by_type_and_tag( 0, etype, &tag, 0, 1, old_tagged );
513  if( MB_SUCCESS != rval && old_tagged.empty() ) return error( rval );
514 
515  // remove entities that were not reordered, unless the tag
516  // is handle type in which case we need to update the data
517  // for all entities, regardless of reordering.
518  std::vector< EntityHandle > newhandles;
519  if( MB_TYPE_HANDLE == tagtype )
520  rval = get_reordered_handles( new_handles, old_tagged, newhandles );
521  else
522  rval = get_new_handles( new_handles, old_tagged, newhandles );CHKERR;
523 
524  if( old_tagged.empty() ) return MB_SUCCESS;
525 
526  // get copy of all tag data
527  std::vector< unsigned char > buffer;
528  std::vector< const void* > pointers;
529  std::vector< int > sizes;
530  // if variable-length tag
531  if( -1 == tagsize )
532  {
533  pointers.resize( old_tagged.size() );
534  sizes.resize( pointers.size() );
535  rval = mMB->tag_get_by_ptr( tag, old_tagged, &pointers[0], &sizes[0] );CHKERR;
536  int total = std::accumulate( sizes.begin(), sizes.end(), 0 );
537  DataType dtype;
538  mMB->tag_get_data_type( tag, dtype );
539  int type_size;
540  switch( dtype )
541  {
542  case MB_TYPE_INTEGER:
543  type_size = sizeof( int );
544  break;
545  case MB_TYPE_DOUBLE:
546  type_size = sizeof( double );
547  break;
548  case MB_TYPE_HANDLE:
549  type_size = sizeof( EntityHandle );
550  break;
551  case MB_TYPE_BIT:
552  type_size = 1;
553  break;
554  case MB_TYPE_OPAQUE:
555  type_size = 1;
556  break;
557  default:
558  return MB_TYPE_OUT_OF_RANGE;
559  }
560  buffer.resize( total * type_size );
561  size_t off = 0;
562  for( size_t j = 0; j < pointers.size(); ++j )
563  {
564  memcpy( &buffer[off], pointers[j], type_size * sizes[j] );
565  pointers[j] = &buffer[off];
566  off += sizes[j] * type_size;
567  }
568  }
569  // if fixed-length tag
570  else
571  {
572  buffer.resize( old_tagged.size() * tagsize );
573  rval = mMB->tag_get_data( tag, old_tagged, &buffer[0] );CHKERR;
574  }
575 
576  // if handle tag, update tag values for reordered handles
577  if( MB_TYPE_HANDLE == tagtype )
578  {
579  assert( !( buffer.size() % sizeof( EntityHandle ) ) );
580  std::vector< unsigned char > buffer2( buffer.size() );
581  rval = get_reordered_handles( new_handles, reinterpret_cast< const EntityHandle* >( &buffer[0] ),
582  reinterpret_cast< EntityHandle* >( &buffer2[0] ),
583  buffer.size() / sizeof( EntityHandle ) );CHKERR;
584  // if var-length tag then do not do swap because pointers[] contains pointers
585  // into old buffer
586  if( -1 == tagsize )
587  memcpy( &buffer[0], &buffer2[0], buffer.size() );
588  else
589  buffer.swap( buffer2 );
590  }
591 
592  // store re-ordered tag data
593  if( -1 == tagsize )
594  {
595  rval = mMB->tag_set_by_ptr( tag, &newhandles[0], newhandles.size(), &pointers[0], &sizes[0] );
596  pointers.clear();
597  sizes.clear();
598  buffer.clear();
599  }
600  else
601  {
602  rval = mMB->tag_set_data( tag, &newhandles[0], newhandles.size(), &buffer[0] );
603  buffer.clear();
604  }
605  CHKERR;
606 
607  // all permutations should be cyclical, but not all permuted
608  // entities necessarily had tag values, so we might need to delete
609  // tags for some entities
610  std::sort( newhandles.begin(), newhandles.end() );
611  std::vector< EntityHandle >::iterator k = newhandles.begin();
612  Range::iterator i = old_tagged.begin();
613  while( i != old_tagged.end() )
614  {
615  while( k != newhandles.end() && *k < *i )
616  ++k;
617  if( k == newhandles.end() ) break;
618 
619  if( *i == *k ) // what old_tagged -= newhandles
620  i = old_tagged.erase( i );
621  else
622  ++i;
623  }
624 
625  if( !old_tagged.empty() )
626  {
627  rval = mMB->tag_delete_data( tag, old_tagged );CHKERR;
628  }
629 
630  return MB_SUCCESS;
631 }
632 
634 {
635  Range sets;
637 
638  std::vector< EntityHandle > old_handles, new_handles;
639  for( Range::iterator i = sets.begin(); i != sets.end(); ++i )
640  {
641  // If set is un-ordered...
642  unsigned opts = 0;
643  mMB->get_meshset_options( *i, opts );
644  if( !( opts & MESHSET_ORDERED ) )
645  {
646  Range contents;
647  rval = mMB->get_entities_by_handle( *i, contents );CHKERR;
648 
649  rval = get_new_handles( nh_tag, contents, new_handles );CHKERR;
650 
651  Range replace;
652  std::sort( new_handles.begin(), new_handles.end() );
653  Range::iterator hint = replace.begin();
654  for( size_t j = 0; j < new_handles.size(); ++j )
655  hint = replace.insert( hint, new_handles[j] );
656  Range common = intersect( contents, replace );
657  contents -= common;
658  replace -= common;
659  assert( contents.size() == replace.size() );
660  if( !contents.empty() )
661  {
662  rval = mMB->remove_entities( *i, contents );CHKERR;
663  rval = mMB->add_entities( *i, replace );
664  }
665  }
666 
667  // If set is ordered...
668  else
669  {
670  // get set contents
671  old_handles.clear();
672  rval = mMB->get_entities_by_handle( *i, old_handles );CHKERR;
673 
674  // get new handles from old contained handles
675  rval = get_reordered_handles( nh_tag, old_handles, new_handles );CHKERR;
676 
677  rval = mMB->clear_meshset( &*i, 1 );CHKERR;
678 
679  rval = mMB->add_entities( *i, &new_handles[0], new_handles.size() );CHKERR;
680  }
681  } // for each set
682 
683  return MB_SUCCESS;
684 }
685 
686 } // namespace moab