gen.cpp

Go to the documentation of this file.

#include <iostream>
#include <iomanip> // for setprecision
#include <limits> // for min/max values
#include <assert.h>
#include <math.h>
#include <time.h>
#include <vector>

#include "moab/GeomTopoTool.hpp"
#include "moab/FileOptions.hpp"
#include "moab/Core.hpp"

#include "meshkit/gen.hpp"
#include "meshkit/zip.hpp"
#include "moab/Skinner.hpp"


const char GEOM_SENSE_2_TAG_NAME[] = "GEOM_SENSE_2";
const char GEOM_SENSE_N_ENTS_TAG_NAME[] = "GEOM_SENSE_N_ENTS";
const char GEOM_SENSE_N_SENSES_TAG_NAME[] = "GEOM_SENSE_N_SENSES"; 

using namespace moab;
namespace gen {

  bool error( const bool error_has_occured, const std::string message ) {
    if(error_has_occured) {
      if("" == message) {
        std::cout << "Error at " << __FILE__ << ":" << __LINE__ 
                  << std::endl;
      } else {
        std::cout << message << std::endl;
      }
      return true;
    } else {
      return false;
    }
  }

void moab_printer(ErrorCode error_code)
{
  if ( error_code == MB_INDEX_OUT_OF_RANGE )
    {
      std::cerr << "ERROR: MB_INDEX_OUT_OF_RANGE" << std::endl;
    }
  if ( error_code == MB_MEMORY_ALLOCATION_FAILED )
    {
      std::cerr << "ERROR: MB_MEMORY_ALLOCATION_FAILED" << std::endl;
    }
  if ( error_code == MB_ENTITY_NOT_FOUND )
    {
      std::cerr << "ERROR: MB_ENTITY_NOT_FOUND" << std::endl;
    }
  if ( error_code == MB_MULTIPLE_ENTITIES_FOUND )
    {
      std::cerr << "ERROR: MB_MULTIPLE_ENTITIES_FOUND" << std::endl;
    }
  if ( error_code == MB_TAG_NOT_FOUND )
    {
      std::cerr << "ERROR: MB_TAG_NOT_FOUND" << std::endl;
    }
  if ( error_code == MB_FILE_DOES_NOT_EXIST )
    {
      std::cerr << "ERROR: MB_FILE_DOES_NOT_EXIST" << std::endl;
    }    
  if ( error_code == MB_FILE_WRITE_ERROR )
    {
      std::cerr << "ERROR: MB_FILE_WRITE_ERROR" << std::endl;
    }    
  if ( error_code == MB_ALREADY_ALLOCATED )
    {
      std::cerr << "ERROR: MB_ALREADY_ALLOCATED" << std::endl;
    }    
  if ( error_code == MB_VARIABLE_DATA_LENGTH )
    {
      std::cerr << "ERROR: MB_VARIABLE_DATA_LENGTH" << std::endl;
    }  
  if ( error_code == MB_INVALID_SIZE )
    {
      std::cerr << "ERROR: MB_INVALID_SIZE" << std::endl;
    }  
  if ( error_code == MB_UNSUPPORTED_OPERATION )
    {
      std::cerr << "ERROR: MB_UNSUPPORTED_OPERATION" << std::endl;
    }  
  if ( error_code == MB_UNHANDLED_OPTION )
    {
      std::cerr << "ERROR: MB_UNHANDLED_OPTION" << std::endl;
    }  
  if ( error_code == MB_FAILURE )
    {
      std::cerr << "ERROR: MB_FAILURE" << std::endl;
    }  
  return;
}





  void print_vertex_cubit( const EntityHandle vertex ) {

    ErrorCode result;
    double coords[3];
    int n_precision = 20;
    result = MBI()->get_coords( &vertex, 1, coords );
    if(gen::error(MB_SUCCESS!=result, "failed to get vertex coords"));
    assert(MB_SUCCESS == result);
    std::cout << "  create vertex " 
              << std::setprecision(n_precision)
              << coords[0] << " " << coords[1] << " " << coords[2]
              << std::endl;
  }

  void print_vertex_coords( const EntityHandle vertex ) {

    ErrorCode result;
    double coords[3];
    result = MBI()->get_coords( &vertex, 1, coords );
    if(MB_SUCCESS!=result) std::cout << "vert=" << vertex << std::endl;
    assert(MB_SUCCESS == result);
    std::cout << "    vertex " << vertex << " coords= (" 
              << coords[0] << "," << coords[1] << "," << coords[2] << ")" 
              << std::endl;
  }

  void print_triangles( const Range tris ) {
    for(Range::const_iterator i=tris.begin(); i!=tris.end(); i++) {
      print_triangle( *i, false );
    }
  }
  // If the edges of the tri are ambiguous, do not print edges!
  void print_triangle( const EntityHandle tri, bool print_edges ) {
    ErrorCode result;
    double area;
    result = triangle_area( tri, area );
    assert(MB_SUCCESS == result);
    std::cout << "    triangle " << tri << " area=" << area << std::endl;
    const EntityHandle *conn;
    int n_verts;
    result = MBI()->get_connectivity( tri, conn, n_verts );
    assert(MB_SUCCESS == result);
    assert(3 == n_verts);
    for(int i=0; i<3; i++) print_vertex_coords( conn[i] );
    
    if(print_edges) {
      Range edges;
      result = MBI()->get_adjacencies( &tri, 1, 1, true, edges );
      if(MB_SUCCESS != result) std::cout << "result=" << result << std::endl;
      assert(MB_SUCCESS == result);
      //std::cout << "      edges: ";
      for(Range::iterator i=edges.begin(); i!=edges.end(); i++) {
        //std::cout << *i << " ";
        print_edge( *i );
      }
      //std::cout << std::endl;
    }
  }

  void print_edge( const EntityHandle edge ) {
    const EntityHandle *conn;
    int n_verts;
    std::cout << "    edge " << edge << std::endl;
    ErrorCode result = MBI()->get_connectivity( edge, conn, n_verts );
    if(gen::error(MB_SUCCESS!=result, "failed to get edge connectivity"));
    assert(MB_SUCCESS == result);
    assert(2 == n_verts);
    print_vertex_coords( conn[0] );   
    print_vertex_coords( conn[1] ); 

    Range tris;
    result = MBI()->get_adjacencies( &edge, 1, 2, false, tris );
    assert(MB_SUCCESS == result);
    std::cout << "     tris: ";
    for(Range::iterator i=tris.begin(); i!=tris.end(); i++) {
      std::cout << *i << " ";
    }
    std::cout << std::endl; 
  }

  void print_range( const Range range ) {
    std::cout << "print range:" << std::endl;
    Range::iterator i;
    for(i=range.begin(); i!=range.end(); i++) {
      std::cout << "    " << *i << std::endl;
    }
  }

  void print_range_of_edges( const Range range ) {
    std::cout << "print range:" << std::endl;
    Range::const_iterator i;
    for(i=range.begin(); i!=range.end(); i++) {
      print_edge( *i );
    }
  }


  void print_vertex_count(const EntityHandle input_meshset) {
  
    // get the range of facets of the surface meshset
    ErrorCode result;
    Range vertices;
    result = MBI()->get_entities_by_type(0, MBVERTEX, vertices);
    if(gen::error(MB_SUCCESS!=result, "failed to get vertex entities from the mesh"));
    assert( MB_SUCCESS == result );
  
    std::cout<< "    " << vertices.size() << " vertices found." << std::endl;
  }

  void print_arcs( const std::vector< std::vector<EntityHandle> > arcs ) {
    for(unsigned int i=0; i<arcs.size(); i++) {
      std::cout << "arc " << i << std::endl;
      print_loop( arcs[i] );
    }
  }

  void print_arc_of_edges( const std::vector<EntityHandle> arc_of_edges ) {
  
    ErrorCode result;
    std::vector<EntityHandle>::const_iterator i;
    double dist = 0;
    for( i=arc_of_edges.begin(); i!=arc_of_edges.end(); i++ ) {
      int n_verts;
      const EntityHandle *conn;
      result = MBI()->get_connectivity( *i, conn, n_verts ); 
      if(gen::error(MB_SUCCESS!=result, "failed to get edge connectivity"));
      assert(MB_SUCCESS == result);
      assert( 2 == n_verts );
      dist += dist_between_verts( conn[0], conn[1] );
      print_vertex_coords( conn[0] );
      print_vertex_coords( conn[1] );
    }
    std::cout << "  dist= " << dist << std::endl;
  }

  void print_loop( const std::vector<EntityHandle> loop_of_verts ) {
   
    std::cout << "  size=" << loop_of_verts.size() << std::endl;
    double dist = 0;
    //std::vector<EntityHandle>::iterator i;
    //for( i=loop_of_verts.begin(); i!=loop_of_verts.end(); i++ ) {
    for(unsigned int i=0; i<loop_of_verts.size(); i++) {
      print_vertex_coords( loop_of_verts[i] );
      if(i != loop_of_verts.size()-1) {
        dist += dist_between_verts( loop_of_verts[i], loop_of_verts[i+1] );
      }
    }
    std::cout << "  dist=" << dist << std::endl;
  }


ErrorCode find_closest_vert( const EntityHandle reference_vert,
                               const std::vector<EntityHandle> arc_of_verts,
                               unsigned &position,
                               const double dist_limit ) {
  ErrorCode rval;
  const bool debug = false;
  double min_dist_sqr = std::numeric_limits<double>::max();
  CartVect ref_coords;
  rval = MBI()->get_coords( &reference_vert, 1, ref_coords.array() );
  if(gen::error(MB_SUCCESS!=rval,"failed to get ref coords")) return rval;
  double length = 0;
  CartVect prev_coords;

  for(unsigned i=0; i<arc_of_verts.size(); ++i) {
    CartVect coords;
    rval = MBI()->get_coords( &arc_of_verts[i], 1, coords.array() );
    if(gen::error(MB_SUCCESS!=rval,"failed to get coords")) return rval;

    // use dist_limit to exit early; avoid checking the entire arc
    if(0!=i) {
      CartVect temp = prev_coords - coords;
      length += temp.length();
      if(length>dist_limit && debug) 
        std::cout << "length=" << length << " dist_limit=" << dist_limit << std::endl;
      if(length > dist_limit) return MB_SUCCESS;
    }
    prev_coords = coords;

    // get distance to ref_vert
    CartVect temp = ref_coords - coords;
    double dist_sqr = temp.length_squared();
    if(dist_sqr < min_dist_sqr) {
      position = i;
      min_dist_sqr = dist_sqr;
      if(debug) std::cout << "min_dist_sqr=" << min_dist_sqr << std::endl;
    }
  }
 
  return MB_SUCCESS;
}



