NestedRefine.cpp

Go to the documentation of this file.

#include "moab/NestedRefine.hpp"
#include "moab/NestedRefineTemplates.hpp"
#include "moab/HalfFacetRep.hpp"
#include "moab/CpuTimer.hpp"
#include "moab/ReadUtilIface.hpp"
#include "Internals.hpp"
#include "MBTagConventions.hpp"
 
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "moab/ParallelMergeMesh.hpp"
#include "moab/Skinner.hpp"
#endif
 
#include <iostream>
#include <cassert>
#include <vector>
#include <limits>
#include <cmath>
 
namespace moab
{
 
NestedRefine::NestedRefine( Core* impl, ParallelComm* comm, EntityHandle rset )
 : mbImpl( impl ), pcomm( comm ), _rset( rset )
{
 ErrorCode error;
 assert( NULL != impl );
 
#ifdef MOAB_HAVE_MPI
 // Get the Parallel Comm instance to prepare all new sets to work in parallel
 // in case the user did not provide any arguments
 if( !comm ) pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 );
#endif
 error = initialize();
 if( error != MB_SUCCESS )
 {
 std::cout << "Error initializing NestedRefine\n" << std::endl;
 exit( 1 );
 }
}
 
NestedRefine::~NestedRefine()
{
#ifdef MOAB_HAVE_AHF
 ahf = NULL;
#else
 delete ahf;
#endif
 delete tm;
}
 
ErrorCode NestedRefine::initialize()
{
 ErrorCode error;
 
 tm = new CpuTimer();
 if( !tm ) return MB_MEMORY_ALLOCATION_FAILED;
 
#ifdef MOAB_HAVE_AHF
 ahf = mbImpl->a_half_facet_rep();
#else
 ahf = new HalfFacetRep( mbImpl, pcomm, _rset, true );
 if( !ahf ) return MB_MEMORY_ALLOCATION_FAILED;
#endif
 
 // Check for mixed entity type
 bool chk_mixed = ahf->check_mixed_entity_type();
 if( chk_mixed ) MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mesh with mixed entity types" );
 
 error = ahf->initialize();MB_CHK_ERR( error );
 error = ahf->get_entity_ranges( _inverts, _inedges, _infaces, _incells );MB_CHK_ERR( error );
 
 // Check for supported entity type
 if( !_incells.empty() )
 {
 EntityType type = mbImpl->type_from_handle( _incells[0] );
 if( type != MBTET && type != MBHEX )
 MB_SET_ERR( MB_FAILURE, "Not supported 3D entity types: MBPRISM, MBPYRAMID, MBKNIFE, MBPOLYHEDRON" );
 
 meshdim = 3;
 elementype = type;
 }
 else if( !_infaces.empty() )
 {
 EntityType type = mbImpl->type_from_handle( _infaces[0] );
 if( type == MBPOLYGON ) MB_SET_ERR( MB_FAILURE, "Not supported 2D entity type: POLYGON" );
 
 meshdim = 2;
 elementype = type;
 }
 else if( !_inedges.empty() )
 {
 meshdim = 1;
 elementype = MBEDGE;
 }
 else
 MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mixed-dimensional or invalid mesh" );
 
 // Initialize std::map to get indices of degrees.
 deg_index[2] = 0;
 deg_index[3] = 1;
 deg_index[5] = 2;
 
 // Set ghost flag to false
 hasghost = false;
 return MB_SUCCESS;
}
 
/************************************************************
 * Interface Functions *
 ************************************************************/
 
