fathomteam-moab-7bde9dfb84a8/src/ReorderTool.cpp

/** \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