  // Return the closest vert and all within tol. This is needed because sometimes
  // the correct vert is not the closest. For example, iter_surf4010 the skin
  // loop has the same point in it twice, at two different locations (center of L).
  // This ensure that both are returned as candidates.
  ErrorCode find_closest_vert( const double tol,
                                 const EntityHandle reference_vert,
                                 const std::vector<EntityHandle> loop_of_verts,
                                 std::vector<unsigned> &positions, 
                                 std::vector<double> &dists) {

    ErrorCode rval;
    positions.clear();
    dists.clear();
    const double TOL_SQR = tol*tol;
    unsigned min_pos;
    double sqr_min_dist = std::numeric_limits<double>::max();
    for(unsigned int i=0; i<loop_of_verts.size(); i++) {
        double sqr_dist = std::numeric_limits<double>::max();
        rval = squared_dist_between_verts(reference_vert, loop_of_verts[i], sqr_dist);
        if(gen::error(MB_SUCCESS!=rval,"could not get dist")) return rval;
        if(sqr_dist < sqr_min_dist) {
          if(sqr_dist >= TOL_SQR) {
            sqr_min_dist = sqr_dist;
            min_pos = i;
          } else {
            sqr_min_dist = TOL_SQR;
            positions.push_back(i);
            dists.push_back( sqrt(sqr_dist) );
          }
        }
    }
 
    if(dists.empty()) {
      dists.push_back( sqrt(sqr_min_dist) );
      positions.push_back( min_pos );
    }

            
   return MB_SUCCESS;
  }

  ErrorCode merge_vertices( Range verts /* in */, const double tol /* in */ ) {

    ErrorCode result;
    const double SQR_TOL = tol*tol;
    // Clean up the created tree, and track verts so that if merged away they are
    // removed from the tree.
    AdaptiveKDTree kdtree(MBI()); //, true, 0, MESHSET_TRACK_OWNER);
    // initialize the KD Tree
    EntityHandle root;
    const char settings[]="MAX_PER_LEAF=6;MAX_DEPTH=50;SPLITS_PER_DIR=1;PLANE_SET=2;MESHSET_FLAGS=0x1;TAG_NAME=0";
    FileOptions fileopts(settings);


    /* Old KDTree settings
    AdaptiveKDTree::Settings settings;
   
    // tells the tree to split leaves with more than 6 entities
    settings.maxEntPerLeaf = 6;
    // tells the tree a maximum depth limit
    settings.maxTreeDepth  = 50;
    // tells the tree how many candidate split planed to consider in each dimension
    settings.candidateSplitsPerDir = 1;
    // tells the tree to use the median vertex coordinate values to set planes
    settings.candidatePlaneSet = AdaptiveKDTree::VERTEX_MEDIAN;
   
    */

    // builds the KD Tree, making the EntityHandle root the root of the tree
    result = kdtree.build_tree( verts, &root, &fileopts);
    assert(MB_SUCCESS == result);
    // create tree iterator to loop over all verts in the tree
    AdaptiveKDTreeIter tree_iter;
    kdtree.get_tree_iterator( root, tree_iter );
 
    //for(unsigned int i=0; i<verts.size(); i++) {
    for(Range::iterator i=verts.begin(); i!=verts.end(); ++i) {
      double from_point[3];
      //EntityHandle vert = *i;
      result = MBI()->get_coords( &(*i), 1, from_point);
      assert(MB_SUCCESS == result);
      std::vector<EntityHandle> leaves_out;
      result = kdtree.distance_search( from_point, tol, leaves_out, root);
      assert(MB_SUCCESS == result);
      for(unsigned int j=0; j<leaves_out.size(); j++) {
        std::vector<EntityHandle> leaf_verts;
        result = MBI()->get_entities_by_type( leaves_out[j], MBVERTEX, leaf_verts);
        assert(MB_SUCCESS == result);
        if(100 < leaf_verts.size()) std::cout << "*i=" << *i << " leaf_verts.size()=" << leaf_verts.size() << std::endl;
        for(unsigned int k=0; k<leaf_verts.size(); k++) {
          if( leaf_verts[k] == *i ) continue; 
          double sqr_dist;
          result = gen::squared_dist_between_verts( *i, leaf_verts[k], sqr_dist);
          assert(MB_SUCCESS == result);

          if(SQR_TOL >= sqr_dist) {
            //  std::cout << "merge_vertices: vert " << leaf_verts[k] << " merged to vert "
            //            << verts[i] << " dist=" << dist << " leaf=" << std::endl;

            // The delete_vert is automatically remove from the tree because it
            // uses tracking meshsets. merge_verts checks for degenerate tris.
            // Update the list of leaf verts to prevent stale handles.
            std::vector<EntityHandle> temp_arc;
            EntityHandle keep_vert   = *i;
            EntityHandle delete_vert = leaf_verts[k];
            result = zip::merge_verts( keep_vert, delete_vert, leaf_verts, temp_arc );
            assert(MB_SUCCESS == result);
            // Erase delete_vert from verts
            // Iterator should remain valid because delete_vert > keep_vert handle.
            verts.erase( delete_vert );
          }
        }
      }
    }   
    return result;
  }

  ErrorCode squared_dist_between_verts( const EntityHandle v0,
                                          const EntityHandle v1,
                                          double &d) {
    ErrorCode result;
    CartVect coords0, coords1;
    result = MBI()->get_coords( &v0, 1, coords0.array() );
    if(MB_SUCCESS != result) {
      std::cout << "dist_between_verts: get_coords on v0=" << v0 << " result=" 
                << result << std::endl; 
      return result;
    }
    result = MBI()->get_coords( &v1, 1, coords1.array() );
    if(MB_SUCCESS != result) {
      std::cout << "dist_between_verts: get_coords on v1=" << v1 << " result=" 
                << result << std::endl; 
      return result;
    }
    const CartVect diff = coords0 - coords1;
    d = diff.length_squared();
    return MB_SUCCESS;
  }

  double dist_between_verts( const CartVect v0, const CartVect v1 ) {
    CartVect v2 = v0 - v1;
    return v2.length();
  }
  ErrorCode dist_between_verts( const EntityHandle v0, const EntityHandle v1, double &d) {
    ErrorCode result;
    CartVect coords0, coords1;
    result = MBI()->get_coords( &v0, 1, coords0.array() );
    if(MB_SUCCESS != result) {
      std::cout << "dist_between_verts: get_coords on v0=" << v0 << " result=" 
                << result << std::endl; 
      return result;
    }
    result = MBI()->get_coords( &v1, 1, coords1.array() );
    if(MB_SUCCESS != result) {
      std::cout << "dist_between_verts: get_coords on v1=" << v1 << " result=" 
                << result << std::endl; 
      return result;
    }
    d = dist_between_verts( coords0, coords1 );
    return MB_SUCCESS;
  }

