AEntityFactory.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.
 *
 */

#include "AEntityFactory.hpp"
#include "Internals.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/Error.hpp"
#include "moab/CN.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "EntitySequence.hpp"
#include "SequenceData.hpp"
#include "SequenceManager.hpp"
#include "RangeSeqIntersectIter.hpp"

#include <cassert>
#include <algorithm>
#include <set>

namespace moab
{

ErrorCode AEntityFactory::get_vertices( EntityHandle h,
                                        const EntityHandle*& vect_out,
                                        int& count_out,
                                        std::vector< EntityHandle >& storage )
{
    ErrorCode result;
    if( MBPOLYHEDRON == TYPE_FROM_HANDLE( h ) )
    {
        storage.clear();
        result    = thisMB->get_adjacencies( &h, 1, 0, false, storage );
        vect_out  = &storage[0];
        count_out = storage.size();
    }
    else
    {
        result = thisMB->get_connectivity( h, vect_out, count_out, false, &storage );
    }
    return result;
}

AEntityFactory::AEntityFactory( Core* mdb )
{
    assert( NULL != mdb );
    thisMB       = mdb;
    mVertElemAdj = false;
}

AEntityFactory::~AEntityFactory()
{
    // clean up all the adjacency information that was created
    EntityType ent_type;

    // iterate through each element type
    for( ent_type = MBVERTEX; ent_type <= MBENTITYSET; ent_type++ )
    {
        TypeSequenceManager::iterator i;
        TypeSequenceManager& seqman = thisMB->sequence_manager()->entity_map( ent_type );
        for( i = seqman.begin(); i != seqman.end(); ++i )
        {
            std::vector< EntityHandle >** adj_list = ( *i )->data()->get_adjacency_data();
            if( !adj_list ) continue;
            adj_list += ( *i )->start_handle() - ( *i )->data()->start_handle();

            for( EntityID j = 0; j < ( *i )->size(); ++j )
            {
                delete adj_list[j];
                adj_list[j] = 0;
            }
        }
    }
}

//! get the elements contained by source_entity, of
//! type target_type, passing back in target_entities; if create_if_missing
//! is true and no entity is found, one is created; if create_adjacency_option
//! is >= 0, adjacencies from entities of that dimension to each target_entity
//! are created (this function uses AEntityFactory::get_element for each element)
ErrorCode AEntityFactory::get_elements( EntityHandle source_entity,
                                        const unsigned int target_dimension,
                                        std::vector< EntityHandle >& target_entities,
                                        const bool create_if_missing,
                                        const int create_adjacency_option )
{
    // check for trivial case first
    const EntityType source_type    = TYPE_FROM_HANDLE( source_entity );
    const unsigned source_dimension = CN::Dimension( source_type );

    if( source_type >= MBENTITYSET || target_dimension < 1 || target_dimension > 3 )
    {
        return MB_TYPE_OUT_OF_RANGE;
    }
    else if( source_dimension == target_dimension )
    {
        target_entities.push_back( source_entity );
        return MB_SUCCESS;
    }

    ErrorCode result;
    if( mVertElemAdj == false )
    {
        result = create_vert_elem_adjacencies();
        if( MB_SUCCESS != result ) return result;
    }

    if( source_dimension == 0 )
    {
        result = get_zero_to_n_elements( source_entity, target_dimension, target_entities, create_if_missing,
                                         create_adjacency_option );
    }
    else if( source_dimension > target_dimension )
    {
        result = get_down_adjacency_elements( source_entity, target_dimension, target_entities, create_if_missing,
                                              create_adjacency_option );
    }
    else  // if(source_dimension < target_dimension)
    {
        result = get_up_adjacency_elements( source_entity, target_dimension, target_entities, create_if_missing,
                                            create_adjacency_option );
    }

    return result;
}

ErrorCode AEntityFactory::get_polyhedron_vertices( const EntityHandle source_entity,
                                                   std::vector< EntityHandle >& target_entities )
{
    // get the connectivity array pointer
    const EntityHandle* connect = NULL;
    int num_connect             = 0;
    ErrorCode result            = thisMB->get_connectivity( source_entity, connect, num_connect );
    if( MB_SUCCESS != result ) return result;

    // now get the union of those polygons' vertices
    result = thisMB->get_adjacencies( connect, num_connect, 0, false, target_entities, Interface::UNION );
    return result;
}

ErrorCode AEntityFactory::get_associated_meshsets( EntityHandle source_entity,
                                                   std::vector< EntityHandle >& target_entities )
{

    ErrorCode result;

    const EntityHandle* adj_vec;
    int num_adj;
    result = get_adjacencies( source_entity, adj_vec, num_adj );
    if( result != MB_SUCCESS || adj_vec == NULL ) return result;

    // find the meshsets in this vector
    DimensionPair dim_pair = CN::TypeDimensionMap[4];
    int dum;
    const EntityHandle* start_ent =
        std::lower_bound( adj_vec, adj_vec + num_adj, CREATE_HANDLE( dim_pair.first, MB_START_ID, dum ) );
    const EntityHandle* end_ent =
        std::lower_bound( start_ent, adj_vec + num_adj, CREATE_HANDLE( dim_pair.second, MB_END_ID, dum ) );

    // copy the the meshsets
    target_entities.insert( target_entities.end(), start_ent, end_ent );

    return result;
}

//! get the element defined by the vertices in vertex_list, of the
//! type target_type, passing back in target_entity; if create_if_missing
//! is true and no entity is found, one is created; if create_adjacency_option
//! is >= 0, adjacencies from entities of that dimension to target_entity
//! are created (only create_adjacency_option=0 is supported right now,
//! so that never creates other ancillary entities)
ErrorCode AEntityFactory::get_element( const EntityHandle* vertex_list,
                                       const int vertex_list_size,
                                       const EntityType target_type,
                                       EntityHandle& target_entity,
                                       const bool create_if_missing,
                                       const EntityHandle source_entity,
                                       const int /*create_adjacency_option*/ )
{

    // look over nodes to see if this entity already exists
    target_entity = 0;
    ErrorCode result;
    const EntityHandle *i_adj, *end_adj;

    // need vertex adjacencies, so create if necessary
    if( mVertElemAdj == false ) create_vert_elem_adjacencies();

    // get the adjacency list
    const EntityHandle* adj_vec;
    int num_adj;
    result = get_adjacencies( vertex_list[0], adj_vec, num_adj );
    if( result != MB_SUCCESS || adj_vec == NULL ) return result;

    // check to see if any of these are equivalent to the vertex list
    int dum;

    // use a fixed-size array, for speed; there should never be more than 5 equivalent entities
    EntityHandle temp_vec[15];
    int temp_vec_size = 0;

    i_adj   = std::lower_bound( adj_vec, adj_vec + num_adj, CREATE_HANDLE( target_type, MB_START_ID, dum ) );
    end_adj = std::lower_bound( i_adj, adj_vec + num_adj, CREATE_HANDLE( target_type, MB_END_ID, dum ) );
    for( ; i_adj != end_adj; ++i_adj )
    {
        if( TYPE_FROM_HANDLE( *i_adj ) != target_type ) continue;

        if( true == entities_equivalent( *i_adj, vertex_list, vertex_list_size, target_type ) )
        {
            temp_vec[temp_vec_size++] = *i_adj;
        }
    }

    if( temp_vec_size == 0 && !create_if_missing ) return result;

    // test for size against fixed-size array
    assert( temp_vec_size <= 15 );

    // test for empty first, 'cuz it's cheap
    if( temp_vec_size == 0 && true == create_if_missing )
    {

        // Create the element with this handle (handle is a return type and should be the last
        // parameter)
        result = thisMB->create_element( target_type, vertex_list, vertex_list_size, target_entity );
    }

    // next most likely is one entity
    else if( temp_vec_size == 1 )
        target_entity = temp_vec[0];

    // least likely, most work - leave for last test
    else
    {
        // multiple entities found - look for direct adjacencies
        if( 0 != source_entity )
        {

            int num_adjs;
            for( dum = 0; dum < temp_vec_size; dum++ )
            {
                result = get_adjacencies( temp_vec[dum], adj_vec, num_adjs );
                if( std::find( adj_vec, ( adj_vec + num_adjs ), source_entity ) != ( adj_vec + num_adjs ) )
                {
                    // found it, return it
                    target_entity = temp_vec[dum];
                    break;
                }
            }

            if( 0 == target_entity &&
                thisMB->dimension_from_handle( source_entity ) > CN::Dimension( target_type ) + 1 )
            {
                // still have multiple entities, and source dimension is two greater than target,
                // so there may not be any explicit adjacencies between the two; look for common
                // entities of the intermediate dimension
                MeshTopoUtil mtu( thisMB );
                int intermed_dim = CN::Dimension( target_type ) + 1;
                for( dum = 0; dum < temp_vec_size; dum++ )
                {
                    if( 0 != mtu.common_entity( temp_vec[dum], source_entity, intermed_dim ) )
                    {
                        target_entity = temp_vec[dum];
                        break;
                    }
                }
            }
        }

        if( target_entity == 0 )
        {
            // if we get here, we didn't find a matching adjacency; just take the first one, but
            // return a non-success result
            target_entity = temp_vec[0];
            result        = MB_MULTIPLE_ENTITIES_FOUND;
        }
    }

    return result;
}

bool AEntityFactory::entities_equivalent( const EntityHandle this_entity,
                                          const EntityHandle* vertex_list,
                                          const int vertex_list_size,
                                          const EntityType target_type )
{
    // compare vertices of this_entity with those in the list, returning true if they
    // represent the same element
    EntityType this_type = TYPE_FROM_HANDLE( this_entity );

    if( this_type != target_type )
        return false;

    else if( this_type == MBVERTEX && ( vertex_list_size > 1 || vertex_list[0] != this_entity ) )
        return false;

    // need to compare the actual vertices
    const EntityHandle* this_vertices = NULL;
    int num_this_vertices             = 0;
    std::vector< EntityHandle > storage;
    thisMB->get_connectivity( this_entity, this_vertices, num_this_vertices, false, &storage );

    // see if we can get one node id to match
    assert( vertex_list_size > 0 );
    int num_corner_verts =
        ( ( this_type == MBPOLYGON || this_type == MBPOLYHEDRON ) ? num_this_vertices
                                                                  : CN::VerticesPerEntity( target_type ) );
    const EntityHandle* iter = std::find( this_vertices, ( this_vertices + num_corner_verts ), vertex_list[0] );
    if( iter == ( this_vertices + num_corner_verts ) ) return false;

    // now lets do connectivity matching
    bool they_match = true;

    // line up our vectors
    int i;
    int offset = iter - this_vertices;

    // first compare forward
    for( i = 1; i < num_corner_verts; ++i )
    {
        if( i >= vertex_list_size )
        {
            they_match = false;
            break;
        }

        if( vertex_list[i] != this_vertices[( offset + i ) % num_corner_verts] )
        {
            they_match = false;
            break;
        }
    }

    if( they_match == true ) return true;

    they_match = true;

    // then compare reverse
    // offset iter to avoid addition inside loop; this just makes sure we don't
    // go off beginning of this_vertices with an index < 0
    offset += num_corner_verts;
    for( i = 1; i < num_corner_verts; i++ )
    {
        if( vertex_list[i] != this_vertices[( offset - i ) % num_corner_verts] )
        {
            they_match = false;
            break;
        }
    }
    return they_match;
}

//! add an adjacency from from_ent to to_ent; if both_ways is true, add one
//! in reverse too
//! NOTE: this function is defined even though we may only be implementing
//! vertex-based up-adjacencies
ErrorCode AEntityFactory::add_adjacency( EntityHandle from_ent, EntityHandle to_ent, const bool both_ways )
{
    EntityType to_type = TYPE_FROM_HANDLE( to_ent );

    if( to_type == MBVERTEX ) return MB_ALREADY_ALLOCATED;

    AdjacencyVector* adj_list_ptr = NULL;
    ErrorCode result              = get_adjacencies( from_ent, adj_list_ptr, true );
    if( MB_SUCCESS != result ) return result;

    // get an iterator to the right spot in this sorted vector
    AdjacencyVector::iterator adj_iter;
    if( !adj_list_ptr->empty() )
    {
        adj_iter = std::lower_bound( adj_list_ptr->begin(), adj_list_ptr->end(), to_ent );

        if( adj_iter == adj_list_ptr->end() || to_ent != *adj_iter )
        {
            adj_list_ptr->insert( adj_iter, to_ent );
        }
    }
    else
        adj_list_ptr->push_back( to_ent );

    // if both_ways is true, recursively call this function
    if( true == both_ways && to_type != MBVERTEX ) result = add_adjacency( to_ent, from_ent, false );

    return result;
}

//! remove an adjacency from from the base_entity.
ErrorCode AEntityFactory::remove_adjacency( EntityHandle base_entity, EntityHandle adj_to_remove )
{
    ErrorCode result;

    if( TYPE_FROM_HANDLE( base_entity ) == MBENTITYSET )
        return thisMB->remove_entities( base_entity, &adj_to_remove, 1 );

    // get the adjacency tag
    AdjacencyVector* adj_list = NULL;
    result                    = get_adjacencies( base_entity, adj_list );
    if( adj_list == NULL || MB_SUCCESS != result ) return result;

    // remove the specified entity from the adjacency list and truncate
    // the list to the new length
    adj_list->erase( std::remove( adj_list->begin(), adj_list->end(), adj_to_remove ), adj_list->end() );

    return result;
}

//! remove all adjacencies from from the base_entity.
ErrorCode AEntityFactory::remove_all_adjacencies( EntityHandle base_entity, const bool delete_adj_list )
{
    ErrorCode result;
    EntityType base_type = TYPE_FROM_HANDLE( base_entity );

    if( base_type == MBENTITYSET ) return thisMB->clear_meshset( &base_entity, 1 );
    const int base_ent_dim = CN::Dimension( base_type );

    // Remove adjacencies from element vertices back to
    // this element.  Also check any elements adjacent
    // to the vertex and of higher dimension than this
    // element for downward adjacencies to this element.
    if( vert_elem_adjacencies() && base_type != MBVERTEX )
    {
        EntityHandle const *connvect = 0, *adjvect = 0;
        int numconn = 0, numadj = 0;
        std::vector< EntityHandle > connstorage;
        result = get_vertices( base_entity, connvect, numconn, connstorage );
        if( MB_SUCCESS != result ) return result;

        for( int i = 0; i < numconn; ++i )
        {
            result = get_adjacencies( connvect[i], adjvect, numadj );
            if( MB_SUCCESS != result ) return result;

            bool remove_this = false;
            for( int j = 0; j < numadj; ++j )
            {
                if( adjvect[j] == base_entity ) remove_this = true;

                if( CN::Dimension( TYPE_FROM_HANDLE( adjvect[j] ) ) != base_ent_dim &&
                    explicitly_adjacent( adjvect[j], base_entity ) )
                    remove_adjacency( adjvect[j], base_entity );
            }

            if( remove_this ) remove_adjacency( connvect[i], base_entity );
        }
    }

    // get the adjacency tag
    AdjacencyVector* adj_list = 0;
    result                    = get_adjacencies( base_entity, adj_list );
    if( MB_SUCCESS != result || !adj_list ) return result;

    // check adjacent entities for references back to this entity
    for( AdjacencyVector::reverse_iterator it = adj_list->rbegin(); it != adj_list->rend(); ++it )
        remove_adjacency( *it, base_entity );

    if( delete_adj_list )
        result = set_adjacency_ptr( base_entity, NULL );
    else
        adj_list->clear();

    return MB_SUCCESS;
}

ErrorCode AEntityFactory::create_vert_elem_adjacencies()
{

    mVertElemAdj = true;

    EntityType ent_type;
    Range::iterator i_range;
    const EntityHandle* connectivity;
    std::vector< EntityHandle > aux_connect;
    int number_nodes;
    ErrorCode result;
    Range handle_range;

    // 1. over all element types, for each element, create vertex-element adjacencies
    for( ent_type = MBEDGE; ent_type != MBENTITYSET; ent_type++ )
    {
        handle_range.clear();

        // get this type of entity
        result = thisMB->get_entities_by_type( 0, ent_type, handle_range );
        if( result != MB_SUCCESS ) return result;

        for( i_range = handle_range.begin(); i_range != handle_range.end(); ++i_range )
        {
            result = get_vertices( *i_range, connectivity, number_nodes, aux_connect );
            if( MB_SUCCESS != result ) return result;

            // add the adjacency
            for( int k = 0; k < number_nodes; k++ )
                if( ( result = add_adjacency( connectivity[k], *i_range ) ) != MB_SUCCESS ) return result;
        }
    }

    return MB_SUCCESS;
}

ErrorCode AEntityFactory::get_adjacencies( EntityHandle entity,
                                           const EntityHandle*& adjacent_entities,
                                           int& num_entities ) const
{
    AdjacencyVector const* vec_ptr = 0;
    ErrorCode result               = get_adjacency_ptr( entity, vec_ptr );
    if( MB_SUCCESS != result || !vec_ptr )
    {
        adjacent_entities = 0;
        num_entities      = 0;
        return result;
    }

    num_entities      = vec_ptr->size();
    adjacent_entities = ( vec_ptr->empty() ) ? NULL : &( ( *vec_ptr )[0] );
    return MB_SUCCESS;
}

ErrorCode AEntityFactory::get_adjacencies( EntityHandle entity, std::vector< EntityHandle >& adjacent_entities ) const
{
    AdjacencyVector const* vec_ptr = 0;
    ErrorCode result               = get_adjacency_ptr( entity, vec_ptr );
    if( MB_SUCCESS != result || !vec_ptr )
    {
        adjacent_entities.clear();
        return result;
    }

    adjacent_entities = *vec_ptr;
    return MB_SUCCESS;
}

ErrorCode AEntityFactory::get_adjacencies( EntityHandle entity, std::vector< EntityHandle >*& adj_vec, bool create )
{
    adj_vec          = 0;
    ErrorCode result = get_adjacency_ptr( entity, adj_vec );
    if( MB_SUCCESS == result && !adj_vec && create )
    {
        adj_vec = new AdjacencyVector;
        result  = set_adjacency_ptr( entity, adj_vec );
        if( MB_SUCCESS != result )
        {
            delete adj_vec;
            adj_vec = 0;
        }
    }
    return result;
}

ErrorCode AEntityFactory::get_adjacencies( const EntityHandle source_entity,
                                           const unsigned int target_dimension,
                                           bool create_if_missing,
                                           std::vector< EntityHandle >& target_entities )
{
    const EntityType source_type    = TYPE_FROM_HANDLE( source_entity );
    const unsigned source_dimension = CN::Dimension( source_type );

    ErrorCode result;
    if( target_dimension == 4 )
    {  // get meshsets 'source' is in
        result = get_associated_meshsets( source_entity, target_entities );
    }
    else if( target_dimension == ( source_type != MBPOLYHEDRON ? 0 : 2 ) )
    {
        std::vector< EntityHandle > tmp_storage;
        const EntityHandle* conn = NULL;
        int len                  = 0;
        result                   = thisMB->get_connectivity( source_entity, conn, len, false, &tmp_storage );
        target_entities.insert( target_entities.end(), conn, conn + len );
    }
    else if( target_dimension == 0 && source_type == MBPOLYHEDRON )
    {
        result = get_polyhedron_vertices( source_entity, target_entities );
    }
    else if( source_dimension == target_dimension )
    {
        target_entities.push_back( source_entity );
        result = MB_SUCCESS;
    }
    else
    {
        if( mVertElemAdj == false )
        {
            result = create_vert_elem_adjacencies();
            if( MB_SUCCESS != result ) return result;
        }

        if( source_dimension == 0 )
        {
            result = get_zero_to_n_elements( source_entity, target_dimension, target_entities, create_if_missing );
        }
        else if( source_dimension > target_dimension )
        {
            result = get_down_adjacency_elements( source_entity, target_dimension, target_entities, create_if_missing );
        }
        else  // if(source_dimension < target_dimension)
        {
            result = get_up_adjacency_elements( source_entity, target_dimension, target_entities, create_if_missing );
        }
    }

    return result;
}

ErrorCode AEntityFactory::notify_create_entity( const EntityHandle entity,
                                                const EntityHandle* node_array,
                                                const int number_nodes )
{
    ErrorCode result = MB_SUCCESS, tmp_result;
    if( vert_elem_adjacencies() )
    {
        // iterate through nodes and add adjacency information
        if( TYPE_FROM_HANDLE( entity ) == MBPOLYHEDRON )
        {
            // polyhedron - get real vertex connectivity
            std::vector< EntityHandle > verts;
            tmp_result = get_adjacencies( entity, 0, false, verts );
            if( MB_SUCCESS != tmp_result ) return tmp_result;
            for( std::vector< EntityHandle >::iterator vit = verts.begin(); vit != verts.end(); ++vit )
            {
                tmp_result = add_adjacency( *vit, entity );
                if( MB_SUCCESS != tmp_result ) result = tmp_result;
            }
        }
        else
        {
            for( unsigned int i = number_nodes; i--; )
            {
                tmp_result = add_adjacency( node_array[i], entity );
                if( MB_SUCCESS != tmp_result ) result = tmp_result;
            }
        }
    }

    return result;
}

ErrorCode AEntityFactory::get_zero_to_n_elements( EntityHandle source_entity,
                                                  const unsigned int target_dimension,
                                                  std::vector< EntityHandle >& target_entities,
                                                  const bool create_if_missing,
                                                  const int /*create_adjacency_option = -1*/ )
{
    AdjacencyVector::iterator start_ent, end_ent;

    // get the adjacency vector
    AdjacencyVector* adj_vec = NULL;
    ErrorCode result         = get_adjacencies( source_entity, adj_vec );
    if( result != MB_SUCCESS || adj_vec == NULL ) return result;

    if( target_dimension < 3 && create_if_missing )
    {
        std::vector< EntityHandle > tmp_ents;

        start_ent = std::lower_bound( adj_vec->begin(), adj_vec->end(),
                                      FIRST_HANDLE( CN::TypeDimensionMap[target_dimension + 1].first ) );

        end_ent = std::lower_bound( start_ent, adj_vec->end(), LAST_HANDLE( CN::TypeDimensionMap[3].second ) );

        std::vector< EntityHandle > elems( start_ent, end_ent );

        // make target_dimension elements from all adjacient higher-dimension elements
        for( start_ent = elems.begin(); start_ent != elems.end(); ++start_ent )
        {
            tmp_ents.clear();
            get_down_adjacency_elements( *start_ent, target_dimension, tmp_ents, create_if_missing, 0 );
        }
    }

    DimensionPair dim_pair = CN::TypeDimensionMap[target_dimension];
    start_ent              = std::lower_bound( adj_vec->begin(), adj_vec->end(), FIRST_HANDLE( dim_pair.first ) );
    end_ent                = std::lower_bound( start_ent, adj_vec->end(), LAST_HANDLE( dim_pair.second ) );
    target_entities.insert( target_entities.end(), start_ent, end_ent );
    return MB_SUCCESS;
}

ErrorCode AEntityFactory::get_down_adjacency_elements( EntityHandle source_entity,
                                                       const unsigned int target_dimension,
                                                       std::vector< EntityHandle >& target_entities,
                                                       const bool create_if_missing,
                                                       const int create_adjacency_option )
{

    EntityType source_type = TYPE_FROM_HANDLE( source_entity );

    if( source_type == MBPOLYHEDRON || source_type == MBPOLYGON )
        return get_down_adjacency_elements_poly( source_entity, target_dimension, target_entities, create_if_missing,
                                                 create_adjacency_option );

    // make this a fixed size to avoid cost of working with STL vectors
    EntityHandle vertex_array[27] = {};
    ErrorCode temp_result;

    const EntityHandle* vertices = NULL;
    int num_verts                = 0;

    // I know there are already vertex adjacencies for this - call
    // another function to get them
    std::vector< EntityHandle > storage;
    ErrorCode result = thisMB->get_connectivity( source_entity, vertices, num_verts, false, &storage );
    if( MB_SUCCESS != result ) return result;

    int has_mid_nodes[4];
    CN::HasMidNodes( source_type, num_verts, has_mid_nodes );

    std::vector< int > index_list;
    int num_sub_ents = CN::NumSubEntities( source_type, target_dimension );

    for( int j = 0; j < num_sub_ents; j++ )
    {
        const CN::ConnMap& cmap = CN::mConnectivityMap[source_type][target_dimension - 1];

        int verts_per_sub = cmap.num_corners_per_sub_element[j];

        // get the corner vertices
        for( int i = 0; i < verts_per_sub; i++ )
            vertex_array[i] = vertices[cmap.conn[j][i]];

        // get the ho nodes for sub-subfacets
        if( has_mid_nodes[1] && target_dimension > 1 )
        {
            // has edge mid-nodes; for each edge, get the right mid-node and put in vertices
            // first get the edge indices
            index_list.clear();
            int int_result = CN::AdjacentSubEntities( source_type, &j, 1, target_dimension, 1, index_list );
            if( 0 != int_result ) return MB_FAILURE;
            for( unsigned int k = 0; k < index_list.size(); k++ )
            {
                int tmp_index = CN::HONodeIndex( source_type, num_verts, 1, index_list[k] );
                if( tmp_index >= (int)num_verts ) return MB_INDEX_OUT_OF_RANGE;

                // put this vertex on the end; reuse verts_per_sub as an index
                vertex_array[verts_per_sub++] = vertices[tmp_index];
            }
        }
        // get the ho nodes for the target dimension
        if( has_mid_nodes[target_dimension] )
        {
            // get the ho node index for this subfacet
            int tmp_index = CN::HONodeIndex( source_type, num_verts, target_dimension, j );
            if( tmp_index >= num_verts ) return MB_INDEX_OUT_OF_RANGE;
            vertex_array[verts_per_sub++] = vertices[tmp_index];
        }

        EntityHandle tmp_target = 0;
        temp_result = get_element( vertex_array, verts_per_sub, cmap.target_type[j], tmp_target, create_if_missing,
                                   source_entity, create_adjacency_option );

        if( temp_result != MB_SUCCESS )
            result = temp_result;
        else if( 0 != tmp_target )
            target_entities.push_back( tmp_target );

        // make sure we're not writing past the end of our fixed-size array
        if( verts_per_sub > 27 ) return MB_INDEX_OUT_OF_RANGE;
    }

    return result;
}

ErrorCode AEntityFactory::get_down_adjacency_elements_poly( EntityHandle source_entity,
                                                            const unsigned int target_dimension,
                                                            std::vector< EntityHandle >& target_entities,
                                                            const bool create_if_missing,
                                                            const int /*create_adjacency_option*/ )
{

    EntityType source_type = TYPE_FROM_HANDLE( source_entity );

    if( !( source_type == MBPOLYHEDRON && target_dimension > 0 && target_dimension < 3 ) &&
        !( source_type == MBPOLYGON && target_dimension == 1 ) )
        return MB_TYPE_OUT_OF_RANGE;

    // make this a fixed size to avoid cost of working with STL vectors
    std::vector< EntityHandle > vertex_array;

    // I know there are already vertex adjacencies for this - call
    // another function to get them
    ErrorCode result = get_adjacencies( source_entity, 0, false, vertex_array );
    if( MB_SUCCESS != result ) return result;

    ErrorCode tmp_result;
    if( source_type == MBPOLYGON )
    {
        result = MB_SUCCESS;
        // put the first vertex on the end so we have a ring
        vertex_array.push_back( *vertex_array.begin() );
        for( unsigned int i = 0; i < vertex_array.size() - 1; i++ )
        {
            Range vrange, adj_edges;
            vrange.insert( vertex_array[i] );
            vrange.insert( vertex_array[i + 1] );
            // account for padded polygons; if the vertices are the same, skip
            if( vrange.size() == 1 ) continue;
            tmp_result = thisMB->get_adjacencies( vrange, 1, false, adj_edges );
            if( MB_SUCCESS != tmp_result ) result = tmp_result;
            if( adj_edges.size() == 1 )
            {
                // single edge - don't check adjacencies
                target_entities.push_back( *adj_edges.begin() );
            }
            else if( adj_edges.size() != 0 )
            {
                // multiple ones - need to check for explicit adjacencies
                unsigned int start_sz = target_entities.size();
                const EntityHandle* explicit_adjs;
                int num_exp;
                for( Range::iterator rit = adj_edges.begin(); rit != adj_edges.end(); ++rit )
                {
                    // TODO check return value
                    this->get_adjacencies( *rit, explicit_adjs, num_exp );
                    if( NULL != explicit_adjs &&
                        std::find( explicit_adjs, explicit_adjs + num_exp, source_entity ) != explicit_adjs + num_exp )
                        target_entities.push_back( *rit );
                }
                if( target_entities.size() == start_sz )
                {
                    result = MB_MULTIPLE_ENTITIES_FOUND;
                    target_entities.push_back( *adj_edges.begin() );
                }
            }
            else
            {
                // we have no adjacent edge yet; we need to create one and also add
                // them to the adjacency of the vertices
                if( create_if_missing )
                {
                    EntityHandle newEdge;
                    EntityHandle v[2] = { vertex_array[i], vertex_array[i + 1] };
                    result            = thisMB->create_element( MBEDGE, v, 2, newEdge );
                    if( MB_SUCCESS != result ) return result;
                    // we also need to add explicit adjacency, so next time we do not
                    // create again (because we do not find the edge if it is not adjacent to the
                    // vertices
                    // if (create_adjacency_option >= 0)
                    //{
                    result = add_adjacency( v[0], newEdge );
                    if( MB_SUCCESS != result ) return result;
                    result = add_adjacency( v[1], newEdge );
                    if( MB_SUCCESS != result ) return result;
                    target_entities.push_back( newEdge );
                    //}
                }
            }
        }
        return result;
    }

    else
    {
        if( target_dimension == 2 )
        {
            result = thisMB->get_connectivity( &source_entity, 1, target_entities );
        }
        else
        {
            std::vector< EntityHandle > dum_vec;
            result = thisMB->get_connectivity( &source_entity, 1, dum_vec );
            if( MB_SUCCESS != result ) return result;
            result = thisMB->get_adjacencies( &dum_vec[0], dum_vec.size(), 1, create_if_missing, target_entities,
                                              Interface::UNION );
            return result;
        }
    }

    return MB_SUCCESS;
}

#if 0
// Do in-place set intersect of two *sorted* containers.
// First container is modifed.  Second is not.
// First container must allow assignment through iterators (in practice,
// must be a container that does not enforce ordering, such
// as std::vector, std::list, or a c-style array)
template <typename T1, typename T2>
static inline T1 intersect( T1 set1_begin, T1 set1_end,
                            T2 set2_begin, T2 set2_end )
{
  T1 set1_write = set1_begin;
  while (set1_begin != set1_end) {
    if (set2_begin == set2_end)
      return set1_write;
    while (*set2_begin < *set1_begin)
      if (++set2_begin == set2_end)
        return set1_write;
    if (!(*set1_begin < *set2_begin)) {
      *set1_write = *set1_begin;
      ++set1_write;
      ++set2_begin;
    }
    ++set1_begin;
  }
  return set1_write;
}


ErrorCode AEntityFactory::get_up_adjacency_elements(
                                   EntityHandle source_entity,
                                   const unsigned int target_dimension,
                                   std::vector<EntityHandle>& target_entities,
                                   const bool create_if_missing,
                                   const int option )
{
  ErrorCode rval;
  const std::vector<EntityHandle> *vtx_adj, *vtx2_adj;
  std::vector<EntityHandle> duplicates;

    // Handle ranges
  const size_t in_size = target_entities.size();
  const EntityType src_type = TYPE_FROM_HANDLE(source_entity);
  DimensionPair target_types = CN::TypeDimensionMap[target_dimension];
  const EntityHandle src_beg_handle = CREATE_HANDLE( src_type, 0 );
  const EntityHandle src_end_handle = CREATE_HANDLE( src_type+1, 0 );
  const EntityHandle tgt_beg_handle = CREATE_HANDLE( target_types.first, 0 );
  const EntityHandle tgt_end_handle = CREATE_HANDLE( target_types.second+1, 0 );

    // get vertices
  assert(TYPE_FROM_HANDLE(source_entity) != MBPOLYHEDRON); // can't go up from a region
  std::vector<EntityHandle> conn_storage;
  const EntityHandle* conn;
  int conn_len;
  rval = thisMB->get_connectivity( source_entity, conn, conn_len, true, &conn_storage );
  if (MB_SUCCESS != rval)
    return rval;

    // shouldn't be here if source entity is not an element
  assert(conn_len > 1);

    // create necessary entities. this only makes sense if there exists of a
    // dimension greater than the target dimension.
  if (create_if_missing && target_dimension < 3 && CN::Dimension(src_type) < 2) {
    for (size_t i = 0; i < conn_len; ++i) {
      rval = get_adjacency_ptr( conn[i], vtx_adj );
      if (MB_SUCCESS != rval)
        return rval;
      assert(vtx_adj != NULL); // should contain at least source_entity

      std::vector<EntityHandle> tmp2, tmp(*vtx_adj); // copy in case adjacency vector is changed
      for (size_t j = 0; j < tmp.size(); ++j) {
        if (CN::Dimension(TYPE_FROM_HANDLE(tmp[j])) <= (int)target_dimension)
          continue;
        if (TYPE_FROM_HANDLE(tmp[j]) == MBENTITYSET)
          break;

        tmp2.clear();
        rval = get_down_adjacency_elements( tmp[j], target_dimension, tmp2, true, option );
        if (MB_SUCCESS != rval)
          return rval;
      }
    }
  }

    // get elements adjacent to first vertex
  rval = get_adjacency_ptr( conn[0], vtx_adj );
  if (MB_SUCCESS != rval)
    return rval;
  assert(vtx_adj != NULL); // should contain at least source_entity
    // get elements adjacent to second vertex
  rval = get_adjacency_ptr( conn[1], vtx2_adj );
  if (MB_SUCCESS != rval)
    return rval;
  assert(vtx2_adj != NULL);

    // Put intersect of all entities except source entity with
    // the same type as the source entity in 'duplicates'
  std::vector<EntityHandle>::const_iterator it1, it2, end1, end2;
  it1 = std::lower_bound( vtx_adj->begin(), vtx_adj->end(), src_beg_handle );
  it2 = std::lower_bound( vtx2_adj->begin(), vtx2_adj->end(), src_beg_handle );
  end1 = std::lower_bound( it1, vtx_adj->end(), src_end_handle );
  end2 = std::lower_bound( it2, vtx2_adj->end(), src_end_handle );
  assert(end1 != it1); // should at least contain source entity
  duplicates.resize( end1 - it1 - 1 );
  std::vector<EntityHandle>::iterator ins = duplicates.begin();
  for (; it1 != end1; ++it1) {
    if (*it1 != source_entity) {
      *ins = *it1;
      ++ins;
    }
  }
  duplicates.erase( intersect( duplicates.begin(), duplicates.end(), it2, end2 ), duplicates.end() );

    // Append to input list any entities of the desired target dimension
  it1 = std::lower_bound( end1, vtx_adj->end(), tgt_beg_handle );
  it2 = std::lower_bound( end2, vtx2_adj->end(), tgt_beg_handle );
  end1 = std::lower_bound( it1, vtx_adj->end(), tgt_end_handle );
  end2 = std::lower_bound( it2, vtx2_adj->end(), tgt_end_handle );
  std::set_intersection( it1, end1, it2, end2, std::back_inserter( target_entities ) );

    // for each additional vertex
  for (int i = 2; i < conn_len; ++i) {
    rval = get_adjacency_ptr( conn[i], vtx_adj );
    if (MB_SUCCESS != rval)
      return rval;
    assert(vtx_adj != NULL); // should contain at least source_entity

    it1 = std::lower_bound( vtx_adj->begin(), vtx_adj->end(), src_beg_handle );
    end1 = std::lower_bound( it1, vtx_adj->end(), src_end_handle );
    duplicates.erase( intersect( duplicates.begin(), duplicates.end(), it1, end1 ), duplicates.end() );

    it1 = std::lower_bound( end1, vtx_adj->end(), tgt_beg_handle );
    end1 = std::lower_bound( it1, vtx_adj->end(), tgt_end_handle );
    target_entities.erase( intersect( target_entities.begin()+in_size, target_entities.end(),
                           it1, end1 ), target_entities.end() );
  }

    // if no duplicates, we're done
  if (duplicates.empty())
    return MB_SUCCESS;

    // Check for explicit adjacencies.  If an explicit adjacency
    // connects candidate target entity to an entity equivalent
    // to the source entity, then assume that source entity is *not*
    // adjacent
  const std::vector<EntityHandle>* adj_ptr;
    // check adjacencies from duplicate entities to candidate targets
  for (size_t i = 0; i < duplicates.size(); ++i) {
    rval = get_adjacency_ptr( duplicates[i], adj_ptr );
    if (MB_SUCCESS != rval)
      return rval;
    if (!adj_ptr)
      continue;

    for (size_t j = 0; j < adj_ptr->size(); ++j) {
      std::vector<EntityHandle>::iterator k =
        std::find( target_entities.begin()+in_size, target_entities.end(), (*adj_ptr)[j] );
      if (k != target_entities.end())
        target_entities.erase(k);
    }
  }

  // If target dimension is 3 and source dimension is 1, also need to
  // check for explicit adjacencies to intermediate faces
  if (CN::Dimension(src_type) > 1 || target_dimension < 3)
    return MB_SUCCESS;

    // Get faces adjacent to each element and check for explict
    // adjacencies from duplicate entities to faces
  for (size_t i = 0; i < duplicates.size(); ++i) {
    rval = get_adjacency_ptr( duplicates[i], adj_ptr );
    if (MB_SUCCESS != rval)
      return rval;
    if (!adj_ptr)
      continue;

    size_t j;
    for (j = 0; j < adj_ptr->size(); ++j) {
      const std::vector<EntityHandle>* adj_ptr2;
      rval = get_adjacency_ptr( (*adj_ptr)[j], adj_ptr2 );
      if (MB_SUCCESS != rval)
        return rval;
      if (!adj_ptr2)
        continue;

      for (size_t k = 0; k < adj_ptr2->size(); ++k) {
        std::vector<EntityHandle>::iterator it;
        it = std::find( target_entities.begin()+in_size, target_entities.end(), (*adj_ptr2)[k] );
        if (it != target_entities.end()) {
          target_entities.erase(it);
          j = adj_ptr->size(); // break outer loop
          break;
        }
      }
    }
  }

  return MB_SUCCESS;
}
#else
ErrorCode AEntityFactory::get_up_adjacency_elements( EntityHandle source_entity,
                                                     const unsigned int target_dimension,
                                                     std::vector< EntityHandle >& target_entities,
                                                     const bool create_if_missing,
                                                     const int /*create_adjacency_option = -1*/ )
{

    EntityType source_type = TYPE_FROM_HANDLE( source_entity );

    const EntityHandle* source_vertices = NULL;
    int num_source_vertices             = 0;
    std::vector< EntityHandle > conn_storage;

    // check to see whether there are any equivalent entities (same verts, different entity);
    // do this by calling get_element with a 0 source_entity, and look for a
    // MB_MULTIPLE_ENTITIES_FOUND return code

    // NOTE: we only want corner vertices here, and for the code below which also uses
    // source_vertices
    ErrorCode result =
        thisMB->get_connectivity( source_entity, source_vertices, num_source_vertices, true, &conn_storage );
    if( MB_SUCCESS != result ) return result;
    EntityHandle temp_entity;
    result = get_element( source_vertices, num_source_vertices, source_type, temp_entity, false, 0 );

    bool equiv_entities = ( result == MB_MULTIPLE_ENTITIES_FOUND ) ? true : false;

    std::vector< EntityHandle > tmp_vec;
    if( !equiv_entities )
    {
        // get elems adjacent to each node
        std::vector< std::vector< EntityHandle > > elems( num_source_vertices );
        int i;
        for( i = 0; i < num_source_vertices; i++ )
        {
            // get elements
            // see comment above pertaining to source_vertices; these are corner vertices only
            get_zero_to_n_elements( source_vertices[i], target_dimension, elems[i], create_if_missing, 0 );
            // sort this element list
            std::sort( elems[i].begin(), elems[i].end() );
        }

        // perform an intersection between all the element lists
        // see comment above pertaining to source_vertices; these are corner vertices only
        for( i = 1; i < num_source_vertices; i++ )
        {
            tmp_vec.clear();

            // intersection between first list and ith list, put result in tmp
            std::set_intersection( elems[0].begin(), elems[0].end(), elems[i].begin(), elems[i].end(),
                                   std::back_insert_iterator< std::vector< EntityHandle > >( tmp_vec ) );
            // tmp has elems[0] contents and elems[0] contents has tmp's contents
            // so that elems[0] always has the intersection of previous operations
            elems[0].swap( tmp_vec );
        }

        // elems[0] contains the intersection, swap with target_entities
        target_entities.insert( target_entities.end(), elems[0].begin(), elems[0].end() );
    }
    else if( source_type == MBPOLYGON )
    {
        // get adjacencies using polyhedra's connectivity vectors
        // first get polyhedra neighboring vertices
        result = thisMB->get_adjacencies( source_vertices, num_source_vertices, 3, false, tmp_vec );
        if( MB_SUCCESS != result ) return result;

        // now filter according to whether each is adjacent to the polygon
        const EntityHandle* connect = NULL;
        int num_connect             = 0;
        std::vector< EntityHandle > storage;
        for( unsigned int i = 0; i < tmp_vec.size(); i++ )
        {
            result = thisMB->get_connectivity( tmp_vec[i], connect, num_connect, false, &storage );
            if( MB_SUCCESS != result ) return result;
            if( std::find( connect, connect + num_connect, source_entity ) != connect + num_connect )
                target_entities.push_back( tmp_vec[i] );
        }
    }

    else
    {
        // else get up-adjacencies directly; code copied from get_zero_to_n_elements

        // get the adjacency vector
        AdjacencyVector* adj_vec = NULL;
        result                   = get_adjacencies( source_entity, adj_vec );

        if( result != MB_SUCCESS )
            return result;
        else if( adj_vec == NULL )
            return MB_SUCCESS;

        DimensionPair dim_pair_dp1 = CN::TypeDimensionMap[CN::Dimension( source_type ) + 1],
                      dim_pair_td  = CN::TypeDimensionMap[target_dimension];
        int dum;

        Range tmp_ents, target_ents;

        // get iterators for start handle of source_dim+1 and target_dim, and end handle
        // of target_dim
        AdjacencyVector::iterator start_ent_dp1 =
                                      std::lower_bound( adj_vec->begin(), adj_vec->end(),
                                                        CREATE_HANDLE( dim_pair_dp1.first, MB_START_ID, dum ) ),

                                  start_ent_td =
                                      std::lower_bound( adj_vec->begin(), adj_vec->end(),
                                                        CREATE_HANDLE( dim_pair_td.first, MB_START_ID, dum ) ),

                                  end_ent_td = std::lower_bound( adj_vec->begin(), adj_vec->end(),
                                                                 CREATE_HANDLE( dim_pair_td.second, MB_END_ID, dum ) );

        // get the adjacencies for source_dim+1 to target_dim-1, and the adjacencies from
        // those to target_dim
        std::copy( start_ent_dp1, start_ent_td, range_inserter( tmp_ents ) );
        result = thisMB->get_adjacencies( tmp_ents, target_dimension, false, target_ents, Interface::UNION );
        if( MB_SUCCESS != result ) return result;

        // now copy the explicit adjacencies to target_dimension
        std::copy( start_ent_td, end_ent_td, range_inserter( target_ents ) );

        // now insert the whole thing into the argument vector
#ifdef MOAB_NO_VECTOR_TEMPLATE_INSERT
        std::copy( target_ents.begin(), target_ents.end(), std::back_inserter( target_entities ) );
#else
        target_entities.insert( target_entities.end(), target_ents.begin(), target_ents.end() );
#endif
    }

    return result;
}
#endif

ErrorCode AEntityFactory::notify_change_connectivity( EntityHandle entity,
                                                      const EntityHandle* old_array,
                                                      const EntityHandle* new_array,
                                                      int number_verts )
{
    EntityType source_type = TYPE_FROM_HANDLE( entity );
    if( source_type == MBPOLYHEDRON ) return MB_NOT_IMPLEMENTED;

    // find out which ones to add and which to remove
    std::vector< EntityHandle > old_verts, new_verts;
    int i;
    for( i = 0; i < number_verts; i++ )
    {
        if( old_array[i] != new_array[i] )
        {
            old_verts.push_back( old_array[i] );
            new_verts.push_back( new_array[i] );
        }
    }

    ErrorCode result;

    if( mVertElemAdj == true )
    {
        // update the vertex-entity adjacencies
        std::vector< EntityHandle >::iterator adj_iter;
        for( adj_iter = old_verts.begin(); adj_iter != old_verts.end(); ++adj_iter )
        {
            if( std::find( new_verts.begin(), new_verts.end(), *adj_iter ) == new_verts.end() )
            {
                result = remove_adjacency( *adj_iter, entity );
                if( MB_SUCCESS != result ) return result;
            }
        }
        for( adj_iter = new_verts.begin(); adj_iter != new_verts.end(); ++adj_iter )
        {
            if( std::find( old_verts.begin(), old_verts.end(), *adj_iter ) == old_verts.end() )
            {
                result = add_adjacency( *adj_iter, entity );
                if( MB_SUCCESS != result ) return result;
            }
        }
    }

    return MB_SUCCESS;
}

//! return true if 2 entities are explicitly adjacent
bool AEntityFactory::explicitly_adjacent( const EntityHandle ent1, const EntityHandle ent2 )
{
    const EntityHandle* explicit_adjs;
    int num_exp;
    get_adjacencies( ent1, explicit_adjs, num_exp );
    if( std::find( explicit_adjs, explicit_adjs + num_exp, ent2 ) != explicit_adjs + num_exp )
        return true;
    else
        return false;
}

ErrorCode AEntityFactory::merge_adjust_adjacencies( EntityHandle entity_to_keep, EntityHandle entity_to_remove )
{
    int ent_dim = CN::Dimension( TYPE_FROM_HANDLE( entity_to_keep ) );
    ErrorCode result;

    // check for newly-formed equivalent entities, and create explicit adjacencies
    // to distinguish them; this must be done before connectivity of higher-dimensional
    // entities is changed below, and only needs to be checked if merging vertices
    if( ent_dim == 0 )
    {
        result = check_equiv_entities( entity_to_keep, entity_to_remove );
        if( MB_SUCCESS != result ) return result;
    }

    // check adjacencies TO removed entity
    for( int dim = 1; dim < ent_dim; dim++ )
    {
        Range adjs;
        result = thisMB->get_adjacencies( &entity_to_remove, 1, dim, false, adjs );
        if( result != MB_SUCCESS ) return result;
        // for any explicit ones, make them adjacent to keeper
        for( Range::iterator rit = adjs.begin(); rit != adjs.end(); ++rit )
        {
            if( this->explicitly_adjacent( *rit, entity_to_remove ) )
            {
                result = this->add_adjacency( *rit, entity_to_keep );
                if( result != MB_SUCCESS ) return result;
            }
        }
    }

    // check adjacencies FROM removed entity
    std::vector< EntityHandle > conn, adjs;
    result = this->get_adjacencies( entity_to_remove, adjs );
    if( result != MB_SUCCESS ) return result;
    // set them all, and if to_entity is a set, add to that one too
    for( unsigned int i = 0; i < adjs.size(); i++ )
    {
        if( TYPE_FROM_HANDLE( adjs[i] ) == MBENTITYSET )
        {
            // result = this->add_adjacency(entity_to_keep, adjs[i]);
            // if(result != MB_SUCCESS) return result;
            // result = thisMB->add_entities(adjs[i], &entity_to_keep, 1);
            // if(result != MB_SUCCESS) return result;
            result = thisMB->replace_entities( adjs[i], &entity_to_remove, &entity_to_keep, 1 );
            if( MB_SUCCESS != result ) return result;
        }
        else if( ent_dim == 0 )
        {
            conn.clear();
            result = thisMB->get_connectivity( &adjs[i], 1, conn );

            if( result == MB_SUCCESS )
            {
                std::replace( conn.begin(), conn.end(), entity_to_remove, entity_to_keep );
                result = thisMB->set_connectivity( adjs[i], &conn[0], conn.size() );
                if( MB_SUCCESS != result ) return result;
            }
            else
                return result;
        }
        else
        {
            result = this->add_adjacency( entity_to_keep, adjs[i] );
            if( result != MB_SUCCESS ) return result;
        }
    }

    return MB_SUCCESS;
}

// check for equivalent entities that may be formed when merging two entities, and
// create explicit adjacencies accordingly
ErrorCode AEntityFactory::check_equiv_entities( EntityHandle entity_to_keep, EntityHandle entity_to_remove )
{
    if( thisMB->dimension_from_handle( entity_to_keep ) > 0 ) return MB_SUCCESS;

    // get all the adjacencies for both entities for all dimensions > 0
    Range adjs_keep, adjs_remove;
    ErrorCode result;

    for( int dim = 1; dim <= 3; dim++ )
    {
        result = thisMB->get_adjacencies( &entity_to_keep, 1, dim, false, adjs_keep, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
        result = thisMB->get_adjacencies( &entity_to_remove, 1, dim, false, adjs_remove, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
    }

    // now look for equiv entities which will be formed
    // algorithm:
    // for each entity adjacent to removed entity:
    EntityHandle two_ents[2];
    for( Range::iterator rit_rm = adjs_remove.begin(); rit_rm != adjs_remove.end(); ++rit_rm )
    {
        two_ents[0] = *rit_rm;

        // - for each entity of same dimension adjacent to kept entity:
        for( Range::iterator rit_kp = adjs_keep.begin(); rit_kp != adjs_keep.end(); ++rit_kp )
        {
            if( TYPE_FROM_HANDLE( *rit_kp ) != TYPE_FROM_HANDLE( *rit_rm ) ) continue;

            Range all_verts;
            two_ents[1] = *rit_kp;
            //   . get union of adjacent vertices to two entities
            result = thisMB->get_adjacencies( two_ents, 2, 0, false, all_verts, Interface::UNION );
            if( MB_SUCCESS != result ) return result;

            assert( all_verts.find( entity_to_keep ) != all_verts.end() &&
                    all_verts.find( entity_to_remove ) != all_verts.end() );

            //   . if # vertices != number of corner vertices + 1, continue
            if( CN::VerticesPerEntity( TYPE_FROM_HANDLE( *rit_rm ) ) + 1 != (int)all_verts.size() ) continue;

            //   . for the two entities adjacent to kept & removed entity:
            result = create_explicit_adjs( *rit_rm );
            if( MB_SUCCESS != result ) return result;
            result = create_explicit_adjs( *rit_kp );
            if( MB_SUCCESS != result ) return result;
            //   . (end for)
        }
        // - (end for)
    }

    return MB_SUCCESS;
}

ErrorCode AEntityFactory::create_explicit_adjs( EntityHandle this_ent )
{
    //     - get adjacent entities of next higher dimension
    Range all_adjs;
    ErrorCode result;
    result = thisMB->get_adjacencies( &this_ent, 1, thisMB->dimension_from_handle( this_ent ) + 1, false, all_adjs,
                                      Interface::UNION );
    if( MB_SUCCESS != result ) return result;

    //     - create explicit adjacency to these entities
    for( Range::iterator rit = all_adjs.begin(); rit != all_adjs.end(); ++rit )
    {
        result = add_adjacency( this_ent, *rit );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

ErrorCode AEntityFactory::get_adjacency_ptr( EntityHandle entity, std::vector< EntityHandle >*& ptr )
{
    ptr = 0;

    EntitySequence* seq;
    ErrorCode rval = thisMB->sequence_manager()->find( entity, seq );
    if( MB_SUCCESS != rval || !seq->data()->get_adjacency_data() ) return rval;

    ptr = seq->data()->get_adjacency_data()[entity - seq->data()->start_handle()];
    return MB_SUCCESS;
}

ErrorCode AEntityFactory::get_adjacency_ptr( EntityHandle entity, const std::vector< EntityHandle >*& ptr ) const
{
    ptr = 0;

    EntitySequence* seq;
    ErrorCode rval = thisMB->sequence_manager()->find( entity, seq );
    if( MB_SUCCESS != rval || !seq->data()->get_adjacency_data() ) return rval;

    ptr = seq->data()->get_adjacency_data()[entity - seq->data()->start_handle()];
    return MB_SUCCESS;
}

ErrorCode AEntityFactory::set_adjacency_ptr( EntityHandle entity, std::vector< EntityHandle >* ptr )
{
    EntitySequence* seq;
    ErrorCode rval = thisMB->sequence_manager()->find( entity, seq );
    if( MB_SUCCESS != rval ) return rval;

    if( !seq->data()->get_adjacency_data() && !seq->data()->allocate_adjacency_data() )
        return MB_MEMORY_ALLOCATION_FAILED;

    const EntityHandle index          = entity - seq->data()->start_handle();
    std::vector< EntityHandle >*& ref = seq->data()->get_adjacency_data()[index];
    delete ref;
    ref = ptr;
    return MB_SUCCESS;
}

void AEntityFactory::get_memory_use( unsigned long long& entity_total, unsigned long long& memory_total )
{
    entity_total = memory_total = 0;

    // iterate through each element type
    SequenceData* prev_data = 0;
    for( EntityType t = MBVERTEX; t != MBENTITYSET; t++ )
    {
        TypeSequenceManager::iterator i;
        TypeSequenceManager& seqman = thisMB->sequence_manager()->entity_map( t );
        for( i = seqman.begin(); i != seqman.end(); ++i )
        {
            if( !( *i )->data()->get_adjacency_data() ) continue;

            if( prev_data != ( *i )->data() )
            {
                prev_data = ( *i )->data();
                memory_total += prev_data->size() * sizeof( AdjacencyVector* );
            }

            const AdjacencyVector* vec;
            for( EntityHandle h = ( *i )->start_handle(); h <= ( *i )->end_handle(); ++h )
            {
                get_adjacency_ptr( h, vec );
                if( vec ) entity_total += vec->capacity() * sizeof( EntityHandle ) + sizeof( AdjacencyVector );
            }
        }
    }

    memory_total += sizeof( *this ) + entity_total;
}

ErrorCode AEntityFactory::get_memory_use( const Range& ents_in,
                                          unsigned long long& min_per_ent,
                                          unsigned long long& amortized )
{
    min_per_ent = amortized = 0;
    SequenceData* prev_data = 0;
    RangeSeqIntersectIter iter( thisMB->sequence_manager() );
    ErrorCode rval = iter.init( ents_in.begin(), ents_in.end() );
    if( MB_SUCCESS != rval ) return rval;

    do
    {
        AdjacencyVector** array = iter.get_sequence()->data()->get_adjacency_data();
        if( !array ) continue;

        EntityID count    = iter.get_end_handle() - iter.get_start_handle() + 1;
        EntityID data_occ = thisMB->sequence_manager()
                                ->entity_map( iter.get_sequence()->type() )
                                .get_occupied_size( iter.get_sequence()->data() );

        if( iter.get_sequence()->data() != prev_data )
        {
            prev_data = iter.get_sequence()->data();
            amortized += sizeof( AdjacencyVector* ) * iter.get_sequence()->data()->size() * count / data_occ;
        }

        array += iter.get_start_handle() - iter.get_sequence()->data()->start_handle();
        for( EntityID i = 0; i < count; ++i )
        {
            if( array[i] ) min_per_ent += sizeof( EntityHandle ) * array[i]->capacity() + sizeof( AdjacencyVector );
        }
    } while( MB_SUCCESS == ( rval = iter.step() ) );

    amortized += min_per_ent;
    return ( rval == MB_FAILURE ) ? MB_SUCCESS : rval;
}

/*!
   calling code is notifying this that an entity is going to be deleted
   from the database
*/
ErrorCode AEntityFactory::notify_delete_entity( EntityHandle entity )
{
    if( TYPE_FROM_HANDLE( entity ) == MBVERTEX )
    {
        std::vector< EntityHandle > adj_entities;
        for( int dim = 1; dim < 4; ++dim )
        {
            ErrorCode rval = get_adjacencies( entity, dim, false, adj_entities );
            if( rval != MB_SUCCESS && rval != MB_ENTITY_NOT_FOUND ) return rval;
            if( !adj_entities.empty() ) return MB_FAILURE;
        }
    }

    // remove any references to this entity from other entities
    return remove_all_adjacencies( entity, true );
}

}  // namespace moab