Skinner.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
#ifdef __MFC_VER
#pragma warning( disable : 4786 )
#endif
 
#ifdef WIN32 /* windows */
#define _USE_MATH_DEFINES // For M_PI
#endif
#include "moab/Skinner.hpp"
#include "moab/Range.hpp"
#include "moab/CN.hpp"
#include <vector>
#include <set>
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "moab/Util.hpp"
#include "Internals.hpp"
#include "MBTagConventions.hpp"
#include "moab/Core.hpp"
#include "AEntityFactory.hpp"
#include "moab/ScdInterface.hpp"
 
#ifdef M_PI
#define SKINNER_PI M_PI
#else
#define SKINNER_PI 3.1415926535897932384626
#endif
 
namespace moab
{
 
Skinner::~Skinner()
{
 // delete the adjacency tag
}
 
ErrorCode Skinner::initialize()
{
 // go through and mark all the target dimension entities
 // that already exist as not deleteable
 // also get the connectivity tags for each type
 // also populate adjacency information
 EntityType type;
 DimensionPair target_ent_types = CN::TypeDimensionMap[mTargetDim];
 
 void* null_ptr = NULL;
 
 ErrorCode result = thisMB->tag_get_handle( "skinner adj", sizeof( void* ), MB_TYPE_OPAQUE, mAdjTag,
 MB_TAG_DENSE | MB_TAG_CREAT, &null_ptr );MB_CHK_ERR( result );
 
 if( mDeletableMBTag == 0 )
 {
 result =
 thisMB->tag_get_handle( "skinner deletable", 1, MB_TYPE_BIT, mDeletableMBTag, MB_TAG_BIT | MB_TAG_CREAT );MB_CHK_ERR( result );
 }
 
 Range entities;
 
 // go through each type at this dimension
 for( type = target_ent_types.first; type <= target_ent_types.second; ++type )
 {
 // get the entities of this type in the MB
 thisMB->get_entities_by_type( 0, type, entities );
 
 // go through each entity of this type in the MB
 // and set its deletable tag to NO
 Range::iterator iter, end_iter;
 end_iter = entities.end();
 for( iter = entities.begin(); iter != end_iter; ++iter )
 {
 unsigned char bit = 0x1;
 result = thisMB->tag_set_data( mDeletableMBTag, &( *iter ), 1, &bit );
 assert( MB_SUCCESS == result );
 // add adjacency information too
 if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) add_adjacency( *iter );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode Skinner::deinitialize()
{
 ErrorCode result;
 if( 0 != mDeletableMBTag )
 {
 result = thisMB->tag_delete( mDeletableMBTag );
 mDeletableMBTag = 0;MB_CHK_ERR( result );
 }
 
 // remove the adjacency tag
 std::vector< std::vector< EntityHandle >* > adj_arr;
 std::vector< std::vector< EntityHandle >* >::iterator i;
 if( 0 != mAdjTag )
 {
 for( EntityType t = MBVERTEX; t != MBMAXTYPE; ++t )
 {
 Range entities;
 result = thisMB->get_entities_by_type_and_tag( 0, t, &mAdjTag, 0, 1, entities );MB_CHK_ERR( result );
 adj_arr.resize( entities.size() );
 result = thisMB->tag_get_data( mAdjTag, entities, &adj_arr[0] );MB_CHK_ERR( result );
 for( i = adj_arr.begin(); i != adj_arr.end(); ++i )
 delete *i;
 }
 
 result = thisMB->tag_delete( mAdjTag );
 mAdjTag = 0;MB_CHK_ERR( result );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode Skinner::add_adjacency( EntityHandle entity )
{
 std::vector< EntityHandle >* adj = NULL;
 const EntityHandle* nodes;
 int num_nodes;
 ErrorCode result = thisMB->get_connectivity( entity, nodes, num_nodes, true );MB_CHK_ERR( result );
 const EntityHandle* iter = std::min_element( nodes, nodes + num_nodes );
 
 if( iter == nodes + num_nodes ) return MB_SUCCESS;
 
 // add this entity to the node
 if( thisMB->tag_get_data( mAdjTag, iter, 1, &adj ) == MB_SUCCESS && adj != NULL )
 {
 adj->push_back( entity );
 }
 // create a new vector and add it
 else
 {
 adj = new std::vector< EntityHandle >;
 adj->push_back( entity );
 result = thisMB->tag_set_data( mAdjTag, iter, 1, &adj );MB_CHK_ERR( result );
 }
 
 return MB_SUCCESS;
}
 
void Skinner::add_adjacency( EntityHandle entity, const EntityHandle* nodes, const int num_nodes )
{
 std::vector< EntityHandle >* adj = NULL;
 const EntityHandle* iter = std::min_element( nodes, nodes + num_nodes );
 
 if( iter == nodes + num_nodes ) return;
 
 // should not be setting adjacency lists in ho-nodes
 assert( TYPE_FROM_HANDLE( entity ) == MBPOLYGON ||
 num_nodes == CN::VerticesPerEntity( TYPE_FROM_HANDLE( entity ) ) );
 
 // add this entity to the node
 if( thisMB->tag_get_data( mAdjTag, iter, 1, &adj ) == MB_SUCCESS && adj != NULL )
 {
 adj->push_back( entity );
 }
 // create a new vector and add it
 else
 {
 adj = new std::vector< EntityHandle >;
 adj->push_back( entity );
 thisMB->tag_set_data( mAdjTag, iter, 1, &adj );
 }
}
 
ErrorCode Skinner::find_geometric_skin( const EntityHandle meshset, Range& forward_target_entities )
{
 // attempts to find whole model skin, using geom topo sets first then
 // normal find_skin function
 bool debug = true;
 
 // look for geom topo sets
 Tag geom_tag;
 ErrorCode result =
 thisMB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
 
 if( MB_SUCCESS != result ) return result;
 
 // get face sets (dimension = 2)
 Range face_sets;
 int two = 2;
 const void* two_ptr = &two;
 result = thisMB->get_entities_by_type_and_tag( meshset, MBENTITYSET, &geom_tag, &two_ptr, 1, face_sets );
 
 Range::iterator it;
 if( MB_SUCCESS != result )
 return result;
 else if( face_sets.empty() )
 return MB_ENTITY_NOT_FOUND;
 
 // ok, we have face sets; use those to determine skin
 Range skin_sets;
 if( debug ) std::cout << "Found " << face_sets.size() << " face sets total..." << std::endl;
 
 for( it = face_sets.begin(); it != face_sets.end(); ++it )
 {
 int num_parents;
 result = thisMB->num_parent_meshsets( *it, &num_parents );
 if( MB_SUCCESS != result )
 return result;
 else if( num_parents == 1 )
 skin_sets.insert( *it );
 }
 
 if( debug ) std::cout << "Found " << skin_sets.size() << " 1-parent face sets..." << std::endl;
 
 if( skin_sets.empty() ) return MB_FAILURE;
 
 // ok, we have the shell; gather up the elements, putting them all in forward for now
 for( it = skin_sets.begin(); it != skin_sets.end(); ++it )
 {
 result = thisMB->get_entities_by_handle( *it, forward_target_entities, true );
 if( MB_SUCCESS != result ) return result;
 }
 
 return result;
}
 
ErrorCode Skinner::find_skin( const EntityHandle meshset,
 const Range& source_entities,
 bool get_vertices,
 Range& output_handles,
 Range* output_reverse_handles,
 bool create_vert_elem_adjs,
 bool create_skin_elements,
 bool look_for_scd )
{
 if( source_entities.empty() ) return MB_SUCCESS;
 
 if( look_for_scd )
 {
 ErrorCode rval = find_skin_scd( source_entities, get_vertices, output_handles, create_skin_elements );
 // if it returns success, it's all scd, and we don't need to do anything more
 if( MB_SUCCESS == rval ) return rval;
 }
 
 Core* this_core = dynamic_cast< Core* >( thisMB );
 if( this_core && create_vert_elem_adjs && !this_core->a_entity_factory()->vert_elem_adjacencies() )
 this_core->a_entity_factory()->create_vert_elem_adjacencies();
 
 return find_skin_vertices( meshset, source_entities, get_vertices ? &output_handles : 0,
 get_vertices ? 0 : &output_handles, output_reverse_handles, create_skin_elements );
}
 
ErrorCode Skinner::find_skin_scd( const Range& source_entities,
 bool get_vertices,
 Range& output_handles,
 bool create_skin_elements )
{
 // get the scd interface and check if it's been initialized
 ScdInterface* scdi = NULL;
 ErrorCode rval = thisMB->query_interface( scdi );
 if( !scdi ) return MB_FAILURE;
 
 // ok, there's scd mesh; see if the entities passed in are all in a scd box
 // a box needs to be wholly included in entities for this to work
 std::vector< ScdBox* > boxes, myboxes;
 Range myrange;
 rval = scdi->find_boxes( boxes );
 if( MB_SUCCESS != rval ) return rval;
 for( std::vector< ScdBox* >::iterator bit = boxes.begin(); bit != boxes.end(); ++bit )
 {
 Range belems( ( *bit )->start_element(), ( *bit )->start_element() + ( *bit )->num_elements() - 1 );
 if( source_entities.contains( belems ) )
 {
 myboxes.push_back( *bit );
 myrange.merge( belems );
 }
 }
 if( myboxes.empty() || myrange.size() != source_entities.size() ) return MB_FAILURE;
 
 // ok, we're all structured; get the skin for each box
 for( std::vector< ScdBox* >::iterator bit = boxes.begin(); bit != boxes.end(); ++bit )
 {
 rval = skin_box( *bit, get_vertices, output_handles, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode Skinner::skin_box( ScdBox* box, bool get_vertices, Range& output_handles, bool create_skin_elements )
{
 HomCoord bmin = box->box_min(), bmax = box->box_max();
 
 // don't support 1d boxes
 if( bmin.j() == bmax.j() && bmin.k() == bmax.k() ) return MB_FAILURE;
 
 int dim = ( bmin.k() == bmax.k() ? 1 : 2 );
 
 ErrorCode rval;
 EntityHandle ent;
 
 // i=min
 for( int k = bmin.k(); k < bmax.k(); k++ )
 {
 for( int j = bmin.j(); j < bmax.j(); j++ )
 {
 ent = 0;
 rval = box->get_adj_edge_or_face( dim, bmin.i(), j, k, 0, ent, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 if( ent ) output_handles.insert( ent );
 }
 }
 // i=max
 for( int k = bmin.k(); k < bmax.k(); k++ )
 {
 for( int j = bmin.j(); j < bmax.j(); j++ )
 {
 ent = 0;
 rval = box->get_adj_edge_or_face( dim, bmax.i(), j, k, 0, ent, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 if( ent ) output_handles.insert( ent );
 }
 }
 // j=min
 for( int k = bmin.k(); k < bmax.k(); k++ )
 {
 for( int i = bmin.i(); i < bmax.i(); i++ )
 {
 ent = 0;
 rval = box->get_adj_edge_or_face( dim, i, bmin.j(), k, 1, ent, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 if( ent ) output_handles.insert( ent );
 }
 }
 // j=max
 for( int k = bmin.k(); k < bmax.k(); k++ )
 {
 for( int i = bmin.i(); i < bmax.i(); i++ )
 {
 ent = 0;
 rval = box->get_adj_edge_or_face( dim, i, bmax.j(), k, 1, ent, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 if( ent ) output_handles.insert( ent );
 }
 }
 // k=min
 for( int j = bmin.j(); j < bmax.j(); j++ )
 {
 for( int i = bmin.i(); i < bmax.i(); i++ )
 {
 ent = 0;
 rval = box->get_adj_edge_or_face( dim, i, j, bmin.k(), 2, ent, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 if( ent ) output_handles.insert( ent );
 }
 }
 // k=max
 for( int j = bmin.j(); j < bmax.j(); j++ )
 {
 for( int i = bmin.i(); i < bmax.i(); i++ )
 {
 ent = 0;
 rval = box->get_adj_edge_or_face( dim, i, j, bmax.k(), 2, ent, create_skin_elements );
 if( MB_SUCCESS != rval ) return rval;
 if( ent ) output_handles.insert( ent );
 }
 }
 
 if( get_vertices )
 {
 Range verts;
 rval = thisMB->get_adjacencies( output_handles, 0, true, verts, Interface::UNION );
 if( MB_SUCCESS != rval ) return rval;
 output_handles.merge( verts );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode Skinner::find_skin_noadj(const Range &source_entities,
 Range &forward_target_entities,
 Range &reverse_target_entities/*,
 bool create_vert_elem_adjs*/)
{
 if( source_entities.empty() ) return MB_FAILURE;
 
 // get our working dimensions
 EntityType type = thisMB->type_from_handle( *( source_entities.begin() ) );
 const int source_dim = CN::Dimension( type );
 mTargetDim = source_dim - 1;
 
 // make sure we can handle the working dimensions
 if( mTargetDim < 0 || source_dim > 3 ) return MB_FAILURE;
 
 initialize();
 
 Range::const_iterator iter, end_iter;
 end_iter = source_entities.end();
 const EntityHandle* conn;
 EntityHandle match;
 
 direction direct;
 ErrorCode result;
 // assume we'll never have more than 32 vertices on a facet (checked
 // with assert later)
 EntityHandle sub_conn[32];
 std::vector< EntityHandle > tmp_conn_vec;
 int num_nodes, num_sub_nodes, num_sides;
 int sub_indices[32]; // Also, assume that no polygon has more than 32 nodes
 // we could increase that, but we will not display it right in visit moab h5m , anyway
 EntityType sub_type;
 
 // for each source entity
 for( iter = source_entities.begin(); iter != end_iter; ++iter )
 {
 // get the connectivity of this entity
 int actual_num_nodes_polygon = 0;
 result = thisMB->get_connectivity( *iter, conn, num_nodes, false, &tmp_conn_vec );
 if( MB_SUCCESS != result ) return result;
 
 type = thisMB->type_from_handle( *iter );
 Range::iterator seek_iter;
 
 // treat separately polygons (also, polyhedra will need special handling)
 if( MBPOLYGON == type )
 {
 // treat padded polygons, if existing; count backwards, see how many of the last nodes
 // are repeated assume connectivity is fine, otherwise we could be in trouble
 actual_num_nodes_polygon = num_nodes;
 while( actual_num_nodes_polygon >= 3 &&
 conn[actual_num_nodes_polygon - 1] == conn[actual_num_nodes_polygon - 2] )
 actual_num_nodes_polygon--;
 num_sides = actual_num_nodes_polygon;
 sub_type = MBEDGE;
 num_sub_nodes = 2;
 }
 else // get connectivity of each n-1 dimension entity
 num_sides = CN::NumSubEntities( type, mTargetDim );
 for( int i = 0; i < num_sides; i++ )
 {
 if( MBPOLYGON == type )
 {
 sub_conn[0] = conn[i];
 sub_conn[1] = conn[i + 1];
 if( i + 1 == actual_num_nodes_polygon ) sub_conn[1] = conn[0];
 }
 else
 {
 CN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, sub_type, num_sub_nodes, sub_indices );
 assert( (size_t)num_sub_nodes <= sizeof( sub_indices ) / sizeof( sub_indices[0] ) );
 for( int j = 0; j < num_sub_nodes; j++ )
 sub_conn[j] = conn[sub_indices[j]];
 }
 
 // see if we can match this connectivity with
 // an existing entity
 find_match( sub_type, sub_conn, num_sub_nodes, match, direct );
 
 // if there is no match, create a new entity
 if( match == 0 )
 {
 EntityHandle tmphndl = 0;
 int indices[MAX_SUB_ENTITY_VERTICES];
 EntityType new_type;
 int num_new_nodes;
 if( MBPOLYGON == type )
 {
 new_type = MBEDGE;
 num_new_nodes = 2;
 }
 else
 {
 CN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
 for( int j = 0; j < num_new_nodes; j++ )
 sub_conn[j] = conn[indices[j]];
 }
 result = thisMB->create_element( new_type, sub_conn, num_new_nodes, tmphndl );
 assert( MB_SUCCESS == result );
 add_adjacency( tmphndl, sub_conn, CN::VerticesPerEntity( new_type ) );
 forward_target_entities.insert( tmphndl );
 }
 // if there is a match, delete the matching entity
 // if we can.
 else
 {
 if( ( seek_iter = forward_target_entities.find( match ) ) != forward_target_entities.end() )
 {
 forward_target_entities.erase( seek_iter );
 remove_adjacency( match );
 if( /*!use_adjs &&*/ entity_deletable( match ) )
 {
 result = thisMB->delete_entities( &match, 1 );
 assert( MB_SUCCESS == result );
 }
 }
 else if( ( seek_iter = reverse_target_entities.find( match ) ) != reverse_target_entities.end() )
 {
 reverse_target_entities.erase( seek_iter );
 remove_adjacency( match );
 if( /*!use_adjs &&*/ entity_deletable( match ) )
 {
 result = thisMB->delete_entities( &match, 1 );
 assert( MB_SUCCESS == result );
 }
 }
 else
 {
 if( direct == FORWARD )
 {
 forward_target_entities.insert( match );
 }
 else
 {
 reverse_target_entities.insert( match );
 }
 }
 }
 }
 }
 
 deinitialize();
 
 return MB_SUCCESS;
}
 
void Skinner::find_match( EntityType type,
 const EntityHandle* conn,
 const int num_nodes,
 EntityHandle& match,
 Skinner::direction& direct )
{
 match = 0;
 
 if( type == MBVERTEX )
 {
 match = *conn;
 direct = FORWARD;
 return;
 }
 
 const EntityHandle* iter = std::min_element( conn, conn + num_nodes );
 
 std::vector< EntityHandle >* adj = NULL;
 
 ErrorCode result = thisMB->tag_get_data( mAdjTag, iter, 1, &adj );
 if( result == MB_FAILURE || adj == NULL )
 {
 return;
 }
 
 std::vector< EntityHandle >::iterator jter, end_jter;
 end_jter = adj->end();
 
 const EntityHandle* tmp;
 int num_verts;
 
 for( jter = adj->begin(); jter != end_jter; ++jter )
 {
 EntityType tmp_type;
 tmp_type = thisMB->type_from_handle( *jter );
 
 if( type != tmp_type ) continue;
 
 result = thisMB->get_connectivity( *jter, tmp, num_verts, false );
 assert( MB_SUCCESS == result && num_verts >= CN::VerticesPerEntity( type ) );
 // FIXME: connectivity_match appears to work only for linear elements,
 // so ignore higher-order nodes.
 if( connectivity_match( conn, tmp, CN::VerticesPerEntity( type ), direct ) )
 {
 match = *jter;
 break;
 }
 }
}
 
bool Skinner::connectivity_match( const EntityHandle* conn1,
 const EntityHandle* conn2,
 const int num_verts,
 Skinner::direction& direct )
{
 const EntityHandle* iter = std::find( conn2, conn2 + num_verts, conn1[0] );
 if( iter == conn2 + num_verts ) return false;
 
 bool they_match = true;
 
 int i;
 unsigned int j = iter - conn2;
 
 // first compare forward
 for( i = 1; i < num_verts; ++i )
 {
 if( conn1[i] != conn2[( j + i ) % num_verts] )
 {
 they_match = false;
 break;
 }
 }
 
 if( they_match == true )
 {
 // need to check for reversed edges here
 direct = ( num_verts == 2 && j ) ? REVERSE : FORWARD;
 return true;
 }
 
 they_match = true;
 
 // then compare reverse
 j += num_verts;
 for( i = 1; i < num_verts; )
 {
 if( conn1[i] != conn2[( j - i ) % num_verts] )
 {
 they_match = false;
 break;
 }
 ++i;
 }
 if( they_match )
 {
 direct = REVERSE;
 }
 return they_match;
}
 
ErrorCode Skinner::remove_adjacency( EntityHandle entity )
{
 std::vector< EntityHandle > nodes, *adj = NULL;
 ErrorCode result = thisMB->get_connectivity( &entity, 1, nodes );MB_CHK_ERR( result );
 std::vector< EntityHandle >::iterator iter = std::min_element( nodes.begin(), nodes.end() );
 
 if( iter == nodes.end() ) return MB_FAILURE;
 
 // remove this entity from the node
 if( thisMB->tag_get_data( mAdjTag, &( *iter ), 1, &adj ) == MB_SUCCESS && adj != NULL )
 {
 iter = std::find( adj->begin(), adj->end(), entity );
 if( iter != adj->end() ) adj->erase( iter );
 }
 
 return result;
}
 
bool Skinner::entity_deletable( EntityHandle entity )
{
 unsigned char deletable = 0;
 ErrorCode result = thisMB->tag_get_data( mDeletableMBTag, &entity, 1, &deletable );
 assert( MB_SUCCESS == result );
 if( MB_SUCCESS == result && deletable == 1 ) return false;
 return true;
}
 
ErrorCode Skinner::classify_2d_boundary( const Range& boundary,
 const Range& bar_elements,
 EntityHandle boundary_edges,
 EntityHandle inferred_edges,
 EntityHandle non_manifold_edges,
 EntityHandle other_edges,
 int& number_boundary_nodes )
{
 Range bedges, iedges, nmedges, oedges;
 ErrorCode result =
 classify_2d_boundary( boundary, bar_elements, bedges, iedges, nmedges, oedges, number_boundary_nodes );MB_CHK_ERR( result );
 
 // now set the input meshsets to the output ranges
 result = thisMB->clear_meshset( &boundary_edges, 1 );MB_CHK_ERR( result );
 result = thisMB->add_entities( boundary_edges, bedges );MB_CHK_ERR( result );
 
 result = thisMB->clear_meshset( &inferred_edges, 1 );MB_CHK_ERR( result );
 result = thisMB->add_entities( inferred_edges, iedges );MB_CHK_ERR( result );
 
 result = thisMB->clear_meshset( &non_manifold_edges, 1 );MB_CHK_ERR( result );
 result = thisMB->add_entities( non_manifold_edges, nmedges );MB_CHK_ERR( result );
 
 result = thisMB->clear_meshset( &other_edges, 1 );MB_CHK_ERR( result );
 result = thisMB->add_entities( other_edges, oedges );MB_CHK_ERR( result );
 
 return MB_SUCCESS;
}
 
ErrorCode Skinner::classify_2d_boundary( const Range& boundary,
 const Range& bar_elements,
 Range& boundary_edges,
 Range& inferred_edges,
 Range& non_manifold_edges,
 Range& other_edges,
 int& number_boundary_nodes )
{
 
 // clear out the edge lists
 
 boundary_edges.clear();
 inferred_edges.clear();
 non_manifold_edges.clear();
 other_edges.clear();
 
 number_boundary_nodes = 0;
 
 // make sure we have something to work with
 if( boundary.empty() )
 {
 return MB_FAILURE;
 }
 
 // get our working dimensions
 EntityType type = thisMB->type_from_handle( *( boundary.begin() ) );
 const int source_dim = CN::Dimension( type );
 
 // make sure we can handle the working dimensions
 if( source_dim != 2 )
 {
 return MB_FAILURE;
 }
 mTargetDim = source_dim - 1;
 
 // initialize
 initialize();
 
 // additional initialization for this routine
 // define a tag for MBEDGE which counts the occurrences of the edge below
 // default should be 0 for existing edges, if any
 
 Tag count_tag;
 int default_count = 0;
 ErrorCode result =
 thisMB->tag_get_handle( 0, 1, MB_TYPE_INTEGER, count_tag, MB_TAG_DENSE | MB_TAG_CREAT, &default_count );MB_CHK_ERR( result );
 
 Range::const_iterator iter, end_iter;
 end_iter = boundary.end();
 
 std::vector< EntityHandle > conn;
 EntityHandle sub_conn[2];
 EntityHandle match;
 
 Range edge_list;
 Range boundary_nodes;
 Skinner::direction direct;
 
 EntityType sub_type;
 int num_edge, num_sub_ent_vert;
 const short* edge_verts;
 
 // now, process each entity in the boundary
 
 for( iter = boundary.begin(); iter != end_iter; ++iter )
 {
 // get the connectivity of this entity
 conn.clear();
 result = thisMB->get_connectivity( &( *iter ), 1, conn, false );
 assert( MB_SUCCESS == result );
 
 // add node handles to boundary_node range
 std::copy( conn.begin(), conn.begin() + CN::VerticesPerEntity( type ), range_inserter( boundary_nodes ) );
 
 type = thisMB->type_from_handle( *iter );
 
 // get connectivity of each n-1 dimension entity (edge in this case)
 const struct CN::ConnMap* conn_map = &( CN::mConnectivityMap[type][0] );
 num_edge = CN::NumSubEntities( type, 1 );
 for( int i = 0; i < num_edge; i++ )
 {
 edge_verts = CN::SubEntityVertexIndices( type, 1, i, sub_type, num_sub_ent_vert );
 assert( sub_type == MBEDGE && num_sub_ent_vert == 2 );
 sub_conn[0] = conn[edge_verts[0]];
 sub_conn[1] = conn[edge_verts[1]];
 int num_sub_nodes = conn_map->num_corners_per_sub_element[i];
 
 // see if we can match this connectivity with
 // an existing entity
 find_match( MBEDGE, sub_conn, num_sub_nodes, match, direct );
 
 // if there is no match, create a new entity
 if( match == 0 )
 {
 EntityHandle tmphndl = 0;
 int indices[MAX_SUB_ENTITY_VERTICES];
 EntityType new_type;
 int num_new_nodes;
 CN::SubEntityNodeIndices( type, conn.size(), 1, i, new_type, num_new_nodes, indices );
 for( int j = 0; j < num_new_nodes; j++ )
 sub_conn[j] = conn[indices[j]];
 
 result = thisMB->create_element( new_type, sub_conn, num_new_nodes, tmphndl );
 assert( MB_SUCCESS == result );
 add_adjacency( tmphndl, sub_conn, num_sub_nodes );
 // target_entities.insert(tmphndl);
 edge_list.insert( tmphndl );
 int count;
 result = thisMB->tag_get_data( count_tag, &tmphndl, 1, &count );
 assert( MB_SUCCESS == result );
 count++;
 result = thisMB->tag_set_data( count_tag, &tmphndl, 1, &count );
 assert( MB_SUCCESS == result );
 }
 else
 {
 // We found a match, we must increment the count on the match
 int count;
 result = thisMB->tag_get_data( count_tag, &match, 1, &count );
 assert( MB_SUCCESS == result );
 count++;
 result = thisMB->tag_set_data( count_tag, &match, 1, &count );
 assert( MB_SUCCESS == result );
 
 // if the entity is not deletable, it was pre-existing in
 // the database. We therefore may need to add it to the
 // edge_list. Since it will not hurt the range, we add
 // whether it was added before or not
 if( !entity_deletable( match ) )
 {
 edge_list.insert( match );
 }
 }
 }
 }
 
 // Any bar elements in the model should be classified separately
 // If the element is in the skin edge_list, then it should be put in
 // the non-manifold edge list. Edges not in the edge_list are stand-alone
 // bars, and we make them simply boundary elements
 
 if( !bar_elements.empty() )
 {
 Range::iterator bar_iter;
 for( iter = bar_elements.begin(); iter != bar_elements.end(); ++iter )
 {
 EntityHandle handle = *iter;
 bar_iter = edge_list.find( handle );
 if( bar_iter != edge_list.end() )
 {
 // it is in the list, erase it and put in non-manifold list
 edge_list.erase( bar_iter );
 non_manifold_edges.insert( handle );
 }
 else
 {
 // not in the edge list, make it a boundary edge
 boundary_edges.insert( handle );
 }
 }
 }
 
 // now all edges should be classified. Go through the edge_list,
 // and put all in the appropriate lists
 
 Range::iterator edge_iter, edge_end_iter;
 edge_end_iter = edge_list.end();
 int count;
 for( edge_iter = edge_list.begin(); edge_iter != edge_end_iter; ++edge_iter )
 {
 // check the count_tag
 result = thisMB->tag_get_data( count_tag, &( *edge_iter ), 1, &count );
 assert( MB_SUCCESS == result );
 if( count == 1 )
 {
 boundary_edges.insert( *edge_iter );
 }
 else if( count == 2 )
 {
 other_edges.insert( *edge_iter );
 }
 else
 {
 non_manifold_edges.insert( *edge_iter );
 }
 }
 
 // find the inferred edges from the other_edge_list
 
 double min_angle_degrees = 20.0;
 find_inferred_edges( const_cast< Range& >( boundary ), other_edges, inferred_edges, min_angle_degrees );
 
 // we now want to remove the inferred_edges from the other_edges
 
 Range temp_range;
 
 std::set_difference( other_edges.begin(), other_edges.end(), inferred_edges.begin(), inferred_edges.end(),
 range_inserter( temp_range ), std::less< EntityHandle >() );
 
 other_edges = temp_range;
 
 // get rid of count tag and deinitialize
 
 result = thisMB->tag_delete( count_tag );
 assert( MB_SUCCESS == result );
 deinitialize();
 
 // set the node count
 number_boundary_nodes = boundary_nodes.size();
 
 return MB_SUCCESS;
}
 
void Skinner::find_inferred_edges( Range& skin_boundary,
 Range& candidate_edges,
 Range& inferred_edges,
 double reference_angle_degrees )
{
 
 // mark all the entities in the skin boundary
 Tag mark_tag;
 ErrorCode result = thisMB->tag_get_handle( 0, 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT );
 assert( MB_SUCCESS == result );
 unsigned char bit = true;
 result = thisMB->tag_clear_data( mark_tag, skin_boundary, &bit );
 assert( MB_SUCCESS == result );
 
 // find the cosine of the reference angle
 
 double reference_cosine = cos( reference_angle_degrees * SKINNER_PI / 180.0 );
 
 // check all candidate edges for an angle greater than the minimum
 
 Range::iterator iter, end_iter = candidate_edges.end();
 std::vector< EntityHandle > adjacencies;
 std::vector< EntityHandle >::iterator adj_iter;
 EntityHandle face[2];
 
 for( iter = candidate_edges.begin(); iter != end_iter; ++iter )
 {
 
 // get the 2D elements connected to this edge
 adjacencies.clear();
 result = thisMB->get_adjacencies( &( *iter ), 1, 2, false, adjacencies );
 if( MB_SUCCESS != result ) continue;
 
 // there should be exactly two, that is why the edge is classified as nonBoundary
 // and manifold
 
 int faces_found = 0;
 for( adj_iter = adjacencies.begin(); adj_iter != adjacencies.end() && faces_found < 2; ++adj_iter )
 {
 // we need to find two of these which are in the skin
 unsigned char is_marked = 0;
 result = thisMB->tag_get_data( mark_tag, &( *adj_iter ), 1, &is_marked );
 assert( MB_SUCCESS == result );
 if( is_marked )
 {
 face[faces_found] = *adj_iter;
 faces_found++;
 }
 }
 
 // assert(faces_found == 2 || faces_found == 0);
 if( 2 != faces_found ) continue;
 
 // see if the two entities have a sufficient angle
 
 if( has_larger_angle( face[0], face[1], reference_cosine ) )
 {
 inferred_edges.insert( *iter );
 }
 }
 
 result = thisMB->tag_delete( mark_tag );
 assert( MB_SUCCESS == result );
}
 
bool Skinner::has_larger_angle( EntityHandle& entity1, EntityHandle& entity2, double reference_angle_cosine )
{
 // compare normals to get angle. We assume that the surface quads
 // which we test here will be approximately planar
 
 double norm[2][3];
 Util::normal( thisMB, entity1, norm[0][0], norm[0][1], norm[0][2] );
 Util::normal( thisMB, entity2, norm[1][0], norm[1][1], norm[1][2] );
 
 double cosine = norm[0][0] * norm[1][0] + norm[0][1] * norm[1][1] + norm[0][2] * norm[1][2];
 
 if( cosine < reference_angle_cosine )
 {
 return true;
 }
 
 return false;
}
 
// get skin entities of prescribed dimension
ErrorCode Skinner::find_skin( const EntityHandle this_set,
 const Range& entities,
 int dim,
 Range& skin_entities,
 bool create_vert_elem_adjs,
 bool create_skin_elements )
{
 Range tmp_skin;
 ErrorCode result =
 find_skin( this_set, entities, ( dim == 0 ), tmp_skin, 0, create_vert_elem_adjs, create_skin_elements );
 if( MB_SUCCESS != result || tmp_skin.empty() ) return result;
 
 if( tmp_skin.all_of_dimension( dim ) )
 {
 if( skin_entities.empty() )
 skin_entities.swap( tmp_skin );
 else
 skin_entities.merge( tmp_skin );
 }
 else
 {
 result = thisMB->get_adjacencies( tmp_skin, dim, create_skin_elements, skin_entities, Interface::UNION );MB_CHK_ERR( result );
 if( this_set ) result = thisMB->add_entities( this_set, skin_entities );
 }
 
 return result;
}
 
ErrorCode Skinner::find_skin_vertices( const EntityHandle this_set,
 const Range& entities,
 Range* skin_verts,
 Range* skin_elems,
 Range* skin_rev_elems,
 bool create_skin_elems,
 bool corners_only )
{
 ErrorCode rval;
 if( entities.empty() ) return MB_SUCCESS;
 
 const int dim = CN::Dimension( TYPE_FROM_HANDLE( entities.front() ) );
 if( dim < 1 || dim > 3 || !entities.all_of_dimension( dim ) ) return MB_TYPE_OUT_OF_RANGE;
 
 // are we skinning all entities
 size_t count = entities.size();
 int num_total;
 rval = thisMB->get_number_entities_by_dimension( this_set, dim, num_total );
 if( MB_SUCCESS != rval ) return rval;
 bool all = ( count == (size_t)num_total );
 
 // Create a bit tag for fast intersection with input entities range.
 // If we're skinning all the entities in the mesh, we really don't
 // need the tag. To save memory, just create it with a default value
 // of one and don't set it. That way MOAB will return 1 for all
 // entities.
 Tag tag;
 char bit = all ? 1 : 0;
 rval = thisMB->tag_get_handle( NULL, 1, MB_TYPE_BIT, tag, MB_TAG_CREAT, &bit );
 if( MB_SUCCESS != rval ) return rval;
 
 // tag all entities in input range
 if( !all )
 {
 std::vector< unsigned char > vect( count, 1 );
 rval = thisMB->tag_set_data( tag, entities, &vect[0] );
 if( MB_SUCCESS != rval )
 {
 thisMB->tag_delete( tag );
 return rval;
 }
 }
 
 switch( dim )
 {
 case 1:
 if( skin_verts )
 rval = find_skin_vertices_1D( tag, entities, *skin_verts );
 else if( skin_elems )
 rval = find_skin_vertices_1D( tag, entities, *skin_elems );
 else
 rval = MB_SUCCESS;
 break;
 case 2:
 rval = find_skin_vertices_2D( this_set, tag, entities, skin_verts, skin_elems, skin_rev_elems,
 create_skin_elems, corners_only );
 break;
 case 3:
 rval = find_skin_vertices_3D( this_set, tag, entities, skin_verts, skin_elems, skin_rev_elems,
 create_skin_elems, corners_only );
 break;
 default:
 rval = MB_TYPE_OUT_OF_RANGE;
 break;
 }
 
 thisMB->tag_delete( tag );
 return rval;
}
 
ErrorCode Skinner::find_skin_vertices_1D( Tag tag, const Range& edges, Range& skin_verts )
{
 // This rather simple algorithm is provided for completeness
 // (not sure how often one really wants the 'skin' of a chain
 // or tangle of edges.)
 //
 // A vertex is on the skin of the edges if it is contained in exactly
 // one of the edges *in the input range*.
 //
 // This function expects the caller to have tagged all edges in the
 // input range with a value of one for the passed bit tag, and all
 // other edges with a value of zero. This allows us to do a faster
 // intersection with the input range and the edges adjacent to a vertex.
 
 ErrorCode rval;
 Range::iterator hint = skin_verts.begin();
 
 // All input entities must be edges.
 if( !edges.all_of_dimension( 1 ) ) return MB_TYPE_OUT_OF_RANGE;
 
 // get all the vertices
 Range verts;
 rval = thisMB->get_connectivity( edges, verts, true );
 if( MB_SUCCESS != rval ) return rval;
 
 // Test how many edges each input vertex is adjacent to.
 std::vector< char > tag_vals;
 std::vector< EntityHandle > adj;
 int n;
 for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
 {
 // get edges adjacent to vertex
 adj.clear();
 rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
 if( MB_SUCCESS != rval ) return rval;
 if( adj.empty() ) continue;
 
 // Intersect adjacent edges with the input list of edges
 tag_vals.resize( adj.size() );
 rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
 if( MB_SUCCESS != rval ) return rval;
#ifdef MOAB_OLD_STD_COUNT
 n = 0;
 std::count( tag_vals.begin(), tag_vals.end(), '\001' );
#else
 n = std::count( tag_vals.begin(), tag_vals.end(), '\001' );
#endif
 // If adjacent to only one input edge, then vertex is on skin
 if( n == 1 )
 {
 hint = skin_verts.insert( hint, *it );
 }
 }
 
 return MB_SUCCESS;
}
 
// A Container for storing a list of element sides adjacent
// to a vertex. The template parameter is the number of
// corners for the side.
template < unsigned CORNERS >
class AdjSides
{
 public:
 /**
 * This struct is used to for a reduced representation of element
 * "sides" adjacent to a give vertex. As such, it
 * a) does not store the initial vertex that all sides are adjacent to
 * b) orders the remaining vertices in a specific way for fast comparison.
 *
 * For edge elements, only the opposite vertex is stored.
 * For triangle elements, only the other two vertices are stored,
 * and they are stored with the smaller of those two handles first.
 * For quad elements, only the other three vertices are stored.
 * They are stored such that the vertex opposite the implicit (not
 * stored) vertex is always in slot 1. The other two vertices
 * (in slots 0 and 2) are ordered such that the handle of the one in
 * slot 0 is smaller than the handle in slot 2.
 *
 * For each side, the adj_elem field is used to store the element that
 * it is a side of as long as the element is considered to be on the skin.
 * The adj_elem field is cleared (set to zero) to indicate that this
 * side is no longer considered to be on the skin (and is the side of
 * more than one element.)
 */
 struct Side
 {
 EntityHandle handles[CORNERS - 1]; //!< side vertices, except for implicit one
 EntityHandle adj_elem; //!< element that this is a side of, or zero
 bool skin() const
 {
 return 0 != adj_elem;
 }
 
 /** construct from connectivity of side
 *\param array The connectivity of the element side.
 *\param idx The index of the implicit vertex (contained
 * in all sides in the list.)
 *\param adj The element that this is a side of.
 */
 Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short ) : adj_elem( adj )
 {
 switch( CORNERS )
 {
 case 3:
 handles[1] = array[( idx + 2 ) % CORNERS];
 // fall through
 case 2:
 if( 3 == CORNERS ) handles[1] = array[( idx + 2 ) % CORNERS];
 if( 2 <= CORNERS ) handles[0] = array[( idx + 1 ) % CORNERS];
 break;
 default:
 assert( false );
 break;
 }
 if( CORNERS == 3 && handles[1] > handles[0] ) std::swap( handles[0], handles[1] );
 }
 
 /** construct from connectivity of parent element
 *\param array The connectivity of the parent element
 *\param idx The index of the implicit vertex (contained
 * in all sides in the list.) This is an index
 * into 'indices', not 'array'.
 *\param adj The element that this is a side of.
 *\param indices The indices into 'array' at which the vertices
 * representing the side occur.
 */
 Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short, const short* indices )
 : adj_elem( adj )
 {
 switch( CORNERS )
 {
 case 3:
 handles[1] = array[indices[( idx + 2 ) % CORNERS]];
 // fall through
 case 2:
 if( 3 == CORNERS ) handles[1] = array[indices[( idx + 2 ) % CORNERS]];
 if( 2 <= CORNERS ) handles[0] = array[indices[( idx + 1 ) % CORNERS]];
 break;
 default:
 assert( false );
 break;
 }
 if( CORNERS == 3 && handles[1] > handles[0] ) std::swap( handles[0], handles[1] );
 }
 
 // Compare two side instances. Relies in the ordering of
 // vertex handles as described above.
 bool operator==( const Side& other ) const
 {
 switch( CORNERS )
 {
 case 3:
 return handles[0] == other.handles[0] && handles[1] == other.handles[1];
 case 2:
 return handles[0] == other.handles[0];
 default:
 assert( false );
 return false;
 }
 }
 };
 
 private:
 std::vector< Side > data; //!< List of sides
 size_t skin_count; //!< Cached count of sides that are skin
 
 public:
 typedef typename std::vector< Side >::iterator iterator;
 typedef typename std::vector< Side >::const_iterator const_iterator;
 const_iterator begin() const
 {
 return data.begin();
 }
 const_iterator end() const
 {
 return data.end();
 }
 
 void clear()
 {
 data.clear();
 skin_count = 0;
 }
 bool empty() const
 {
 return data.empty();
 }
 
 AdjSides() : skin_count( 0 ) {}
 
 size_t num_skin() const
 {
 return skin_count;
 }
 
 /** \brief insert side, specifying side connectivity
 *
 * Either insert a new side, or if the side is already in the
 * list, mark it as not on the skin.
 *
 *\param handles The connectivity of the element side.
 *\param skip_idx The index of the implicit vertex (contained
 * in all sides in the list.)
 *\param adj_elem The element that this is a side of.
 *\param elem_side Which side of adj_elem are we storing
 * (CN side number.)
 */
 void insert( const EntityHandle* handles, int skip_idx, EntityHandle adj_elem, unsigned short elem_side )
 {
 Side side( handles, skip_idx, adj_elem, elem_side );
 iterator p = std::find( data.begin(), data.end(), side );
 if( p == data.end() )
 {
 data.push_back( side );
 ++skin_count; // not in list yet, so skin side (so far)
 }
 else if( p->adj_elem )
 {
 p->adj_elem = 0; // mark as not on skin
 --skin_count; // decrement cached count of skin elements
 }
 }
 
 /** \brief insert side, specifying list of indices into parent element
 * connectivity.
 *
 * Either insert a new side, or if the side is already in the
 * list, mark it as not on the skin.
 *
 *\param handles The connectivity of the parent element
 *\param skip_idx The index of the implicit vertex (contained
 * in all sides in the list.) This is an index
 * into 'indices', not 'handles'.
 *\param adj_elem The element that this is a side of (parent handle).
 *\param indices The indices into 'handles' at which the vertices
 * representing the side occur.
 *\param elem_side Which side of adj_elem are we storing
 * (CN side number.)
 */
 void insert( const EntityHandle* handles,
 int skip_idx,
 EntityHandle adj_elem,
 unsigned short elem_side,
 const short* indices )
 {
 Side side( handles, skip_idx, adj_elem, elem_side, indices );
 iterator p = std::find( data.begin(), data.end(), side );
 if( p == data.end() )
 {
 data.push_back( side );
 ++skin_count; // not in list yet, so skin side (so far)
 }
 else if( p->adj_elem )
 {
 p->adj_elem = 0; // mark as not on skin
 --skin_count; // decrement cached count of skin elements
 }
 }
 
 /**\brief Search list for a given side, and if found, mark as not skin.
 *
 *\param other Connectivity of side
 *\param skip_index Index in 'other' at which implicit vertex occurs.
 *\param elem_out If return value is true, the element that the side is a
 * side of. If return value is false, not set.
 *\return true if found and marked not-skin, false if not found.
 *
 * Given the connectivity of some existing element, check if it occurs
 * in the list. If it does, clear the "is skin" state of the side so
 * that we know that we don't need to later create the side element.
 */
 bool find_and_unmark( const EntityHandle* other, int skip_index, EntityHandle& elem_out )
 {
 Side s( other, skip_index, 0, 0 );
 iterator p = std::find( data.begin(), data.end(), s );
 if( p == data.end() || !p->adj_elem )
 return false;
 else
 {
 elem_out = p->adj_elem;
 p->adj_elem = 0; // clear "is skin" state for side
 --skin_count; // decrement cached count of skin sides
 return true;
 }
 }
};
 
/** construct from connectivity of side
 *\param array The connectivity of the element side.
 *\param idx The index of the implicit vertex (contained
 * in all sides in the list.)
 *\param adj The element that this is a side of.
 */
template <>
AdjSides< 4 >::Side::Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short ) : adj_elem( adj )
{
 const unsigned int CORNERS = 4;
 handles[2] = array[( idx + 3 ) % CORNERS];
 handles[1] = array[( idx + 2 ) % CORNERS];
 handles[0] = array[( idx + 1 ) % CORNERS];
 if( handles[2] > handles[0] ) std::swap( handles[0], handles[2] );
}
 
/** construct from connectivity of parent element
 *\param array The connectivity of the parent element
 *\param idx The index of the implicit vertex (contained
 * in all sides in the list.) This is an index
 * into 'indices', not 'array'.
 *\param adj The element that this is a side of.
 *\param indices The indices into 'array' at which the vertices
 * representing the side occur.
 */
template <>
AdjSides< 4 >::Side::Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short, const short* indices )
 : adj_elem( adj )
{
 const unsigned int CORNERS = 4;
 handles[2] = array[indices[( idx + 3 ) % CORNERS]];
 handles[1] = array[indices[( idx + 2 ) % CORNERS]];
 handles[0] = array[indices[( idx + 1 ) % CORNERS]];
 if( handles[2] > handles[0] ) std::swap( handles[0], handles[2] );
}
 
// Compare two side instances. Relies in the ordering of
// vertex handles as described above.
template <>
bool AdjSides< 4 >::Side::operator==( const Side& other ) const
{
 return handles[0] == other.handles[0] && handles[1] == other.handles[1] && handles[2] == other.handles[2];
}
 
// Utility function used by find_skin_vertices_2D and
// find_skin_vertices_3D to create elements representing
// the skin side of a higher-dimension element if one
// does not already exist.
//
// Some arguments may seem redundant, but they are used
// to create the correct order of element when the input
// element contains higher-order nodes.
//
// This function always creates elements that have a "forward"
// orientation with respect to the parent element (have
// nodes ordered the same as CN returns for the "side").
//
// elem - The higher-dimension element for which to create
// a lower-dim element representing the side.
// side_type - The EntityType of the desired side.
// side_conn - The connectivity of the new side.
ErrorCode Skinner::create_side( const EntityHandle this_set,
 EntityHandle elem,
 EntityType side_type,
 const EntityHandle* side_conn,
 EntityHandle& side_elem )
{
 const int max_side = 9;
 const EntityHandle* conn;
 int len, side_len, side, sense, offset, indices[max_side];
 ErrorCode rval;
 EntityType type = TYPE_FROM_HANDLE( elem ), tmp_type;
 const int ncorner = CN::VerticesPerEntity( side_type );
 const int d = CN::Dimension( side_type );
 std::vector< EntityHandle > storage;
 
 // Get the connectivity of the parent element
 rval = thisMB->get_connectivity( elem, conn, len, false, &storage );
 if( MB_SUCCESS != rval ) return rval;
 
 // treat separately MBPOLYGON; we want to create the edge in the
 // forward sense always ; so figure out the sense first, then get out
 if( MBPOLYGON == type && 1 == d && MBEDGE == side_type )
 {
 // first find the first vertex in the conn list
 int i = 0;
 for( i = 0; i < len; i++ )
 {
 if( conn[i] == side_conn[0] ) break;
 }
 if( len == i ) return MB_FAILURE; // not found, big error
 // now, what if the polygon is padded?
 // the previous index is fine always. but the next one could be trouble :(
 int prevIndex = ( i + len - 1 ) % len;
 int nextIndex = ( i + 1 ) % len;
 // if the next index actually point to the same node, as current, it means it is padded
 if( conn[nextIndex] == conn[i] )
 {
 // it really means we are at the end of proper nodes, the last nodes are repeated, so it
 // should be the first node
 nextIndex = 0; // this is the first node!
 }
 EntityHandle conn2[2] = { side_conn[0], side_conn[1] };
 if( conn[prevIndex] == side_conn[1] )
 {
 // reverse, so the edge will be forward
 conn2[0] = side_conn[1];
 conn2[1] = side_conn[0];
 }
 else if( conn[nextIndex] != side_conn[1] )
 return MB_FAILURE; // it is not adjacent to the polygon
 
 rval = thisMB->create_element( MBEDGE, conn2, 2, side_elem );MB_CHK_ERR( rval );
 if( this_set )
 {
 rval = thisMB->add_entities( this_set, &side_elem, 1 );MB_CHK_ERR( rval );
 }
 return MB_SUCCESS;
 }
 // Find which side we are creating and get indices of all nodes
 // (including higher-order, if any.)
 CN::SideNumber( type, conn, side_conn, ncorner, d, side, sense, offset );
 CN::SubEntityNodeIndices( type, len, d, side, tmp_type, side_len, indices );
 assert( side_len <= max_side );
 assert( side_type == tmp_type );
 
 // NOTE: re-create conn array even when no higher-order nodes
 // because we want it to always be forward with respect
 // to the side ordering.
 EntityHandle side_conn_full[max_side];
 for( int i = 0; i < side_len; ++i )
 side_conn_full[i] = conn[indices[i]];
 
 rval = thisMB->create_element( side_type, side_conn_full, side_len, side_elem );MB_CHK_ERR( rval );
 if( this_set )
 {
 rval = thisMB->add_entities( this_set, &side_elem, 1 );MB_CHK_ERR( rval );
 }
 return MB_SUCCESS;
 ;
}
 
// Test if an edge is reversed with respect CN's ordering
// for the "side" of a face.
bool Skinner::edge_reversed( EntityHandle face, const EntityHandle* edge_ends )
{
 const EntityHandle* conn;
 int len, idx;
 ErrorCode rval = thisMB->get_connectivity( face, conn, len, true );
 if( MB_SUCCESS != rval )
 {
 assert( false );
 return false;
 }
 idx = std::find( conn, conn + len, edge_ends[0] ) - conn;
 if( idx == len )
 {
 assert( false );
 return false;
 }
 return ( edge_ends[1] == conn[( idx + len - 1 ) % len] );
}
 
// Test if a 2D element representing the side or face of a
// volume element is reversed with respect to the CN node
// ordering for the corresponding region element side.
bool Skinner::face_reversed( EntityHandle region, const EntityHandle* face_corners, EntityType face_type )
{
 const EntityHandle* conn;
 int len, side, sense, offset;
 ErrorCode rval = thisMB->get_connectivity( region, conn, len, true );
 if( MB_SUCCESS != rval )
 {
 assert( false );
 return false;
 }
 short r = CN::SideNumber( TYPE_FROM_HANDLE( region ), conn, face_corners, CN::VerticesPerEntity( face_type ),
 CN::Dimension( face_type ), side, sense, offset );
 assert( 0 == r );
 return ( !r && sense == -1 );
}
 
ErrorCode Skinner::find_skin_vertices_2D( const EntityHandle this_set,
 Tag tag,
 const Range& faces,
 Range* skin_verts,
 Range* skin_edges,
 Range* reversed_edges,
 bool create_edges,
 bool corners_only )
{
 // This function iterates over all the vertices contained in the
 // input face list. For each such vertex, it then iterates over
 // all of the sides of the face elements adjacent to the vertex.
 // If an adjacent side is the side of only one of the input
 // faces, then that side is on the skin.
 //
 // This algorithm will visit each skin vertex exactly once. It
 // will visit each skin side once for each vertex in the side.
 //
 // This function expects the caller to have created the passed bit
 // tag and set it to one only for the faces in the passed range. This
 // tag is used to do a fast intersection of the faces adjacent to a
 // vertex with the faces in the input range (discard any for which the
 // tag is not set to one.)
 
 ErrorCode rval;
 std::vector< EntityHandle >::iterator i, j;
 Range::iterator hint;
 if( skin_verts ) hint = skin_verts->begin();
 std::vector< EntityHandle > storage;
 const EntityHandle* conn;
 int len;
 bool find_edges = skin_edges || create_edges;
 bool printed_nonconformal_ho_warning = false;
 EntityHandle face;
 
 if( !faces.all_of_dimension( 2 ) ) return MB_TYPE_OUT_OF_RANGE;
 
 // get all the vertices
 Range verts;
 rval = thisMB->get_connectivity( faces, verts, true );
 if( MB_SUCCESS != rval ) return rval;
 
 std::vector< char > tag_vals;
 std::vector< EntityHandle > adj;
 AdjSides< 2 > adj_edges;
 for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
 {
 bool higher_order = false;
 
 // get all adjacent faces
 adj.clear();
 rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
 if( MB_SUCCESS != rval ) return rval;
 if( adj.empty() ) continue;
 
 // remove those not in the input list (intersect with input list)
 i = j = adj.begin();
 tag_vals.resize( adj.size() );
 rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
 if( MB_SUCCESS != rval ) return rval;
 // remove non-tagged entries
 i = j = adj.begin();
 for( ; i != adj.end(); ++i )
 if( tag_vals[i - adj.begin()] ) *( j++ ) = *i;
 adj.erase( j, adj.end() );
 
 // For each adjacent face, check the edges adjacent to the current vertex
 adj_edges.clear(); // other vertex for adjacent edges
 for( i = adj.begin(); i != adj.end(); ++i )
 {
 rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
 if( MB_SUCCESS != rval ) return rval;
 
 // For a single face element adjacent to this vertex, there
 // will be exactly two sides (edges) adjacent to the vertex.
 // Find the other vertex for each of the two edges.
 
 EntityHandle prev, next; // vertices of two adjacent edge-sides
 const int idx = std::find( conn, conn + len, *it ) - conn;
 assert( idx != len );
 
 if( TYPE_FROM_HANDLE( *i ) == MBTRI && len > 3 )
 {
 len = 3;
 higher_order = true;
 if( idx > 2 )
 { // skip higher-order nodes for now
 if( !printed_nonconformal_ho_warning )
 {
 printed_nonconformal_ho_warning = true;
 std::cerr << "Non-conformal higher-order mesh detected in skinner: "
 << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
 << "some elements and a higher-order node in others" << std::endl;
 }
 continue;
 }
 }
 else if( TYPE_FROM_HANDLE( *i ) == MBQUAD && len > 4 )
 {
 len = 4;
 higher_order = true;
 if( idx > 3 )
 { // skip higher-order nodes for now
 if( !printed_nonconformal_ho_warning )
 {
 printed_nonconformal_ho_warning = true;
 std::cerr << "Non-conformal higher-order mesh detected in skinner: "
 << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
 << "some elements and a higher-order node in others" << std::endl;
 }
 continue;
 }
 }
 
 // so it must be a MBPOLYGON
 const int prev_idx = ( idx + len - 1 ) % len; // this should be fine, always, even for padded case
 prev = conn[prev_idx];
 next = conn[( idx + 1 ) % len];
 if( next == conn[idx] ) // it must be the padded case, so roll to the beginning
 next = conn[0];
 
 // Insert sides (edges) in our list of candidate skin sides
 adj_edges.insert( &prev, 1, *i, prev_idx );
 adj_edges.insert( &next, 1, *i, idx );
 }
 
 // If vertex is not on skin, advance to next vertex.
 // adj_edges handled checking for duplicates on insertion.
 // If every candidate skin edge occurred more than once (was
 // not in fact on the skin), then we're done with this vertex.
 if( 0 == adj_edges.num_skin() ) continue;
 
 // If user requested Range of *vertices* on the skin...
 if( skin_verts )
 {
 // Put skin vertex in output list
 hint = skin_verts->insert( hint, *it );
 
 // Add mid edge nodes to vertex list
 if( !corners_only && higher_order )
 {
 for( AdjSides< 2 >::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p )
 {
 if( p->skin() )
 {
 face = p->adj_elem;
 EntityType type = TYPE_FROM_HANDLE( face );
 
 rval = thisMB->get_connectivity( face, conn, len, false );
 if( MB_SUCCESS != rval ) return rval;
 if( !CN::HasMidEdgeNodes( type, len ) ) continue;
 
 EntityHandle ec[2] = { *it, p->handles[0] };
 int side, sense, offset;
 CN::SideNumber( type, conn, ec, 2, 1, side, sense, offset );
 offset = CN::HONodeIndex( type, len, 1, side );
 assert( offset >= 0 && offset < len );
 skin_verts->insert( conn[offset] );
 }
 }
 }
 }
 
 // If user requested Range of *edges* on the skin...
 if( find_edges )
 {
 // Search list of existing adjacent edges for any that are on the skin
 adj.clear();
 rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
 if( MB_SUCCESS != rval ) return rval;
 for( i = adj.begin(); i != adj.end(); ++i )
 {
 rval = thisMB->get_connectivity( *i, conn, len, true );
 if( MB_SUCCESS != rval ) return rval;
 
 // bool equality expression within find_and_unmark call
 // will be evaluate to the index of *it in the conn array.
 //
 // Note that the order of the terms in the if statement is important.
 // We want to unmark any existing skin edges even if we aren't
 // returning them. Otherwise we'll end up creating duplicates
 // if create_edges is true and skin_edges is not.
 if( adj_edges.find_and_unmark( conn, ( conn[1] == *it ), face ) && skin_edges )
 {
 if( reversed_edges && edge_reversed( face, conn ) )
 reversed_edges->insert( *i );
 else
 skin_edges->insert( *i );
 }
 }
 }
 
 // If the user requested that we create new edges for sides
 // on the skin for which there is no existing edge, and there
 // are still skin sides for which no corresponding edge was found...
 if( create_edges && adj_edges.num_skin() )
 {
 // Create any skin edges that don't exist
 for( AdjSides< 2 >::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p )
 {
 if( p->skin() )
 {
 EntityHandle edge, ec[] = { *it, p->handles[0] };
 rval = create_side( this_set, p->adj_elem, MBEDGE, ec, edge );
 if( MB_SUCCESS != rval ) return rval;
 if( skin_edges ) skin_edges->insert( edge );
 }
 }
 }
 
 } // end for each vertex
 
 return MB_SUCCESS;
}
 
ErrorCode Skinner::find_skin_vertices_3D( const EntityHandle this_set,
 Tag tag,
 const Range& entities,
 Range* skin_verts,
 Range* skin_faces,
 Range* reversed_faces,
 bool create_faces,
 bool corners_only )
{
 // This function iterates over all the vertices contained in the
 // input vol elem list. For each such vertex, it then iterates over
 // all of the sides of the vol elements adjacent to the vertex.
 // If an adjacent side is the side of only one of the input
 // elements, then that side is on the skin.
 //
 // This algorithm will visit each skin vertex exactly once. It
 // will visit each skin side once for each vertex in the side.
 //
 // This function expects the caller to have created the passed bit
 // tag and set it to one only for the elements in the passed range. This
 // tag is used to do a fast intersection of the elements adjacent to a
 // vertex with the elements in the input range (discard any for which the
 // tag is not set to one.)
 //
 // For each vertex, iterate over each adjacent element. Construct
 // lists of the sides of each adjacent element that contain the vertex.
 //
 // A list of three-vertex sides is kept for all triangular faces,
 // included three-vertex faces of type MBPOLYGON. Putting polygons
 // in the same list ensures that we find polyhedron and non-polyhedron
 // elements that are adjacent.
 //
 // A list of four-vertex sides is kept for all quadrilateral faces,
 // including four-vertex faces of type MBPOLYGON.
 //
 // Sides with more than four vertices must have an explicit MBPOLYGON
 // element representing them because MBPOLYHEDRON connectivity is a
 // list of faces rather than vertices. So the third list (vertices>=5),
 // need contain only the handle of the face rather than the vertex handles.
 
 ErrorCode rval;
 std::vector< EntityHandle >::iterator i, j;
 Range::iterator hint;
 if( skin_verts ) hint = skin_verts->begin();
 std::vector< EntityHandle > storage, storage2; // temp storage for conn lists
 const EntityHandle *conn, *conn2;
 int len, len2;
 bool find_faces = skin_faces || create_faces;
 int clen, side, sense, offset, indices[9];
 EntityType face_type;
 EntityHandle elem;
 bool printed_nonconformal_ho_warning = false;
 
 if( !entities.all_of_dimension( 3 ) ) return MB_TYPE_OUT_OF_RANGE;
 
 // get all the vertices
 Range verts;
 rval = thisMB->get_connectivity( entities, verts, true );
 if( MB_SUCCESS != rval ) return rval;
 // if there are polyhedra in the input list, need to make another
 // call to get vertices from faces
 if( !verts.all_of_dimension( 0 ) )
 {
 Range::iterator it = verts.upper_bound( MBVERTEX );
 Range pfaces;
 pfaces.merge( it, verts.end() );
 verts.erase( it, verts.end() );
 rval = thisMB->get_connectivity( pfaces, verts, true );
 if( MB_SUCCESS != rval ) return rval;
 assert( verts.all_of_dimension( 0 ) );
 }
 
 AdjSides< 4 > adj_quads; // 4-node sides adjacent to a vertex
 AdjSides< 3 > adj_tris; // 3-node sides adjacent to a vertex
 AdjSides< 2 > adj_poly; // n-node sides (n>5) adjacent to vertex
 // (must have an explicit polygon, so store
 // polygon handle rather than vertices.)
 std::vector< char > tag_vals;
 std::vector< EntityHandle > adj;
 for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
 {
 bool higher_order = false;
 
 // get all adjacent elements
 adj.clear();
 rval = thisMB->get_adjacencies( &*it, 1, 3, false, adj );
 if( MB_SUCCESS != rval ) return rval;
 if( adj.empty() ) continue;
 
 // remove those not tagged (intersect with input range)
 i = j = adj.begin();
 tag_vals.resize( adj.size() );
 rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
 if( MB_SUCCESS != rval ) return rval;
 for( ; i != adj.end(); ++i )
 if( tag_vals[i - adj.begin()] ) *( j++ ) = *i;
 adj.erase( j, adj.end() );
 
 // Build lists of sides of 3D element adjacent to the current vertex
 adj_quads.clear(); // store three other vertices for each adjacent quad face
 adj_tris.clear(); // store two other vertices for each adjacent tri face
 adj_poly.clear(); // store handle of each adjacent polygonal face
 int idx;
 for( i = adj.begin(); i != adj.end(); ++i )
 {
 const EntityType type = TYPE_FROM_HANDLE( *i );
 
 // Special case for POLYHEDRA
 if( type == MBPOLYHEDRON )
 {
 rval = thisMB->get_connectivity( *i, conn, len );
 if( MB_SUCCESS != rval ) return rval;
 for( int k = 0; k < len; ++k )
 {
 rval = thisMB->get_connectivity( conn[k], conn2, len2, true, &storage2 );
 if( MB_SUCCESS != rval ) return rval;
 idx = std::find( conn2, conn2 + len2, *it ) - conn2;
 if( idx == len2 ) // vertex not in this face
 continue;
 
 // Treat 3- and 4-vertex faces specially, so that
 // if the mesh contains both elements and polyhedra,
 // we don't miss one type adjacent to the other.
 switch( len2 )
 {
 case 3:
 adj_tris.insert( conn2, idx, *i, k );
 break;
 case 4:
 adj_quads.insert( conn2, idx, *i, k );
 break;
 default:
 adj_poly.insert( conn + k, 1, *i, k );
 break;
 }
 }
 }
 else
 {
 rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
 if( MB_SUCCESS != rval ) return rval;
 
 idx = std::find( conn, conn + len, *it ) - conn;
 assert( idx != len );
 
 if( len > CN::VerticesPerEntity( type ) )
 {
 higher_order = true;
 // skip higher-order nodes for now
 if( idx >= CN::VerticesPerEntity( type ) )
 {
 if( !printed_nonconformal_ho_warning )
 {
 printed_nonconformal_ho_warning = true;
 std::cerr << "Non-conformal higher-order mesh detected in skinner: "
 << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
 << "some elements and a higher-order node in others" << std::endl;
 }
 continue;
 }
 }
 
 // For each side of the element...
 const int num_faces = CN::NumSubEntities( type, 2 );
 for( int f = 0; f < num_faces; ++f )
 {
 int num_vtx;
 const short* face_indices = CN::SubEntityVertexIndices( type, 2, f, face_type, num_vtx );
 const short face_idx = std::find( face_indices, face_indices + num_vtx, (short)idx ) - face_indices;
 // skip sides that do not contain vertex from outer loop
 if( face_idx == num_vtx ) continue; // current vertex not in this face
 
 assert( num_vtx <= 4 ); // polyhedra handled above
 switch( face_type )
 {
 case MBTRI:
 adj_tris.insert( conn, face_idx, *i, f, face_indices );
 break;
 case MBQUAD:
 adj_quads.insert( conn, face_idx, *i, f, face_indices );
 break;
 default:
 return MB_TYPE_OUT_OF_RANGE;
 }
 }
 }
 } // end for (adj[3])
 
 // If vertex is not on skin, advance to next vertex
 if( 0 == ( adj_tris.num_skin() + adj_quads.num_skin() + adj_poly.num_skin() ) ) continue;
 
 // If user requested that skin *vertices* be passed back...
 if( skin_verts )
 {
 // Put skin vertex in output list
 hint = skin_verts->insert( hint, *it );
 
 // Add mid-edge and mid-face nodes to vertex list
 if( !corners_only && higher_order )
 {
 for( AdjSides< 3 >::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t )
 {
 if( t->skin() )
 {
 elem = t->adj_elem;
 EntityType type = TYPE_FROM_HANDLE( elem );
 
 rval = thisMB->get_connectivity( elem, conn, len, false );
 if( MB_SUCCESS != rval ) return rval;
 if( !CN::HasMidNodes( type, len ) ) continue;
 
 EntityHandle ec[3] = { *it, t->handles[0], t->handles[1] };
 CN::SideNumber( type, conn, ec, 3, 2, side, sense, offset );
 CN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
 assert( MBTRI == face_type );
 for( int k = 3; k < clen; ++k )
 skin_verts->insert( conn[indices[k]] );
 }
 }
 for( AdjSides< 4 >::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q )
 {
 if( q->skin() )
 {
 elem = q->adj_elem;
 EntityType type = TYPE_FROM_HANDLE( elem );
 
 rval = thisMB->get_connectivity( elem, conn, len, false );
 if( MB_SUCCESS != rval ) return rval;
 if( !CN::HasMidNodes( type, len ) ) continue;
 
 EntityHandle ec[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
 CN::SideNumber( type, conn, ec, 4, 2, side, sense, offset );
 CN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
 assert( MBQUAD == face_type );
 for( int k = 4; k < clen; ++k )
 skin_verts->insert( conn[indices[k]] );
 }
 }
 }
 }
 
 // If user requested that we pass back the list of 2D elements
 // representing the skin of the mesh...
 if( find_faces )
 {
 // Search list of adjacent faces for any that are on the skin
 adj.clear();
 rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
 if( MB_SUCCESS != rval ) return rval;
 
 for( i = adj.begin(); i != adj.end(); ++i )
 {
 rval = thisMB->get_connectivity( *i, conn, len, true );
 if( MB_SUCCESS != rval ) return rval;
 const int idx2 = std::find( conn, conn + len, *it ) - conn;
 if( idx2 >= len )
 {
 assert( printed_nonconformal_ho_warning );
 continue;
 }
 
 // Note that the order of the terms in the if statements below
 // is important. We want to unmark any existing skin faces even
 // if we aren't returning them. Otherwise we'll end up creating
 // duplicates if create_faces is true.
 if( 3 == len )
 {
 if( adj_tris.find_and_unmark( conn, idx2, elem ) && skin_faces )
 {
 if( reversed_faces && face_reversed( elem, conn, MBTRI ) )
 reversed_faces->insert( *i );
 else
 skin_faces->insert( *i );
 }
 }
 else if( 4 == len )
 {
 if( adj_quads.find_and_unmark( conn, idx2, elem ) && skin_faces )
 {
 if( reversed_faces && face_reversed( elem, conn, MBQUAD ) )
 reversed_faces->insert( *i );
 else
 skin_faces->insert( *i );
 }
 }
 else
 {
 if( adj_poly.find_and_unmark( &*i, 1, elem ) && skin_faces ) skin_faces->insert( *i );
 }
 }
 }
 
 // If user does not want use to create new faces representing
 // sides for which there is currently no explicit element,
 // skip the remaining code and advance the outer loop to the
 // next vertex.
 if( !create_faces ) continue;
 
 // Polyhedra always have explicitly defined faces, so
 // there is no way we could need to create such a face.
 assert( 0 == adj_poly.num_skin() );
 
 // Create any skin tris that don't exist
 if( adj_tris.num_skin() )
 {
 for( AdjSides< 3 >::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t )
 {
 if( t->skin() )
 {
 EntityHandle tri, c[3] = { *it, t->handles[0], t->handles[1] };
 rval = create_side( this_set, t->adj_elem, MBTRI, c, tri );
 if( MB_SUCCESS != rval ) return rval;
 if( skin_faces ) skin_faces->insert( tri );
 }
 }
 }
 
 // Create any skin quads that don't exist
 if( adj_quads.num_skin() )
 {
 for( AdjSides< 4 >::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q )
 {
 if( q->skin() )
 {
 EntityHandle quad, c[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
 rval = create_side( this_set, q->adj_elem, MBQUAD, c, quad );
 if( MB_SUCCESS != rval ) return rval;
 if( skin_faces ) skin_faces->insert( quad );
 }
 }
 }
 } // end for each vertex
 
 return MB_SUCCESS;
}
 
} // namespace moab