  double dist_between_verts( double coords0[], double coords1[] ) {
    return sqrt( (coords0[0]-coords1[0])*(coords0[0]-coords1[0]) +
                 (coords0[1]-coords1[1])*(coords0[1]-coords1[1]) +
                 (coords0[2]-coords1[2])*(coords0[2]-coords1[2]) );
  }
  double dist_between_verts( EntityHandle vert0, EntityHandle vert1 ) {
    double coords0[3], coords1[3];
    ErrorCode result;
    result = MBI()->get_coords( &vert0, 1, coords0 );
    if(MB_SUCCESS!=result) std::cout << "result=" << result << " vert="
                                     << vert0 << std::endl;
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &vert1, 1, coords1 );
    if(MB_SUCCESS!=result) std::cout << "result=" << result << " vert="
                                     << vert1 << std::endl;
    assert(MB_SUCCESS == result);
    return dist_between_verts( coords0, coords1 );
  }

  // Return the length of the curve defined by MBEDGEs or ordered MBVERTEXs.
  double length( std::vector<EntityHandle> edges ) {
    if(edges.empty()) return 0;

    ErrorCode result;
    std::vector<EntityHandle>::iterator i;
    double dist = 0;
    EntityType type = MBI()->type_from_handle( edges[0] );

    // if vector has both edges and verts, only use edges
    // NOTE: The curve sets from ReadCGM do not contain duplicate endpoints for loops!
    EntityType end_type = MBI()->type_from_handle( edges.back() );
    if(type != end_type) {
      for(std::vector<EntityHandle>::iterator i=edges.begin(); i!=edges.end(); i++) {
        if(MBVERTEX == MBI()->type_from_handle( *i )) {
          i = edges.erase(i) - 1;
        }
      }
    }

    // determine if vector defines an arc by edges of verts
    type = MBI()->type_from_handle( edges[0] ); 
    if        (MBEDGE == type) {
      if(edges.empty()) return 0.0; 
      for( i=edges.begin(); i!=edges.end(); i++ ) {
        int n_verts;
        const EntityHandle *conn;
        result = MBI()->get_connectivity( *i, conn, n_verts );
        if( MB_SUCCESS!=result ) std::cout << "result=" << result << std::endl;
        assert(MB_SUCCESS == result);
        assert( 2 == n_verts );
        if(conn[0] == conn[1]) continue;
        dist += dist_between_verts( conn[0], conn[1] );
        //std::cout << "length: " << dist << std::endl;
      }
    } else if (MBVERTEX == type) {
      if(2 > edges.size()) return 0.0;
      EntityHandle front_vert = edges.front();
      for( i=edges.begin()+1; i!=edges.end(); i++) {
        dist += dist_between_verts( front_vert, *i );
        front_vert = *i;
      }
    } else return MB_FAILURE;

    return dist;
  }

  // Given a vertex and vector of edges, return the number of edges adjacent to the vertex.
  unsigned int n_adj_edges( EntityHandle vert, Range edges ) {
    ErrorCode result;
    Range adj_edges;
    result = MBI()->get_adjacencies( &vert, 1, 1, false, adj_edges );
    gen::error(MB_SUCCESS!=result, "could not get edges adjacent to the vertex");
    assert(MB_SUCCESS == result);
    //adj_edges = adj_edges.intersect(edges);
    adj_edges = intersect( adj_edges, edges );
    return adj_edges.size();
  }



  // Return true if the edges share a vertex. Does not check for coincident edges.
  bool edges_adjacent( EntityHandle edge0, EntityHandle edge1 ) {
    ErrorCode result;
    Range verts0, verts1;
    result = MBI()->get_adjacencies( &edge0, 1, 0, false, verts0 );
    gen::error(MB_SUCCESS!=result, "could not get edge0 adjacencies");
    assert( MB_SUCCESS == result );
    assert( 2 == verts0.size() );
    result = MBI()->get_adjacencies( &edge1, 1, 0, false, verts1 );
    gen::error(MB_SUCCESS!=result, "could not get edge1 adjacencies");
    assert( MB_SUCCESS == result );
    assert( 2 == verts1.size() );
    if      ( verts0.front() == verts1.front() ) return true;
    else if ( verts0.front() == verts1.back()  ) return true;
    else if ( verts0.back()  == verts1.back()  ) return true;
    else if ( verts0.back()  == verts1.front() ) return true;
    else                                         return false;
  }

  // get the direction unit vector from one vertex to another vertex
  ErrorCode get_direction( const EntityHandle from_vert, const EntityHandle to_vert,
                             CartVect &dir ) {
    // double d[3];
    ErrorCode result;
    CartVect coords0, coords1;
    result = MBI()->get_coords( &from_vert, 1, coords0.array() );
    assert(MB_SUCCESS==result);
    result = MBI()->get_coords( &to_vert, 1, coords1.array() );
    assert(MB_SUCCESS==result);
    dir = coords1 - coords0;
    if(0 == dir.length()) {
      CartVect zero_vector( 0.0 );
      dir = zero_vector;
      std::cout << "direction vector has 0 magnitude" << std::endl;
      return MB_SUCCESS;
    }
    dir.normalize();
    return result;
  }    
 
  // from http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry1
  double edge_point_dist( const CartVect a, const CartVect b, const CartVect c ) {
    CartVect ab, bc, ba, ac;
    ab = b - a;
    bc = c - b;
    ba = a - b;
    ac = c - a;
  
    // find the magnitude of the cross product and test the line
    CartVect cross_product = ab*ac;
    double dist = cross_product.length() / dist_between_verts(a,b);
 
    // test endpoint1
    if (ab%bc > 0) {
      //std::cout << "edge_point_dist=" << dist_between_verts(b,c) 
      //<< " at endpt1" << std::endl;
      return dist_between_verts(b,c);
    }

    // test endpoint0
    if (ba%ac > 0) {
      //std::cout << "edge_point_dist=" << dist_between_verts(a,c) 
      //<< " at endpt0" << std::endl;
      return dist_between_verts(a,c);
    }
 
    //std::cout << "edge_point_dist=" << fabs(dist) << " at middle" 
    //<< std::endl;      
    return fabs(dist);
  }
  double edge_point_dist( const EntityHandle endpt0, const EntityHandle endpt1,
                          const EntityHandle pt ) {
    ErrorCode result;
    CartVect a, b, c;
    result = MBI()->get_coords( &endpt0, 1, a.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS==result);
    result = MBI()->get_coords( &endpt1, 1, b.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS==result);
    result = MBI()->get_coords( &pt,     1, c.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS==result);
    return edge_point_dist( a, b, c);
  }
  double edge_point_dist( const EntityHandle edge, const EntityHandle pt ) {
    ErrorCode result;
    const EntityHandle *conn;
    int n_verts;     
    result = MBI()->get_connectivity( edge, conn, n_verts );
    gen::error(MB_SUCCESS!=result, "could not get edge connectivity");
    assert(MB_SUCCESS==result);
    assert( 2 == n_verts );
    return edge_point_dist( conn[0], conn[1], pt );
  }
  /*
  ErrorCode  point_curve_min_dist( const std::vector<EntityHandle> curve, // of verts
                                     const EntityHandle pt,
                                     double &min_dist,
                                     const double max_dist_along_curve ) { 
  
    min_dist = std::numeric_limits<double>::max();
    double cumulative_dist = 0;
    bool last_edge = false;
  
    // it is a curve of verts or a curve of edges?
    EntityType type = MBI()->type_from_handle( curve.front() );
    std::vector<EntityHandle>::const_iterator i;
    if(MBVERTEX == type) {
      EntityHandle front_vert = curve.front();
      for( i=curve.begin()+1; i!=curve.end(); i++) {
        // if we are using verts, do not explicitly create the edge in MOAB
        cumulative_dist += gen::dist_between_verts( front_vert, *i );
        //std::cout << "  point_curve_min_dist: cumulative_dist="
        //        << cumulative_dist << std::endl;
        double d = edge_point_dist( front_vert, *i, pt );
        //std::cout << "  point_curve_min_dist: d=" << d << std::endl;        
        if(d < min_dist) {
          min_dist = d;
          //std::cout << "min_dist=" << min_dist << std::endl;
          //print_vertex_coords( front_vert );
          //print_vertex_coords( *i );
        }
        // check one edge past the point after max_dist_along_curve
        if(last_edge) return MB_SUCCESS;
        if(max_dist_along_curve<cumulative_dist) last_edge = true;
      
        front_vert = *i;
      }
    } //else if(MBEDGE == type) {
      //for( i=curve.begin(); i!=curve.end(); i++) {
//      double d = edge_point_dist( *i, pt );
        //if(d < min_dist) min_dist = d;
      //}
     // } 
       else return MB_FAILURE;
    //std::cout << "point_curve_min_dist=" << min_dist << " curve.size()=" << curve.size() << std::endl;    
    return MB_SUCCESS;
  }

  ErrorCode  point_curve_min_dist( const std::vector<EntityHandle> curve, // of verts
                                     const EntityHandle pt,
                                     double &min_dist ) { 
    const double max_dist_along_curve = std::numeric_limits<double>::max();
    return point_curve_min_dist( curve, pt, min_dist, max_dist_along_curve ); 
  }
*/
  double triangle_area( const CartVect a, const CartVect b,
                        const CartVect c) {
    CartVect d = c - a;
    CartVect e = c - b;
    CartVect f = d*e;
    return 0.5*f.length();
  }
  ErrorCode triangle_area( const EntityHandle conn[], double &area ) {
    CartVect coords[3];
    ErrorCode result = MBI()->get_coords( conn, 3, coords[0].array() );
    gen::error(MB_SUCCESS!=result, "could not get triangle vertex coordinates");
    assert(MB_SUCCESS == result);
    area = triangle_area( coords[0], coords[1], coords[2] );
    return result;
  }
  ErrorCode triangle_area( const EntityHandle tri, double &area ) {
    ErrorCode result;
    const EntityHandle *conn;
    int n_verts;
    result = MBI()->get_connectivity( tri, conn, n_verts );
    gen::error(MB_SUCCESS!=result, "could not get trangle vertices");
    assert(MB_SUCCESS == result);
    assert(3 == n_verts);
    
    result = triangle_area( conn, area );
    gen::error(MB_SUCCESS!=result, "could not get triangle area");
    assert(MB_SUCCESS == result);
    return result;
  }
  double triangle_area( const Range tris ) {
    double a, area = 0;
    ErrorCode result;
    for(Range::iterator i=tris.begin(); i!=tris.end(); i++) {
      result = triangle_area( *i, a);
      gen::error(MB_SUCCESS!=result, "could not get triangle area");
      assert(MB_SUCCESS == result);
      area += a;
    }
    return area;
  }

  bool triangle_degenerate( const EntityHandle tri ) {
    ErrorCode result;
    const EntityHandle *conn;
    int n_verts;
    result = MBI()->get_connectivity( tri, conn, n_verts );
    gen::error(MB_SUCCESS!=result, "could not get triangle vertices");
    assert(MB_SUCCESS == result);
    assert(3 == n_verts);
    return triangle_degenerate( conn[0], conn[1], conn[2] );
  }

  bool triangle_degenerate( const EntityHandle v0, const EntityHandle v1,
                            const EntityHandle v2 ) {
    if(v0==v1 || v1==v2 || v2==v0) return true;
    return false;
  }

  ErrorCode triangle_normals( const Range tris, std::vector<CartVect> &normals ) {
    ErrorCode result;
    normals.clear();
    for(Range::const_iterator i=tris.begin(); i!=tris.end(); i++) {
      CartVect normal;
      result = triangle_normal( *i, normal );
      gen::error(MB_SUCCESS!=result, "could not get triangle normal vector");
      assert(MB_SUCCESS==result || MB_ENTITY_NOT_FOUND==result);
      normals.push_back( normal );
    }
    // if we've gotten here, then we have succeeded
    result = MB_SUCCESS;
    return result;
  }

  ErrorCode triangle_normal( const EntityHandle tri, CartVect &normal) {
    ErrorCode result;
    const EntityHandle *conn;
    int n_verts;
    result = MBI()->get_connectivity( tri, conn, n_verts );
    if(MB_ENTITY_NOT_FOUND == result) {
      std::cout << "triangle_normal: triangle not found" << std::endl;
      CartVect zero_vector( 0.0 );
      normal = zero_vector;
      return result;
    }else if(MB_SUCCESS != result) {
      return result;
    } else {
      assert(3 == n_verts);
      return triangle_normal( conn[0], conn[1], conn[2], normal );
    }
  }

  ErrorCode triangle_normal( const EntityHandle v0, const EntityHandle v1,
                               const EntityHandle v2, CartVect &normal ) {

    // if tri is degenerate return 0,0,0
    if( triangle_degenerate(v0, v1, v2) ) {
      CartVect zero_vector( 0.0 );
      normal = zero_vector;
      std::cout << "  normal=" << normal << std::endl;
      return MB_SUCCESS;
    }

    EntityHandle conn[3];
    conn[0] = v0;
    conn[1] = v1;
    conn[2] = v2;
    ErrorCode result;
    CartVect coords[3];
    result = MBI()->get_coords( conn, 3, coords[0].array() );
    gen::error(MB_SUCCESS!=result, "could not get coordinates of the triangle vertices");
    assert(MB_SUCCESS == result);
    return triangle_normal( coords[0], coords[1], coords[2], normal );
  }

  ErrorCode triangle_normal( const CartVect coords0, const CartVect coords1,
                               const CartVect coords2, CartVect &normal ) {
    CartVect edge0, edge1;
    edge0 = coords1-coords0;
    edge1 = coords2-coords0;
    normal = edge0*edge1;

    // do not normalize if magnitude is zero (avoid nans)
    if(0 == normal.length()) return MB_SUCCESS;

    normal.normalize();
    //if(debug) std::cout << "  normal=" << normal << std::endl;
    return MB_SUCCESS;
  }
    


  // Distance between a point and line. The line is defined by two verts.
  // We are using a line and not a line segment!
  // http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
  ErrorCode line_point_dist( const EntityHandle line_pt1, const EntityHandle line_pt2,
                               const EntityHandle pt0, double &dist ) {
    ErrorCode result;
    CartVect x0, x1, x2;
    result = MBI()->get_coords( &line_pt1, 1, x1.array() );
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &line_pt2, 1, x2.array() );
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt0, 1, x0.array() );
    assert(MB_SUCCESS == result);

    dist = ( ((x0-x1)*(x0-x2)).length() ) / ( (x2-x1).length() );
    return result;
  }

  // Project the point onto the line. Not the line segment!
  ErrorCode point_line_projection( const EntityHandle line_pt1,
                                     const EntityHandle line_pt2,
                                     const EntityHandle pt0 ) {
    CartVect projected_coords;
    double parameter;
    ErrorCode result = point_line_projection( line_pt1, line_pt2,
                                                pt0, projected_coords,
                                                parameter );
    assert(MB_SUCCESS == result);
    result = MBI()->set_coords( &pt0, 1, projected_coords.array() );   
    assert(MB_SUCCESS == result);
    return result;
  }    
  ErrorCode point_line_projection( const EntityHandle line_pt1,
                                     const EntityHandle line_pt2,
                                     const EntityHandle pt0,
                                     CartVect &projected_coords,
                                     double &parameter  ) {

    ErrorCode result;
    CartVect coords[3];
    result = MBI()->get_coords( &line_pt1, 1, coords[1].array() );
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &line_pt2, 1, coords[2].array() );           
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt0, 1, coords[0].array() );          
    assert(MB_SUCCESS == result);

    // project the t_joint between the endpts                  
    // http://en.wikipedia.org/wiki/Vector_projection           
    CartVect a = coords[0] - coords[1];
    CartVect b = coords[2] - coords[1];
    parameter    = (a%b)/(b%b);
    CartVect c = parameter*b;
    projected_coords = c     + coords[1];      
    return result;
  }
  ErrorCode point_line_projection( const EntityHandle line_pt1,
                                     const EntityHandle line_pt2,
                                     const EntityHandle pt0,
                                     double &dist_along_edge  ) {

    ErrorCode result;
    CartVect coords[3];
    result = MBI()->get_coords( &line_pt1, 1, coords[1].array() );
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &line_pt2, 1, coords[2].array() );           
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt0, 1, coords[0].array() );          
    assert(MB_SUCCESS == result);

    // project the t_joint between the endpts                  
    // http://en.wikipedia.org/wiki/Vector_projection           
    CartVect a = coords[0] - coords[1];
    CartVect b = coords[2] - coords[1];
    dist_along_edge = a%b / b.length();      
    return result;
  }
  

  double area2( const EntityHandle pt_a, const EntityHandle pt_b,
                const EntityHandle pt_c, const CartVect plane_normal ) {
    //std::cout << "area2: a=" << pt_a << " b=" << pt_b << " c=" << pt_c << std::endl;
    ErrorCode result;
    CartVect a, b, c;
    result = MBI()->get_coords( &pt_a, 1, a.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt_b, 1, b.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt_c, 1, c.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS == result);
    CartVect d = b - a;
    CartVect e = c - a;
    // project onto a plane defined by the plane's normal vector

    return (d*e)%plane_normal;
  }

  // Is point c to the left of line ab?
  bool left( const EntityHandle a, const EntityHandle b,
             const EntityHandle c, const CartVect n ) {
    double area_2 = area2(a,b,c,n);
    //std::cout << "left: a=" << a << " b=" << b << " c=" << c
    //          << " area2=" << area_2 << std::endl;
    if(area_2 > 0) return true;
    else return false;
  } 

  // Is point c to the left of line ab or collinear?
  bool left_on( const EntityHandle a, const EntityHandle b,
                const EntityHandle c, const CartVect n ) {
    double area_2 = area2(a,b,c,n);
    //std::cout << "left_on: a=" << a << " b=" << b << " c=" << c
    //          << " area2=" << area_2 << std::endl;
    if(area_2 >= 0) return true;
    else return false;
  } 

  // Are pts a,b,c collinear?
  bool collinear( const EntityHandle a, const EntityHandle b,
                  const EntityHandle c, const CartVect n ) {
    double area_2 = area2(a,b,c,n);
    //std::cout << "collinear: a=" << a << " b=" << b << " c=" << c
    //          << " area2=" << area_2 << std::endl;
    if( area_2 ==0) return true;
    else return false;
  } 

  // Exclusive or: T iff exactly one argument is true
  bool logical_xor( const bool x, const bool y ) {
    return (x || y) && !(x && y);
  }

  bool intersect_prop( const EntityHandle a, const EntityHandle b,
                       const EntityHandle c, const EntityHandle d,
                       const CartVect n ) {
    if( collinear(a,b,c,n) ||
        collinear(a,b,d,n) ||
        collinear(c,d,a,n) ||
        collinear(c,d,b,n) ) {
      return false;
    } else {
      return logical_xor(left(a,b,c,n), left(a,b,d,n)) && 
        logical_xor(left(c,d,a,n), left(c,d,b,n));
    }
  }
      
  bool between( const EntityHandle pt_a, const EntityHandle pt_b,
                const EntityHandle pt_c, const CartVect n) {
    if( !collinear(pt_a,pt_b,pt_c,n) ) return false;

    ErrorCode result;
    CartVect a, b, c;
    result = MBI()->get_coords( &pt_a, 1, a.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt_b, 1, b.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS == result);
    result = MBI()->get_coords( &pt_c, 1, c.array() );
    gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
    assert(MB_SUCCESS == result);

    // if ab not vertical, check betweenness on x; else on y.
    if(a[0] != b[0]) {
      return ((a[0] <= c[0]) && (c[0] <= b[0])) || ((a[0] >= c[0]) && (c[0] >= b[0]));
    }else if(a[1] != b[1]) {
      return ((a[1] <= c[1]) && (c[1] <= b[1])) || ((a[1] >= c[1]) && (c[1] >= b[1]));
    } else {
      return ((a[2] <= c[2]) && (c[2] <= b[2])) || ((a[2] >= c[2]) && (c[2] >= b[2]));
    }
  }

  bool intersect( const EntityHandle a, const EntityHandle b,
                  const EntityHandle c, const EntityHandle d,
                  const CartVect n ) {
    if(intersect_prop(a,b,c,d,n)) return true;
    else if( between(a,b,c,n) ||
             between(a,b,d,n) ||
             between(c,d,a,n) ||
             between(c,d,b,n) ) return true;
    else return false;
  }

  // verts is an ordered polygon of verts
  bool diagonalie( const EntityHandle a, const EntityHandle b,
                   const CartVect n,
                   const std::vector<EntityHandle> verts ) {
    for(unsigned int i=0; i<verts.size(); i++) {
      EntityHandle c = verts[i];
      EntityHandle c1;
      if(verts.size()-1 == i) c1 = verts[0];
      else                    c1 = verts[i+1];

      if( (c != a) && (c1 != a) &&
          (c != b) && (c1 != b) &&
          intersect( a, b, c, c1, n ) ) {
        //std::cout << "diagonalie a=" << a << " b=" << b << " c="
        //        << c << " c1=" << c1 << " result=";
        //std::cout << "false" << std::endl;
        return false;
      }
    }
    //std::cout << "diagonalie a=" << a << " b=" << b << " result=";
    //std::cout << "true" << std::endl;
    return true;
  }

  // verts is an ordered polygon of verts
  bool in_cone( const EntityHandle a, const EntityHandle b,
                const CartVect n,
                const std::vector<EntityHandle> verts ) {
    std::vector<EntityHandle>::const_iterator a_iter;
    a_iter = find( verts.begin(), verts.end(), a );
    EntityHandle a0, a1;
    // a0 is before a
    if(verts.begin() == a_iter) a0 = verts[verts.size()-1];
    else a0 = *(a_iter-1);
    // a1 is after a
    if(verts.end()-1 == a_iter) a1 = verts[0];
    else a1 = *(a_iter+1);

    //std::cout << "in_cone: a=" << a << " b=" << b << " a0="
    //        << a0 << " a1=" << a1 << std::endl;
    // if a is a convex vertex
    if(left_on(a,a1,a0,n)) return left(a,b,a0,n) && left(b,a,a1,n);

    // else a is reflex
    else return !(left_on(a,b,a1,n) && left_on(b,a,a0,n));
  }

  bool diagonal( const EntityHandle a, const EntityHandle b,
                 const CartVect n,
                 const std::vector<EntityHandle> verts ) {
    bool result = in_cone(a,b,n,verts) && in_cone(b,a,n,verts) && diagonalie(a,b,n,verts);
    //std::cout << "diagonal a=" << a << " b=" << b << " result="
    //          << result << std::endl;
    return result;
  }


  // Determine if each vertex is an ear. Input an ordered polygon of verts.
  ErrorCode ear_init( const std::vector<EntityHandle> verts,
                        const CartVect n, // plane normal vector
                        std::vector<bool> &is_ear ) {
    if(verts.size() != is_ear.size()) return MB_FAILURE;
    for(unsigned int i=0; i<verts.size(); i++) {
      EntityHandle prev, next;
      if(0 == i) prev = verts.back();
      else prev = verts[i-1];
      if(verts.size()-1 == i) next = verts[0];
      else next = verts[i+1];
      is_ear[i] = diagonal(prev,next,n,verts);
      //std::cout << "is_ear[" << i << "]=" << is_ear[i] << std::endl;
    }
    return MB_SUCCESS;
  }



  // Input an ordered polygon of verts and a normal vector of the plane
  // that the polygon is mostly in. The vector is required for orientation.
  ErrorCode ear_clip_polygon( std::vector<EntityHandle> verts,
                                CartVect n,
                                Range &new_tris ) {

    // initialize the status of ears
    //std::cout << "begin ear clipping----------------------" << std::endl;
    //for(unsigned int i=0; i<verts.size(); i++) {
    //  print_vertex_cubit( verts[i] );
    //}

    //print_loop( verts );
    ErrorCode result;
    std::vector<bool> is_ear( verts.size() );
    result = ear_init( verts, n, is_ear );
    assert(MB_SUCCESS == result);
 
    // if the polygon intersects itself the algorithm will not stop
    int counter = 0; 
    int n_initial_verts = verts.size();

    while(3 < verts.size()) {
      for(unsigned int i=0; i<verts.size(); i++) {
        if(is_ear[i]) {
          EntityHandle v0, v1, v2, v3, v4;
          if(0 == i)      v0 = verts[verts.size()-2];
          else if(1 == i) v0 = verts[verts.size()-1];
          else            v0 = verts[i-2]; 
          if(0 == i) v1 = verts[verts.size()-1];
          else       v1 = verts[i-1];
          v2 = verts[i];
          if(verts.size()-1 == i) v3 = verts[0];
          else                    v3 = verts[i+1];
          if(verts.size()-2 == i)       v4 = verts[0];
          else if (verts.size()-1 == i) v4 = verts[1];
          else                          v4 = verts[i+2];

          //std::cout << "ear_clip_polygon: triangle=" << std::endl;
          //print_vertex_coords( v1 ); 
          //print_vertex_coords( v2 ); 
          //print_vertex_coords( v3 );
          EntityHandle new_tri;
          EntityHandle conn[3] = {v1,v2,v3};
          result = MBI()->create_element( MBTRI, conn, 3, new_tri );
          assert(MB_SUCCESS == result);
          new_tris.insert( new_tri );

          // update ear status
          if(0 == i) is_ear[verts.size()-1] = diagonal( v0, v3, n, verts );
          else       is_ear[i-1]            = diagonal( v0, v3, n, verts );
          if(verts.size()-1 == i) is_ear[0] = diagonal( v1, v4, n, verts );
          else                    is_ear[i+1] = diagonal( v1, v4, n, verts );

          // cut off the ear at i
          verts.erase( verts.begin()+i );
          is_ear.erase( is_ear.begin()+i );
          //i--;
          break;
        }
      }

      // If the polygon has intersecting edges this loop will continue until it
      // hits this return.
      //std::cout << "counter=" << counter << " verts.size()=" << verts.size() << std::endl;
      if(counter > n_initial_verts) {
        //std::cout << "ear_clip_polygon: no ears found" << std::endl;
        result = MBI()->delete_entities( new_tris );
        assert(MB_SUCCESS == result);
        new_tris.clear();
        return MB_FAILURE;
      }
      counter++;
    } 
    //std::cout << "ear_clip_polygon: triangle=" << std::endl;
    //print_vertex_coords( verts[0] ); 
    //print_vertex_coords( verts[1] ); 
    //print_vertex_coords( verts[2] );
    EntityHandle new_tri;
    EntityHandle conn[3] = {verts[0],verts[1],verts[2]};
    result = MBI()->create_element( MBTRI, conn, 3, new_tri );
    assert(MB_SUCCESS == result);
    new_tris.insert( new_tri );
    
    return result; 
  }

  int geom_id_by_handle( const EntityHandle set ) {
    ErrorCode result;
    Tag id_tag;
    //result = MBI()->tag_create( GLOBAL_ID_TAG_NAME, sizeof(int), MB_TAG_DENSE,
    //                                MB_TYPE_INTEGER, id_tag, 0, true );           
    result = MBI()->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER,id_tag,MB_TAG_DENSE);
    gen::error(MB_SUCCESS!=result && MB_ALREADY_ALLOCATED != result, "could not get the tag handle");
    assert(MB_SUCCESS==result || MB_ALREADY_ALLOCATED==result);
    int id;
    result = MBI()->tag_get_data( id_tag, &set, 1, &id );                  
    //assert(MB_SUCCESS == result);
    return id;
  }
  
  ErrorCode save_normals( Range tris, Tag normal_tag ) {
    std::vector<CartVect> normals(tris.size());
    ErrorCode result;
    result = triangle_normals( tris, normals );
    gen::error(MB_SUCCESS!=result, "could not get triangle normals");
    assert(MB_SUCCESS == result);

    result = MBI()->tag_set_data(normal_tag, tris, &normals[0]);
    gen::error(MB_SUCCESS!=result, "could not get the tag set data");
    assert(MB_SUCCESS == result);
    return result;
  }

  ErrorCode flip(const EntityHandle tri, const EntityHandle vert0,
                   const EntityHandle vert2, const EntityHandle surf_set) {

    // get the triangles in the surface. The tri and adj_tri must be in the surface.
    Range surf_tris;
    ErrorCode result = MBI()->get_entities_by_type( surf_set, MBTRI, surf_tris);
    assert(MB_SUCCESS == result);

    // get the triangle across the edge that will be flipped
    Range adj_tri;
    EntityHandle edge[2] = {vert0, vert2};
    result = MBI()->get_adjacencies( edge, 2, 2, false, adj_tri );
    assert(MB_SUCCESS == result);
    adj_tri = intersect(adj_tri, surf_tris);
    assert(2 == adj_tri.size());
    adj_tri.erase( tri );
    print_triangle( adj_tri.front(), false );
    //result = MBI()->list_entity( adj_tri.front() );
    //assert(MB_SUCCESS == result);

    // get the remaining tri vert
    Range tri_verts;
    result = MBI()->get_adjacencies( &tri, 1, 0, false, tri_verts );
    assert(MB_SUCCESS == result);
    assert(3 == tri_verts.size());
    tri_verts.erase(vert0);
    tri_verts.erase(vert2);
    assert(1 == tri_verts.size());
    EntityHandle vert1 = tri_verts.front();

    // get the remaining adj_tri vert
    Range adj_tri_verts;
    result = MBI()->get_adjacencies( &adj_tri.front(), 1, 0, false, adj_tri_verts );
    assert(MB_SUCCESS == result);
    assert(3 == adj_tri_verts.size());
    adj_tri_verts.erase(vert0);
    adj_tri_verts.erase(vert2);
    assert(1 == adj_tri_verts.size());
    EntityHandle vert3 = adj_tri_verts.front();

    // original tri_conn    = {vert0, vert1, vert2}
    // original adj_tri_conn= {vert2, vert3, vert0}
 
    // set the new connectivity
    EntityHandle tri_conn[3] = {vert0, vert1, vert3};
    result = MBI()->set_connectivity( tri, tri_conn, 3 );
    assert(MB_SUCCESS == result);
    EntityHandle adj_tri_conn[3] = {vert1, vert2, vert3};
    result = MBI()->set_connectivity( adj_tri.front(), adj_tri_conn, 3 );
    assert(MB_SUCCESS == result);
    print_triangle( tri, false );
    print_triangle( adj_tri.front(), false );
    return result;
  }

  ErrorCode ordered_verts_from_ordered_edges( const std::vector<EntityHandle> ordered_edges,
                                                std::vector<EntityHandle> &ordered_verts ) {
    ErrorCode result = MB_SUCCESS;
    ordered_verts.clear();
    ordered_verts.reserve(ordered_edges.size()+1);

    // Save the back of the previous edge to check for continuity.
    EntityHandle previous_back_vert = 0;
    
    for(std::vector<EntityHandle>::const_iterator i=ordered_edges.begin();
        i!=ordered_edges.end(); i++) {
      const EntityHandle *conn;
      int n_verts;
      result = MBI()->get_connectivity( *i, conn, n_verts);
      gen::error(MB_SUCCESS!=result, "could not get edge connectivity");
      assert(MB_SUCCESS == result);
      assert(2 == n_verts);
      if(ordered_edges.begin() == i) {
        ordered_verts.push_back(conn[0]);
      } else {
        gen::error(previous_back_vert!=conn[0], "edges are not adjacent");
        assert(previous_back_vert == conn[0]);
      }
      ordered_verts.push_back(conn[1]);
      previous_back_vert = conn[1];
    }
    return result;
  }

  /* Find the distance between two arcs. Assume that their endpoints are somewhat
     close together. */
  ErrorCode dist_between_arcs( bool debug,
                                 const std::vector<EntityHandle> arc0,
                                 const std::vector<EntityHandle> arc1,
                                 double &dist ) {
    dist = 0;

    // Special Case: arcs have no verts.
    if( arc0.empty() || arc1.empty() ) {
      std::cout << "arc has no vertices" << std::endl;
      return MB_FAILURE;
    }

    //print_loop(arc0);
    //print_loop(arc1);

    // for simplicity, put arcs into the same structure
    std::vector<EntityHandle> arcs[2] = {arc0, arc1};

    // Special Case: Remove duplicate vert handles
    for(unsigned int i=0; i<2; ++i) {
      if( 2>arcs[i].size() ) continue;
      for(std::vector<EntityHandle>::iterator j=arcs[i].begin()+1; j!=arcs[i].end(); ++j) {
        if(*j == *(j-1)) {
          if(debug) {
            gen::print_loop( arcs[i] );
            std::cout << "dist_between_arcs: found duplicate vert handle in arc" << std::endl;
          }
          j = arcs[i].erase(j) - 1;
        }
      }
    }
    
    // get the coords in one call per arc. For speed, do not ask MOAB again for coords.
    ErrorCode result;
    std::vector<CartVect> coords[2];
    for(unsigned int i=0; i<2; i++) {
      coords[i].resize( arcs[i].size() );
      result = MBI()->get_coords( &arcs[i][0], arcs[i].size(), coords[i][0].array());
      gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
      assert(MB_SUCCESS == result);
    }

    // Special case: arc has 1 vert or a length of zero
    bool point_arc_exists = false;
    unsigned int point_arc_index;
    for(unsigned int i=0; i<2; ++i) {
      if( 1==arcs[i].size() ) {
        point_arc_exists = true;
        point_arc_index  = i;
        break;
      } else if ( 0==length(arcs[i]) ) {
        point_arc_exists = true;
        point_arc_index  = i;
        break;
      }
    }
    // If the special case exists, we can still get an average distance
    if(point_arc_exists) {
      unsigned int other_arc_index = (0==point_arc_index) ? 1 : 0;
      // Both are point arcs
      if(1==arcs[other_arc_index].size()) {
        dist = dist_between_verts( arcs[point_arc_index].front(), arcs[other_arc_index].front());
        return MB_SUCCESS;
        // The other arc has more than one point
      } else {
        double area = 0.0;
        for(unsigned int i=0; i<arcs[other_arc_index].size()-1; ++i) {
          area += fabs(triangle_area( coords[other_arc_index][i], 
                                      coords[other_arc_index][i+1], 
                                      coords[point_arc_index].front() ));
        }
        dist = area / gen::length(arcs[other_arc_index]);
        return MB_SUCCESS;
      }
    }

    // get the arc length and parametric coordinates. The parametrize the arcs by
    // length from 0 to 1.
    double arc_len[2] = {0.0};
    std::vector<double> params[2];
    for(unsigned int i=0; i<2; i++) {
      params[i].resize( arcs[i].size() );
      params[i][0] = 0.0;
      for(unsigned int j=0; j<coords[i].size()-1; j++) {
        double d = dist_between_verts( coords[i][j], coords[i][j+1] );
        arc_len[i] += d;
        params[i][1+j] = arc_len[i];
      }
      for(unsigned int j=0; j<coords[i].size(); j++) {
        params[i][j] /= arc_len[i];
        //std::cout << "params[" << i << "][" << j << "]=" << params[i][j] << std::endl;
      }
    }

    // Merge the two sets of parameters into one ordered set without duplicates.
    std::vector<double> mgd_params;
    mgd_params.reserve( params[0].size() + params[1].size() );
    unsigned int a=0, b=0;
    while(a<params[0].size() && b<params[1].size()) {
      if       (params[0][a] < params[1][b]) {
        mgd_params.push_back( params[0][a] );
        a++;
      } else if(params[0][a] > params[1][b]) {
        mgd_params.push_back( params[1][b] );
        b++;
      } else {
        mgd_params.push_back( params[0][a] );
        a++;
        b++;
      }
    }

    for(unsigned int i=0; i<mgd_params.size(); i++) {
      //std::cout << "mgd_params[" << i << "]=" << mgd_params[i] << std::endl;
    }

    // Insert new points to match the other arcs points, by parameter.
    for(unsigned int i=0; i<2; i++) {
      for(unsigned int j=0; j<mgd_params.size(); j++) {
        //std::cout << "params[" << i << "][" << j << "]=" << params[i][j] 
        //          << " mgd_params[" << j << "]=" << mgd_params[j] << std::endl;
        if(params[i][j] > mgd_params[j]) {
          double ratio = (mgd_params[j]-params[i][j-1]) / (params[i][j]-params[i][j-1]);
          //std::cout << "j=" << j << " ratio=" << ratio << std::endl;
          CartVect pt = coords[i][j-1] + ratio*(coords[i][j]-coords[i][j-1]);
          coords[i].insert( coords[i].begin()+j, pt);
          params[i].insert( params[i].begin()+j, mgd_params[j]); 
        }
      }
    }


    // Each arc should now have the same number of points
    if(coords[0].size() != coords[1].size()) {
      for(unsigned int i=0; i<2; i++) {
        for(unsigned int j=0; j<coords[i].size(); j++) {
          std::cout << "coords[" << i << "][" << j << "]=" << coords[i][j] << std::endl;
        }
      }
    }
    assert( coords[0].size() == coords[1].size() );

    // Get the area between arcs. Use absolute value to prevent cancelling.
    double area = 0;
    for(unsigned int i=0; i<coords[0].size()-1; i++) {
      area += fabs(triangle_area( coords[0][i], coords[0][i+1], coords[1][i+1] ));
      //std::cout << "area0=" << area << std::endl;
      area += fabs(triangle_area( coords[0][i], coords[1][i+1], coords[1][i] ));
      //std::cout << "area1=" << area << std::endl;
    }

    // Divide the area by the average length to get the average distance between arcs.
    dist = fabs(2*area / (arc_len[0] + arc_len[1] ));
    //std::cout << "dist_between_arcs=" << dist << std::endl;
    return MB_SUCCESS;
  }


    // qsort struct comparision function        
    // If the first handle is the same, compare the second                                        
    int compare_edge(const void *a, const void *b) {                                 
      struct edge *ia = (struct edge *)a;           
      struct edge *ib = (struct edge *)b;
      if(ia->v0 == ib->v0) {        
        return (int)(ia->v1 - ib->v1);                   
      } else {
        return (int)(ia->v0 - ib->v0);                   
      }
      /* float comparison: returns negative if b > a           
      and positive if a > b. We multiplied result by 100.0        
      to preserve decimal fraction */                
    }  

  // WARNING: This skinner goes 10x faster by assuming that no edges already exist
  // in the MOAB instance. Otherwise checking to see if an edge exists before
  // creating a new one if very slow. This is partly the reason that Skinner is
  // very slow.
  ErrorCode find_skin( Range tris, const int dim,
 //                         std::vector<std::vector<EntityHandle> > &skin_edges,
                         Range &skin_edges,
                         const bool temp_bool ) {    

    const bool local_debug = false;
    //Skinner tool(MBI());
    //Range skin_verts;
    //return tool.find_skin( tris, dim, skin_edges, temp_bool );
    //return tool.find_skin_vertices( tris, skin_verts, &skin_edges, true );
  
    if(1 != dim) return MB_NOT_IMPLEMENTED;
    if(MBTRI != MBI()->type_from_handle(tris.front())) return MB_NOT_IMPLEMENTED;

    skin_edges.clear();
    if(tris.empty()) return MB_ENTITY_NOT_FOUND;

    // This implementation gets some of its speed due to not checking for edges
    ErrorCode result;
    int n_edges;
    result = MBI()->get_number_entities_by_type( 0, MBEDGE, n_edges );
    assert(MB_SUCCESS == result);
    if(0 != n_edges) {
      Range temp_edges;
      result = MBI()->get_entities_by_type( 0, MBEDGE, temp_edges);
      assert(MB_SUCCESS == result);
      result = MBI()->list_entities( temp_edges );
      assert(MB_SUCCESS == result);
    }      
    assert(0 == n_edges);

    // Get connectivity. Do not create edges.
    edge *edges = new edge[3*tris.size()];
    int n_verts;
    int ii = 0;
    for(Range::iterator i=tris.begin(); i!=tris.end(); i++) {
      const EntityHandle *conn;
      result = MBI()->get_connectivity( *i, conn, n_verts );
      assert(MB_SUCCESS == result);
      assert(3 == n_verts);
      // shouldn't be degenerate
      assert(conn[0] != conn[1]);
      assert(conn[1] != conn[2]);
      assert(conn[2] != conn[0]);
      // make edges
      edges[3*ii+0].v0 = conn[0];
      edges[3*ii+0].v1 = conn[1];
      edges[3*ii+1].v0 = conn[1];
      edges[3*ii+1].v1 = conn[2];
      edges[3*ii+2].v0 = conn[2];
      edges[3*ii+2].v1 = conn[0];
      ii++;
    }

    // Change the first handle to be lowest
    for(unsigned int i=0; i<3*tris.size(); ++i) {
      if(edges[i].v0 > edges[i].v1) {
        EntityHandle temp = edges[i].v0;
        edges[i].v0 = edges[i].v1;
        edges[i].v1 = temp;
      }
    }

    // Sort by first handle, then second handle. Do not sort the extra edge on the
    // back.
    qsort(edges, 3*tris.size(), sizeof(struct edge), compare_edge);    

    // Go through array, saving edges that are not paired.
    for(unsigned int i=0; i<3*tris.size(); i++) {
      // If the last edge has not been paired, create it. This avoids overrunning
      // the edges array with i+1.
      if(3*tris.size()-1 == i) {
        const EntityHandle conn[2] = {edges[i].v0, edges[i].v1};
        EntityHandle edge;
        result = MBI()->create_element( MBEDGE, conn, 2, edge );
        assert(MB_SUCCESS == result);
        skin_edges.insert(edge);
    
      // If a match exists, skip ahead
      } else if(edges[i].v0==edges[i+1].v0 && edges[i].v1==edges[i+1].v1) {
        i++;
        // test to make sure surface is manifold
        if(3*tris.size() != i+1) { // avoid overrunning array
          while( edges[i].v0==edges[i+1].v0 && edges[i].v1==edges[i+1].v1 ) {
            if(local_debug) {
              std::cout << "find_skin WARNING: non-manifold edge" << std::endl;
              MBI()->list_entity( edges[i].v0 );
              MBI()->list_entity( edges[i].v1 );
            }
            ++i;
          }
        }
        // otherwise a skin edge has been found
      } else {
        const EntityHandle conn[2] = {edges[i].v0, edges[i].v1};
        EntityHandle edge;
        result = MBI()->create_element( MBEDGE, conn, 2, edge );
        if(gen::error(MB_SUCCESS!=result, "could not create edge element")) return result;
        skin_edges.insert( edge );
      } 
    }
    delete[] edges;
    return MB_SUCCESS;
  }
  /*  ErrorCode find_skin( Range tris, const int dim, Range &skin_edges, const bool temp ) {
    std::vector<std::vector<EntityHandle> > skin_edges_vctr;
    ErrorCode result = find_skin( tris, dim, skin_edges_vctr, temp );
    assert(MB_SUCCESS == result);
    for(std::vector<std::vector<EntityHandle> >::const_iterator i=skin_edges_vctr.begin();
        i!=skin_edges_vctr.end(); i++) {
      EntityHandle edge;
      result = MBI()->create_element( MBEDGE, &(*i)[0], 2, edge );
      if(MB_SUCCESS != result) std::cout << "result=" << result << std::endl;
      assert(MB_SUCCESS == result);
      skin_edges.insert( edge );
    }
    return MB_SUCCESS;
    }*/

