ReorderTool.cpp

Go to the documentation of this file.

/** \file ReorderTool.cpp
 * \author Jason Kraftcheck
 * \date 2011-05-23
 */
 
#include "moab/ReorderTool.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
 
#include "SequenceManager.hpp"
#include "TypeSequenceManager.hpp"
#include "EntitySequence.hpp"
 
#include <algorithm>
#include <numeric>
#include <set>
#include <iostream>
 
namespace moab
{
 
// no-op function as a convenient spot to set a breakpoint
static inline ErrorCode error( ErrorCode val )
{
 return val;
}
 
#define CHKERR \
 if( MB_SUCCESS != rval ) return error( rval )
 
#define UNRECOVERABLE( ERRCODE ) \
 do \
 { \
 if( MB_SUCCESS != ( ERRCODE ) ) \
 { \
 error( ( ERRCODE ) ); \
 std::cerr << "Unreconverable error during mesh reorder." << std::endl \
 << "Error Code " << ( ERRCODE ) << " at " << __FILE__ << ":" << __LINE__ << std::endl; \
 std::cerr.flush(); \
 abort(); \
 } \
 } while( false )
 
static ErrorCode check_tag_type( Interface* moab, Tag tag, DataType exp_type, int exp_size )
{
 ErrorCode rval;
 DataType act_type;
 int act_size;
 
 rval = moab->tag_get_data_type( tag, act_type );CHKERR;
 
 rval = moab->tag_get_bytes( tag, act_size );CHKERR;
 
 if( act_type != exp_type || act_size != exp_size ) return MB_TYPE_OUT_OF_RANGE;
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::handle_order_from_int_tag( Tag tag, int skip_value, Tag& new_handles )
{
 ErrorCode rval;
 
 // check that input tag handles are what we expect
 rval = check_tag_type( mMB, tag, MB_TYPE_INTEGER, sizeof( int ) );CHKERR;
 EntityHandle zero = 0;
 rval = mMB->tag_get_handle( 0, 1, MB_TYPE_HANDLE, new_handles, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_EXCL, &zero );CHKERR;
 
 // Can only reorder within same type/connectivity (or vertex/num_coords)
 // groupings. Call helper function for each such grouping.
 for( EntityType t = MBVERTEX; t < MBENTITYSET; ++t )
 {
 // Get list of all connectivity lengths (or vertex dimensions)
 // that exist for type t.
 TypeSequenceManager& seqs = mMB->sequence_manager()->entity_map( t );
 TypeSequenceManager::iterator i;
 std::set< int > values;
 for( i = seqs.begin(); i != seqs.end(); ++i )
 {
 EntitySequence* seq = *i;
 // 0 values per entity implies structured data, which
 // we cannot reorder. Skip those.
 if( t == MBVERTEX || 0 < seq->values_per_entity() ) values.insert( seq->values_per_entity() );
 }
 
 // Call helper function for each (type,size) tuple.
 std::set< int >::iterator j;
 for( j = values.begin(); j != values.end(); ++j )
 {
 rval = handle_order_from_int_tag( t, *j, tag, skip_value, new_handles );
 if( MB_SUCCESS != rval )
 {
 mMB->tag_delete( new_handles );
 return error( rval );
 }
 }
 }
 
 return MB_SUCCESS;
}
 
void ReorderTool::get_entities( EntityType t, int vals_per_ent, Range& entities )
{
 Range::iterator hint = entities.begin();
 ;
 TypeSequenceManager& seqs = mMB->sequence_manager()->entity_map( t );
 TypeSequenceManager::iterator s;
 for( s = seqs.begin(); s != seqs.end(); ++s )
 {
 EntitySequence* seq = *s;
 if( seq->values_per_entity() == vals_per_ent )
 hint = entities.insert( hint, seq->start_handle(), seq->end_handle() );
 }
}
 
ErrorCode ReorderTool::handle_order_from_int_tag( EntityType t,
 int vals_per_ent,
 Tag tag,
 int skip_value,
 Tag new_handles )
{
 ErrorCode rval;
 
 // check that input tag handles are what we expect
 rval = check_tag_type( mMB, tag, MB_TYPE_INTEGER, sizeof( int ) );CHKERR;
 rval = check_tag_type( mMB, new_handles, MB_TYPE_HANDLE, sizeof( EntityHandle ) );CHKERR;
 
 // get entities to re-order
 Range entities;
 get_entities( t, vals_per_ent, entities );
 
 // get entity order values
 std::vector< int > sortvals( entities.size() );
 rval = mMB->tag_get_data( tag, entities, &sortvals[0] );CHKERR;
 
 // remove any entities for which value is skip_value
 size_t r = 0, w = 0;
 for( Range::iterator i = entities.begin(); i != entities.end(); ++r )
 {
 if( sortvals[r] == skip_value )
 i = entities.erase( i );
 else
 {
 sortvals[w++] = sortvals[r];
 ++i;
 }
 }
 sortvals.resize( w );
 
 // sort
 std::sort( sortvals.begin(), sortvals.end() );
 // Convert to unique list and offsets for each value
 // When done, sortvals will contain unique, sortvals and
 // offsets will contain, for each unique value in sortvals,
 // the number of values that occured in the non-unique list
 // before the first instance of that value.
 std::vector< size_t > offsets;
 offsets.push_back( 0 );
 offsets.push_back( 1 );
 for( w = 0, r = 1; r < sortvals.size(); ++r )
 {
 if( sortvals[r] == sortvals[w] )
 {
 ++offsets.back();
 }
 else
 {
 ++w;
 sortvals[w] = sortvals[r];
 offsets.push_back( offsets.back() + 1 );
 }
 }
 ++w;
 assert( w + 1 == offsets.size() );
 sortvals.resize( w );
 
 // Tag each entity with its new handle
 for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
 {
 int val;
 rval = mMB->tag_get_data( tag, &*i, 1, &val );CHKERR;
 w = std::lower_bound( sortvals.begin(), sortvals.end(), val ) - sortvals.begin();
 assert( w < sortvals.size() );
 size_t offset = offsets[w];
 ++offsets[w];
 // should maybe copy range into array to avoid possibly n^2 behavior here
 EntityHandle h = *( entities.begin() + offset );
 rval = mMB->tag_set_data( new_handles, &*i, 1, &h );CHKERR;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::handle_order_from_sets_and_adj( const Range& sets, Tag& handle_tag )
{
 ErrorCode rval;
 
 if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
 
 Tag order_tag;
 const int negone = -1;
 rval = mMB->tag_get_handle( 0, 1, MB_TYPE_INTEGER, order_tag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_EXCL, &negone );
 if( MB_SUCCESS != rval )
 {
 mMB->tag_delete( handle_tag );
 handle_tag = 0;
 return error( rval );
 }
 
 std::vector< std::vector< EntityHandle >* > data;
 rval = int_order_from_sets_and_adj( sets, order_tag, negone, data );
 for( size_t i = 0; i < data.size(); ++i )
 delete data[i];
 if( MB_SUCCESS != rval )
 {
 mMB->tag_delete( order_tag );
 return error( rval );
 }
 
 rval = handle_order_from_int_tag( order_tag, negone, handle_tag );
 if( MB_SUCCESS != rval )
 {
 mMB->tag_delete( order_tag );
 return error( rval );
 }
 
 rval = mMB->tag_delete( order_tag );
 if( MB_SUCCESS != rval ) return error( rval );
 
 return MB_SUCCESS;
}
 
// Compare function to use for a map keyed on pointers to sorted vectors
struct CompSortedVect
{
 bool operator()( const std::vector< EntityHandle >* v1, const std::vector< EntityHandle >* v2 ) const
 {
 std::vector< EntityHandle >::const_iterator i1, i2;
 for( i1 = v1->begin(), i2 = v2->begin(); i1 != v1->end(); ++i1, ++i2 )
 {
 if( i2 == v2->end() || *i1 > *i2 )
 return false;
 else if( *i1 < *i2 )
 return true;
 }
 return i2 != v2->end();
 }
};
 
ErrorCode ReorderTool::int_order_from_sets_and_adj( const Range& sets,
 Tag order_tag,
 int skip_val,
 std::vector< std::vector< EntityHandle >* >& revMap )
{
 ErrorCode rval;
 
 if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
 
 rval = check_tag_type( mMB, order_tag, MB_TYPE_INTEGER, sizeof( int ) );CHKERR;
 
 // Compare function to use for a map keyed on pointers to sorted vectors
 CompSortedVect lessthan;
 // Map from sorted list of handles to index
 std::map< std::vector< EntityHandle >*, int, CompSortedVect > forMap;
 std::map< std::vector< EntityHandle >*, int, CompSortedVect >::iterator fiter, fiter2;
 std::vector< EntityHandle > sharing; // tmp storage for entity
 
 // for each set
 for( Range::iterator s = sets.begin(); s != sets.end(); ++s )
 {
 
 // gather up all entities and adjacencies
 Range tmp, ents, adj[4]; // indexed by dimension
 for( int dim = 0; dim < 4; ++dim )
 {
 rval = mMB->get_entities_by_dimension( *s, dim, tmp );CHKERR;
 for( int ldim = 0; ldim < dim; ++ldim )
 {
 rval = mMB->get_adjacencies( tmp, ldim, false, adj[ldim], Interface::UNION );CHKERR;
 }
 for( int udim = dim + 1; udim <= 3; ++udim )
 {
 rval = mMB->get_adjacencies( tmp, udim, false, adj[udim], Interface::UNION );CHKERR;
 }
 ents.merge( tmp );
 tmp.clear();
 }
 for( int dim = 0; dim < 4; ++dim )
 {
 ents.merge( adj[dim] );
 adj[dim].clear();
 }
 
 // process each entity
 for( Range::iterator e = ents.begin(); e != ents.end(); ++e )
 {
 int val;
 rval = mMB->tag_get_data( order_tag, &*e, 1, &val );CHKERR;
 
 // If this entity is already in one or more of the sets (either
 // directly or through adjacency) then get the existing list of
 // sets and append this set handle (we are iterating over sets
 // in sorted order, so appending should aways maintain a sorted
 // list.)
 sharing.clear();
 if( val != skip_val )
 {
 sharing = *revMap[val];
 assert( std::lower_bound( sharing.begin(), sharing.end(), *s ) == sharing.end() );
 }
 sharing.push_back( *s );
 
 // Check if new sharing list already exists in forward map
 fiter = forMap.lower_bound( &sharing );
 if( fiter == forMap.end() || lessthan( fiter->first, &sharing ) )
 {
 // Add new sharing list to forward and reverse maps.
 std::vector< EntityHandle >* newvec = new std::vector< EntityHandle >;
 newvec->swap( sharing );
 if( (int)revMap.size() == skip_val ) revMap.push_back( 0 );
 fiter2 =
 forMap.insert( fiter, std::pair< std::vector< EntityHandle >*, int >( newvec, revMap.size() ) );
 assert( fiter2 != fiter );
 fiter = fiter2;
 revMap.push_back( newvec );
 }
 
 // Update index on entity
 val = fiter->second;
 rval = mMB->tag_set_data( order_tag, &*e, 1, &val );CHKERR;
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::get_reordered_handles( Tag tag,
 const Range& old_handles,
 std::vector< EntityHandle >& new_handles )
{
 new_handles.resize( old_handles.size() );
 ErrorCode rval = mMB->tag_get_data( tag, old_handles, ( new_handles.empty() ) ? NULL : &new_handles[0] );CHKERR;
 
 Range::const_iterator it1 = old_handles.begin();
 std::vector< EntityHandle >::iterator it2 = new_handles.begin();
 for( ; it1 != old_handles.end(); ++it1, ++it2 )
 if( 0 == *it2 ) *it2 = *it1;
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::get_reordered_handles( Tag tag,
 const std::vector< EntityHandle >& old_handles,
 std::vector< EntityHandle >& new_handles )
{
 new_handles.resize( old_handles.size() );
 return get_reordered_handles( tag, &old_handles[0], &new_handles[0], old_handles.size() );
}
 
ErrorCode ReorderTool::get_reordered_handles( Tag tag,
 const EntityHandle* old_handles,
 EntityHandle* new_handles,
 size_t num_handles )
{
 ErrorCode rval = mMB->tag_get_data( tag, old_handles, num_handles, new_handles );CHKERR;
 
 for( size_t i = 0; i < num_handles; ++i )
 if( 0 == new_handles[i] ) new_handles[i] = old_handles[i];
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::get_new_handles( Tag tag, Range& old_handles, std::vector< EntityHandle >& newhandles )
{
 // get new handles for tagged entities
 newhandles.resize( old_handles.size() );
 ErrorCode rval = mMB->tag_get_data( tag, old_handles, ( newhandles.empty() ) ? NULL : &newhandles[0] );CHKERR;
 
 // remove entities that were not reordered
 Range::iterator i = old_handles.begin();
 size_t w = 0;
 for( size_t r = 0; r < newhandles.size(); ++r )
 {
 if( 0 != newhandles[r] )
 {
 newhandles[w] = newhandles[r];
 ++w;
 ++i;
 }
 else
 {
 i = old_handles.erase( i );
 }
 }
 newhandles.resize( w );
 assert( newhandles.size() == old_handles.size() );
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::reorder_entities( Tag new_handles )
{
 ErrorCode rval;
 
 rval = check_tag_type( mMB, new_handles, MB_TYPE_HANDLE, sizeof( EntityHandle ) );CHKERR;
 EntityHandle defval;
 rval = mMB->tag_get_default_value( new_handles, &defval );CHKERR;
 if( 0 != defval ) return error( MB_INDEX_OUT_OF_RANGE );
 
 // Can only reorder within same type/connectivity (or vertex/num_coords)
 // groupings.
 for( EntityType t = MBVERTEX; t < MBENTITYSET; ++t )
 {
 // Get list of all connectivity lengths (or vertex dimensions)
 // that exist for type t.
 TypeSequenceManager& seqs = mMB->sequence_manager()->entity_map( t );
 TypeSequenceManager::iterator i;
 std::set< int > values;
 for( i = seqs.begin(); i != seqs.end(); ++i )
 {
 EntitySequence* seq = *i;
 // 0 values per entity implies structured data, which
 // we cannot reorder. Skip those.
 if( t == MBVERTEX || 0 < seq->values_per_entity() ) values.insert( seq->values_per_entity() );
 }
 
 // reorder primary data for each (type,size) tuple.
 std::set< int >::iterator j;
 for( j = values.begin(); j != values.end(); ++j )
 {
 Range entities;
 get_entities( t, *j, entities );
 std::vector< EntityHandle > handles;
 rval = get_reordered_handles( new_handles, entities, handles );
 UNRECOVERABLE( rval );
 
 if( t == MBVERTEX )
 {
 std::vector< double > coords( entities.size() * 3 );
 rval = mMB->get_coords( entities, &coords[0] );
 UNRECOVERABLE( rval );
 rval = mMB->set_coords( &handles[0], handles.size(), &coords[0] );
 UNRECOVERABLE( rval );
 }
 else
 {
 std::vector< EntityHandle > conn;
 conn.reserve( entities.size() * *j );
 std::vector< EntityHandle > old_handles;
 old_handles.resize( entities.size() );
 std::copy( entities.begin(), entities.end(), old_handles.begin() );
 rval = mMB->get_connectivity( &old_handles[0], old_handles.size(), conn, false );
 UNRECOVERABLE( rval );
 old_handles.clear();
 old_handles = conn;
 rval = get_reordered_handles( new_handles, old_handles, conn );
 UNRECOVERABLE( rval );
 for( unsigned int h = 0; h < handles.size(); ++h )
 {
 rval = mMB->set_connectivity( handles[h], &conn[h * *j], *j );
 UNRECOVERABLE( rval );
 }
 }
 }
 }
 
 // now update tag data
 std::vector< Tag > tag_handles;
 mMB->tag_get_tags( tag_handles );
 for( size_t i = 0; i < tag_handles.size(); ++i )
 {
 Tag tag = tag_handles[i];
 if( tag == new_handles ) // don't mess up mapping from old to new handles
 continue;
 
 for( EntityType t = MBVERTEX; t <= MBENTITYSET; ++t )
 {
 rval = reorder_tag_data( t, new_handles, tag );
 UNRECOVERABLE( rval );
 }
 }
 
 rval = update_set_contents( new_handles );
 UNRECOVERABLE( rval );
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::reorder_tag_data( EntityType etype, Tag new_handles, Tag tag )
{
 ErrorCode rval;
 
 int tagsize;
 DataType tagtype;
 rval = mMB->tag_get_data_type( tag, tagtype );
 if( MB_SUCCESS != rval ) return error( rval );
 if( MB_TYPE_BIT == tagtype )
 tagsize = 1;
 else
 {
 rval = mMB->tag_get_bytes( tag, tagsize );
 if( MB_VARIABLE_DATA_LENGTH == rval )
 tagsize = -1;
 else if( MB_SUCCESS != rval )
 return error( rval );
 }
 
 // we don't re-order set handles, so we don't care about sets
 // unless the tag contains handles that need to be updated
 if( MBENTITYSET == etype && MB_TYPE_HANDLE != tagtype ) return MB_SUCCESS;
 
 // get tagged entities
 Range old_tagged;
 rval = mMB->get_entities_by_type_and_tag( 0, etype, &tag, 0, 1, old_tagged );
 if( MB_SUCCESS != rval && old_tagged.empty() ) return error( rval );
 
 // remove entities that were not reordered, unless the tag
 // is handle type in which case we need to update the data
 // for all entities, regardless of reordering.
 std::vector< EntityHandle > newhandles;
 if( MB_TYPE_HANDLE == tagtype )
 rval = get_reordered_handles( new_handles, old_tagged, newhandles );
 else
 rval = get_new_handles( new_handles, old_tagged, newhandles );CHKERR;
 
 if( old_tagged.empty() ) return MB_SUCCESS;
 
 // get copy of all tag data
 std::vector< unsigned char > buffer;
 std::vector< const void* > pointers;
 std::vector< int > sizes;
 // if variable-length tag
 if( -1 == tagsize )
 {
 pointers.resize( old_tagged.size() );
 sizes.resize( pointers.size() );
 rval = mMB->tag_get_by_ptr( tag, old_tagged, &pointers[0], &sizes[0] );CHKERR;
 int total = std::accumulate( sizes.begin(), sizes.end(), 0 );
 DataType dtype;
 mMB->tag_get_data_type( tag, dtype );
 int type_size;
 switch( dtype )
 {
 case MB_TYPE_INTEGER:
 type_size = sizeof( int );
 break;
 case MB_TYPE_DOUBLE:
 type_size = sizeof( double );
 break;
 case MB_TYPE_HANDLE:
 type_size = sizeof( EntityHandle );
 break;
 case MB_TYPE_BIT:
 type_size = 1;
 break;
 case MB_TYPE_OPAQUE:
 type_size = 1;
 break;
 default:
 return MB_TYPE_OUT_OF_RANGE;
 }
 buffer.resize( total * type_size );
 size_t off = 0;
 for( size_t j = 0; j < pointers.size(); ++j )
 {
 memcpy( &buffer[off], pointers[j], type_size * sizes[j] );
 pointers[j] = &buffer[off];
 off += sizes[j] * type_size;
 }
 }
 // if fixed-length tag
 else
 {
 buffer.resize( old_tagged.size() * tagsize );
 rval = mMB->tag_get_data( tag, old_tagged, &buffer[0] );CHKERR;
 }
 
 // if handle tag, update tag values for reordered handles
 if( MB_TYPE_HANDLE == tagtype )
 {
 assert( !( buffer.size() % sizeof( EntityHandle ) ) );
 std::vector< unsigned char > buffer2( buffer.size() );
 rval = get_reordered_handles( new_handles, reinterpret_cast< const EntityHandle* >( &buffer[0] ),
 reinterpret_cast< EntityHandle* >( &buffer2[0] ),
 buffer.size() / sizeof( EntityHandle ) );CHKERR;
 // if var-length tag then do not do swap because pointers[] contains pointers
 // into old buffer
 if( -1 == tagsize )
 memcpy( &buffer[0], &buffer2[0], buffer.size() );
 else
 buffer.swap( buffer2 );
 }
 
 // store re-ordered tag data
 if( -1 == tagsize )
 {
 rval = mMB->tag_set_by_ptr( tag, &newhandles[0], newhandles.size(), &pointers[0], &sizes[0] );
 pointers.clear();
 sizes.clear();
 buffer.clear();
 }
 else
 {
 rval = mMB->tag_set_data( tag, &newhandles[0], newhandles.size(), &buffer[0] );
 buffer.clear();
 }
 CHKERR;
 
 // all permutations should be cyclical, but not all permuted
 // entities necessarily had tag values, so we might need to delete
 // tags for some entities
 std::sort( newhandles.begin(), newhandles.end() );
 std::vector< EntityHandle >::iterator k = newhandles.begin();
 Range::iterator i = old_tagged.begin();
 while( i != old_tagged.end() )
 {
 while( k != newhandles.end() && *k < *i )
 ++k;
 if( k == newhandles.end() ) break;
 
 if( *i == *k ) // what old_tagged -= newhandles
 i = old_tagged.erase( i );
 else
 ++i;
 }
 
 if( !old_tagged.empty() )
 {
 rval = mMB->tag_delete_data( tag, old_tagged );CHKERR;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode ReorderTool::update_set_contents( Tag nh_tag )
{
 Range sets;
 ErrorCode rval = mMB->get_entities_by_type( 0, MBENTITYSET, sets );CHKERR;
 
 std::vector< EntityHandle > old_handles, new_handles;
 for( Range::iterator i = sets.begin(); i != sets.end(); ++i )
 {
 // If set is un-ordered...
 unsigned opts = 0;
 mMB->get_meshset_options( *i, opts );
 if( !( opts & MESHSET_ORDERED ) )
 {
 Range contents;
 rval = mMB->get_entities_by_handle( *i, contents );CHKERR;
 
 rval = get_new_handles( nh_tag, contents, new_handles );CHKERR;
 
 Range replace;
 std::sort( new_handles.begin(), new_handles.end() );
 Range::iterator hint = replace.begin();
 for( size_t j = 0; j < new_handles.size(); ++j )
 hint = replace.insert( hint, new_handles[j] );
 Range common = intersect( contents, replace );
 contents -= common;
 replace -= common;
 assert( contents.size() == replace.size() );
 if( !contents.empty() )
 {
 rval = mMB->remove_entities( *i, contents );CHKERR;
 rval = mMB->add_entities( *i, replace );
 }
 }
 
 // If set is ordered...
 else
 {
 // get set contents
 old_handles.clear();
 rval = mMB->get_entities_by_handle( *i, old_handles );CHKERR;
 
 // get new handles from old contained handles
 rval = get_reordered_handles( nh_tag, old_handles, new_handles );CHKERR;
 
 rval = mMB->clear_meshset( &*i, 1 );CHKERR;
 
 rval = mMB->add_entities( *i, &new_handles[0], new_handles.size() );CHKERR;
 }
 } // for each set
 
 return MB_SUCCESS;
}
 
} // namespace moab