ErrorCode NestedRefine::generate_mesh_hierarchy( int num_level,
 int* level_degrees,
 std::vector< EntityHandle >& level_sets,
 bool optimize )
{
 assert( num_level > 0 );
 nlevels = num_level;
 
 ErrorCode error;
 std::vector< moab::EntityHandle > hmsets( num_level );
 
 if( meshdim <= 2 )
 {
 for( int i = 0; i < num_level; i++ )
 {
 assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) || ( level_degrees[i] == 5 ) );
 level_dsequence[i] = level_degrees[i];
 }
 }
 else
 {
 for( int i = 0; i < num_level; i++ )
 {
 assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) );
 level_dsequence[i] = level_degrees[i];
 }
 }
 
 error = generate_hm( level_degrees, num_level, &hmsets[0], optimize );MB_CHK_ERR( error );
 
 // copy the entity handles
 level_sets.resize( num_level + 1 );
 level_sets[0] = _rset;
 for( int i = 0; i < num_level; i++ )
 level_sets[i + 1] = hmsets[i];
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn )
{
 ErrorCode error;
 EntityType type = mbImpl->type_from_handle( ent );
 EntityHandle start_ent;
 if( !conn.empty() ) conn.clear();
 if( level > 0 )
 {
 if( type == MBEDGE )
 {
 conn.reserve( 2 );
 start_ent = level_mesh[level - 1].start_edge;
 EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
 conn.push_back( level_mesh[level - 1].edge_conn[2 * offset] );
 conn.push_back( level_mesh[level - 1].edge_conn[2 * offset + 1] );
 }
 else if( type == MBTRI || type == MBQUAD )
 {
 int num_corners = ahf->lConnMap2D[type - 2].num_verts_in_face;
 conn.reserve( num_corners );
 start_ent = level_mesh[level - 1].start_face;
 EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
 
 for( int i = 0; i < num_corners; i++ )
 conn.push_back( level_mesh[level - 1].face_conn[num_corners * offset + i] );
 }
 else if( type == MBTET || type == MBHEX )
 {
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 int num_corners = ahf->lConnMap3D[index].num_verts_in_cell;
 conn.reserve( num_corners );
 start_ent = level_mesh[level - 1].start_cell;
 EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
 for( int i = 0; i < num_corners; i++ )
 conn.push_back( level_mesh[level - 1].cell_conn[num_corners * offset + i] );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting connectivity for an unsupported entity type" );
 }
 else
 {
 error = mbImpl->get_connectivity( &ent, 1, conn );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords )
{
 if( level > 0 )
 {
 EntityID vstart = ID_FROM_HANDLE( level_mesh[level - 1].start_vertex );
 for( int i = 0; i < num_verts; i++ )
 {
 const EntityHandle& vid = verts[i];
 EntityID offset = ID_FROM_HANDLE( vid ) - vstart;
 coords[3 * i] = level_mesh[level - 1].coordinates[0][offset];
 coords[3 * i + 1] = level_mesh[level - 1].coordinates[1][offset];
 coords[3 * i + 2] = level_mesh[level - 1].coordinates[2][offset];
 }
 }
 else
 {
 ErrorCode error;
 error = mbImpl->get_coords( verts, num_verts, coords );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_adjacencies( const EntityHandle source_entity,
 const unsigned int target_dimension,
 std::vector< EntityHandle >& target_entities )
 
{
 ErrorCode error;
 error = ahf->get_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent )
{
 assert( ( child_level > 0 ) && ( child_level > parent_level ) );
 EntityType type = mbImpl->type_from_handle( child );
 assert( type != MBVERTEX );
 
 int child_index;
 if( type == MBEDGE )
 child_index = child - level_mesh[child_level - 1].start_edge;
 else if( type == MBTRI || type == MBQUAD )
 child_index = child - level_mesh[child_level - 1].start_face;
 else if( type == MBTET || type == MBHEX )
 child_index = child - level_mesh[child_level - 1].start_cell;
 else
 MB_SET_ERR( MB_FAILURE, "Requesting parent for unsupported entity type" );
 
 int parent_index;
 int l = child_level - parent_level;
 for( int i = 0; i < l; i++ )
 {
 int d = get_index_from_degree( level_dsequence[child_level - i - 1] );
 int nch = refTemplates[type - 1][d].total_new_ents;
 child_index = child_index / nch;
 }
 parent_index = child_index;
 
 if( type == MBEDGE )
 {
 if( parent_level > 0 )
 *parent = level_mesh[parent_level - 1].start_edge + parent_index;
 else
 *parent = _inedges[parent_index];
 }
 else if( type == MBTRI || type == MBQUAD )
 {
 if( parent_level > 0 )
 *parent = level_mesh[parent_level - 1].start_face + parent_index;
 else
 *parent = _infaces[parent_index];
 }
 else if( type == MBTET || type == MBHEX )
 {
 if( parent_level > 0 )
 *parent = level_mesh[parent_level - 1].start_cell + parent_index;
 else
 *parent = _incells[parent_index];
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::parent_to_child( EntityHandle parent,
 int parent_level,
 int child_level,
 std::vector< EntityHandle >& children )
{
 assert( ( child_level > 0 ) && ( child_level > parent_level ) );
 EntityType type = mbImpl->type_from_handle( parent );
 assert( type != MBVERTEX );
 
 int parent_index;
 if( type == MBEDGE )
 {
 if( parent_level > 0 )
 parent_index = parent - level_mesh[parent_level - 1].start_edge;
 else
 parent_index = _inedges.index( parent );
 }
 else if( type == MBTRI || type == MBQUAD )
 {
 if( parent_level > 0 )
 parent_index = parent - level_mesh[parent_level - 1].start_face;
 else
 parent_index = _infaces.index( parent );
 }
 else if( type == MBTET || type == MBHEX )
 {
 if( parent_level > 0 )
 parent_index = parent - level_mesh[parent_level - 1].start_cell;
 else
 parent_index = _incells.index( parent );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting children for unsupported entity type" );
 
 int start, end;
 start = end = parent_index;
 for( int i = parent_level; i < child_level; i++ )
 {
 int d = get_index_from_degree( level_dsequence[i] );
 int nch = refTemplates[type - 1][d].total_new_ents;
 start = start * nch;
 end = end * nch + nch - 1;
 }
 
 int num_child = end - start;
 children.reserve( num_child );
 
 for( int i = start; i <= end; i++ )
 {
 EntityHandle child;
 if( type == MBEDGE )
 child = level_mesh[child_level - 1].start_edge + i;
 else if( type == MBTRI || type == MBQUAD )
 child = level_mesh[child_level - 1].start_face + i;
 else if( type == MBTET || type == MBHEX )
 child = level_mesh[child_level - 1].start_cell + i;
 
 children.push_back( child );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::vertex_to_entities_up( EntityHandle vertex,
 int vert_level,
 int parent_level,
 std::vector< EntityHandle >& incident_entities )
{
 assert( vert_level > parent_level );
 ErrorCode error;
 
 // Step 1: Get the incident entities at the current level
 std::vector< EntityHandle > inents;
 if( meshdim == 1 )
 {
 error = ahf->get_up_adjacencies_1d( vertex, inents );MB_CHK_ERR( error );
 }
 else if( meshdim == 2 )
 {
 error = ahf->get_up_adjacencies_vert_2d( vertex, inents );MB_CHK_ERR( error );
 }
 else if( meshdim == 3 )
 {
 error = ahf->get_up_adjacencies_vert_3d( vertex, inents );MB_CHK_ERR( error );
 }
 
 // Step 2: Loop over all the incident entities at the current level and gather their parents
 for( int i = 0; i < (int)inents.size(); i++ )
 {
 EntityHandle ent = inents[i];
 EntityHandle parent;
 error = child_to_parent( ent, vert_level, parent_level, &parent );MB_CHK_ERR( error );
 incident_entities.push_back( parent );
 }
 
 // Step 3: Sort and remove duplicates
 std::sort( incident_entities.begin(), incident_entities.end() );
 incident_entities.erase( std::unique( incident_entities.begin(), incident_entities.end() ),
 incident_entities.end() );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::vertex_to_entities_down( EntityHandle vertex,
 int vert_level,
 int child_level,
 std::vector< EntityHandle >& incident_entities )
{
 assert( vert_level < child_level );
 ErrorCode error;
 
 // Step 1: Get the incident entities at the current level
 std::vector< EntityHandle > inents;
 if( meshdim == 1 )
 {
 error = ahf->get_up_adjacencies_1d( vertex, inents );MB_CHK_ERR( error );
 }
 else if( meshdim == 2 )
 {
 error = ahf->get_up_adjacencies_vert_2d( vertex, inents );MB_CHK_ERR( error );
 }
 else if( meshdim == 3 )
 {
 error = ahf->get_up_adjacencies_vert_3d( vertex, inents );MB_CHK_ERR( error );
 }
 
 // Step 2: Loop over all the incident entities at the current level and gather their parents
 std::vector< EntityHandle > childs;
 for( int i = 0; i < (int)inents.size(); i++ )
 {
 childs.clear();
 EntityHandle ent = inents[i];
 error = parent_to_child( ent, vert_level, child_level, childs );MB_CHK_ERR( error );
 for( int j = 0; j < (int)childs.size(); j++ )
 incident_entities.push_back( childs[j] );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_vertex_duplicates( EntityHandle vertex, int level, EntityHandle& dupvertex )
{
 if( ( vertex - *_inverts.begin() ) > _inverts.size() )
 MB_SET_ERR( MB_FAILURE, "Requesting duplicates for non-coarse vertices" );
 
 dupvertex = level_mesh[level - 1].start_vertex + ( vertex - *_inverts.begin() );
 
 return MB_SUCCESS;
}
 
bool NestedRefine::is_entity_on_boundary( const EntityHandle& entity )
{
 bool is_border = false;
 EntityType type = mbImpl->type_from_handle( entity );
 
 if( type == MBVERTEX )
 is_border = is_vertex_on_boundary( entity );
 else if( type == MBEDGE )
 is_border = is_edge_on_boundary( entity );
 else if( type == MBTRI || type == MBQUAD )
 is_border = is_face_on_boundary( entity );
 else if( type == MBTET || type == MBHEX )
 is_border = is_cell_on_boundary( entity );
 else
 MB_SET_ERR( MB_FAILURE, "Requesting boundary information for unsupported entity type" );
 
 return is_border;
}
 
ErrorCode NestedRefine::exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers )
{
 ErrorCode error;
 
 if( hasghost ) return MB_SUCCESS;
 
 hasghost = true;
#ifdef MOAB_HAVE_MPI
 error = pcomm->exchange_ghost_cells( meshdim, 0, num_glayers, 0, true, false );MB_CHK_ERR( error );
 {
 Range empty_range;
 error = pcomm->exchange_tags( GLOBAL_ID_TAG_NAME, empty_range );MB_CHK_ERR( error );
 // error = pcomm->assign_global_ids(lsets[i], 0, 1, false, true, false);MB_CHK_ERR(error);
 }
#else
 MB_SET_ERR( MB_FAILURE, "Requesting ghost layers for a serial mesh" );
#endif
 
 Range* lverts = new Range[lsets.size()];
 Range* lents = new Range[lsets.size()];
 for( size_t i = 0; i < lsets.size(); i++ )
 {
 error = mbImpl->get_entities_by_dimension( lsets[i], meshdim, lents[i] );MB_CHK_ERR( error );
 error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
 
 for( int gl = 0; gl < num_glayers; gl++ )
 {
 error = mbImpl->get_adjacencies( lverts[i], meshdim, false, lents[i], Interface::UNION );MB_CHK_ERR( error );
 error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
 }
 }
 for( size_t i = 0; i < lsets.size(); i++ )
 {
 error = mbImpl->add_entities( lsets[i], lverts[i] );MB_CHK_ERR( error );
 error = mbImpl->add_entities( lsets[i], lents[i] );MB_CHK_ERR( error );
 }
 
 delete[] lverts;
 delete[] lents;
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_special_tags( int level, EntityHandle& lset )
{
 assert( level > 0 && level < nlevels + 1 );
 
 ErrorCode error;
 std::vector< Tag > mtags( 3 );
 
 error = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[0] );MB_CHK_ERR( error );
 error = mbImpl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[1] );MB_CHK_ERR( error );
 error = mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[2] );MB_CHK_ERR( error );
 
 for( int i = 0; i < 3; i++ )
 {
 // Gather sets of a particular tag
 Range sets;
 error = mbImpl->get_entities_by_type_and_tag( _rset, MBENTITYSET, &mtags[i], NULL, 1, sets );MB_CHK_ERR( error );
 
 // Loop over all sets, gather entities in each set and add their children at all levels to
 // the set
 Range set_ents;
 Range::iterator set_it;
 std::vector< EntityHandle > childs;
 
 for( set_it = sets.begin(); set_it != sets.end(); ++set_it )
 {
 // Get the entities in the set, recursively
 set_ents.clear();
 childs.clear();
 error = mbImpl->get_entities_by_handle( *set_it, set_ents, true );MB_CHK_ERR( error );
 
 // Gather child entities at the input level
 for( Range::iterator sit = set_ents.begin(); sit != set_ents.end(); sit++ )
 {
 EntityType type = mbImpl->type_from_handle( *sit );
 if( type == MBVERTEX )
 {
 Range conn;
 std::vector< EntityHandle > cents;
 error = vertex_to_entities_down( *sit, 0, level, cents );MB_CHK_ERR( error );
 error = mbImpl->get_connectivity( &cents[0], (int)cents.size(), conn, true );MB_CHK_ERR( error );
 childs.insert( childs.end(), cents.begin(), cents.end() );
 }
 else
 {
 error = parent_to_child( *sit, 0, level, childs );MB_CHK_ERR( error );
 }
 
 std::sort( childs.begin(), childs.end() );
 childs.erase( std::unique( childs.begin(), childs.end() ), childs.end() );
 
 // Add child entities to tagged sets
 error = mbImpl->add_entities( *set_it, &childs[0], childs.size() );MB_CHK_ERR( error );
 }
 
 // Remove the coarse entities
 error = mbImpl->remove_entities( *set_it, set_ents );MB_CHK_ERR( error );
 
 // Add
 error = mbImpl->add_entities( lset, &( *set_it ), 1 );MB_CHK_ERR( error );
 }
 }
 return MB_SUCCESS;
}
 
/***********************************************
 * Basic functionalities: generate HM *
 ***********************************************/
 
ErrorCode NestedRefine::estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] )
{
 ErrorCode error;
 
 // Obtain the size of input mesh.
 int nverts_prev, nedges_prev, nfaces_prev, ncells_prev;
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nedges_prev = level_mesh[cur_level - 1].num_edges;
 nfaces_prev = level_mesh[cur_level - 1].num_faces;
 ncells_prev = level_mesh[cur_level - 1].num_cells;
 }
 else
 {
 nverts_prev = _inverts.size();
 nedges_prev = _inedges.size();
 nfaces_prev = _infaces.size();
 ncells_prev = _incells.size();
 }
 
 // Estimate mesh size of current level mesh.
 int nedges = 0, nfaces = 0;
 error = count_subentities( set, cur_level - 1, &nedges, &nfaces );MB_CHK_ERR( error );
 
 int d = get_index_from_degree( level_degree );
 int nverts = refTemplates[MBEDGE - 1][d].nv_edge * nedges;
 hmest[0] = nverts_prev + nverts;
 hmest[1] = nedges_prev * refTemplates[MBEDGE - 1][d].total_new_ents;
 hmest[2] = 0;
 hmest[3] = 0;
 
 int findex, cindex;
 if( nfaces_prev != 0 )
 {
 EntityHandle start_face;
 if( cur_level )
 start_face = level_mesh[cur_level - 1].start_face;
 else
 start_face = *_infaces.begin();
 findex = mbImpl->type_from_handle( start_face ) - 1;
 hmest[2] = nfaces_prev * refTemplates[findex][d].total_new_ents;
 
 if( meshdim == 2 ) hmest[0] += refTemplates[findex][d].nv_face * nfaces_prev;
 
 if( meshdim == 3 ) hmest[1] += nfaces_prev * intFacEdg[findex - 1][d].nie;
 }
 
 if( ncells_prev != 0 )
 {
 cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
 hmest[3] = ncells_prev * refTemplates[cindex][d].total_new_ents;
 
 hmest[0] += refTemplates[cindex][d].nv_face * nfaces;
 hmest[0] += refTemplates[cindex][d].nv_cell * ncells_prev;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] )
{
 // Obtain chunks of memory for the current level. Add them to a particular meshset.
 EntityHandle set_handle;
 ErrorCode error = mbImpl->create_meshset( MESHSET_SET, set_handle );MB_CHK_SET_ERR( error, "Cannot create mesh for the current level" );
 *set = set_handle;
 
 ReadUtilIface* read_iface;
 error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
 
 // Vertices
 error = read_iface->get_node_coords( 3, estL[0], 0, level_mesh[cur_level].start_vertex,
 level_mesh[cur_level].coordinates );MB_CHK_ERR( error );
 level_mesh[cur_level].num_verts = estL[0];
 
 Range newverts( level_mesh[cur_level].start_vertex, level_mesh[cur_level].start_vertex + estL[0] - 1 );
 error = mbImpl->add_entities( *set, newverts );MB_CHK_ERR( error );
 level_mesh[cur_level].verts = newverts;
 
 Tag gidtag;
 error = mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, gidtag );MB_CHK_ERR( error );
 error = read_iface->assign_ids( gidtag, newverts, level_mesh[cur_level].start_vertex );MB_CHK_ERR( error );
 
 // Edges
 if( estL[1] )
 {
 error = read_iface->get_element_connect( estL[1], 2, MBEDGE, 0, level_mesh[cur_level].start_edge,
 level_mesh[cur_level].edge_conn );MB_CHK_ERR( error );
 level_mesh[cur_level].num_edges = estL[1];
 
 Range newedges( level_mesh[cur_level].start_edge, level_mesh[cur_level].start_edge + estL[1] - 1 );
 error = mbImpl->add_entities( *set, newedges );MB_CHK_ERR( error );
 level_mesh[cur_level].edges = newedges;
 }
 else
 level_mesh[cur_level].num_edges = 0;
 
 // Faces
 if( estL[2] )
 {
 EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
 int nvpf = ahf->lConnMap2D[type - 2].num_verts_in_face;
 error = read_iface->get_element_connect( estL[2], nvpf, type, 0, level_mesh[cur_level].start_face,
 level_mesh[cur_level].face_conn );MB_CHK_ERR( error );
 level_mesh[cur_level].num_faces = estL[2];
 
 Range newfaces( level_mesh[cur_level].start_face, level_mesh[cur_level].start_face + estL[2] - 1 );
 error = mbImpl->add_entities( *set, newfaces );MB_CHK_ERR( error );
 level_mesh[cur_level].faces = newfaces;
 }
 else
 level_mesh[cur_level].num_faces = 0;
 
 // Cells
 if( estL[3] )
 {
 EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
 error = read_iface->get_element_connect( estL[3], nvpc, type, 0, level_mesh[cur_level].start_cell,
 level_mesh[cur_level].cell_conn );MB_CHK_ERR( error );
 level_mesh[cur_level].num_cells = estL[3];
 
 Range newcells( level_mesh[cur_level].start_cell, level_mesh[cur_level].start_cell + estL[3] - 1 );
 error = mbImpl->add_entities( *set, newcells );MB_CHK_ERR( error );
 level_mesh[cur_level].cells = newcells;
 }
 else
 level_mesh[cur_level].num_cells = 0;
 
 // Resize the ahf maps
 error = ahf->resize_hf_maps( level_mesh[cur_level].start_vertex, level_mesh[cur_level].num_verts,
 level_mesh[cur_level].start_edge, level_mesh[cur_level].num_edges,
 level_mesh[cur_level].start_face, level_mesh[cur_level].num_faces,
 level_mesh[cur_level].start_cell, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
 
 error = ahf->update_entity_ranges( *set );MB_CHK_ERR( error );
 
 // If the mesh type changes, then update the member variable in ahf to use the applicable
 // adjacency matrix
 MESHTYPE nwmesh = ahf->get_mesh_type( level_mesh[cur_level].num_verts, level_mesh[cur_level].num_edges,
 level_mesh[cur_level].num_faces, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
 if( ahf->thismeshtype != nwmesh ) ahf->thismeshtype = nwmesh;
 
 return MB_SUCCESS;
}
 
/**********************************
 * Hierarchical Mesh Generation *
 * *********************************/
 
ErrorCode NestedRefine::generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize )
{
 ErrorCode error;
 
 Tag gidtag;
 error = mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, gidtag );MB_CHK_ERR( error );
 
 nlevels = num_level;
 
 timeall.tm_total = 0;
 timeall.tm_refine = 0;
 timeall.tm_resolve = 0;
 
 for( int l = 0; l < num_level; l++ )
 {
 double tstart;
 tstart = tm->time_elapsed();
 
 // Estimate storage
 int hmest[4] = { 0, 0, 0, 0 };
 EntityHandle set;
 if( l )
 set = hm_set[l - 1];
 else
 set = _rset;
 error = estimate_hm_storage( set, level_degrees[l], l, hmest );MB_CHK_ERR( error );
 
 // Create arrays for storing the current level
 error = create_hm_storage_single_level( &hm_set[l], l, hmest );MB_CHK_ERR( error );
 
 // Copy the old vertices along with their coordinates
 error = copy_vertices_from_prev_level( l );MB_CHK_ERR( error );
 
 // Create the new entities and new vertices
 error = construct_hm_entities( l, level_degrees[l] );MB_CHK_ERR( error );
 
 timeall.tm_refine += tm->time_elapsed() - tstart;
 
 // Go into parallel communication
 if( !optimize )
 {
#ifdef MOAB_HAVE_MPI
 if( pcomm && ( pcomm->size() > 1 ) )
 {
 double tpstart = tm->time_elapsed();
 error = resolve_shared_ents_parmerge( l, hm_set[l] );MB_CHK_ERR( error );
 timeall.tm_resolve += tm->time_elapsed() - tpstart;
 }
#endif
 }
 }
 
 if( optimize )
 {
#ifdef MOAB_HAVE_MPI
 if( pcomm && ( pcomm->size() > 1 ) )
 {
 double tpstart = tm->time_elapsed();
 error = resolve_shared_ents_opt( hm_set, nlevels );MB_CHK_ERR( error );
 timeall.tm_resolve = tm->time_elapsed() - tpstart;
 }
#endif
 }
 timeall.tm_total = timeall.tm_refine + timeall.tm_resolve;
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::construct_hm_entities( int cur_level, int deg )
{
 ErrorCode error;
 
 // Generate mesh for current level by refining previous level.
 if( ahf->thismeshtype == CURVE )
 {
 error = construct_hm_1D( cur_level, deg );MB_CHK_ERR( error );
 }
 else if( ahf->thismeshtype == SURFACE || ahf->thismeshtype == SURFACE_MIXED )
 {
 error = construct_hm_2D( cur_level, deg );MB_CHK_ERR( error );
 }
 else
 {
 error = construct_hm_3D( cur_level, deg );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::construct_hm_1D( int cur_level, int deg )
{
 ErrorCode error;
 int nverts_prev, nents_prev;
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_edges;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _inedges.size();
 }
 
 int d = get_index_from_degree( deg );
 int vtotal = 2 + refTemplates[0][d].total_new_verts;
 std::vector< EntityHandle > vbuffer( vtotal );
 
 std::vector< EntityHandle > conn;
 int count_nents = 0;
 int count_verts = nverts_prev;
 
 // Step 1: Create the subentities via refinement of the previous mesh
 for( int eid = 0; eid < nents_prev; eid++ )
 {
 conn.clear();
 
 // EntityHandle of the working edge
 EntityHandle edge;
 if( cur_level )
 edge = level_mesh[cur_level - 1].start_edge + eid;
 else
 edge = _inedges[eid]; // Makes the assumption initial mesh is contiguous in memory
 
 error = get_connectivity( edge, cur_level, conn );MB_CHK_ERR( error );
 
 // Add the vertex handles to vbuffer for the current level for the working edge
 
 // Since the old vertices are copied first, their local indices do not change as new levels
 // are added. Clearly the local indices of the new vertices introduced in the current level
 // is still the same when the old vertices are copied. Thus, there is no need to explicitly
 // store another map between the old and duplicates in the subsequent levels. The second
 // part in the following sum is the local index in the previous level.
 
 // Add the corners to the vbuffer first.
 
 for( int i = 0; i < (int)conn.size(); i++ )
 {
 if( cur_level )
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
 else
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
 }
 
 // Adding rest of the entityhandles to working buffer for vertices.
 int num_new_verts = refTemplates[0][d].total_new_verts;
 
 for( int i = 0; i < num_new_verts; i++ )
 {
 vbuffer[i + 2] = level_mesh[cur_level].start_vertex + count_verts;
 count_verts += 1;
 }
 
 // Use the template to obtain the subentities
 int id1, id2;
 int etotal = refTemplates[0][d].total_new_ents;
 std::vector< EntityHandle > ent_buffer( etotal );
 
 for( int i = 0; i < etotal; i++ )
 {
 id1 = refTemplates[0][d].ents_conn[i][0];
 id2 = refTemplates[0][d].ents_conn[i][1];
 level_mesh[cur_level].edge_conn[2 * ( count_nents )] = vbuffer[id1];
 level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
 ent_buffer[i] = level_mesh[cur_level].start_edge + count_nents;
 count_nents += 1;
 };
 
 error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
 
 // Compute the coordinates of the new vertices: Linear interpolation
 int idx;
 double xi;
 for( int i = 0; i < num_new_verts; i++ )
 {
 xi = refTemplates[0][d].vert_nat_coord[i][0];
 idx = vbuffer[i + 2] - level_mesh[cur_level].start_vertex; // index of new vertex in current level
 if( cur_level )
 {
 id1 = conn[0] - level_mesh[cur_level - 1].start_vertex; // index of old end vertices in current level
 id2 = conn[1] - level_mesh[cur_level - 1].start_vertex;
 }
 else
 {
 id1 = _inverts.index( conn[0] );
 id2 = _inverts.index( conn[1] );
 }
 
 level_mesh[cur_level].coordinates[0][idx] =
 ( 1 - xi ) * level_mesh[cur_level].coordinates[0][id1] + xi * level_mesh[cur_level].coordinates[0][id2];
 level_mesh[cur_level].coordinates[1][idx] =
 ( 1 - xi ) * level_mesh[cur_level].coordinates[1][id1] + xi * level_mesh[cur_level].coordinates[1][id2];
 level_mesh[cur_level].coordinates[2][idx] =
 ( 1 - xi ) * level_mesh[cur_level].coordinates[2][id1] + xi * level_mesh[cur_level].coordinates[2][id2];
 }
 }
 
 error = update_global_ahf( MBEDGE, cur_level, deg );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::construct_hm_1D( int cur_level,
 int deg,
 EntityType type,
 std::vector< EntityHandle >& trackverts )
{
 ErrorCode error;
 
 int nedges_prev;
 if( cur_level )
 nedges_prev = level_mesh[cur_level - 1].num_edges;
 else
 nedges_prev = _inedges.size();
 
 int d = get_index_from_degree( deg );
 int nve = refTemplates[0][d].nv_edge;
 int vtotal = 2 + refTemplates[0][d].total_new_verts;
 int etotal = refTemplates[0][d].total_new_ents;
 int ne = 0, dim = 0, index = 0;
 if( type == MBTRI || type == MBQUAD )
 {
 index = type - 2;
 ne = ahf->lConnMap2D[index].num_verts_in_face;
 dim = 2;
 }
 else if( type == MBTET || type == MBHEX )
 {
 index = ahf->get_index_in_lmap( *( _incells.begin() ) );
 ne = ahf->lConnMap3D[index].num_edges_in_cell;
 dim = 3;
 }
 
 std::vector< EntityHandle > vbuffer( vtotal );
 std::vector< EntityHandle > ent_buffer( etotal );
 
 std::vector< EntityHandle > adjents, econn, fconn;
 std::vector< int > leids;
 int count_nents = 0;
 
 // Loop over all the edges and gather the vertices to be used for refinement
 for( int eid = 0; eid < nedges_prev; eid++ )
 {
 adjents.clear();
 leids.clear();
 econn.clear();
 fconn.clear();
 for( int i = 0; i < vtotal; i++ )
 vbuffer[i] = 0;
 for( int i = 0; i < etotal; i++ )
 ent_buffer[i] = 0;
 
 EntityHandle edge;
 if( cur_level )
 edge = level_mesh[cur_level - 1].start_edge + eid;
 else
 edge = _inedges[eid];
 
 error = get_connectivity( edge, cur_level, econn );MB_CHK_ERR( error );
 
 for( int i = 0; i < (int)econn.size(); i++ )
 {
 if( cur_level )
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - level_mesh[cur_level - 1].start_vertex );
 else
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - *_inverts.begin() );
 }
 
 int fid = -1, lid = -1, idx1 = -1, idx2 = -1;
 
 if( dim == 2 )
 {
 error = ahf->get_up_adjacencies_2d( edge, adjents, &leids );MB_CHK_ERR( error );
 if( cur_level )
 fid = adjents[0] - level_mesh[cur_level - 1].start_face;
 else
 fid = _infaces.index( adjents[0] );
 
 lid = leids[0];
 idx1 = lid;
 idx2 = ahf->lConnMap2D[index].next[lid];
 }
 else if( dim == 3 )
 {
 error = ahf->get_up_adjacencies_edg_3d( edge, adjents, &leids );MB_CHK_ERR( error );
 if( cur_level )
 fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
 else
 fid = _incells.index( adjents[0] );
 
 lid = leids[0];
 idx1 = ahf->lConnMap3D[index].e2v[lid][0];
 idx2 = ahf->lConnMap3D[index].e2v[lid][1];
 }
 
 error = get_connectivity( adjents[0], cur_level, fconn );MB_CHK_ERR( error );
 
 bool orient = false;
 if( ( fconn[idx1] == econn[0] ) && ( fconn[idx2] == econn[1] ) ) orient = true;
 
 if( orient )
 {
 for( int j = 0; j < nve; j++ )
 vbuffer[j + 2] = trackverts[fid * ne * nve + nve * lid + j];
 }
 else
 {
 for( int j = 0; j < nve; j++ )
 vbuffer[( nve - j - 1 ) + 2] = trackverts[fid * ne * nve + nve * lid + j];
 }
 
 // Use the template to obtain the subentities
 int id1, id2;
 
 for( int i = 0; i < etotal; i++ )
 {
 id1 = refTemplates[0][d].ents_conn[i][0];
 id2 = refTemplates[0][d].ents_conn[i][1];
 level_mesh[cur_level].edge_conn[2 * ( count_nents )] = vbuffer[id1];
 level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
 ent_buffer[i] = level_mesh[cur_level].start_edge + count_nents;
 
 count_nents += 1;
 };
 
 error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
 }
 
 error = update_global_ahf_1D_sub( cur_level, deg );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::construct_hm_2D( int cur_level, int deg )
{
 ErrorCode error;
 int nverts_prev, nents_prev;
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_faces;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _infaces.size();
 }
 
 // Create some book-keeping arrays over the old mesh to avoid introducing duplicate vertices and
 // calculating vertices more than once.
 EntityType ftype = mbImpl->type_from_handle( *_infaces.begin() );
 int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
 int findex = ftype - 1;
 
 int d = get_index_from_degree( deg );
 int tnv = refTemplates[findex][d].total_new_verts;
 int vtotal = nepf + tnv;
 std::vector< EntityHandle > vbuffer( vtotal );
 int etotal = refTemplates[findex][d].total_new_ents;
 std::vector< EntityHandle > ent_buffer( etotal );
 
 int nve = refTemplates[findex][d].nv_edge;
 std::vector< EntityHandle > trackvertsF( nents_prev * nepf * nve, 0 );
 int cur_nverts = level_mesh[cur_level].num_verts;
 std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
 
 int count_nverts = nverts_prev;
 int count_nents = 0;
 std::vector< EntityHandle > conn, cur_conn;
 
 // Step 1: Create the subentities via refinement of the previous mesh
 for( int fid = 0; fid < nents_prev; fid++ )
 {
 conn.clear();
 cur_conn.clear();
 for( int i = 0; i < vtotal; i++ )
 vbuffer[i] = 0;
 for( int i = 0; i < etotal; i++ )
 ent_buffer[i] = 0;
 
 // EntityHandle of the working face
 EntityHandle face;
 if( cur_level )
 face = level_mesh[cur_level - 1].start_face + fid;
 else
 face = _infaces[fid];
 
 error = get_connectivity( face, cur_level, conn );MB_CHK_ERR( error );
 
 // Step 1: Add vertices from the current level for the working face that will be used for
 // subdivision.
 // Add the corners to vbuffer
 for( int i = 0; i < (int)conn.size(); i++ )
 {
 if( cur_level )
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
 else
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
 
 cur_conn.push_back( vbuffer[i] );
 }
 
 // Gather vertices already added to tracking array due to refinement of the sibling faces
 
 for( int i = 0; i < nepf; i++ )
 {
 for( int j = 0; j < nve; j++ )
 {
 int id = refTemplates[findex][d].vert_on_edges[i][j];
 vbuffer[id] = trackvertsF[fid * nve * nepf + nve * i + j];
 }
 }
 
 // Add the remaining vertex handles to vbuffer for the current level for the working face
 for( int i = 0; i < tnv; i++ )
 {
 if( !vbuffer[i + nepf] )
 {
 vbuffer[i + nepf] = level_mesh[cur_level].start_vertex + count_nverts;
 count_nverts += 1;
 }
 }
 
 // Step 2: Create the subentities using the template and the vbuffer
 int idx;
 for( int i = 0; i < etotal; i++ )
 {
 for( int k = 0; k < nepf; k++ )
 {
 idx = refTemplates[findex][d].ents_conn[i][k];
 level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
 }
 ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
 count_nents += 1;
 }
 
 // Step 3: Update the local AHF maps
 error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
 
 // Step 4: Add the new vertices to the tracking array
 int id;
 
 for( int i = 0; i < nepf; i++ )
 {
 // Add the vertices to trackvertsF for fid
 for( int j = 0; j < nve; j++ )
 {
 id = refTemplates[findex][d].vert_on_edges[i][j];
 trackvertsF[fid * nepf * nve + nve * i + j] = vbuffer[id];
 }
 
 std::vector< EntityHandle > sibfids;
 std::vector< int > sibleids;
 std::vector< int > siborient;
 
 // Add the vertices to trackvertsF for siblings of fid, if any.
 error = ahf->get_up_adjacencies_2d( face, i, false, sibfids, &sibleids, &siborient );MB_CHK_ERR( error );
 
 if( !sibfids.size() ) continue;
 
 for( int s = 0; s < (int)sibfids.size(); s++ )
 {
 int sibid;
 if( cur_level )
 sibid = sibfids[s] - level_mesh[cur_level - 1].start_face;
 else
 sibid = sibfids[s] - *_infaces.begin();
 
 if( siborient[s] ) // Same half-edge direction as the current half-edge
 {
 for( int j = 0; j < nve; j++ )
 {
 id = refTemplates[findex][d].vert_on_edges[i][j];
 trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
 }
 }
 else
 {
 for( int j = 0; j < nve; j++ )
 {
 id = refTemplates[findex][d].vert_on_edges[i][nve - j - 1];
 trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
 }
 }
 }
 }
 
 // Step 5: Compute the coordinates of the new vertices, avoids computing more than once via
 // the flag_verts array.
 std::vector< double > corner_coords( nepf * 3 );
 error = get_coordinates( &cur_conn[0], nepf, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
 
 error = compute_coordinates( cur_level, deg, ftype, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
 nverts_prev );MB_CHK_ERR( error );
 }
 
 // Step 6: Update the global maps
 error = update_global_ahf( ftype, cur_level, deg );MB_CHK_ERR( error );
 
 // Step 7: If edges exists, refine them.
 if( !_inedges.empty() )
 {
 error = construct_hm_1D( cur_level, deg, ftype, trackvertsF );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::construct_hm_2D( int cur_level,
 int deg,
 EntityType type,
 std::vector< EntityHandle >& trackvertsE,
 std::vector< EntityHandle >& trackvertsF )
{
 ErrorCode error;
 
 EntityType ftype = MBTRI;
 if( type == MBHEX ) ftype = MBQUAD;
 
 int d = get_index_from_degree( deg );
 int findex = ftype - 1;
 int cidx = ahf->get_index_in_lmap( *( _incells.begin() ) );
 
 int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
 int nepc = ahf->lConnMap3D[cidx].num_edges_in_cell;
 int nfpc = ahf->lConnMap3D[cidx].num_faces_in_cell;
 
 int tnv = refTemplates[findex][d].total_new_verts;
 int nve = refTemplates[findex][d].nv_edge;
 int nvf = refTemplates[findex][d].nv_face;
 int vtotal = nepf + tnv;
 int etotal = refTemplates[findex][d].total_new_ents;
 
 std::vector< EntityHandle > vbuffer( vtotal );
 std::vector< EntityHandle > ent_buffer( etotal );
 
 std::vector< EntityHandle > adjents, fconn, cconn;
 std::vector< int > leids;
 int count_nents = 0;
 
 int nents_prev, ecount;
 if( cur_level )
 {
 nents_prev = level_mesh[cur_level - 1].num_faces;
 ecount = level_mesh[cur_level - 1].num_edges * refTemplates[MBEDGE - 1][d].total_new_ents;
 ;
 }
 else
 {
 nents_prev = _infaces.size();
 ecount = _inedges.size() * refTemplates[MBEDGE - 1][d].total_new_ents;
 ;
 }
 
 // Step 1: Create the subentities via refinement of the previous mesh
 for( int it = 0; it < nents_prev; it++ )
 {
 fconn.clear();
 cconn.clear();
 adjents.clear();
 leids.clear();
 for( int i = 0; i < vtotal; i++ )
 vbuffer[i] = 0;
 for( int i = 0; i < etotal; i++ )
 ent_buffer[i] = 0;
 
 // EntityHandle of the working face
 EntityHandle face;
 if( cur_level )
 face = level_mesh[cur_level - 1].start_face + it;
 else
 face = _infaces[it];
 
 error = get_connectivity( face, cur_level, fconn );MB_CHK_ERR( error );
 
 // Add the new handles for old connectivity in the buffer
 for( int i = 0; i < (int)fconn.size(); i++ )
 {
 if( cur_level )
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - level_mesh[cur_level - 1].start_vertex );
 else
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - *_inverts.begin() );
 }
 
 // Add handles for vertices on edges and faces from the already refined cell
 int fid, lid;
 error = ahf->get_up_adjacencies_face_3d( face, adjents, &leids );MB_CHK_ERR( error );
 
 if( cur_level )
 fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
 else
 fid = _incells.index( adjents[0] );
 
 lid = leids[0];
 
 error = get_connectivity( adjents[0], cur_level, cconn );MB_CHK_ERR( error );
 
 // Find the orientation w.r.t the half-face and then add vertices properly.
 std::vector< EntityHandle > fac_conn( nepf );
 std::vector< EntityHandle > lfac_conn( nepf );
 for( int j = 0; j < nepf; j++ )
 {
 fac_conn[j] = fconn[j];
 int id = ahf->lConnMap3D[cidx].hf2v[lid][j];
 lfac_conn[j] = cconn[id];
 }
 
 std::vector< int > le_idx, indices;
 
 error = reorder_indices( deg, &fac_conn[0], &lfac_conn[0], nepf, le_idx, indices );MB_CHK_ERR( error );
 
 // Add the existing vertices on edges of the already refined cell to the vbuffer
 for( int j = 0; j < nepf; j++ )
 {
 int id = le_idx[j]; // Corresponding local edge
 int idx = ahf->lConnMap3D[cidx].f2leid[lid][id]; // Local edge in the cell
 
 // Get the orientation of the local edge of the face wrt the corresponding local edge in
 // the cell
 bool eorient = false;
 int fnext = ahf->lConnMap2D[ftype - 2].next[j];
 int idx1 = ahf->lConnMap3D[cidx].e2v[idx][0];
 int idx2 = ahf->lConnMap3D[cidx].e2v[idx][1];
 if( ( fconn[j] == cconn[idx1] ) && ( fconn[fnext] == cconn[idx2] ) ) eorient = true;
 
 if( eorient )
 {
 for( int k = 0; k < nve; k++ )
 {
 int ind = refTemplates[findex][d].vert_on_edges[j][k];
 vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
 }
 }
 else
 {
 for( int k = 0; k < nve; k++ )
 {
 int ind = refTemplates[findex][d].vert_on_edges[j][nve - k - 1];
 vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
 }
 }
 }
 
 // Add the existing vertices on the face of the refine cell to vbuffer
 if( nvf )
 {
 for( int k = 0; k < nvf; k++ )
 {
 int ind = refTemplates[findex][d].vert_on_faces[0][k];
 vbuffer[ind] = trackvertsF[fid * nfpc * nvf + nvf * lid + indices[k] - 1];
 }
 }
 
 // Create the subentities using the template and the vbuffer
 for( int i = 0; i < etotal; i++ )
 {
 for( int k = 0; k < nepf; k++ )
 {
 int idx = refTemplates[findex][d].ents_conn[i][k];
 level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
 }
 ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
 count_nents += 1;
 }
 
 error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
 
 // Create the interior edges
 int id1, id2;
 
 int ne = intFacEdg[ftype - 2][d].nie;
 for( int i = 0; i < ne; i++ )
 {
 id1 = intFacEdg[ftype - 2][d].ieconn[i][0];
 id2 = intFacEdg[ftype - 2][d].ieconn[i][1];
 level_mesh[cur_level].edge_conn[2 * ( ecount )] = vbuffer[id1];
 level_mesh[cur_level].edge_conn[2 * ( ecount ) + 1] = vbuffer[id2];
 ecount += 1;
 }
 }
 
 // Step 6: Update the global maps
 error = update_global_ahf_2D_sub( cur_level, deg );MB_CHK_ERR( error );
 
 // Step 7: Update the hf-maps for the edges
 error = update_ahf_1D( cur_level );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::construct_hm_3D( int cur_level, int deg )
{
 ErrorCode error;
 EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
 if( type == MBTET )
 {
 error = subdivide_tets( cur_level, deg );MB_CHK_ERR( error );
 }
 else
 {
 error = subdivide_cells( type, cur_level, deg );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::subdivide_cells( EntityType type, int cur_level, int deg )
{
 ErrorCode error;
 int nverts_prev, nents_prev;
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_cells;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _incells.size();
 }
 
 // Create some book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
 int cindex = type - 1;
 int d = get_index_from_degree( deg );
 int ne = refTemplates[cindex][d].nv_edge;
 int nvf = refTemplates[cindex][d].nv_face;
 int nvtotal = refTemplates[cindex][d].total_new_verts;
 
 int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
 int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
 int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
 int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
 
 int vtotal = nvpc + nvtotal;
 std::vector< EntityHandle > vbuffer( vtotal );
 
 std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
 std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
 
 int cur_nverts = level_mesh[cur_level].num_verts;
 std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
 
 int count_nverts = nverts_prev;
 int count_ents = 0;
 std::vector< EntityHandle > conn, cur_conn;
 
 // Step 1: Create the subentities via refinement of the previous mesh
 for( int cid = 0; cid < nents_prev; cid++ )
 {
 conn.clear();
 cur_conn.clear();
 for( int i = 0; i < vtotal; i++ )
 vbuffer[i] = 0;
 
 // EntityHandle of the working cell
 EntityHandle cell;
 if( cur_level )
 cell = level_mesh[cur_level - 1].start_cell + cid;
 else
 cell = _incells[cid];
 
 error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
 
 // Step 1: Add vertices from the current level for the working face that will be used for
 // subdivision.
 // Add the corners to vbuffer
 for( int i = 0; i < (int)conn.size(); i++ )
 {
 if( cur_level )
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
 else
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
 
 cur_conn.push_back( vbuffer[i] );
 }
 
 // Gather vertices already added to tracking array due to refinement of the sibling cells
 for( int i = 0; i < nepc; i++ )
 {
 for( int j = 0; j < ne; j++ )
 {
 int idx = refTemplates[cindex][d].vert_on_edges[i][j];
 vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
 }
 }
 
 // Add remaining new vertex handles
 for( int i = 0; i < nfpc; i++ )
 {
 for( int j = 0; j < nvf; j++ )
 {
 int idx = refTemplates[cindex][d].vert_on_faces[i][j];
 vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
 }
 }
 
 // Add the remaining vertex handles to vbuffer for the current level for the working cell
 for( int i = 0; i < nvtotal; i++ )
 {
 if( !vbuffer[i + nvpc] )
 {
 vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
 count_nverts += 1;
 }
 }
 
 // Step 2: Use the template to obtain the subentities. The coordinates and local ahf maps
 // are also constructed. Connectivity of the children
 int etotal = refTemplates[type - 1][d].total_new_ents;
 std::vector< EntityHandle > ent_buffer( etotal );
 
 for( int i = 0; i < etotal; i++ )
 {
 for( int k = 0; k < nvpc; k++ )
 {
 int idx = refTemplates[type - 1][d].ents_conn[i][k];
 level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
 }
 ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
 count_ents += 1;
 }
 
 // Step 3: Update local ahf maps
 error = update_local_ahf( deg, type, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
 
 // Step 4: Update tracking information
 error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
 
 // Step 5: Coordinates of the new vertices
 std::vector< double > corner_coords( nvpc * 3 );
 error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
 
 error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
 nverts_prev );MB_CHK_ERR( error );
 }
 
 // error = ahf->print_tags(3);
 
 // Step 6: Update the global maps
 error = update_global_ahf( type, cur_level, deg );MB_CHK_ERR( error );
 
 // Step 7: If edges exists, refine them as well.
 if( level_mesh[cur_level].num_edges != 0 )
 {
 error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
 }
 
 // Step 8: If faces exists, refine them as well.
 if( !_infaces.empty() )
 {
 error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
 }
 
 // error = ahf->print_tags(3);
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::subdivide_tets( int cur_level, int deg )
{
 ErrorCode error;
 int nverts_prev, nents_prev;
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_cells;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _incells.size();
 }
 
 EntityType type = MBTET;
 int cindex = type - 1;
 int d = get_index_from_degree( deg );
 int ne = refTemplates[cindex][d].nv_edge;
 int nvf = refTemplates[cindex][d].nv_face;
 int nvtotal = refTemplates[cindex][d].total_new_verts;
 
 int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
 int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
 int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
 int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
 
 // Create vertex buffer
 int vtotal = nvpc + nvtotal;
 std::vector< EntityHandle > vbuffer( vtotal );
 
 // Create book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
 std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
 std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
 
 int cur_nverts = level_mesh[cur_level].num_verts;
 std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
 std::vector< int > cell_patterns( nents_prev, 0 );
 
 int count_nverts = nverts_prev;
 int count_ents = 0;
 std::vector< EntityHandle > conn, cur_conn;
 
 // Step 1: Create the subentities via refinement of the previous mesh
 for( int cid = 0; cid < nents_prev; cid++ )
 {
 conn.clear();
 cur_conn.clear();
 for( int i = 0; i < vtotal; i++ )
 vbuffer[i] = 0;
 
 // EntityHandle of the working cell
 EntityHandle cell;
 if( cur_level )
 cell = level_mesh[cur_level - 1].start_cell + cid;
 else
 cell = _incells[cid];
 
 error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
 
 // Step 1: Add vertices from the current level for the working face that will be used for
 // subdivision.
 // Add the corners to vbuffer
 for( int i = 0; i < (int)conn.size(); i++ )
 {
 if( cur_level )
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
 else
 vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
 
 cur_conn.push_back( vbuffer[i] );
 }
 
 // Gather vertices already added to tracking array due to refinement of the sibling cells
 for( int i = 0; i < nepc; i++ )
 {
 for( int j = 0; j < ne; j++ )
 {
 int idx = refTemplates[cindex][d].vert_on_edges[i][j];
 vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
 }
 }
 
 // Add remaining new vertex handles
 for( int i = 0; i < nfpc; i++ )
 {
 for( int j = 0; j < nvf; j++ )
 {
 int idx = refTemplates[cindex][d].vert_on_faces[i][j];
 vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
 }
 }
 
 // Add the remaining vertex handles to vbuffer for the current level for the working cell
 for( int i = 0; i < nvtotal; i++ )
 {
 if( !vbuffer[i + nvpc] )
 {
 vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
 count_nverts += 1;
 }
 }
 
 // Step 2: Coordinates of the new vertices
 std::vector< double > corner_coords( nvpc * 3 );
 error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
 
 error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
 nverts_prev );MB_CHK_ERR( error );
 
 // Step 3: Choose the tet refine pattern to be used for this tet
 int diag = find_shortest_diagonal_octahedron( cur_level, deg, &vbuffer[0] );
 int pat_id = diag + 2;
 cell_patterns[cid] = pat_id;
 
 // Step 4: Use the template to obtain the subentities. The coordinates and local ahf maps
 // are also constructed. Connectivity of the children
 int etotal = refTemplates[pat_id][d].total_new_ents;
 std::vector< EntityHandle > ent_buffer( etotal );
 
 for( int i = 0; i < etotal; i++ )
 {
 for( int k = 0; k < nvpc; k++ )
 {
 int idx = refTemplates[pat_id][d].ents_conn[i][k];
 level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
 }
 ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
 count_ents += 1;
 }
 
 // Step 5: Update local ahf maps
 error = update_local_ahf( deg, MBTET, pat_id, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
 
 // Step 6: Update tracking information
 error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
 }
 
 // Step 7: Update the global maps
 // error = update_global_ahf(cur_level, deg, cell_patterns); MB_CHK_ERR(error);
 error = update_global_ahf( type, cur_level, deg, &cell_patterns );MB_CHK_ERR( error );
 
 // Step 8: If edges exists, refine them as well.
 if( level_mesh[cur_level].num_edges != 0 )
 {
 error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
 }
 
 // Step 9: If faces exists, refine them as well.
 if( !_infaces.empty() )
 {
 error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::compute_coordinates( int cur_level,
 int deg,
 EntityType type,
 EntityHandle* vbuffer,
 int vtotal,
 double* corner_coords,
 std::vector< int >& vflag,
 int nverts_prev )
{
 EntityHandle vstart = level_mesh[cur_level].start_vertex;
 int d = get_index_from_degree( deg );
 
 if( type == MBTRI )
 {
 double xi, eta, N[3];
 int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
 
 for( int i = 3; i < vtotal; i++ )
 {
 if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
 
 xi = refTemplates[findex][d].vert_nat_coord[i - 3][0];
 eta = refTemplates[findex][d].vert_nat_coord[i - 3][1];
 N[0] = 1 - xi - eta;
 N[1] = xi;
 N[2] = eta;
 
 double x = 0, y = 0, z = 0;
 for( int j = 0; j < 3; j++ )
 {
 x += N[j] * corner_coords[3 * j];
 y += N[j] * corner_coords[3 * j + 1];
 z += N[j] * corner_coords[3 * j + 2];
 }
 
 level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
 level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
 level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
 vflag[vbuffer[i] - vstart - nverts_prev] = 1;
 }
 }
 else if( type == MBQUAD )
 {
 double xi, eta, N[4];
 int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
 
 for( int i = 4; i < vtotal; i++ )
 {
 if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
 
 xi = refTemplates[findex][d].vert_nat_coord[i - 4][0];
 eta = refTemplates[findex][d].vert_nat_coord[i - 4][1];
 N[0] = ( 1 - xi ) * ( 1 - eta ) / 4;
 N[1] = ( 1 + xi ) * ( 1 - eta ) / 4;
 N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
 
 double x = 0, y = 0, z = 0;
 for( int j = 0; j < 4; j++ )
 {
 x += N[j] * corner_coords[3 * j];
 y += N[j] * corner_coords[3 * j + 1];
 z += N[j] * corner_coords[3 * j + 2];
 }
 
 level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
 level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
 level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
 vflag[vbuffer[i] - vstart - nverts_prev] = 1;
 }
 }
 else if( type == MBTET )
 {
 double xi, eta, mu, N[4];
 int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
 
 for( int i = 4; i < vtotal; i++ )
 {
 if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
 
 xi = refTemplates[cindex][d].vert_nat_coord[i - 4][0];
 eta = refTemplates[cindex][d].vert_nat_coord[i - 4][1];
 mu = refTemplates[cindex][d].vert_nat_coord[i - 4][2];
 
 N[0] = 1 - xi - eta - mu;
 N[1] = xi;
 N[2] = eta, N[3] = mu;
 
 double x = 0, y = 0, z = 0;
 for( int j = 0; j < 4; j++ )
 {
 x += N[j] * corner_coords[3 * j];
 y += N[j] * corner_coords[3 * j + 1];
 z += N[j] * corner_coords[3 * j + 2];
 }
 
 level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
 level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
 level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
 vflag[vbuffer[i] - vstart - nverts_prev] = 1;
 }
 }
 else if( type == MBPRISM )
 {
 double xi, eta, mu, N[6];
 int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
 
 for( int i = 6; i < vtotal; i++ )
 {
 if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
 
 xi = refTemplates[cindex][d].vert_nat_coord[i - 6][0];
 eta = refTemplates[cindex][d].vert_nat_coord[i - 6][1];
 mu = refTemplates[cindex][d].vert_nat_coord[i - 6][2];
 
 N[0] = ( 1 - xi - eta ) * ( 1 - mu ), N[1] = xi * ( 1 - mu ), N[2] = eta * ( 1 - mu ),
 N[3] = ( 1 - xi - eta ) * ( 1 + mu ), N[4] = xi * ( 1 + mu ), N[5] = eta * ( 1 + mu );
 
 double x = 0, y = 0, z = 0;
 for( int j = 0; j < 6; j++ )
 {
 x += N[j] * corner_coords[3 * j];
 y += N[j] * corner_coords[3 * j + 1];
 z += N[j] * corner_coords[3 * j + 2];
 }
 
 level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
 level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
 level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
 vflag[vbuffer[i] - vstart - nverts_prev] = 1;
 }
 }
 else if( type == MBHEX )
 {
 double xi, eta, mu, N[8];
 double d1, d2, d3, s1, s2, s3;
 int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
 
 for( int i = 8; i < vtotal; i++ )
 {
 
 if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
 
 xi = refTemplates[cindex][d].vert_nat_coord[i - 8][0];
 eta = refTemplates[cindex][d].vert_nat_coord[i - 8][1];
 mu = refTemplates[cindex][d].vert_nat_coord[i - 8][2];
 
 d1 = 1 - xi;
 d2 = 1 - eta;
 d3 = 1 - mu;
 s1 = 1 + xi;
 s2 = 1 + eta;
 s3 = 1 + mu;
 N[0] = ( d1 * d2 * d3 ) / 8;
 N[1] = ( s1 * d2 * d3 ) / 8;
 N[2] = ( s1 * s2 * d3 ) / 8;
 N[3] = ( d1 * s2 * d3 ) / 8;
 N[4] = ( d1 * d2 * s3 ) / 8;
 N[5] = ( s1 * d2 * s3 ) / 8;
 N[6] = ( s1 * s2 * s3 ) / 8;
 N[7] = ( d1 * s2 * s3 ) / 8;
 
 double x = 0, y = 0, z = 0;
 for( int j = 0; j < 8; j++ )
 {
 x += N[j] * corner_coords[3 * j];
 y += N[j] * corner_coords[3 * j + 1];
 z += N[j] * corner_coords[3 * j + 2];
 }
 
 level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
 level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
 level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
 
 vflag[vbuffer[i] - vstart - nverts_prev] = 1;
 }
 }
 return MB_SUCCESS;
}
/**********************************
 * Parallel Communication *
 * ********************************/
 
#ifdef MOAB_HAVE_MPI
ErrorCode NestedRefine::resolve_shared_ents_parmerge( int level, EntityHandle levelset )
{
 // TEMP: Add the adjacencies for MOAB-native DS
 // NOTE (VSM): This is expensive since it creates a doubly
 // redundant copy of the adjacency data in both MOAB-native
 // and AHF. Need to fix this with AHF optimized branch.
 ErrorCode error;
 ReadUtilIface* read_iface;
 error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
 if( level_mesh[level].num_edges != 0 )
 {
 error = read_iface->update_adjacencies( level_mesh[level].start_edge, level_mesh[level].num_edges, 2,
 level_mesh[level].edge_conn );MB_CHK_ERR( error );
 }
 if( level_mesh[level].num_faces != 0 )
 {
 EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
 int nvpf = ahf->lConnMap2D[type - 2].num_verts_in_face;
 error = read_iface->update_adjacencies( level_mesh[level].start_face, level_mesh[level].num_faces, nvpf,
 level_mesh[level].face_conn );MB_CHK_ERR( error );
 }
 if( level_mesh[level].num_cells != 0 )
 {
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
 error = read_iface->update_adjacencies( level_mesh[level].start_cell, level_mesh[level].num_cells, nvpc,
 level_mesh[level].cell_conn );MB_CHK_ERR( error );
 }
 
 if( pcomm->size() > 1 )
 {
 
 // get all entities on the rootset
 moab::Range vtxs, edgs, facs, elms;
 error = mbImpl->get_entities_by_dimension( levelset, 0, vtxs, false );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( levelset, 1, edgs, false );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( levelset, 2, facs, false );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( levelset, 3, elms, false );MB_CHK_ERR( error );
 
 // set the parallel partition tag data
 moab::Tag part_tag;
 int partid = pcomm->rank(), dum_id = -1;
 error = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
 error = mbImpl->tag_set_data( part_tag, &levelset, 1, &partid );MB_CHK_ERR( error );
 //
 // Now that we have the local piece of the mesh refined consistently,
 // call parallel merge instead of resolved_shared to stitch together the meshes
 // and to handle parallel communication of remote proc/entity-handle pairs for
 // shared + non-owned interfaces entities.
 //
 // TODO: This needs to be replaced by the following scheme in the future as
 // an optimization step.
 // > Assign global IDs consistently for shared entities so that parallel
 // > resolve shared ents can happen out of the box. This is the fastest option.
 
 ParallelMergeMesh pm( pcomm, 1e-08 );
 error = pm.merge( levelset, true );MB_CHK_ERR( error );
 
 //
 // Parallel Communication complete - all entities resolved
 //
 }
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels )
{
 assert( pcomm->size() > 1 );
 
 ErrorCode error;
 
 // Step 1A: Pre-processing: setting the parallel partition tag data
 for( int i = 0; i < num_levels; i++ )
 {
 Tag part_tag;
 int partid = pcomm->rank(), dum_id = -1;
 error = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
 error = mbImpl->tag_set_data( part_tag, &hm_set[i], 1, &partid );MB_CHK_ERR( error );
 }
 
 // Step 1B: Pre-processing: gather all shared entities and list entities shared with each
 // sharing processor
 
 // All shared processors
 std::set< unsigned int > shprocs;
 error = pcomm->get_comm_procs( shprocs );MB_CHK_ERR( error );
 
 std::vector< int > sharedprocs;
 for( std::set< unsigned int >::iterator it = shprocs.begin(); it != shprocs.end(); it++ )
 sharedprocs.push_back( *it );
 int nprocs = sharedprocs.size();
 
 // Create buffer variables storing the entities to be sent and received
 std::vector< std::vector< int > > nsharedEntsperproc( nprocs );
 std::vector< std::vector< EntityHandle > > localBuffs( nprocs );
 std::vector< std::vector< EntityHandle > > remlocalBuffs( nprocs );
 std::vector< std::vector< EntityHandle > > remoteBuffs( nprocs );
 
 int i;
 Range sharedentities;
 
 for( i = 0; i < nprocs; i++ )
 {
 // List of shared entities at the coarsest level
 sharedentities.clear();
 error = pcomm->get_shared_entities( sharedprocs[i], sharedentities, -1, true );MB_CHK_ERR( error );
 
 // Get the list shared edges and vertices that are not part of the shared edges
 Range allEnts;
 error = collect_shared_entities_by_dimension( sharedentities, allEnts );MB_CHK_ERR( error );
 
 Range V0, E0, F0;
 V0 = allEnts.subset_by_dimension( 0 );
 E0 = allEnts.subset_by_dimension( 1 );
 F0 = allEnts.subset_by_dimension( 2 );
 
 // Step 2A: Prepare msg to be sent:
 //
 // FList = <F, FC, FE> where F, FC and FE are vectors containing the faces, their
 // connectivities and edges. F = <F_0, F_1, F_2, ..., F_L> where
 // F_0 is the list of shared faces at the coarsest level,
 // F_i are the list of children faces at subsequent levels.
 // FC vector contains the connectivities of F_i's and FE vector contains the bounding edges
 // of F_i's.
 //
 // EList = <E, EC> where E and EC are vectors containing the edges and their connectivities.
 // E = <E_0, E_1, ..., E_L> where
 // E_0 is the list of shared edges at the coarsest level and
 // E_i are the list of children edges at subsequent levels.
 // The EC vector contains the connectivities of E_i's.
 //
 // VList = <V> = <V_0, V_1, ...., V_L> where V_0 are shared vertices at the coarsest level
 // and V_i are the duplicates in the subsequent levels.
 
 std::vector< EntityHandle > locFList, remFList;
 std::vector< EntityHandle > locEList, remEList;
 std::vector< EntityHandle > locVList, remVList;
 
 // collect faces
 if( !F0.empty() )
 {
 error = collect_FList( sharedprocs[i], F0, locFList, remFList );MB_CHK_ERR( error );
 }
 
 // collect edges
 if( !E0.empty() )
 {
 error = collect_EList( sharedprocs[i], E0, locEList, remEList );MB_CHK_ERR( error );
 }
 
 // collect vertices
 if( !V0.empty() )
 {
 error = collect_VList( sharedprocs[i], V0, locVList, remVList );MB_CHK_ERR( error );
 }
 
 // Step 2B: Add data to NR local buffer to be sent
 std::vector< int > msgsz;
 msgsz.push_back( F0.size() );
 msgsz.push_back( locFList.size() );
 msgsz.push_back( E0.size() );
 msgsz.push_back( locEList.size() );
 msgsz.push_back( V0.size() );
 msgsz.push_back( locVList.size() );
 nsharedEntsperproc[i].insert( nsharedEntsperproc[i].end(), msgsz.begin(), msgsz.end() );
 
 if( !F0.empty() )
 {
 localBuffs[i].insert( localBuffs[i].end(), locFList.begin(), locFList.end() );
 remlocalBuffs[i].insert( remlocalBuffs[i].end(), remFList.begin(), remFList.end() );
 }
 if( !E0.empty() )
 {
 localBuffs[i].insert( localBuffs[i].end(), locEList.begin(), locEList.end() );
 remlocalBuffs[i].insert( remlocalBuffs[i].end(), remEList.begin(), remEList.end() );
 }
 if( !V0.empty() )
 {
 localBuffs[i].insert( localBuffs[i].end(), locVList.begin(), locVList.end() );
 remlocalBuffs[i].insert( remlocalBuffs[i].end(), remVList.begin(), remVList.end() );
 }
 }
 
 // Step 3: Send and receive the remote collection of child ents
 error = pcomm->send_recv_entities( sharedprocs, nsharedEntsperproc, remlocalBuffs, remoteBuffs );MB_CHK_ERR( error );
 
 // Step 5: Resolve shared child entities and update parallel tags
 std::multimap< EntityHandle, int > rprocs;
 std::multimap< EntityHandle, EntityHandle > rhandles;
 
 error = decipher_remote_handles( sharedprocs, nsharedEntsperproc, localBuffs, remoteBuffs, rprocs, rhandles );MB_CHK_ERR( error );
 
 // Step 6: Update pcomm tags
 error = update_parallel_tags( rprocs, rhandles );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts )
{
 ErrorCode error;
 
 Range F0, E0, V0, E0all, V0all;
 std::vector< EntityHandle > ents;
 
 F0 = sharedEnts.subset_by_dimension( 2 );
 E0all = sharedEnts.subset_by_dimension( 1 );
 V0all = sharedEnts.subset_by_dimension( 0 );
 
 if( !F0.empty() )
 {
 Range edges, verts;
 
 for( Range::iterator it = F0.begin(); it != F0.end(); it++ )
 {
 ents.clear();
 error = ahf->get_adjacencies( *it, 1, ents );MB_CHK_ERR( error );
 std::copy( ents.begin(), ents.end(), range_inserter( edges ) );
 ents.clear();
 error = mbImpl->get_connectivity( &( *it ), 0, ents );MB_CHK_ERR( error );
 std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
 }
 
 E0 = subtract( E0all, edges );
 if( !E0.empty() )
 {
 for( Range::iterator it = E0.begin(); it != E0.end(); it++ )
 {
 ents.clear();
 error = mbImpl->get_connectivity( &( *it ), 0, ents );MB_CHK_ERR( error );
 std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
 }
 }
 V0 = subtract( V0all, verts );
 }
 else if( !E0all.empty() )
 {
 Range verts;
 for( Range::iterator it = E0all.begin(); it != E0all.end(); it++ )
 {
 ents.clear();
 error = mbImpl->get_connectivity( &( *it ), 1, ents );MB_CHK_ERR( error );
 std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
 }
 E0 = E0all;
 V0 = subtract( V0all, verts );
 }
 else if( !V0all.empty() )
 {
 V0 = V0all;
 }
 else
 MB_SET_ERR( MB_FAILURE, "Trying to pack unsupported sub-entities for shared interface entities" );
 
 if( !F0.empty() ) std::copy( F0.begin(), F0.end(), range_inserter( allEnts ) );
 if( !E0.empty() ) std::copy( E0.begin(), E0.end(), range_inserter( allEnts ) );
 if( !V0.empty() ) std::copy( V0.begin(), V0.end(), range_inserter( allEnts ) );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::collect_FList( int to_proc,
 Range faces,
 std::vector< EntityHandle >& FList,
 std::vector< EntityHandle >& RList )
{
 ErrorCode error;
 
 FList.clear();
 std::vector< EntityHandle > F, FC, FE, lF, lFC, lFE, rF, rFC, rFE;
 std::vector< EntityHandle > childEnts, conn, fedges;
 
 for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
 {
 EntityHandle face = *it;
 conn.clear();
 fedges.clear();
 error = mbImpl->get_connectivity( &face, 1, conn );MB_CHK_ERR( error );
 error = ahf->get_face_edges( *it, fedges );MB_CHK_ERR( error );
 
 // add local handles
 lF.push_back( *it );
 lFC.insert( lFC.end(), conn.begin(), conn.end() );
 lFE.insert( lFE.end(), fedges.begin(), fedges.end() );
 
 // replace local handles with remote handles on to_proc
 EntityHandle rval = 0;
 error = pcomm->get_remote_handles( &face, &rval, 1, to_proc );MB_CHK_ERR( error );
 rF.push_back( rval );
 
 for( int i = 0; i < (int)conn.size(); i++ )
 {
 error = pcomm->get_remote_handles( &conn[i], &rval, 1, to_proc );MB_CHK_ERR( error );
 rFC.push_back( rval );
 
 error = pcomm->get_remote_handles( &fedges[i], &rval, 1, to_proc );MB_CHK_ERR( error );
 rFE.push_back( rval );
 }
 }
 
 for( int l = 0; l < nlevels; l++ )
 {
 for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
 {
 childEnts.clear();
 error = parent_to_child( *it, 0, l + 1, childEnts );MB_CHK_ERR( error );
 
 for( int i = 0; i < (int)childEnts.size(); i++ )
 {
 conn.clear();
 fedges.clear();
 error = mbImpl->get_connectivity( &childEnts[i], 1, conn );MB_CHK_ERR( error );
 error = ahf->get_face_edges( childEnts[i], fedges );MB_CHK_ERR( error );
 
 F.push_back( childEnts[i] );
 FC.insert( FC.end(), conn.begin(), conn.end() );
 FE.insert( FE.end(), fedges.begin(), fedges.end() );
 }
 }
 }
 
 FList.insert( FList.end(), lF.begin(), lF.end() );
 FList.insert( FList.end(), F.begin(), F.end() );
 FList.insert( FList.end(), lFC.begin(), lFC.end() );
 FList.insert( FList.end(), FC.begin(), FC.end() );
 FList.insert( FList.end(), lFE.begin(), lFE.end() );
 FList.insert( FList.end(), FE.begin(), FE.end() );
 
 RList.insert( RList.end(), rF.begin(), rF.end() );
 RList.insert( RList.end(), F.begin(), F.end() );
 RList.insert( RList.end(), rFC.begin(), rFC.end() );
 RList.insert( RList.end(), FC.begin(), FC.end() );
 RList.insert( RList.end(), rFE.begin(), rFE.end() );
 RList.insert( RList.end(), FE.begin(), FE.end() );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::collect_EList( int to_proc,
 Range edges,
 std::vector< EntityHandle >& EList,
 std::vector< EntityHandle >& RList )
{
 ErrorCode error;
 EList.clear();
 std::vector< EntityHandle > E, EC, lE, lEC, rE, rEC;
 std::vector< EntityHandle > childEnts, conn;
 
 // Add the edges and their connectivities at the coarsest level first.
 for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
 {
 EntityHandle edg = *it;
 conn.clear();
 error = mbImpl->get_connectivity( &edg, 1, conn );MB_CHK_ERR( error );
 
 // add local handles
 lE.push_back( edg );
 lEC.insert( lEC.end(), conn.begin(), conn.end() );
 
 // replace local handles with remote handle on to_proc
 EntityHandle rval = 0;
 error = pcomm->get_remote_handles( &edg, &rval, 1, to_proc );MB_CHK_ERR( error );
 rE.push_back( rval );
 error = pcomm->get_remote_handles( &conn[0], &rval, 1, to_proc );MB_CHK_ERR( error );
 rEC.push_back( rval );
 error = pcomm->get_remote_handles( &conn[1], &rval, 1, to_proc );MB_CHK_ERR( error );
 rEC.push_back( rval );
 }
 
 // Add the edges and their connectivities at subsequent levels.
 for( int l = 0; l < nlevels; l++ )
 {
 for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
 {
 childEnts.clear();
 error = parent_to_child( *it, 0, l + 1, childEnts );MB_CHK_ERR( error );
 
 for( int i = 0; i < (int)childEnts.size(); i++ )
 {
 conn.clear();
 error = mbImpl->get_connectivity( &childEnts[i], 1, conn );MB_CHK_ERR( error );
 E.push_back( childEnts[i] );
 EC.insert( EC.end(), conn.begin(), conn.end() );
 }
 }
 }
 
 EList.insert( EList.end(), lE.begin(), lE.end() );
 EList.insert( EList.end(), E.begin(), E.end() );
 EList.insert( EList.end(), lEC.begin(), lEC.end() );
 EList.insert( EList.end(), EC.begin(), EC.end() );
 
 RList.insert( RList.end(), rE.begin(), rE.end() );
 RList.insert( RList.end(), E.begin(), E.end() );
 RList.insert( RList.end(), rEC.begin(), rEC.end() );
 RList.insert( RList.end(), EC.begin(), EC.end() );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::collect_VList( int to_proc,
 Range verts,
 std::vector< EntityHandle >& VList,
 std::vector< EntityHandle >& RList )
{
 ErrorCode error;
 std::vector< EntityHandle > V, lV, rV;
 VList.clear();
 
 // Add the vertices at the coarsest level first.
 for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
 {
 EntityHandle v = *it;
 
 lV.push_back( v );
 EntityHandle rval = 0;
 error = pcomm->get_remote_handles( &v, &rval, 1, to_proc );MB_CHK_ERR( error );
 
 rV.push_back( rval );
 }
 
 // Add the vertices at the subsequent levels .
 for( int l = 0; l < nlevels; l++ )
 {
 for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
 {
 EntityHandle dupvert = 0;
 error = get_vertex_duplicates( *it, l + 1, dupvert );MB_CHK_ERR( error );
 V.push_back( dupvert );
 }
 }
 
 // local vertex handles at the coarsest level
 VList.insert( VList.end(), lV.begin(), lV.end() );
 VList.insert( VList.end(), V.begin(), V.end() );
 
 // remote vertex handles at the coarsest level
 RList.insert( RList.end(), rV.begin(), rV.end() );
 RList.insert( RList.end(), V.begin(), V.end() );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::decipher_remote_handles( std::vector< int >& sharedprocs,
 std::vector< std::vector< int > >& auxinfo,
 std::vector< std::vector< EntityHandle > >& localbuffers,
 std::vector< std::vector< EntityHandle > >& remotebuffers,
 std::multimap< EntityHandle, int >& remProcs,
 std::multimap< EntityHandle, EntityHandle >& remHandles )
{
 ErrorCode error;
 
 int i;
 int nprocs = sharedprocs.size();
 
 for( i = 0; i < nprocs; i++ )
 {
 std::vector< int > msgsz;
 for( int j = 0; j < (int)auxinfo[i].size(); j++ )
 {
 msgsz.push_back( auxinfo[i][j] );
 }
 
 if( msgsz[0] != 0 ) // Faces
 {
 // Get the local and remote face handles from the buffers
 std::vector< EntityHandle > LFList, RFList;
 LFList.insert( LFList.end(), localbuffers[i].begin(), localbuffers[i].begin() + msgsz[1] );
 RFList.insert( RFList.end(), remotebuffers[i].begin(), remotebuffers[i].begin() + msgsz[1] );
 
 error = decipher_remote_handles_face( sharedprocs[i], msgsz[0], LFList, RFList, remProcs, remHandles );MB_CHK_ERR( error );
 
 if( msgsz[2] != 0 ) // Edges
 {
 std::vector< EntityHandle > LEList, REList;
 LEList.insert( LEList.end(), localbuffers[i].begin() + msgsz[1] + 1,
 localbuffers[i].begin() + msgsz[1] + msgsz[3] );
 REList.insert( REList.end(), remotebuffers[i].begin() + msgsz[1] + 1,
 remotebuffers[i].begin() + msgsz[1] + msgsz[3] );
 
 error = decipher_remote_handles_edge( sharedprocs[i], msgsz[2], LEList, REList, remProcs, remHandles );MB_CHK_ERR( error );
 
 if( msgsz[4] != 0 ) // Vertices
 {
 // Get the local and remote face handles from the buffers
 std::vector< EntityHandle > LVList, RVList;
 LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[1] + msgsz[3] + 1,
 localbuffers[i].begin() + msgsz[1] + msgsz[3] + msgsz[5] );
 RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[1] + msgsz[3] + 1,
 remotebuffers[i].begin() + msgsz[1] + msgsz[3] + msgsz[5] );
 
 error = decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs,
 remHandles );MB_CHK_ERR( error );
 }
 }
 else if( msgsz[4] != 0 ) // Vertices
 {
 // Get the local and remote face handles from the buffers
 std::vector< EntityHandle > LVList, RVList;
 LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[1] + 1,
 localbuffers[i].begin() + msgsz[1] + msgsz[5] );
 RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[1] + 1,
 remotebuffers[i].begin() + msgsz[1] + msgsz[5] );
 
 error =
 decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
 }
 }
 
 else if( msgsz[2] != 0 ) // Edges
 {
 // Get the local and remote face handles from the buffers
 std::vector< EntityHandle > LEList, REList;
 LEList.insert( LEList.end(), localbuffers[i].begin(), localbuffers[i].begin() + msgsz[3] );
 REList.insert( REList.end(), remotebuffers[i].begin(), remotebuffers[i].begin() + msgsz[3] );
 
 error = decipher_remote_handles_edge( sharedprocs[i], msgsz[2], LEList, REList, remProcs, remHandles );MB_CHK_ERR( error );
 
 if( msgsz[4] != 0 ) // Vertices
 {
 // Get the local and remote face handles from the buffers
 std::vector< EntityHandle > LVList, RVList;
 LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[3] + 1,
 localbuffers[i].begin() + msgsz[3] + msgsz[5] );
 RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[3] + 1,
 remotebuffers[i].begin() + msgsz[3] + msgsz[5] );
 
 error =
 decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
 }
 }
 
 else if( msgsz[4] != 0 ) // Vertices
 {
 // Get the local and remote face handles from the buffers
 std::vector< EntityHandle > LVList, RVList;
 LVList.insert( LVList.end(), localbuffers[i].begin(), localbuffers[i].end() );
 RVList.insert( RVList.end(), remotebuffers[i].begin(), remotebuffers[i].end() );
 
 error = decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Trying to decipher entities other than verts, edges, faces" );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::decipher_remote_handles_face( int shared_proc,
 int numfaces,
 std::vector< EntityHandle >& localFaceList,
 std::vector< EntityHandle >& remFaceList,
 std::multimap< EntityHandle, int >& remProcs,
 std::multimap< EntityHandle, EntityHandle >& remHandles )
{
 ErrorCode error;
 
 for( int i = 0; i < numfaces; i++ )
 {
 // find local and remote handles of the coarsest face
 EntityHandle Lface = localFaceList[i];
 int Rface_idx =
 ( std::find( remFaceList.begin(), remFaceList.begin() + numfaces - 1, Lface ) ) - remFaceList.begin();
 
 // get connectivities of the local and remote coarsest faces
 std::vector< EntityHandle > Lface_conn, Rface_conn;
 error = get_data_from_buff( 2, 1, 0, i, numfaces, localFaceList, Lface_conn );MB_CHK_ERR( error );
 error = get_data_from_buff( 2, 1, 0, Rface_idx, numfaces, remFaceList, Rface_conn );MB_CHK_ERR( error );
 
 // find the combination difference between local and remote coarsest face
 std::vector< int > cmap;
 int comb = 0;
 int nvF = (int)Lface_conn.size();
 error = reorder_indices( &Lface_conn[0], &Rface_conn[0], nvF, &cmap[0], comb );MB_CHK_ERR( error );
 
 // go into loop over all levels
 std::vector< EntityHandle > lchildents, lparents;
 std::vector< EntityHandle > lcents, rcents;
 int lidx, ridx;
 
 for( int l = 0; l < nlevels; l++ )
 {
 lchildents.clear();
 lparents.clear();
 lcents.clear();
 rcents.clear();
 
 // obtain children at the current level
 error = get_data_from_buff( 2, 0, l + 1, i, numfaces, localFaceList, lchildents );MB_CHK_ERR( error );
 
 // obtain parents at the previous level
 if( l == 0 )
 {
 lparents.push_back( Lface );
 }
 else
 {
 error = get_data_from_buff( 2, 0, l, i, numfaces, localFaceList, lparents );MB_CHK_ERR( error );
 }
 
 //#children at the previous level and the current level
 EntityType ftype = mbImpl->type_from_handle( Lface );
 int d = get_index_from_degree( level_dsequence[l] );
 int nch = refTemplates[ftype - 1][d].total_new_ents;
 std::vector< int > fmap;
 error = reorder_indices( level_dsequence[l], nvF, comb, &fmap[0] );MB_CHK_ERR( error );
 
 // loop over all the lparents
 for( int j = 0; j < (int)lparents.size(); j++ )
 {
 // list local childrent at the current level
 lidx = std::find( localFaceList.begin(), localFaceList.end(), lparents[j] ) - localFaceList.begin();
 error = get_data_from_buff( 2, 0, l + 1, lidx, numfaces, localFaceList, lcents );MB_CHK_ERR( error );
 
 // find the corresponding remote of lparent and its children
 EntityHandle rparent = 0;
 error = check_for_parallelinfo( lparents[j], shared_proc, remHandles, remProcs, rparent );MB_CHK_ERR( error );
 ridx = std::find( remFaceList.begin(), remFaceList.end(), rparent ) - remFaceList.begin();
 error = get_data_from_buff( 2, 0, l + 1, ridx, numfaces, remFaceList, rcents );MB_CHK_ERR( error );
 
 // match up local face with remote handles according to cmap
 std::vector< EntityHandle > lconn, rconn, ledg, redg;
 for( int k = 0; k < nch; k++ )
 {
 lconn.clear();
 rconn.clear();
 ledg.clear();
 redg.clear();
 
 // matching the local children face handles with their remotes
 bool found = check_for_parallelinfo( lcents[k], shared_proc, remProcs );
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lcents[k], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lcents[k], rcents[fmap[k]] ) );
 }
 
 // find indices of the matched child face
 lidx = std::find( localFaceList.begin(), localFaceList.end(), lcents[k] ) - localFaceList.begin();
 ridx = std::find( remFaceList.begin(), remFaceList.end(), rcents[fmap[k]] ) - remFaceList.begin();
 
 // find bounding edges and connectivity of the matched child face
 error = get_data_from_buff( 2, 2, l + 1, lidx, numfaces, localFaceList, ledg );MB_CHK_ERR( error );
 error = get_data_from_buff( 2, 1, l + 1, lidx, numfaces, localFaceList, lconn );MB_CHK_ERR( error );
 error = get_data_from_buff( 2, 2, l + 1, ridx, numfaces, remFaceList, redg );MB_CHK_ERR( error );
 error = get_data_from_buff( 2, 1, l + 1, ridx, numfaces, remFaceList, rconn );MB_CHK_ERR( error );
 
 // now match the handles of the bounding edges and the vertices using
 // combination difference
 for( int m = 0; m < (int)ledg.size(); m++ )
 {
 found = check_for_parallelinfo( ledg[m], shared_proc, remProcs );
 
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( ledg[m], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( ledg[m], redg[cmap[m]] ) );
 }
 
 found = check_for_parallelinfo( lconn[m], shared_proc, remProcs );
 
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lconn[m], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[m], rconn[cmap[m]] ) );
 }
 }
 }
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::decipher_remote_handles_edge( int shared_proc,
 int numedges,
 std::vector< EntityHandle >& localEdgeList,
 std::vector< EntityHandle >& remEdgeList,
 std::multimap< EntityHandle, int >& remProcs,
 std::multimap< EntityHandle, EntityHandle >& remHandles )
{
 ErrorCode error;
 
 for( int i = 0; i < numedges; i++ )
 {
 EntityHandle Ledge = localEdgeList[i];
 int Redge_idx =
 ( std::find( remEdgeList.begin(), remEdgeList.begin() + numedges - 1, Ledge ) ) - remEdgeList.begin();
 
 std::vector< EntityHandle > Ledge_conn, Redge_conn;
 error = get_data_from_buff( 1, 1, 0, i, numedges, localEdgeList, Ledge_conn );MB_CHK_ERR( error );
 error = get_data_from_buff( 1, 1, 0, Redge_idx, numedges, remEdgeList, Redge_conn );MB_CHK_ERR( error );
 
 bool orient = true;
 if( ( Ledge_conn[0] == Redge_conn[1] ) && ( Ledge_conn[1] == Redge_conn[0] ) ) orient = false;
 
 if( orient ) assert( ( Ledge_conn[0] == Redge_conn[0] ) && ( Ledge_conn[1] == Redge_conn[1] ) );
 
 std::vector< EntityHandle > lchildEdgs, rchildEdgs, lconn, rconn;
 for( int l = 0; l < nlevels; l++ )
 {
 lchildEdgs.clear();
 rchildEdgs.clear();
 error = get_data_from_buff( 1, 0, l + 1, i, numedges, localEdgeList, lchildEdgs );MB_CHK_ERR( error );
 error = get_data_from_buff( 1, 0, l + 1, Redge_idx, numedges, remEdgeList, rchildEdgs );MB_CHK_ERR( error );
 
 int nchd = lchildEdgs.size();
 if( orient )
 {
 for( int j = 0; j < nchd; j++ )
 {
 // match entityhandles of child edges
 bool found = check_for_parallelinfo( lchildEdgs[j], shared_proc, remProcs );
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lchildEdgs[j], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lchildEdgs[j], rchildEdgs[j] ) );
 }
 
 // match entityhandles of child vertices
 lconn.clear();
 rconn.clear();
 
 int lidx =
 std::find( localEdgeList.begin(), localEdgeList.end(), lchildEdgs[j] ) - localEdgeList.begin();
 int ridx = std::find( remEdgeList.begin(), remEdgeList.end(), rchildEdgs[j] ) - remEdgeList.begin();
 
 error = get_data_from_buff( 1, 1, l + 1, lidx, numedges, localEdgeList, lconn );MB_CHK_ERR( error );
 error = get_data_from_buff( 1, 1, l + 1, ridx, numedges, remEdgeList, rconn );MB_CHK_ERR( error );
 
 found = check_for_parallelinfo( lconn[0], shared_proc, remProcs );
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lconn[0], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[0], rconn[0] ) );
 }
 found = check_for_parallelinfo( lconn[1], shared_proc, remProcs );
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lconn[1], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[1], rconn[0] ) );
 }
 }
 }
 
 else
 {
 for( int j = 0; j < nchd; j++ )
 {
 // match entityhandles of child edges
 bool found = check_for_parallelinfo( lchildEdgs[j], shared_proc, remProcs );
 
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lchildEdgs[j], shared_proc ) );
 remHandles.insert(
 std::pair< EntityHandle, EntityHandle >( lchildEdgs[j], rchildEdgs[nchd - j - 1] ) );
 }
 
 // match entityhandles of child vertices
 lconn.clear();
 rconn.clear();
 
 int lidx =
 std::find( localEdgeList.begin(), localEdgeList.end(), lchildEdgs[j] ) - localEdgeList.begin();
 int ridx = std::find( remEdgeList.begin(), remEdgeList.end(), rchildEdgs[nchd - j - 1] ) -
 remEdgeList.begin();
 
 error = get_data_from_buff( 1, 1, l + 1, lidx, numedges, localEdgeList, lconn );MB_CHK_ERR( error );
 error = get_data_from_buff( 1, 1, l + 1, ridx, numedges, remEdgeList, rconn );MB_CHK_ERR( error );
 found = check_for_parallelinfo( lconn[0], shared_proc, remProcs );
 
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lconn[0], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[0], rconn[1] ) );
 }
 
 found = check_for_parallelinfo( lconn[1], shared_proc, remProcs );
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lconn[1], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[1], rconn[1] ) );
 }
 }
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::decipher_remote_handles_vertex( int shared_proc,
 int numverts,
 std::vector< EntityHandle >& localVertexList,
 std::vector< EntityHandle >& remVertexList,
 std::multimap< EntityHandle, int >& remProcs,
 std::multimap< EntityHandle, EntityHandle >& remHandles )
{
 // LVList = <V> where V = <LV0, LV1, ..,LVL>
 // RVList = <V> where V = <LV0', RV1, ..,RVL>, LV0' is the local handles of coarsest vertices
 // but in the order in which they appear on the remote/shared proc
 
 ErrorCode error;
 
 for( int i = 0; i < numverts; i++ )
 {
 EntityHandle v = localVertexList[i];
 int Rvert_idx =
 ( std::find( remVertexList.begin(), remVertexList.begin() + numverts - 1, v ) ) - remVertexList.begin();
 
 std::vector< EntityHandle > lverts, rverts;
 for( int l = 0; l < nlevels; l++ )
 {
 error = get_data_from_buff( 0, 0, l + 1, i, numverts, localVertexList, lverts );MB_CHK_ERR( error );
 error = get_data_from_buff( 0, 0, l + 1, Rvert_idx, numverts, remVertexList, rverts );MB_CHK_ERR( error );
 
 bool found = check_for_parallelinfo( lverts[0], shared_proc, remProcs );
 if( !found )
 {
 remProcs.insert( std::pair< EntityHandle, int >( lverts[0], shared_proc ) );
 remHandles.insert( std::pair< EntityHandle, EntityHandle >( lverts[0], rverts[0] ) );
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
 std::multimap< EntityHandle, EntityHandle >& remHandles )
{
 ErrorCode error;
 
 std::vector< int > rprocs;
 std::vector< EntityHandle > rhandles;
 
 std::multimap< EntityHandle, int >::iterator it;
 std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_procs;
 std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
 std::multimap< EntityHandle, EntityHandle >::iterator >
 it_handles;
 
 it = remProcs.begin();
 while( it != remProcs.end() )
 {
 rprocs.clear();
 rhandles.clear();
 
 EntityHandle entity = it->first;
 it_procs = remProcs.equal_range( entity );
 it_handles = remHandles.equal_range( entity );
 
 for( std::multimap< EntityHandle, int >::iterator pit = it_procs.first; pit != it_procs.second; pit++ )
 rprocs.push_back( pit->second );
 for( std::multimap< EntityHandle, EntityHandle >::iterator pit = it_handles.first; pit != it_handles.second;
 pit++ )
 rhandles.push_back( pit->second );
 
 error = pcomm->update_remote_data( entity, rprocs, rhandles );MB_CHK_ERR( error );
 
 it = remProcs.upper_bound( it->first );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_data_from_buff( int listtype,
 int datatype,
 int level,
 int entity_index,
 int nentities,
 std::vector< EntityHandle >& buffer,
 std::vector< EntityHandle >& data )
{
 /**
 * listtype = 2, 1, 0 for FList (faces), EList (edge) and VList (vertices), respectively
 * datatype = 0(entityhandles of entities at level 0 and their children upto nlevels, in the
 *case of vertices its the duplicate vertices at subsequent levels), = 1(connectivity for the
 *above), = 2(subentities, only case is edges from FList) entity_index = integer index of the
 *entity in the buffer nentities = number of entities at the coarsest level i.e., |F0| or |E0|
 *or |V0|
 *
 * Two types of queries:
 * 1) Given an entity, find its children at a particular level.
 * 2) Given any entity at any level , find its connectivity (or bounding edges)
 *
 **/
 
 data.clear();
 
 if( listtype == 2 ) // FList
 {
 // FList = <F, FC, FE> where F = <F0, F1, ..,FL>, FC = <FC0, FC1, .., FCL> and FE = <FE0,
 // FE1, .., FEL>
 if( datatype == 0 ) // requesting child entities
 {
 EntityType ftype = mbImpl->type_from_handle( buffer[0] );
 
 int start, end, toadd, prev;
 start = end = entity_index;
 toadd = prev = nentities;
 for( int i = 0; i < level; i++ )
 {
 int d = get_index_from_degree( level_dsequence[i] );
 int nch = refTemplates[ftype - 1][d].total_new_ents;
 start = start * nch;
 end = end * nch + nch - 1;
 if( i < level - 1 )
 {
 prev = prev * nch;
 toadd += prev;
 }
 }
 
 start += toadd;
 end += toadd;
 
 int num_child = end - start + 1;
 data.reserve( num_child );
 
 for( int i = start; i <= end; i++ )
 {
 EntityHandle child = buffer[i];
 data.push_back( child );
 }
 }
 else if( datatype == 1 ) // requesting connectivity of an entity
 {
 EntityType ftype = mbImpl->type_from_handle( buffer[0] );
 int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
 
 int toadd = nentities, prev = nentities;
 for( int i = 0; i < nlevels; i++ )
 {
 int d = get_index_from_degree( level_dsequence[i] );
 int nch = refTemplates[ftype - 1][d].total_new_ents;
 prev = prev * nch;
 toadd += prev;
 }
 
 for( int i = 0; i < nepf; i++ )
 data.push_back( buffer[toadd + nepf * entity_index + i] );
 }
 else if( datatype == 2 ) // requesting bounding edges of an entity
 {
 EntityType ftype = mbImpl->type_from_handle( buffer[0] );
 int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
 
 int toadd = nentities, prev = nentities;
 for( int i = 0; i < nlevels; i++ )
 {
 int d = get_index_from_degree( level_dsequence[i] );
 int nch = refTemplates[ftype - 1][d].total_new_ents;
 prev = prev * nch;
 toadd += prev;
 }
 toadd += toadd * nepf;
 
 for( int i = 0; i < nepf; i++ )
 data.push_back( buffer[toadd + nepf * entity_index + i] );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for faces" );
 }
 else if( listtype == 1 ) // EList
 {
 // EList = <E, EC> where E = <E0, E1, ..,EL> and EC = <EC0, EC1, .., ECL>
 if( datatype == 0 ) // requesting child entities
 {
 int start, end, toadd, prev;
 start = end = entity_index;
 toadd = prev = nentities;
 for( int i = 0; i < level; i++ )
 {
 int d = get_index_from_degree( level_dsequence[i] );
 int nch = refTemplates[MBEDGE - 1][d].total_new_ents;
 start = start * nch;
 end = end * nch + nch - 1;
 
 if( i < level - 1 )
 {
 prev = prev * nch;
 toadd += prev;
 }
 }
 
 start += toadd;
 end += toadd;
 
 int num_child = end - start + 1;
 data.reserve( num_child );
 
 for( int i = start; i <= end; i++ )
 {
 EntityHandle child = buffer[i];
 data.push_back( child );
 }
 }
 else if( datatype == 1 ) // requesting connectivities of child entities
 {
 int toadd = nentities, prev = nentities;
 for( int i = 0; i < nlevels; i++ )
 {
 int d = get_index_from_degree( level_dsequence[i] );
 int nch = refTemplates[MBEDGE - 1][d].total_new_ents;
 prev = prev * nch;
 toadd += prev;
 }
 
 data.push_back( buffer[toadd + 2 * entity_index] );
 data.push_back( buffer[toadd + 2 * entity_index + 1] );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for edges" );
 }
 else if( listtype == 0 ) // VList
 {
 // VList = <V> where V = <V0, V1, ..,VL>
 if( datatype == 0 )
 {
 int idx = level * nentities + entity_index;
 data.push_back( buffer[idx] );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for vertices" );
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer" );
 
 return MB_SUCCESS;
}
 
bool NestedRefine::check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs )
{
 bool found = false;
 
 std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_hes;
 it_hes = remProcs.equal_range( entity );
 
 for( std::multimap< EntityHandle, int >::iterator it = it_hes.first; it != it_hes.second; ++it )
 {
 if( it->second == proc )
 {
 found = true;
 break;
 }
 }
 
 return found;
}
 
ErrorCode NestedRefine::check_for_parallelinfo( EntityHandle entity,
 int proc,
 std::multimap< EntityHandle, EntityHandle >& remHandles,
 std::multimap< EntityHandle, int >& remProcs,
 EntityHandle& rhandle )
{
 // shared procs for given entity
 std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_ps;
 it_ps = remProcs.equal_range( entity );
 
 // shared remote handles for given entity
 std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
 std::multimap< EntityHandle, EntityHandle >::iterator >
 it_hs;
 it_hs = remHandles.equal_range( entity );
 
 std::multimap< EntityHandle, int >::iterator itp;
 std::multimap< EntityHandle, EntityHandle >::iterator ith;
 
 for( itp = it_ps.first, ith = it_hs.first; itp != it_ps.second; ++itp, ++ith )
 {
 if( itp->second == proc )
 {
 rhandle = ith->second;
 break;
 }
 }
 
 return MB_SUCCESS;
}
 
#endif
 
/**********************************
 * Update AHF maps *
 * ********************************/
 
ErrorCode NestedRefine::update_local_ahf( int deg,
 EntityType type,
 int pat_id,
 EntityHandle* vbuffer,
 EntityHandle* ent_buffer,
 int etotal )
{
 ErrorCode error;
 int nhf = 0, nv = 0, total_new_verts = 0;
 int d = get_index_from_degree( deg );
 
 // Get the number of half-facets
 if( type == MBEDGE )
 {
 nhf = 2;
 nv = 2;
 total_new_verts = refTemplates[0][d].total_new_verts;
 }
 else if( type == MBTRI || type == MBQUAD )
 {
 nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
 nv = nhf;
 total_new_verts = refTemplates[pat_id][d].total_new_verts;
 }
 else if( type == MBTET || type == MBHEX )
 {
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 nhf = ahf->lConnMap3D[index].num_faces_in_cell;
 nv = ahf->lConnMap3D[index].num_verts_in_cell;
 total_new_verts = refTemplates[pat_id][d].total_new_verts;
 }
 
 std::vector< EntityHandle > ent;
 std::vector< int > lid;
 
 // Update the vertex to half-facet map
 for( int i = 0; i < total_new_verts; i++ )
 {
 ent.clear();
 lid.clear();
 EntityHandle vid = vbuffer[i + nv];
 error = ahf->get_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
 
 if( ent[0] ) continue;
 
 int id = refTemplates[pat_id][d].v2hf[i + nv][0] - 1;
 ent[0] = ent_buffer[id];
 lid[0] = refTemplates[pat_id][d].v2hf[i + nv][1];
 
 error = ahf->set_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
 }
 
 // Update the sibling half-facet map
 for( int i = 0; i < etotal; i++ )
 {
 std::vector< EntityHandle > sib_entids( nhf );
 std::vector< int > sib_lids( nhf );
 
 error = ahf->get_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 for( int l = 0; l < nhf; l++ )
 {
 if( sib_entids[l] ) continue;
 
 // Fill out the sibling values
 int id = refTemplates[pat_id][d].ents_opphfs[i][2 * l];
 if( id )
 {
 sib_entids[l] = ent_buffer[id - 1];
 sib_lids[l] = refTemplates[pat_id][d].ents_opphfs[i][2 * l + 1];
 }
 else
 {
 sib_entids[l] = 0;
 sib_lids[l] = 0;
 }
 }
 
 error = ahf->set_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 for( int l = 0; l < nhf; l++ )
 {
 if( sib_entids[l] )
 {
 
 EntityHandle set_entid = ent_buffer[i];
 int set_lid = l;
 
 error = ahf->set_sibling_map( type, sib_entids[l], sib_lids[l], set_entid, set_lid );MB_CHK_ERR( error );
 }
 }
 }
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_local_ahf( int deg,
 EntityType type,
 EntityHandle* vbuffer,
 EntityHandle* ent_buffer,
 int etotal )
{
 ErrorCode error;
 assert( type != MBTET );
 error = update_local_ahf( deg, type, type - 1, vbuffer, ent_buffer, etotal );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids )
{
 
 ErrorCode error;
 
 // Get the number of half-facets and number of children of each type
 if( type == MBEDGE )
 {
 assert( pattern_ids == NULL );
 error = update_global_ahf_1D( cur_level, deg );MB_CHK_ERR( error );
 }
 else if( type == MBTRI || type == MBQUAD )
 {
 assert( pattern_ids == NULL );
 error = update_global_ahf_2D( cur_level, deg );MB_CHK_ERR( error );
 }
 else if( type == MBHEX )
 {
 assert( pattern_ids == NULL );
 error = update_global_ahf_3D( cur_level, deg );MB_CHK_ERR( error );
 }
 else if( type == MBTET )
 {
 assert( pattern_ids != NULL );
 error = update_global_ahf_3D( cur_level, deg, pattern_ids );MB_CHK_ERR( error );
 }
 else
 MB_SET_ERR( MB_NOT_IMPLEMENTED, "Requesting AHF update for an unsupported mesh entity type" );
 
 return MB_SUCCESS;
}
 
/*ErrorCode NestedRefine::update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids)
{
 ErrorCode error;
 error = update_global_ahf_3D(cur_level, deg, pattern_ids); MB_CHK_ERR(error);
 
 return MB_SUCCESS;
}*/
 
ErrorCode NestedRefine::update_global_ahf_1D( int cur_level, int deg )
{
 ErrorCode error;
 int d = get_index_from_degree( deg );
 int nhf, nchilds, nverts_prev, nents_prev;
 nhf = 2;
 nchilds = refTemplates[0][d].total_new_ents;
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_edges;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _inedges.size();
 }
 
 std::vector< EntityHandle > inci_ent, child_ents;
 std::vector< int > inci_lid, child_lids;
 
 // Update the vertex to half-facet maps for duplicate vertices
 for( int i = 0; i < nverts_prev; i++ )
 {
 inci_ent.clear();
 inci_lid.clear();
 child_ents.clear();
 child_lids.clear();
 
 // Vertex id in the previous mesh and the current one
 EntityHandle vid;
 if( cur_level )
 vid = level_mesh[cur_level - 1].start_vertex + i;
 else
 vid = _inverts[i];
 EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
 
 // Get the incident half-vert in the previous mesh
 error = ahf->get_incident_map( MBEDGE, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
 
 // Obtain the corresponding incident child in the current mesh
 int lvid = get_local_vid( vid, inci_ent[0], cur_level - 1 );
 if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
 int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
 
 int pid;
 if( cur_level )
 pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
 else
 pid = inci_ent[0] - *_inedges.begin();
 
 int ind = nchilds * pid;
 
 child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
 child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
 
 error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
 }
 
 // Update the sibling half-facet maps across entities
 for( int i = 0; i < nents_prev; i++ )
 {
 EntityHandle ent;
 if( cur_level )
 ent = level_mesh[cur_level - 1].start_edge + i;
 else
 ent = _inedges[i];
 
 std::vector< EntityHandle > sib_entids( nhf );
 std::vector< int > sib_lids( nhf );
 
 error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 int id, idx;
 
 for( int l = 0; l < nhf; l++ )
 {
 if( !sib_entids[l] ) continue;
 
 // Find the child incident on the half-facet
 id = refTemplates[0][d].ents_on_pent[l][1] - 1;
 idx = nchilds * i;
 EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
 int ch_lid = l;
 
 // Find the sibling of the child
 std::vector< EntityHandle > sib_childs( nhf );
 std::vector< int > sib_chlids( nhf );
 
 error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
 
 // If the sibling already exists, dont do anything
 if( sib_childs[ch_lid] ) continue;
 
 // Get the correponding child of the sibling of the current parent
 int psib;
 if( cur_level )
 psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
 else
 psib = sib_entids[l] - *_inedges.begin();
 
 int plid = sib_lids[l];
 
 id = refTemplates[0][d].ents_on_pent[plid][1] - 1;
 idx = nchilds * psib;
 
 EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
 int psib_chlid = plid;
 
 // Set the siblings
 sib_childs[ch_lid] = psib_child;
 sib_chlids[ch_lid] = psib_chlid;
 
 error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_global_ahf_1D_sub( int cur_level, int deg )
{
 ErrorCode error;
 int d = get_index_from_degree( deg );
 int nhf, nchilds, nents_prev;
 nhf = 2;
 nchilds = refTemplates[0][d].total_new_ents;
 if( cur_level )
 {
 nents_prev = level_mesh[cur_level - 1].num_edges;
 }
 else
 {
 nents_prev = _inedges.size();
 }
 
 // Update the sibling half-facet maps across entities
 
 std::vector< EntityHandle > conn;
 for( int i = 0; i < nents_prev; i++ )
 {
 EntityHandle ent;
 if( cur_level )
 ent = level_mesh[cur_level - 1].start_edge + i;
 else
 ent = _inedges[i];
 
 // Set incident hv maps
 conn.clear();
 error = get_connectivity( ent, cur_level, conn );MB_CHK_ERR( error );
 
 std::vector< EntityHandle > inci_ent, child_ents;
 std::vector< int > inci_lid, child_lids;
 for( int j = 0; j < 2; j++ )
 {
 inci_ent.clear();
 inci_lid.clear();
 child_ents.clear();
 child_lids.clear();
 
 // Get the entityhandle of the vertex from previous level in the current level
 EntityHandle cur_vid;
 if( cur_level )
 cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - level_mesh[cur_level - 1].start_vertex );
 else
 cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - *_inverts.begin() );
 
 // Obtain the incident half-facet. If exists, then no need to assign another
 error = ahf->get_incident_map( MBEDGE, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
 if( inci_ent[0] != 0 ) continue;
 
 // Get the incident half-facet on the old vertex
 error = ahf->get_incident_map( MBEDGE, conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
 
 // Obtain the corresponding incident child in the current mesh
 int lvid = get_local_vid( conn[j], inci_ent[0], cur_level - 1 );
 if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
 int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
 
 int pid;
 if( cur_level )
 pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
 else
 pid = inci_ent[0] - *_inedges.begin();
 
 int ind = nchilds * pid;
 
 child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
 child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
 
 error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
 }
 
 std::vector< EntityHandle > sib_entids( nhf );
 std::vector< int > sib_lids( nhf );
 
 error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 int id, idx;
 
 for( int l = 0; l < nhf; l++ )
 {
 if( !sib_entids[l] ) continue;
 
 // Find the child incident on the half-facet
 id = refTemplates[0][d].ents_on_pent[l][1] - 1;
 idx = nchilds * i;
 EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
 int ch_lid = l;
 
 // Find the sibling of the child
 std::vector< EntityHandle > sib_childs( nhf );
 std::vector< int > sib_chlids( nhf );
 
 error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
 
 // If the sibling already exists, dont do anything
 if( sib_childs[ch_lid] ) continue;
 
 // Get the correponding child of the sibling of the current parent
 int psib;
 if( cur_level )
 psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
 else
 psib = sib_entids[l] - *_inedges.begin();
 
 int plid = sib_lids[l];
 
 id = refTemplates[0][d].ents_on_pent[plid][1] - 1;
 idx = nchilds * psib;
 
 EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
 int psib_chlid = plid;
 
 // Set the siblings
 sib_childs[ch_lid] = psib_child;
 sib_chlids[ch_lid] = psib_chlid;
 
 error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_ahf_1D( int cur_level )
{
 ErrorCode error;
 error = ahf->determine_sibling_halfverts( level_mesh[cur_level].verts, level_mesh[cur_level].edges );MB_CHK_ERR( error );
 
 error = ahf->determine_incident_halfverts( level_mesh[cur_level].edges );MB_CHK_ERR( error );
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_global_ahf_2D( int cur_level, int deg )
{
 ErrorCode error;
 
 EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
 int nhf, nchilds, nverts_prev, nents_prev;
 
 nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
 int d = get_index_from_degree( deg );
 nchilds = refTemplates[type - 1][d].total_new_ents;
 
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_faces;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _infaces.size();
 }
 
 std::vector< EntityHandle > inci_ent, child_ents;
 std::vector< int > inci_lid, child_lids;
 
 // Update the vertex to half-edge maps for old/duplicate vertices
 for( int i = 0; i < nverts_prev; i++ )
 {
 inci_ent.clear();
 inci_lid.clear();
 child_ents.clear();
 child_lids.clear();
 
 // Vertex id in the previous mesh
 EntityHandle vid;
 if( cur_level )
 vid = level_mesh[cur_level - 1].start_vertex + i;
 else
 vid = _inverts[i];
 EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
 
 // Get the incident half-vert in the previous mesh
 error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
 
 // Obtain the corresponding incident child in the current mesh
 for( int j = 0; j < (int)inci_ent.size(); j++ )
 {
 int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
 if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
 int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
 
 int pid;
 if( cur_level )
 pid = inci_ent[j] - level_mesh[cur_level - 1].start_face;
 else
 pid = inci_ent[j] - *_infaces.begin();
 
 int ind = nchilds * pid;
 
 child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
 child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
 }
 error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
 }
 
 EntityHandle fedge[2];
 
 // Update the sibling half-facet maps across entities
 for( int i = 0; i < nents_prev; i++ )
 {
 EntityHandle ent;
 if( cur_level )
 ent = level_mesh[cur_level - 1].start_face + i;
 else
 ent = _infaces[i];
 
 std::vector< EntityHandle > fid_conn;
 error = get_connectivity( ent, cur_level, fid_conn );
 if( MB_SUCCESS != error ) return error;
 
 std::vector< EntityHandle > sib_entids( nhf );
 std::vector< int > sib_lids( nhf );
 
 error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 int id, idx;
 
 for( int l = 0; l < nhf; l++ )
 {
 if( !sib_entids[l] ) continue;
 
 int nidx = ahf->lConnMap2D[type - 2].next[l];
 fedge[0] = fid_conn[l];
 fedge[1] = fid_conn[nidx];
 
 EntityHandle sfid = sib_entids[l];
 int slid = sib_lids[l];
 
 std::vector< EntityHandle > conn;
 error = get_connectivity( sfid, cur_level, conn );
 if( MB_SUCCESS != error ) return error;
 
 bool orient = true;
 nidx = ahf->lConnMap2D[type - 2].next[slid];
 if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
 
 if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
 
 // Find the childrens incident on the half-facet
 int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
 idx = nchilds * i;
 
 // Loop over all the incident childrens
 for( int k = 0; k < nch; k++ )
 {
 id = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
 EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
 int child_lid = l;
 
 // Find the sibling of the child
 EntityHandle child_sibent;
 int child_siblid;
 error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
 
 if( child_sibent != 0 ) continue;
 
 // Get the correponding child of the sibling of the current parent
 int psib;
 if( cur_level )
 psib = sfid - level_mesh[cur_level - 1].start_face;
 else
 psib = sfid - *_infaces.begin();
 
 int plid = slid;
 
 if( orient )
 id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
 else
 id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
 
 int sidx = nchilds * psib;
 
 EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
 int psib_chlid = plid;
 
 // Set the siblings
 error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_global_ahf_2D_sub( int cur_level, int deg )
{
 ErrorCode error;
 int d = get_index_from_degree( deg );
 EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
 int nhf, nchilds, nents_prev;
 nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
 nchilds = refTemplates[type - 1][d].total_new_ents;
 
 if( cur_level )
 nents_prev = level_mesh[cur_level - 1].num_faces;
 else
 nents_prev = _infaces.size();
 
 EntityHandle fedge[2];
 
 // Update the sibling half-facet maps across entities
 for( int i = 0; i < nents_prev; i++ )
 {
 EntityHandle ent;
 if( cur_level )
 ent = level_mesh[cur_level - 1].start_face + i;
 else
 ent = _infaces[i];
 
 std::vector< EntityHandle > fid_conn;
 error = get_connectivity( ent, cur_level, fid_conn );
 if( MB_SUCCESS != error ) return error;
 
 std::vector< EntityHandle > inci_ent, child_ents;
 std::vector< int > inci_lid, child_lids;
 
 // Set incident half-edges
 for( int j = 0; j < nhf; j++ )
 {
 inci_ent.clear();
 inci_lid.clear();
 child_ents.clear();
 child_lids.clear();
 EntityHandle cur_vid;
 if( cur_level )
 cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - level_mesh[cur_level - 1].start_vertex );
 else
 cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - *_inverts.begin() );
 
 // Obtain the incident half-facet. If exists, then no need to assign another
 error = ahf->get_incident_map( type, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
 if( inci_ent[0] != 0 ) continue;
 
 // Get the incident half-facet on the old vertex
 error = ahf->get_incident_map( type, fid_conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
 
 // Obtain the corresponding incident child in the current mesh
 for( int k = 0; k < (int)inci_ent.size(); k++ )
 {
 int lvid = get_local_vid( fid_conn[j], inci_ent[k], cur_level - 1 );
 if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
 int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
 
 int pid;
 if( cur_level )
 pid = inci_ent[k] - level_mesh[cur_level - 1].start_face;
 else
 pid = inci_ent[k] - *_infaces.begin();
 
 int ind = nchilds * pid;
 
 child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
 child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
 }
 
 error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
 }
 
 // Set sibling half-edges
 std::vector< EntityHandle > sib_entids( nhf );
 std::vector< int > sib_lids( nhf );
 
 error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 int id, idx;
 
 for( int l = 0; l < nhf; l++ )
 {
 if( !sib_entids[l] ) continue;
 
 int nidx = ahf->lConnMap2D[type - 2].next[l];
 fedge[0] = fid_conn[l];
 fedge[1] = fid_conn[nidx];
 
 EntityHandle sfid = sib_entids[l];
 int slid = sib_lids[l];
 
 std::vector< EntityHandle > conn;
 error = get_connectivity( sfid, cur_level, conn );MB_CHK_ERR( error );
 
 assert( (int)conn.size() > nidx && (int)conn.size() > slid );
 
 bool orient = true;
 nidx = ahf->lConnMap2D[type - 2].next[slid];
 if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
 
 if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
 
 // Find the childrens incident on the half-facet
 int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
 idx = nchilds * i;
 
 // Loop over all the incident childrens
 for( int k = 0; k < nch; k++ )
 {
 id = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
 EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
 int child_lid = l;
 
 // Find the sibling of the child
 EntityHandle child_sibent;
 int child_siblid;
 error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
 
 if( child_sibent != 0 ) continue;
 
 // Get the correponding child of the sibling of the current parent
 int psib;
 if( cur_level )
 psib = sfid - level_mesh[cur_level - 1].start_face;
 else
 psib = sfid - *_infaces.begin();
 
 int plid = slid;
 
 if( orient )
 id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
 else
 id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
 
 int sidx = nchilds * psib;
 
 EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
 int psib_chlid = plid;
 
 // Set the siblings
 error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids )
{
 ErrorCode error;
 int nvpc, ne, nhf, nchilds, nverts_prev, nents_prev;
 
 EntityType type = mbImpl->type_from_handle( *_incells.begin() );
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 int d = get_index_from_degree( deg );
 
 nhf = ahf->lConnMap3D[index].num_faces_in_cell;
 ne = ahf->lConnMap3D[index].num_edges_in_cell;
 nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
 nchilds = refTemplates[type - 1][d].total_new_ents;
 
 if( cur_level )
 {
 nverts_prev = level_mesh[cur_level - 1].num_verts;
 nents_prev = level_mesh[cur_level - 1].num_cells;
 }
 else
 {
 nverts_prev = _inverts.size();
 nents_prev = _incells.size();
 }
 
 std::vector< EntityHandle > inci_ent, child_ents;
 std::vector< int > inci_lid, child_lids;
 
 // Step 1: Update the V2HF maps for old/duplicate vertices
 for( int i = 0; i < nverts_prev; i++ )
 {
 inci_ent.clear();
 inci_lid.clear();
 child_ents.clear();
 child_lids.clear();
 
 // Vertex id in the previous mesh
 EntityHandle vid;
 if( cur_level )
 vid = level_mesh[cur_level - 1].start_vertex + i;
 else
 vid = _inverts[i];
 EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
 
 // Get the incident half-vert in the previous mesh
 error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
 
 // Obtain the corresponding incident child in the current mesh
 for( int j = 0; j < (int)inci_ent.size(); j++ )
 {
 int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
 if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
 int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
 
 int pid;
 if( cur_level )
 pid = inci_ent[j] - level_mesh[cur_level - 1].start_cell;
 else
 pid = inci_ent[j] - *_incells.begin();
 
 int ind = nchilds * pid;
 
 // EntityHandle child_ent = level_mesh[cur_level].start_cell + ind+chid ;
 // int child_lid = refTemplates[type-1][d].v2hf[lvid][1];
 child_ents.push_back( level_mesh[cur_level].start_cell + ind + chid );
 child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
 }
 
 error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
 }
 
 // error = ahf->determine_incident_halffaces( level_mesh[cur_level].cells);MB_CHK_ERR(error);
 
 // Step 2: Update SIBHFS maps
 for( int i = 0; i < nents_prev; i++ )
 {
 EntityHandle ent;
 if( cur_level )
 ent = level_mesh[cur_level - 1].start_cell + i;
 else
 ent = _incells[i];
 
 std::vector< EntityHandle > sib_entids( nhf );
 std::vector< int > sib_lids( nhf );
 
 error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
 
 int id, idx;
 
 for( int l = 0; l < nhf; l++ )
 {
 
 if( !sib_entids[l] ) continue;
 
 // Get the number of children incident on this half-face
 int pat_id;
 if( type == MBTET )
 pat_id = ( *pattern_ids )[i];
 else
 pat_id = type - 1;
 int nch = refTemplates[pat_id][d].ents_on_pent[l][0];
 
 // Get the order of children indices incident on this half-face
 std::vector< int > id_sib( nch );
 for( int k = 0; k < nch; k++ )
 id_sib[k] = 0;
 
 error = reorder_indices( cur_level, deg, ent, l, sib_entids[l], sib_lids[l], 1, &id_sib[0] );MB_CHK_ERR( error );
 
 // Get the parent index of the sibling cell
 int psib;
 if( cur_level )
 psib = sib_entids[l] - level_mesh[cur_level - 1].start_cell;
 else
 psib = sib_entids[l] - *_incells.begin();
 
 int plid = sib_lids[l];
 int sidx = nchilds * psib;
 int sibpat_id;
 if( type == MBTET )
 sibpat_id = ( *pattern_ids )[psib];
 else
 sibpat_id = type - 1;
 
 // Loop over all the childs incident on the working half-face
 idx = nchilds * i;
 
 for( int k = 0; k < nch; k++ )
 {
 id = refTemplates[pat_id][d].ents_on_pent[l][k + 1] - 1;
 EntityHandle child_ent = level_mesh[cur_level].start_cell + idx + id;
 int child_lid = l;
 
 // Find the sibling of the working child
 EntityHandle child_sibent;
 int child_siblid;
 error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
 
 if( child_sibent != 0 ) continue;
 
 // Get the correponding child of the sibling of the current parent
 // We have already computed the order the children on incident corresponding to the
 // working half-face
 id = refTemplates[sibpat_id][d].ents_on_pent[plid][id_sib[k]] - 1;
 
 EntityHandle psib_child = level_mesh[cur_level].start_cell + sidx + id;
 int psib_chlid = plid;
 
 // Set the siblings of children incident on current half-face
 error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
 
 // Set the sibling of the sibling of the children to the children
 error = ahf->set_sibling_map( type, psib_child, psib_chlid, child_ent, child_lid );MB_CHK_ERR( error );
 }
 }
 
 // Loop over edges to check if there are any non-manifold edges. If there are then the v2hfs
 // map should be updated for the new vertices on it.
 const EntityHandle* conn;
 error = mbImpl->get_connectivity( ent, conn, nvpc );MB_CHK_ERR( error );
 
 int nv = refTemplates[type - 1][d].nv_edge; //#verts on each edge
 for( int l = 0; l < ne; l++ )
 {
 id = ahf->lConnMap3D[index].e2v[l][0];
 EntityHandle v_start = conn[id];
 id = ahf->lConnMap3D[index].e2v[l][1];
 EntityHandle v_end = conn[id];
 
 bool visited = false;
 
 std::vector< EntityHandle > inci_ent1, inci_ent2;
 std::vector< int > inci_lid1, inci_lid2;
 error = ahf->get_incident_map( type, v_start, inci_ent1, inci_lid1 );MB_CHK_ERR( error );
 error = ahf->get_incident_map( type, v_end, inci_ent2, inci_lid2 );MB_CHK_ERR( error );
 
 if( inci_ent1.size() > 1 && inci_ent2.size() > 1 )
 {
 std::vector< EntityHandle > cell_comps;
 std::vector< int > leid_comps;
 
 error = ahf->get_up_adjacencies_edg_3d_comp( ent, l, cell_comps, &leid_comps );MB_CHK_ERR( error );
 
 int ncomps = cell_comps.size();
 std::vector< EntityHandle > edgverts;
 std::vector< EntityHandle > compchildents( nv * ncomps );
 std::vector< int > compchildlfids( nv * ncomps );
 
 for( int s = 0; s < nv * ncomps; s++ )
 {
 compchildents[s] = 0;
 compchildlfids[s] = 0;
 }
 
 for( int j = 0; j < (int)cell_comps.size(); j++ )
 {
 int ind;
 if( cur_level )
 ind = level_mesh[cur_level - 1].cells.index( cell_comps[j] );
 else
 ind = _incells.index( cell_comps[j] );
 
 for( int k = 0; k < nv; k++ )
 {
 int chid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k] - 1;
 int lfid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 1];
 int lvid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 2];
 
 EntityHandle childcell = level_mesh[cur_level].start_cell + ind * nchilds + chid;
 
 const EntityHandle* econn;
 error = mbImpl->get_connectivity( childcell, econn, nvpc );MB_CHK_ERR( error );
 
 EntityHandle vert = econn[lvid];
 
 if( ahf->check_nonmanifold_vertices( type, vert ) )
 {
 visited = true;
 break;
 }
 
 if( edgverts.empty() )
 {
 edgverts.push_back( vert );
 compchildents[0] = childcell;
 compchildlfids[0] = lfid;
 }
 else
 {
 std::vector< EntityHandle >::iterator it;
 it = find( edgverts.begin(), edgverts.end(), vert );
 int indx = it - edgverts.begin();
 
 if( it == edgverts.end() )
 {
 edgverts.push_back( vert );
 compchildents[k * ncomps] = childcell;
 compchildlfids[k * ncomps] = lfid;
 }
 else
 {
 compchildents[indx * ncomps + j] = childcell;
 compchildlfids[indx * ncomps + j] = lfid;
 }
 }
 }
 }
 
 if( visited )
 {
 break;
 }
 
 // Set the incident half-facet map
 for( int k = 0; k < nv; k++ )
 {
 std::vector< EntityHandle > set_childents;
 std::vector< int > set_childlfids;
 for( int j = 0; j < ncomps; j++ )
 {
 set_childents.push_back( compchildents[k * ncomps + j] );
 set_childlfids.push_back( compchildlfids[k * ncomps + j] );
 }
 
 error = ahf->set_incident_map( type, edgverts[k], set_childents, set_childlfids );MB_CHK_ERR( error );
 }
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_lid_inci_child( EntityType type,
 int deg,
 int lfid,
 int leid,
 std::vector< int >& child_ids,
 std::vector< int >& child_lvids )
{
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 int d = get_index_from_degree( deg );
 
 // int lv0 = ahf->lConnMap3D[index].e2v[leid][0];
 // int lv1 = ahf->lConnMap3D[index].e2v[leid][1];
 int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
 
 int nv = refTemplates[type - 1][d].nv_edge;
 int nch = refTemplates[type - 1][d].ents_on_pent[lfid][0];
 
 for( int i = 0; i < nch; i++ )
 {
 int id = refTemplates[type - 1][d].ents_on_pent[lfid][i + 1] - 1;
 for( int j = 0; j < nvpc; j++ )
 {
 int lv = refTemplates[type - 1][d].ents_conn[id][j];
 for( int k = 0; k < nv; k++ )
 {
 if( lv == refTemplates[type - 1][d].vert_on_edges[leid][k] )
 {
 child_ids.push_back( id );
 child_lvids.push_back( j );
 }
 }
 }
 }
 
 return MB_SUCCESS;
}
 
/* **********************************
 * * Boundary Functions *
 ************************************/
 
bool NestedRefine::is_vertex_on_boundary( const EntityHandle& vertex )
{
 ErrorCode error;
 EntityHandle sibents[27];
 int siblids[27];
 std::vector< EntityHandle > ent;
 std::vector< int > lid;
 
 int nhf;
 if( elementype == MBEDGE )
 nhf = 2;
 else if( ( elementype == MBTRI ) || ( elementype == MBQUAD ) )
 nhf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
 else if( ( elementype == MBTET ) || ( elementype == MBHEX ) )
 {
 int idx = ahf->get_index_in_lmap( *_incells.begin() );
 nhf = ahf->lConnMap3D[idx].num_faces_in_cell;
 }
 else
 MB_SET_ERR( MB_FAILURE, "Requesting vertex boundary information for an unsupported entity type" );
 
 error = ahf->get_incident_map( elementype, vertex, ent, lid );MB_CHK_ERR( error );
 error = ahf->get_sibling_map( elementype, ent[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
 
 return ( sibents[lid[0]] == 0 );
}
 
bool NestedRefine::is_edge_on_boundary( const EntityHandle& entity )
{
 ErrorCode error;
 bool is_border = false;
 if( meshdim == 1 ) // The edge has a vertex on the boundary in the curve mesh
 {
 EntityHandle sibents[2];
 int siblids[2];
 error = ahf->get_sibling_map( MBEDGE, entity, &sibents[0], &siblids[0], 2 );MB_CHK_ERR( error );
 for( int i = 0; i < 2; i++ )
 {
 if( sibents[i] == 0 )
 {
 is_border = true;
 break;
 }
 }
 }
 else if( meshdim == 2 ) // The edge is on the boundary of the 2d mesh
 {
 std::vector< EntityHandle > adjents;
 error = ahf->get_up_adjacencies_2d( entity, adjents );MB_CHK_ERR( error );
 if( adjents.size() == 1 ) is_border = true;
 }
 else if( meshdim == 3 ) // The edge is part of a face on the boundary of the 3d mesh
 {
 std::vector< EntityHandle > adjents;
 std::vector< int > leids;
 error = ahf->get_up_adjacencies_edg_3d( entity, adjents, &leids );MB_CHK_ERR( error );
 assert( !adjents.empty() );
 
 int index = ahf->get_index_in_lmap( adjents[0] );
 int nhf = ahf->lConnMap3D[index].num_faces_in_cell;
 
 for( int i = 0; i < (int)adjents.size(); i++ )
 {
 EntityHandle sibents[6];
 int siblids[6];
 error = ahf->get_sibling_map( elementype, adjents[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
 for( int k = 0; k < 2; k++ )
 {
 int hf = ahf->lConnMap3D[index].e2hf[leids[0]][k];
 if( sibents[hf] == 0 )
 {
 is_border = true;
 break;
 }
 }
 }
 }
 return is_border;
}
 
bool NestedRefine::is_face_on_boundary( const EntityHandle& entity )
{
 ErrorCode error;
 bool is_border = false;
 
 if( meshdim == 1 )
 MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a face entity type on a curve mesh" );
 else if( meshdim == 2 ) // The face has a local edge on the boundary of the 2d mesh
 {
 EntityHandle sibents[4];
 int siblids[4];
 int nepf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
 error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nepf );MB_CHK_ERR( error );
 
 for( int i = 0; i < nepf; i++ )
 {
 if( sibents[i] == 0 )
 {
 is_border = true;
 break;
 }
 }
 }
 else if( meshdim == 3 ) // The face lies on the boundary of the 3d mesh
 {
 std::vector< EntityHandle > adjents;
 error = ahf->get_up_adjacencies_face_3d( entity, adjents );MB_CHK_ERR( error );
 if( adjents.size() == 1 ) is_border = true;
 }
 return is_border;
}
 
bool NestedRefine::is_cell_on_boundary( const EntityHandle& entity )
{
 if( meshdim != 3 )
 MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a cell entity type on a curve or surface mesh" );
 
 bool is_border = false;
 int index = ahf->get_index_in_lmap( *_incells.begin() );
 int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
 EntityHandle sibents[6];
 int siblids[6];
 
 ErrorCode error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nfpc );MB_CHK_ERR( error );
 
 for( int i = 0; i < nfpc; i++ )
 {
 if( sibents[i] == 0 )
 {
 is_border = true;
 break;
 }
 }
 return is_border;
}
 
/* **********************************
 * Helper Functions *
 ************************************/
ErrorCode NestedRefine::copy_vertices_from_prev_level( int cur_level )
{
 ErrorCode error;
 
 if( cur_level )
 {
 int nverts_prev = level_mesh[cur_level - 1].num_verts;
 for( int i = 0; i < nverts_prev; i++ )
 {
 level_mesh[cur_level].coordinates[0][i] = level_mesh[cur_level - 1].coordinates[0][i];
 level_mesh[cur_level].coordinates[1][i] = level_mesh[cur_level - 1].coordinates[1][i];
 level_mesh[cur_level].coordinates[2][i] = level_mesh[cur_level - 1].coordinates[2][i];
 }
 }
 else // Copy the vertices from the input mesh
 {
 int nverts_in = _inverts.size();
 std::vector< double > vcoords( 3 * nverts_in );
 error = mbImpl->get_coords( _inverts, &vcoords[0] );MB_CHK_ERR( error );
 
 for( int i = 0; i < nverts_in; i++ )
 {
 level_mesh[cur_level].coordinates[0][i] = vcoords[3 * i];
 level_mesh[cur_level].coordinates[1][i] = vcoords[3 * i + 1];
 level_mesh[cur_level].coordinates[2][i] = vcoords[3 * i + 2];
 }
 }
 return MB_SUCCESS;
 // To add: Map from old vertices to new duplicates: NOT NEEDED
}
 
ErrorCode NestedRefine::update_tracking_verts( EntityHandle cid,
 int cur_level,
 int deg,
 std::vector< EntityHandle >& trackvertsC_edg,
 std::vector< EntityHandle >& trackvertsC_face,
 EntityHandle* vbuffer )
{
 // The vertices in the vbuffer are added to appropriate edges and faces of cells that are
 // incident on the working cell.
 ErrorCode error;
 
 EntityHandle cstart_prev;
 if( cur_level )
 cstart_prev = level_mesh[cur_level - 1].start_cell;
 else
 cstart_prev = *_incells.begin();
 
 EntityType cell_type = mbImpl->type_from_handle( cstart_prev );
 int cindex = cell_type - 1;
 int d = get_index_from_degree( deg );
 
 int nve = refTemplates[cindex][d].nv_edge;
 int nvf = refTemplates[cindex][d].nv_face;
 
 int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
 int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
 int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
 
 // Step 1: Add the vertices on an edge of the working cell to tracking array of incident cells.
 for( int i = 0; i < nepc; i++ )
 {
 // Add the vertices to edges of the current cell
 for( int j = 0; j < nve; j++ )
 {
 int id = refTemplates[cindex][d].vert_on_edges[i][j];
 int idx = cid - cstart_prev;
 int aid = idx * nve * nepc + nve * i + j;
 
 if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
 }
 
 // Obtain all the incident cells
 std::vector< EntityHandle > inc_cids;
 std::vector< int > inc_leids, inc_orient;
 
 error = ahf->get_up_adjacencies_edg_3d( cid, i, inc_cids, &inc_leids, &inc_orient );MB_CHK_ERR( error );
 
 if( inc_cids.size() == 1 ) continue;
 
 // Add the vertices to the edges of the incident cells
 for( int k = 0; k < (int)inc_cids.size(); k++ )
 {
 if( inc_cids[k] == cid ) continue;
 
 int idx = inc_cids[k] - cstart_prev;
 
 if( inc_orient[k] ) // Same edge direction as the current edge
 {
 for( int j = 0; j < nve; j++ )
 {
 int id = refTemplates[cindex][d].vert_on_edges[i][j];
 int aid = idx * nve * nepc + nve * inc_leids[k] + j;
 
 if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
 }
 }
 else
 {
 for( int j = 0; j < nve; j++ )
 {
 int id = refTemplates[cindex][d].vert_on_edges[i][nve - j - 1];
 int aid = idx * nve * nepc + nve * inc_leids[k] + j;
 
 if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
 }
 }
 }
 }
 
 // Step 2: Add the vertices on a face of the working cell to tracking array of incident cells.
 if( nvf )
 {
 
 for( int i = 0; i < nfpc; i++ )
 {
 // Add vertices to the tracking array of vertices on faces for the current cell
 std::vector< EntityHandle > face_vbuf( nvf, 0 );
 for( int j = 0; j < nvf; j++ )
 {
 int id = refTemplates[cindex][d].vert_on_faces[i][j];
 int idx = cid - cstart_prev;
 int aid = idx * nvf * nfpc + nvf * i + j;
 
 if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = vbuffer[id];
 
 face_vbuf[j] = vbuffer[id];
 }
 
 // Obtain all the incident cells
 std::vector< EntityHandle > sib_cids;
 std::vector< int > sib_lfids;
 error = ahf->get_up_adjacencies_face_3d( cid, i, sib_cids, &sib_lfids );MB_CHK_ERR( error );
 
 if( sib_cids.size() == 1 ) continue;
 
 // Reorder the vertex local ids incident on the half-face
 std::vector< int > id_sib( nvf );
 for( int k = 0; k < nvf; k++ )
 id_sib[k] = 0;
 
 error = reorder_indices( cur_level, deg, sib_cids[1], sib_lfids[1], cid, i, 0, &id_sib[0] );MB_CHK_ERR( error );
 
 // Add vertices to the tracking array of vertices on faces for the sibling cell of the
 // current cell
 for( int j = 0; j < nvf; j++ )
 {
 int idx = sib_cids[1] - cstart_prev;
 int aid = idx * nvf * nfpc + nvf * sib_lfids[1] + j;
 
 if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = face_vbuf[id_sib[j] - 1];
 }
 }
 }
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::reorder_indices( int cur_level,
 int deg,
 EntityHandle cell,
 int lfid,
 EntityHandle sib_cell,
 int sib_lfid,
 int index,
 int* id_sib )
{
 // Reorders the indices of either vertices or children cell local ids to match with order of the
 // given cell and a local face. index = 0 : vertices,
 // = 1 : face
 
 assert( deg == 2 || deg == 3 );
 
 ErrorCode error;
 int idx = ahf->get_index_in_lmap( *_incells.begin() );
 int nvF = ahf->lConnMap3D[idx].hf2v_num[lfid];
 int nco = permutation[nvF - 3].num_comb;
 
 if( !index && ( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) ) )
 {
 id_sib[0] = 1;
 }
 else
 {
 // Get connectivity of the cell and its sibling cell
 std::vector< EntityHandle > conn, sib_conn;
 error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
 
 error = get_connectivity( sib_cell, cur_level, sib_conn );MB_CHK_ERR( error );
 
 // Get the connectivity of the local face in the cell and its sibling
 std::vector< EntityHandle > lface( nvF );
 std::vector< EntityHandle > lface_sib( nvF );
 for( int i = 0; i < nvF; i++ )
 {
 int id = ahf->lConnMap3D[idx].hf2v[lfid][i];
 lface[i] = conn[id];
 
 id = ahf->lConnMap3D[idx].hf2v[sib_lfid][i];
 lface_sib[i] = sib_conn[id];
 }
 
 // Find the combination
 int c = 0;
 for( int i = 0; i < nco; i++ )
 {
 int count = 0;
 for( int j = 0; j < nvF; j++ )
 {
 int id = permutation[nvF - 3].comb[i][j];
 if( lface[j] == lface_sib[id] ) count += 1;
 }
 
 if( count == nvF )
 {
 c = i;
 break;
 }
 }
 
 if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
 
 // Get the ordered indices
 if( ( ( !index ) && ( nvF == 4 ) && ( deg == 3 ) ) || ( deg == 2 ) )
 {
 for( int i = 0; i < 4; i++ )
 id_sib[i] = permutation[nvF - 3].porder2[c][i];
 }
 else
 {
 for( int i = 0; i < 9; i++ )
 id_sib[i] = permutation[nvF - 3].porder3[c][i];
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::reorder_indices( int deg,
 EntityHandle* face1_conn,
 EntityHandle* face2_conn,
 int nvF,
 std::vector< int >& lemap,
 std::vector< int >& vidx,
 int* leorient )
{
 // Given the connectivities of two faces, get the permuted indices w.r.t first face.
 // Step 1: First find the orientation
 int nco = permutation[nvF - 3].num_comb;
 int c = 0;
 for( int i = 0; i < nco; i++ )
 {
 int count = 0;
 for( int j = 0; j < nvF; j++ )
 {
 int id = permutation[nvF - 3].comb[i][j];
 if( face1_conn[j] == face2_conn[id] ) count += 1;
 }
 
 if( count == nvF )
 {
 c = i;
 break;
 }
 }
 
 if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
 
 // Add the corresponding local edges
 lemap.reserve( nvF );
 for( int i = 0; i < nvF; i++ )
 {
 lemap.push_back( permutation[nvF - 3].lemap[c][i] );
 }
 if( leorient ) leorient[0] = permutation[nvF - 3].orient[c];
 
 if( nvF == 3 && deg == 2 ) return MB_SUCCESS;
 
 if( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) )
 {
 vidx.push_back( 1 );
 }
 else if( nvF == 4 && deg == 3 )
 {
 for( int i = 0; i < 4; i++ )
 vidx.push_back( permutation[nvF - 3].porder2[c][i] );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::reorder_indices( int deg, int nvF, int comb, int* childfid_map )
{
 // Given connectivities of two faces and a degree, get the permuted indices of the children
 // faces w.r.t first face.
 
 assert( deg == 2 || deg == 3 );
 
 // Get the ordered indices
 if( deg == 2 )
 {
 for( int i = 0; i < 4; i++ )
 childfid_map[i] = permutation[nvF - 3].porder2[comb][i];
 }
 else
 {
 for( int i = 0; i < 9; i++ )
 childfid_map[i] = permutation[nvF - 3].porder3[comb][i];
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::reorder_indices( EntityHandle* face1_conn,
 EntityHandle* face2_conn,
 int nvF,
 int* conn_map,
 int& comb,
 int* orient )
{
 // Given connectivities of two faces and a degree, get the permuted indices of the children
 // faces w.r.t first face.
 
 // Step 1: First find the combination
 int nco = permutation[nvF - 3].num_comb;
 int c = 0;
 for( int i = 0; i < nco; i++ )
 {
 int count = 0;
 for( int j = 0; j < nvF; j++ )
 {
 int id = permutation[nvF - 3].comb[i][j];
 if( face1_conn[j] == face2_conn[id] ) count += 1;
 }
 
 if( count == nvF )
 {
 c = i;
 break;
 }
 }
 
 if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
 
 comb = c;
 
 if( orient ) orient[0] = permutation[nvF - 3].orient[c];
 
 for( int j = 0; j < nvF; j++ )
 {
 conn_map[j] = permutation[nvF - 3].comb[c][j];
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces )
{
 ErrorCode error;
 
 if( cur_level >= 0 )
 {
 Range edges, faces, cells;
 
 error = mbImpl->get_entities_by_dimension( set, 1, edges );MB_CHK_ERR( error );
 
 error = mbImpl->get_entities_by_dimension( set, 2, faces );MB_CHK_ERR( error );
 
 error = mbImpl->get_entities_by_dimension( set, 3, cells );MB_CHK_ERR( error );
 
 error = ahf->count_subentities( edges, faces, cells, nedges, nfaces );MB_CHK_ERR( error );
 }
 else
 {
 error = ahf->count_subentities( _inedges, _infaces, _incells, nedges, nfaces );MB_CHK_ERR( error );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords )
{
 int lid[6] = { 0, 0, 0, 0, 0, 0 };
 
 if( deg == 2 )
 {
 lid[0] = 5;
 lid[1] = 8;
 lid[2] = 9;
 lid[3] = 6;
 lid[4] = 4;
 lid[5] = 7;
 }
 else if( deg == 3 )
 {
 lid[0] = 19;
 lid[1] = 16;
 lid[2] = 18;
 lid[3] = 9;
 lid[4] = 4;
 lid[5] = 10;
 }
 
 EntityHandle vstart = level_mesh[cur_level].start_vertex;
 
 for( int i = 0; i < 6; i++ )
 {
 EntityHandle vid = vbuffer[lid[i]];
 ocoords[3 * i] = level_mesh[cur_level].coordinates[0][vid - vstart];
 ocoords[3 * i + 1] = level_mesh[cur_level].coordinates[1][vid - vstart];
 ocoords[3 * i + 2] = level_mesh[cur_level].coordinates[2][vid - vstart];
 }
 
 return MB_SUCCESS;
}
 
int NestedRefine::find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer )
{
 ErrorCode error;
 double coords[18];
 error = get_octahedron_corner_coords( cur_level, deg, vbuffer, coords );
 if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in obtaining octahedron corner coordinates" );
 
 int diag_map[6] = { 1, 3, 2, 4, 5, 0 };
 double length = std::numeric_limits< double >::max();
 
 int diag = 0;
 double x, y, z;
 x = y = z = 0;
 
 for( int d = 0; d < 3; d++ )
 {
 int id1 = diag_map[2 * d];
 int id2 = diag_map[2 * d + 1];
 x = coords[3 * id1] - coords[3 * id2];
 y = coords[3 * id1 + 1] - coords[3 * id2 + 1];
 z = coords[3 * id1 + 2] - coords[3 * id2 + 2];
 double dlen = sqrt( x * x + y * y + z * z );
 if( dlen < length )
 {
 length = dlen;
 diag = d + 1;
 }
 }
 
 return diag;
}
 
int NestedRefine::get_local_vid( EntityHandle vid, EntityHandle ent, int level )
{
 ErrorCode error;
 // Given a vertex, find its local id in the given entity
 std::vector< EntityHandle > conn;
 
 error = get_connectivity( ent, level + 1, conn );
 if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in getting connectivity of the requested entity" );
 
 int lid = -1;
 for( int i = 0; i < (int)conn.size(); i++ )
 {
 if( conn[i] == vid )
 {
 lid = i;
 break;
 }
 }
 if( lid < 0 ) MB_SET_ERR( MB_FAILURE, "Error in getting local vertex id in the given entity" );
 return lid;
}
 
int NestedRefine::get_index_from_degree( int degree )
{
 int d = deg_index.find( degree )->second;
 return d;
}
 
/*
ErrorCode NestedRefine::print_maps_1D(int level)
{
 ErrorCode error;
 int nv, ne;
 nv = level_mesh[level].num_verts;
 ne = level_mesh[level].num_edges;
 
 EntityHandle start_edge = level_mesh[level].start_edge;
 
 //V2HV
 std::cout<<"<V2HV_EID, V2HV_LVID>"<<std::endl;
 for (int i=0; i<nv; i++)
 {
 EntityHandle eid=0;
 int lvid=0;
 EntityHandle vid = level_mesh[level].start_vertex+i;
 error = ahf->get_incident_map(MBEDGE, vid, eid, lvid); MB_CHK_ERR(error);
 
 std::cout<<"For vertex = "<<vid<<"::Incident halfvertex "<<eid<<" "<<lvid<<std::endl;
 }
 
 //SIBHVS
 std::cout<<"start_edge = "<<start_edge<<std::endl;
 std::cout<<"<SIBHVS_EID,SIBHVS_LVID>"<<std::endl;
 for (int i=0; i<ne; i++)
 {
 EntityHandle ent = start_edge+i;
 
 EntityHandle eid[2]; int lvid[2];
 error = ahf->get_sibling_map(MBEDGE, ent, &eid[0], &lvid[0], 2); MB_CHK_ERR(error);
 std::cout<<"<"<<eid[0]<<","<<lvid[0]<<">"<<" "<<"<"<<eid[1]<<","<<lvid[1]<<">"<<std::endl;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::print_maps_2D(int level, EntityType type)
{
 ErrorCode error;
 int nv, nf;
 nv = level_mesh[level].num_verts;
 nf = level_mesh[level].num_faces;
 
 EntityHandle start_face = level_mesh[level].start_face;
 
 //V2HV
 std::cout<<"<V2HE_FID, V2HE_LEID>"<<std::endl;
 for (int i=0; i<nv; i++)
 {
 EntityHandle fid=0;
 int leid=0;
 EntityHandle vid = level_mesh[level].start_vertex+i;
 error = ahf->get_incident_map(type, vid, fid, leid); MB_CHK_ERR(error);
 
 std::cout<<"For vertex = "<<vid<<"::Incident halfedge "<<fid<<" "<<leid<<std::endl;
 }
 
 //SIBHES
 std::cout<<"start_face = "<<start_face<<std::endl;
 std::cout<<"<SIBHES_FID,SIBHES_LEID>"<<std::endl;
 EntityType ftype = mbImpl->type_from_handle(*_infaces.begin());
 int nepf = ahf->lConnMap2D[ftype-2].num_verts_in_face;
 
 EntityHandle *fid = new EntityHandle[nepf];
 int *leid = new int[nepf];
 
 for (int i=0; i<nf; i++)
 {
 for (int j=0; j<nepf; j++)
 {
 fid[j] = 0;
 leid[j] = 0;
 }
 
 EntityHandle ent = start_face+i;
 error = ahf->get_sibling_map(type, ent, fid, leid, nepf); MB_CHK_ERR(error);
 
 for (int j=0; j<nepf; j++){
 std::cout<<"<"<<fid[j]<<","<<leid[j]<<">"<<" ";
 }
 std::cout<<std::endl;
 }
 
 delete [] fid;
 delete [] leid;
 
 return MB_SUCCESS;
}
 
ErrorCode NestedRefine::print_maps_3D(int level, EntityType type)
{
 ErrorCode error;
 int nv, nc;
 nv = level_mesh[level].num_verts;
 nc = level_mesh[level].num_cells;
 EntityHandle start_cell = level_mesh[level].start_cell;
 
 //V2HF
 std::cout<<"<V2HF_CID, V2HF_LFID>"<<std::endl;
 for (int i=0; i<nv; i++)
 {
 EntityHandle cid=0;
 int lfid=0;
 EntityHandle vid = level_mesh[level].start_vertex+i;
 error = ahf->get_incident_map(type, vid, cid, lfid); MB_CHK_ERR(error);
 
 std::cout<<"For vertex = "<<vid<<"::Incident halfface "<<cid<<" "<<lfid<<std::endl;
 }
 
 //SIBHFS
 std::cout<<"start_cell = "<<start_cell<<std::endl;
 std::cout<<"<SIBHFS_CID,SIBHFS_LFID>"<<std::endl;
 int index = ahf->get_index_in_lmap(start_cell);
 int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
 
 EntityHandle *cid = new EntityHandle[nfpc];
 int *lfid = new int[nfpc];
 for (int i=0; i<nc; i++)
 {
 for (int k=0; k<nfpc; k++)
 {
 cid[k] = 0;
 lfid[k] = 0;
 }
 
 EntityHandle ent = start_cell+i;
 error = ahf->get_sibling_map(type, ent, cid, lfid, nfpc); MB_CHK_ERR(error);
 
 for (int j=0; j<nfpc; j++){
 std::cout<<"<"<<cid[j]<<","<<lfid[j]<<">"<<" ";
 }
 std::cout<<std::endl;
 }
 
 delete [] cid;
 delete [] lfid;
 
 return MB_SUCCESS;
}
 
*/
} // namespace moab