// calculate volume of polyhedron
// Copied from DagMC, without index_by_handle. The dagmc function will
// segfault if build_indices is not first called. For sealing there is
// no need to build_indices.
  ErrorCode measure_volume( const EntityHandle volume, double& result, bool debug, bool verbose )
{
  ErrorCode rval;
  std::vector<EntityHandle> surfaces, surf_volumes;
  result = 0.0;
  
 

  // get surfaces from volume
  rval = MBI()->get_child_meshsets( volume, surfaces );
  if (MB_SUCCESS != rval)
    {
      return rval;
    }
  if(debug) std::cout << "in measure_volume 1" << std::endl;

    // get surface senses
  std::vector<int> senses( surfaces.size() );


  if (rval != MB_SUCCESS)
    {
      std::cout << "Could not measure volume" << std::endl;
      std::cout << "This can happen for 2 reasons, there are no volumes" << std::endl;
      std::cout << "or the pointer to the Moab instance is poor" << std::endl;
      exit(rval);
    }
  if(debug) std::cout << surfaces.size() << " " << result << std::endl;

  rval = surface_sense( volume, surfaces.size(), &surfaces[0], &senses[0] );
  

  if(debug) std::cout << "in measure_volume 2" << std::endl;

  if (MB_SUCCESS != rval)
    {
      std::cerr << "ERROR: Surface-Volume relative sense not available. "
                << "Cannot calculate volume." << std::endl;
      return rval;
    }


  for (unsigned i = 0; i < surfaces.size(); ++i) 
    {
      // skip non-manifold surfaces
      if (!senses[i])
        continue;
    
      // get triangles in surface
      Range triangles;
      rval = MBI()->get_entities_by_dimension( surfaces[i], 2, triangles );
      if (MB_SUCCESS != rval)
        {
          return rval;
        }
    if (!triangles.all_of_type(MBTRI)) 
      {
        std::cout << "WARNING: Surface " << geom_id_by_handle(surfaces[i])
                  << " contains non-triangle elements. Volume calculation may be incorrect."
                  << std::endl;
        triangles.clear();
        rval = MBI()->get_entities_by_type( surfaces[i], MBTRI, triangles );
        if (MB_SUCCESS != rval)
          {
            return rval;
          }
    }
    
      // calculate signed volume beneath surface (x 6.0)
    double surf_sum = 0.0;
    const EntityHandle *conn;
    int len;
    CartVect coords[3];
    for (Range::iterator j = triangles.begin(); j != triangles.end(); ++j) {
      rval = MBI()->get_connectivity( *j, conn, len, true );
      if (MB_SUCCESS != rval) return rval;
      assert(3 == len);
      rval = MBI()->get_coords( conn, 3, coords[0].array() );
      if (MB_SUCCESS != rval) return rval;
    
      coords[1] -= coords[0];
      coords[2] -= coords[0];
      surf_sum += (coords[0] % (coords[1] * coords[2]));
    }
    result += senses[i] * surf_sum;
  }
  
  result /= 6.0;
  if(debug)  std::cout << result << std::endl;
  return MB_SUCCESS;
}

  ErrorCode get_signed_volume( const EntityHandle surf_set, double &signed_volume) {
    ErrorCode rval;
    Range tris;
    rval = MBI()->get_entities_by_type( surf_set, MBTRI, tris );
    if(MB_SUCCESS != rval) return rval;
    signed_volume = 0.0;
    const EntityHandle *conn;
    int len;
    CartVect coords[3];
    for (Range::iterator j = tris.begin(); j != tris.end(); ++j) {
      rval = MBI()->get_connectivity( *j, conn, len, true );           
      if (MB_SUCCESS != rval) return rval;
      assert(3 == len);                    
      rval = MBI()->get_coords( conn, 3, coords[0].array() );              
      if (MB_SUCCESS != rval) return rval;
                            
      coords[1] -= coords[0];                              
      coords[2] -= coords[0];                     
      signed_volume += (coords[0] % (coords[1] * coords[2]));       
    }                                     
    return MB_SUCCESS;
  }                 

ErrorCode measure( const EntityHandle set, const Tag geom_tag, double &size, bool debug,  bool verbose ) {
    ErrorCode result;
    int dim;
    result = MBI()->tag_get_data( geom_tag, &set, 1, &dim );                  
    assert(MB_SUCCESS == result);
    if(0 == dim) {
      std::cout << "measure: cannot measure vertex" << std::endl;
      return MB_FAILURE;

    } else if(1 == dim) {
      std::vector<EntityHandle> vctr;
      result = arc::get_meshset( set, vctr );
      assert(MB_SUCCESS == result);
      size = length( vctr );

    } else if(2 == dim) {
      Range tris;
      result = MBI()->get_entities_by_type( set, MBTRI, tris );
      assert(MB_SUCCESS == result);
      size = triangle_area( tris );

    } else if(3 == dim) {
      
      result = measure_volume( set, size, debug, verbose );
      
      if(MB_SUCCESS != result) {
        std::cout << "result=" << result << " vol_id=" 
                  << gen::geom_id_by_handle(set) << std::endl;
      }
    } else {
      std::cout << "measure: incorrect dimension" << std::endl;
      return MB_FAILURE;
    }
    return MB_SUCCESS;
  }

  // From CGMA/builds/dbg/include/CubitDefines
  ErrorCode get_curve_surf_sense( const EntityHandle surf_set, const EntityHandle curve_set,
                                    int &sense, bool debug ) {
    std::vector<EntityHandle> surfs;
    std::vector<int> senses;
    ErrorCode rval;
    GeomTopoTool gt( MBI(), false);
    rval = gt.get_senses( curve_set, surfs, senses );
    if(gen::error(MB_SUCCESS!=rval,"failed to get_senses")) return rval;
    
    unsigned counter = 0;
    int edim;
    for(unsigned i=0; i<surfs.size(); i++) {
      edim=gt.dimension(surfs[i]);
      if( edim == -1)
      { 
        surfs.erase(surfs.begin()+i);
        senses.erase(senses.begin()+i);
        i=0;
      }
     if(surf_set==surfs[i]) {
        sense = senses[i];
        
        ++counter;
      }
    }

   if(debug)
   {
    for(unsigned int index=0; index < surfs.size() ; index++)
    {
     std::cout << "surf = " << geom_id_by_handle(surfs[index]) << std::endl;
     std::cout << "sense = " << senses[index] << std::endl;
    }
   }
    // special case: parent surface does not exist
    if(gen::error(0==counter,"failed to find a surf in sense list")) return MB_FAILURE;

   

     // special case: ambiguous
    if(1<counter) sense = SENSE_UNKNOWN;
    return MB_SUCCESS;
  }

  ErrorCode surface_sense( EntityHandle volume,
                           int num_surfaces,
                           const EntityHandle* surfaces,
                           int* senses_out )
  {
    std::vector<EntityHandle> surf_volumes( 2*num_surfaces );
    Tag senseTag = get_tag( "GEOM_SENSE_2", 2, MB_TAG_SPARSE, MB_TYPE_HANDLE, NULL, false );
    ErrorCode rval = MBI()->tag_get_data( senseTag , surfaces, num_surfaces, &surf_volumes[0] );
    if (MB_SUCCESS != rval)
      {
        return rval;
      }
  
    const EntityHandle* end = surfaces + num_surfaces;
    std::vector<EntityHandle>::const_iterator surf_vols = surf_volumes.begin();
    while (surfaces != end) 
      {
        EntityHandle forward = *surf_vols; ++surf_vols;
        EntityHandle reverse = *surf_vols; ++surf_vols;
        if (volume == forward) 
          {
            *senses_out = (volume != reverse); // zero if both, otherwise 1
          }
        else if (volume == reverse)
          {
            *senses_out = SENSE_UNKNOWN;
          }
        else 
          {
            return MB_ENTITY_NOT_FOUND;
          }
    
        ++surfaces;
        ++senses_out;
      }
  
    return MB_SUCCESS;
  }

  ErrorCode surface_sense( EntityHandle volume,
                             EntityHandle surface,
                             int& sense_out )
  {
    // get sense of surfaces wrt volumes
    EntityHandle surf_volumes[2];
    Tag senseTag = get_tag( "GEOM_SENSE_2", 2, MB_TAG_SPARSE, MB_TYPE_HANDLE, NULL, false );
    ErrorCode rval = MBI()->tag_get_data( senseTag , &surface, 1, surf_volumes );
    if (MB_SUCCESS != rval)
      {
        return rval;
      }
  
    if (surf_volumes[0] == volume)
      {
        sense_out = (surf_volumes[1] != volume); // zero if both, otherwise 1
      }
    else if (surf_volumes[1] == volume)
      {
        sense_out = SENSE_UNKNOWN;
      }
    else
      {
        return MB_ENTITY_NOT_FOUND;
      }
  
  return MB_SUCCESS;
 }

 Tag get_tag( const char* name, int size, TagType store,
                     DataType type, const void* def_value,
                     bool create_if_missing)
 {
   Tag retval = 0;
   unsigned flags = store|MB_TAG_CREAT;
   if (!create_if_missing)
     {
       flags |= MB_TAG_EXCL;
     }
   ErrorCode result = MBI()->tag_get_handle(name, size, type, retval, flags, def_value);
   if (create_if_missing && MB_SUCCESS != result)
     {
       std::cerr << "Couldn't find nor create tag named " << name << std::endl;
     }

   return retval;
 }



ErrorCode delete_surface( EntityHandle surf, Tag geom_tag, Range tris, int id, bool debug, bool verbose ) {

  ErrorCode result;

  //measure area of the surface  
        double area;
        result = gen::measure( surf, geom_tag, area, debug, verbose );
        if(gen::error(MB_SUCCESS!=result,"could not measure area")) return result;
        assert(MB_SUCCESS == result);

  //remove triagngles from the surface
        result = MBI()->remove_entities( surf, tris);                          
        if(gen::error(MB_SUCCESS!=result,"could not remove tris")) return result;
        assert(MB_SUCCESS == result);
  //print information about the deleted surface if requested by the user
        if(debug) std::cout << "  deleted surface " << id
                  << ", area=" << area << " cm^2, n_facets=" << tris.size() << std::endl;

  // delete triangles from mesh data
        result = MBI()->delete_entities( tris );                               
        if(gen::error(MB_SUCCESS!=result,"could not delete tris")) return result;
        assert(MB_SUCCESS == result);

  //remove the sense data for the surface from its child curves
        result = remove_surf_sense_data( surf, debug);
        if(gen::error(MB_SUCCESS!=result, "could not remove surface's sense data")) return result;

  //remove the surface set itself
        result = MBI()->delete_entities( &(surf), 1);
        if(gen::error(MB_SUCCESS!=result,"could not delete surface set")) return result;
        assert(MB_SUCCESS == result);


  return MB_SUCCESS;
 }


ErrorCode remove_surf_sense_data(EntityHandle del_surf, bool debug) {
 
  ErrorCode result;
  GeomTopoTool gt(MBI(), false);
    int edim = gt.dimension(del_surf);

    if(gen::error(edim!=2,"could not remove sense data: entity is of the wrong dimension")) return MB_FAILURE;

 // get the curves of the surface
        Range del_surf_curves;
        result = MBI() -> get_child_meshsets( del_surf, del_surf_curves);
        if(gen::error(MB_SUCCESS!=result,"could not get the curves of the surface to delete")) return result;
        if (debug) std::cout << "got the curves" << std::endl;
       
  

        if (debug)
        {
          std::cout << "number of curves to the deleted surface = " << del_surf_curves.size() << std::endl;
          for(unsigned int index =0 ; index < del_surf_curves.size() ; index++)
          {
            std::cout << "deleted surface's child curve id " << index << " = " << gen::geom_id_by_handle(del_surf_curves[index]) << std::endl;
          }        
        }
        //get the sense data for each curve

        //get sense_tag handles from MOAB
        Tag senseEnts, senseSenses;
        unsigned flags = MB_TAG_SPARSE;
     
        //get tag for the entities with sense data associated with a given moab entity
        result = MBI()-> tag_get_handle(GEOM_SENSE_N_ENTS_TAG_NAME, 0, MB_TYPE_HANDLE, senseEnts, flags);
        if(gen::error(MB_SUCCESS!=result, "could not get senseEnts tag")) return result;
        
        //get tag for the sense data associated with the senseEnts entities for a given moab entity
        result = MBI()-> tag_get_handle(GEOM_SENSE_N_SENSES_TAG_NAME, 0, MB_TYPE_INTEGER, senseSenses, flags);
        if(gen::error(MB_SUCCESS!=result,"could not get senseSenses tag")) return result;
        
        //initialize vectors for entities and sense data
        std::vector<EntityHandle> surfaces;
        std::vector<int> senses;
        for(Range::iterator i=del_surf_curves.begin(); i!=del_surf_curves.end(); i++ )
        {
       
        result = gt.get_senses(*i, surfaces, senses);
        if(gen::error(MB_SUCCESS!=result, "could not get the senses for the del_surf_curve")) return result;
          // if the surface to be deleted (del_surf) exists in the sense data (which it should), then remove it
          for(unsigned int index = 0; index < senses.size() ; index++)
          {
            if(surfaces[index]==del_surf)
            {
             surfaces.erase(surfaces.begin() + index);
             senses.erase(senses.begin() +index);
             index=-1;
            }
          }
          //remove existing sense entity data for the curve
          result= MBI()-> tag_delete_data( senseEnts, &*i, 1);
          if(gen::error(MB_SUCCESS!=result, "could not delete sense entity data")) return result;
       
          //remove existing sense data for the curve
          result = MBI()-> tag_delete_data(senseSenses, &*i, 1);
          if(gen::error(MB_SUCCESS!=result, "could not delete sense data")) return result;

          //reset the sense data for each curve 
          result = gt.set_senses( *i, surfaces, senses);
          if(gen::error(MB_SUCCESS!=result, "could not update sense data for surface deletion")) return result;

        }

 return MB_SUCCESS;
 } 

 ErrorCode combine_merged_curve_senses( std::vector<EntityHandle> &curves, Tag merge_tag, bool debug) {

  ErrorCode result;
  
  for(std::vector<EntityHandle>::iterator j=curves.begin(); j!=curves.end(); j++) {


  EntityHandle merged_curve;
  result = MBI() -> tag_get_data( merge_tag, &(*j), 1, &merged_curve);
  if(gen::error(MB_SUCCESS!=result && MB_TAG_NOT_FOUND!=result, "could not get the merge_tag data of the curve")) return result;

  if(MB_SUCCESS==result) { // we have found a merged curve pairing
    // add the senses from the curve_to_delete to curve_to keep
    // create vectors for the senses and surfaces of each curve
    std::vector<EntityHandle> curve_to_keep_surfs, curve_to_delete_surfs, combined_surfs;
    std::vector<int> curve_to_keep_senses, curve_to_delete_senses, combined_senses;

    //initialize GeomTopoTool.cpp instance in MOAB
    GeomTopoTool gt(MBI(), false);
    // get senses of the iterator curve and place them in the curve_to_delete vectors
    result = gt.get_senses( *j, curve_to_delete_surfs, curve_to_delete_senses);
    if(gen::error(MB_SUCCESS!=result, "could not get the surfs/senses of the curve to delete")) return result;
    // get surfaces/senses of the merged_curve and place them in the curve_to_keep vectors
    result = gt.get_senses( merged_curve, curve_to_keep_surfs, curve_to_keep_senses);
    if(gen::error(MB_SUCCESS!=result, "could not get the surfs/senses of the curve to delete")) return result;
    
    if(debug){
          std::cout << "curve to keep id = " << gen::geom_id_by_handle(merged_curve) << std::endl;
          std::cout << "curve to delete id = " << gen::geom_id_by_handle(*j) << std::endl;
          for(unsigned int index=0; index < curve_to_keep_surfs.size(); index++)
          {
           std::cout << "curve_to_keep_surf " << index << " id = " << gen::geom_id_by_handle(curve_to_keep_surfs[index]) << std::endl;
           std::cout << "curve_to_keep_sense " << index << " = " << curve_to_keep_senses[index] << std::endl;
           
          }
          for(unsigned int index=0; index < curve_to_keep_surfs.size(); index++)
          {
          
          std::cout << "curve_to_delete_surf " << index << " id = " << gen::geom_id_by_handle(curve_to_delete_surfs[index]) << std::endl;
          std::cout << "curve_to_delete_sense " << index << " = " << curve_to_delete_senses[index] << std::endl;
          }
    } // end of debug st.

    // combine the surface and sense data for both curves into the same vector 

    combined_surfs.insert( combined_surfs.end(), curve_to_keep_surfs.begin(), 
                                                       curve_to_keep_surfs.end() );
    combined_surfs.insert( combined_surfs.end(), curve_to_delete_surfs.begin(), 
                                                       curve_to_delete_surfs.end() );
    combined_senses.insert(combined_senses.end(),curve_to_keep_senses.begin(),
                                                       curve_to_keep_senses.end() );
    combined_senses.insert(combined_senses.end(),curve_to_delete_senses.begin(),
                                                       curve_to_delete_senses.end() );

    if(debug){
      std::cout << combined_surfs.size() << std::endl;
      std::cout << combined_senses.size() << std::endl;
      for(unsigned int index=0; index < combined_senses.size(); index++)
          {

           std::cout << "combined_surfs{"<< index << "] = " << gen::geom_id_by_handle(combined_surfs[index]) << std::endl;
           std::cout << "combined_sense["<< index << "] = " << combined_senses[index] << std::endl;
          }
    } // end debug st. 

    result = gt.set_senses(merged_curve, combined_surfs, combined_senses);
    if(gen::error(MB_SUCCESS!=result && MB_MULTIPLE_ENTITIES_FOUND!=result,"failed to set senses: "));


    
    // add the duplicate curve_to_keep to the vector of curves
    *j = merged_curve;

   } //end merge_tag result if st.


  } //end curves loop




  return MB_SUCCESS;
  }

ErrorCode get_sealing_mesh_tags( double &facet_tol,
                           double &sme_resabs_tol,
                           Tag &geom_tag,
                           Tag &id_tag,
                           Tag &normal_tag,
                           Tag &merge_tag,
                           Tag &faceting_tol_tag,
                           Tag &geometry_resabs_tag,
                           Tag &size_tag,
                           Tag &orig_curve_tag) {

  ErrorCode result;

    result = MBI()->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1,
                                MB_TYPE_INTEGER, geom_tag, MB_TAG_DENSE|MB_TAG_CREAT );
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    result = MBI()->tag_get_handle( GLOBAL_ID_TAG_NAME, 1,
                                MB_TYPE_INTEGER, id_tag, MB_TAG_DENSE|MB_TAG_CREAT);
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    result = MBI()->tag_get_handle( "NORMAL", sizeof(CartVect), MB_TYPE_OPAQUE,
        normal_tag, MB_TAG_DENSE|MB_TAG_CREAT);
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    result = MBI()->tag_get_handle( "MERGE", 1, MB_TYPE_HANDLE,
        merge_tag, MB_TAG_SPARSE|MB_TAG_CREAT );
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      } 
    result = MBI()->tag_get_handle( "FACETING_TOL", 1, MB_TYPE_DOUBLE,
        faceting_tol_tag , MB_TAG_SPARSE|MB_TAG_CREAT );
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    result = MBI()->tag_get_handle( "GEOMETRY_RESABS", 1,     MB_TYPE_DOUBLE,
                             geometry_resabs_tag, MB_TAG_SPARSE|MB_TAG_CREAT  );
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    result = MBI()->tag_get_handle( "GEOM_SIZE", 1, MB_TYPE_DOUBLE,
                                    size_tag, MB_TAG_DENSE|MB_TAG_CREAT  );
    assert( (MB_SUCCESS == result) );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    int true_int = 1;    
    result = MBI()->tag_get_handle( "ORIG_CURVE", 1,
                                MB_TYPE_INTEGER, orig_curve_tag, MB_TAG_DENSE|MB_TAG_CREAT, &true_int );
    assert( MB_SUCCESS == result );
    if ( result != MB_SUCCESS )
      {
        moab_printer(result);
      }
    // PROBLEM: MOAB is not consistent with file_set behavior. The tag may not be
    // on the file_set.
    Range file_set;
    result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &faceting_tol_tag,
                                                  NULL, 1, file_set );

    if(gen::error(MB_SUCCESS!=result,"could not get faceting_tol_tag"))
      {
        return result;
      }

    gen::error(file_set.empty(),"file set not found");

    if(gen::error(1!=file_set.size(),"Refacet with newer version of ReadCGM.")) 
      {
        return MB_FAILURE;
      }

    result = MBI()->tag_get_data( faceting_tol_tag, &file_set.front(), 1,  
                                  &facet_tol );
    assert(MB_SUCCESS == result);
    result = MBI()->tag_get_data( geometry_resabs_tag, &file_set.front(), 1,  
                                  &sme_resabs_tol );
    if(MB_SUCCESS != result)
      {
        std::cout <<  "absolute tolerance could not be read from file" << std::endl;
      }



  return MB_SUCCESS;

  }

  ErrorCode get_geometry_meshsets( Range geometry_sets[], Tag geom_tag, bool verbose) {

    ErrorCode result;

    // get all geometry sets
    for(unsigned dim=0; dim<4; dim++) 
      {
        void *val[] = {&dim};
        result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag,
                                                    val, 1, geometry_sets[dim] );
        assert(MB_SUCCESS == result);

        // make sure that sets TRACK membership and curves are ordered
        // MESHSET_TRACK_OWNER=0x1, MESHSET_SET=0x2, MESHSET_ORDERED=0x4      

        for(Range::iterator i=geometry_sets[dim].begin(); i!=geometry_sets[dim].end(); i++)
          {
            unsigned int options;
            result = MBI()->get_meshset_options(*i, options );
            assert(MB_SUCCESS == result);
    
            // if options are wrong change them
            if(dim==1) 
              {
                if( !(MESHSET_TRACK_OWNER&options) || !(MESHSET_ORDERED&options) ) 
                  {
                    result = MBI()->set_meshset_options(*i, MESHSET_TRACK_OWNER|MESHSET_ORDERED);
                    assert(MB_SUCCESS == result);
                  }
              } 
            else 
              {
                if( !(MESHSET_TRACK_OWNER&options) ) 
                  {        
                    result = MBI()->set_meshset_options(*i, MESHSET_TRACK_OWNER);
                    assert(MB_SUCCESS == result);
                  }
              }
          }
      }
 
    return result;

    }

  ErrorCode check_for_geometry_sets(Tag geom_tag, bool verbose){

    ErrorCode result;
        // go get all geometry sets
        Range geometry_sets[4];
        result = get_geometry_meshsets( geometry_sets, geom_tag, false);
        if(gen::error(MB_SUCCESS!=result,"could not get the geometry meshsets")) return result;

        //make sure they're there
        for(unsigned dim=0; dim<4; dim++){
 
          if(geometry_sets[dim].size() == 0) return MB_FAILURE;

        }
 
    return MB_SUCCESS;
  }


} //EOL

Interface *MBI()
{
    static Core instance;
    return &instance;
}