DualTool.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

#include "moab/DualTool.hpp"
#include "moab/Range.hpp"
// using Core for call to check_adjacencies
#include "moab/Core.hpp"
#include "Internals.hpp"
#include "MBTagConventions.hpp"
#include "moab/Skinner.hpp"
#include "moab/Core.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "AEntityFactory.hpp"
#include "moab/CN.hpp"
#include <string>
#include <algorithm>
#include <iostream>
#include <sstream>
#include <cassert>

#define RR \
    if( MB_SUCCESS != result ) return result
#define SWAP( a, b )                    \
    {                                   \
        EntityHandle tmp_ent = a;       \
        ( a )                = b;       \
        ( b )                = tmp_ent; \
    }

namespace moab
{

bool debug    = false;
bool debug_ap = false;

//! tag name for dual surfaces
const char* DualTool::DUAL_SURFACE_TAG_NAME = "DUAL_SURFACE";

//! tag name for dual curves
const char* DualTool::DUAL_CURVE_TAG_NAME = "DUAL_CURVE";

//! tag name for dual cells
const char* DualTool::IS_DUAL_CELL_TAG_NAME = "__IS_DUAL_CELL";

//! tag name for dual entities
const char* DualTool::DUAL_ENTITY_TAG_NAME = "__DUAL_ENTITY";

//! tag name for extra dual entities
const char* DualTool::EXTRA_DUAL_ENTITY_TAG_NAME = "__EXTRA_DUAL_ENTITY";

//! tag name for graphics point
const char* DualTool::DUAL_GRAPHICS_POINT_TAG_NAME = "__DUAL_GRAPHICS_POINT";

// const int DualTool::GP_SIZE = 20;

DualTool::DualTool( Interface* impl ) : mbImpl( impl )
{
    EntityHandle dum_handle = 0;
    ErrorCode result;

    result = mbImpl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dualSurfaceTag,
                                     MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
    assert( MB_SUCCESS == result );

    result = mbImpl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dualCurveTag, MB_TAG_SPARSE | MB_TAG_CREAT,
                                     &dum_handle );
    assert( MB_SUCCESS == result );

    unsigned int dummy = 0;
    result             = mbImpl->tag_get_handle( IS_DUAL_CELL_TAG_NAME, 1, MB_TYPE_INTEGER, isDualCellTag,
                                                 MB_TAG_SPARSE | MB_TAG_CREAT, &dummy );
    assert( MB_SUCCESS == result );

    result = mbImpl->tag_get_handle( DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, dualEntityTag,
                                     MB_TAG_DENSE | MB_TAG_CREAT, &dum_handle );
    assert( MB_SUCCESS == result );

    result = mbImpl->tag_get_handle( EXTRA_DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, extraDualEntityTag,
                                     MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
    assert( MB_SUCCESS == result );

    static const char dum_name[CATEGORY_TAG_SIZE] = { 0 };
    result = mbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
                                     MB_TAG_SPARSE | MB_TAG_CREAT, dum_name );
    assert( MB_SUCCESS == result );

    DualTool::GraphicsPoint dum_pt( 0.0, 0.0, 0.0, -1 );
    result = mbImpl->tag_get_handle( DUAL_GRAPHICS_POINT_TAG_NAME, sizeof( DualTool::GraphicsPoint ), MB_TYPE_DOUBLE,
                                     dualGraphicsPointTag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_BYTES, &dum_pt );
    assert( MB_SUCCESS == result );

    globalIdTag = mbImpl->globalId_tag();

    if( MB_SUCCESS == result )
    {
    }  // empty statement to get rid of warning.

    maxHexId = -1;
}

DualTool::~DualTool() {}

//! construct the dual entities for the entire mesh
ErrorCode DualTool::construct_dual( EntityHandle* entities, const int num_entities )
{
    // allocate a dual entity for each primal entity in the mesh, starting
    // with highest dimension and working downward; do each dimension in a separate code
    // block, since they're all handled slightly differently

    Range regions, faces, edges, vertices;
    ErrorCode result;

    if( NULL == entities || 0 == num_entities )
    {

        // first, construct all the aentities, since they're currently needed to
        // compute the dual
        result = mbImpl->get_entities_by_dimension( 0, 0, vertices );
        if( MB_SUCCESS != result ) return result;

        result = MeshTopoUtil( mbImpl ).construct_aentities( vertices );
        if( MB_SUCCESS != result ) return result;

        // get all edges, faces and regions now, so we don't need to filter out dual
        // entities later

        result = mbImpl->get_entities_by_dimension( 0, 1, edges );
        if( MB_SUCCESS != result ) return result;
        result = mbImpl->get_entities_by_dimension( 0, 2, faces );
        if( MB_SUCCESS != result ) return result;
        result = mbImpl->get_entities_by_dimension( 0, 3, regions );
        if( MB_SUCCESS != result ) return result;

        // get the max global id for hexes, we'll need for modification ops
        std::vector< int > gid_vec( regions.size() );
        result = mbImpl->tag_get_data( globalId_tag(), regions, &gid_vec[0] );
        if( MB_SUCCESS != result ) return result;
        maxHexId = -1;
        Range::iterator rit;
        unsigned int i;
        for( rit = regions.begin(), i = 0; rit != regions.end(); ++rit, i++ )
        {
            if( gid_vec[i] > maxHexId && mbImpl->type_from_handle( *rit ) == MBHEX ) maxHexId = gid_vec[i];
        }
    }
    else
    {
        // get entities of various dimensions adjacent to these
        result = mbImpl->get_adjacencies( entities, num_entities, 0, true, vertices, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
        result = mbImpl->get_adjacencies( entities, num_entities, 1, true, edges, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
        result = mbImpl->get_adjacencies( entities, num_entities, 2, true, faces, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
        result = mbImpl->get_adjacencies( entities, num_entities, 3, true, regions, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
    }

    Range dual_verts;
    result = construct_dual_vertices( regions, dual_verts );
    if( MB_SUCCESS != result || dual_verts.size() != regions.size() ) return result;
    if( debug ) std::cout << "Constructed " << dual_verts.size() << " dual vertices." << std::endl;

    // don't really need dual edges, but construct 'em anyway
    Range dual_edges;
    result = construct_dual_edges( faces, dual_edges );
    if( MB_SUCCESS != result || dual_edges.size() != faces.size() ) return result;
    if( debug ) std::cout << "Constructed " << dual_edges.size() << " dual edges." << std::endl;

    // construct dual faces
    Range dual_faces;
    result = construct_dual_faces( edges, dual_faces );
    if( MB_SUCCESS != result || dual_faces.size() != edges.size() ) return result;
    if( debug ) std::cout << "Constructed " << dual_faces.size() << " dual faces." << std::endl;

    // construct dual cells
    Range dual_cells;
    result = construct_dual_cells( vertices, dual_cells );
    if( MB_SUCCESS != result || dual_cells.size() != vertices.size() ) return result;
    if( debug ) std::cout << "Constructed " << dual_cells.size() << " dual cells." << std::endl;

    return MB_SUCCESS;
}

ErrorCode DualTool::construct_dual_vertices( const Range& all_regions, Range& dual_ents )
{
    if( all_regions.empty() ) return MB_SUCCESS;

    // make sure they're all regions
    assert( 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.begin() ) ) &&
            3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.rbegin() ) ) );

    Range::const_iterator rit;
    EntityHandle dual_ent;
    ErrorCode tmp_result = MB_SUCCESS;
    ErrorCode result     = MB_SUCCESS;

    for( rit = all_regions.begin(); rit != all_regions.end(); ++rit )
    {
        if( tmp_result != MB_SUCCESS ) result = tmp_result;

        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        if( MB_SUCCESS == tmp_result && 0 != dual_ent )
        {
            dual_ents.insert( dual_ent );
            continue;
        }
        else if( MB_SUCCESS != tmp_result )
            continue;

        tmp_result = construct_dual_vertex( *rit, dual_ent, false, true );
        if( MB_SUCCESS != tmp_result ) continue;

        // save it in the list of new dual ents
        dual_ents.insert( dual_ent );
    }

    return result;
}

ErrorCode DualTool::construct_dual_vertex( EntityHandle entity,
                                           EntityHandle& dual_ent,
                                           const bool extra,
                                           const bool add_graphics_pt )
{
    // no dual entity; construct one; first need the avg coordinates
    unsigned int is_dual = 0x1;
    double avg_pos[3];
    ErrorCode result = MeshTopoUtil( mbImpl ).get_average_position( entity, avg_pos );
    if( MB_SUCCESS != result ) return result;

    // now construct the new dual entity
    result = mbImpl->create_vertex( avg_pos, dual_ent );
    if( MB_SUCCESS != result ) return result;

    // tag it indicating it's a dual entity
    result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
    if( MB_SUCCESS != result ) return result;

    // tag the primal entity with its dual entity and vica versa
    if( extra )
        result = mbImpl->tag_set_data( extraDualEntity_tag(), &( entity ), 1, &dual_ent );
    else
        result = mbImpl->tag_set_data( dualEntity_tag(), &( entity ), 1, &dual_ent );
    if( MB_SUCCESS != result ) return result;

    result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( entity ) );
    if( MB_SUCCESS != result ) return result;

    if( add_graphics_pt )
        // put a graphics point on that vertex too
        result = add_graphics_point( dual_ent, avg_pos );

    return result;
}

ErrorCode DualTool::add_graphics_point( EntityHandle entity, double* avg_pos )
{
    // add a graphics pt, placed at the same position as the vertex
    double my_pos[3];
    ErrorCode result;

    if( NULL == avg_pos )
    {
        result = MeshTopoUtil( mbImpl ).get_average_position( entity, my_pos );
        if( MB_SUCCESS != result ) return result;
    }
    else
        for( int i = 0; i < 3; i++ )
            my_pos[i] = avg_pos[i];

    DualTool::GraphicsPoint dum_pt( my_pos, -1 );
    result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), &entity, 1, &dum_pt );
    return result;
}

ErrorCode DualTool::construct_dual_edges( const Range& all_faces, Range& dual_ents )
{
    if( all_faces.empty() ) return MB_SUCCESS;

    // make sure they're all faces
    assert( 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.begin() ) ) &&
            2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.rbegin() ) ) );

    Range::const_iterator rit;
    EntityHandle dual_ent;
    unsigned int is_dual = 0x1;
    ErrorCode tmp_result = MB_SUCCESS;
    ErrorCode result     = MB_SUCCESS;

    for( rit = all_faces.begin(); rit != all_faces.end(); ++rit )
    {
        if( tmp_result != MB_SUCCESS ) result = tmp_result;

        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        if( MB_SUCCESS == tmp_result && 0 != dual_ent )
        {
            dual_ents.insert( dual_ent );
            continue;
        }

        // no dual entity; construct one; get the bounding regions
        std::vector< EntityHandle > out_ents;
        tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, out_ents );
        if( MB_SUCCESS != tmp_result || out_ents.empty() ) continue;

        // get the dual vertices
        std::vector< EntityHandle > dual_verts( out_ents.size() );
        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &out_ents[0], out_ents.size(), &dual_verts[0] );
        if( MB_SUCCESS != tmp_result ) continue;
        assert( dual_verts.size() <= 2 );

        double avg_pos[3];
        bool bdy_face = ( dual_verts.size() == 1 ? true : false );
        if( bdy_face )
        {
            // boundary face - make a dual vertex at the face center and put in list
            tmp_result = construct_dual_vertex( *rit, dual_ent, true, true );

            // put it on vertex list
            dual_verts.push_back( dual_ent );
        }

        assert( dual_verts.size() == 2 );

        // now create the dual edge
        tmp_result = mbImpl->create_element( MBEDGE, &dual_verts[0], 2, dual_ent );
        if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;

        // save it in the list of new dual ents
        dual_ents.insert( dual_ent );

        // tag the primal entity with its dual entity and vica versa
        tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        if( MB_SUCCESS != tmp_result ) continue;

        tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
        if( MB_SUCCESS != tmp_result ) continue;

        // tag the edge indicating it's a dual entity
        tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
        if( MB_SUCCESS != tmp_result ) continue;

        // add a graphics point to the edge; position depends on whether it's a
        // bdy face (mid-pt of dual edge) or not (mid-pt of primal face)
        if( bdy_face )
            tmp_result = add_graphics_point( dual_ent );
        else
        {
            // get the face's position
            tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
            if( MB_SUCCESS != tmp_result ) continue;
            tmp_result = add_graphics_point( dual_ent, avg_pos );
        }
        if( MB_SUCCESS != tmp_result ) continue;
    }

    return result;
}

ErrorCode DualTool::construct_dual_faces( const Range& all_edges, Range& dual_ents )
{
    if( all_edges.empty() ) return MB_SUCCESS;

    // make sure they're all edges
    assert( 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.begin() ) ) &&
            1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.rbegin() ) ) );

    Range::const_iterator rit;
    EntityHandle dual_ent;
    unsigned int is_dual = 0x1;
    ErrorCode tmp_result = MB_SUCCESS;
    ErrorCode result     = MB_SUCCESS;
    Range equiv_edges;
#define TRC                        \
    if( MB_SUCCESS != tmp_result ) \
    {                              \
        result = tmp_result;       \
        continue;                  \
    }
    for( rit = all_edges.begin(); rit != all_edges.end(); ++rit )
    {

        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        if( MB_SUCCESS == tmp_result && 0 != dual_ent )
        {
            dual_ents.insert( dual_ent );
            continue;
        }

        // no dual entity; construct one; get the dual vertices bounding the edge in radial order,
        // then construct the dual face
        std::vector< EntityHandle > rad_dverts;
        bool bdy_edge;
        tmp_result = get_radial_dverts( *rit, rad_dverts, bdy_edge );
        TRC if( rad_dverts.empty() ) continue;

        tmp_result = mbImpl->create_element( MBPOLYGON, &rad_dverts[0], rad_dverts.size(), dual_ent );
        TRC

            // tag it indicating it's a dual entity, and tag primal/dual with dual/primal
            tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
        TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
        TRC

            // save it in the list of new dual ents
            dual_ents.insert( dual_ent );

        // add a graphics point to the cell; position depends on whether it's a
        // bdy cell (mid-pt of cell's vertices) or not (mid-pt of primal edge)
        double avg_pos[3];
        tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
        TRC if( bdy_edge )
        {

            // add a new dual edge betw last 2 verts
            EntityHandle new_edge;
            tmp_result     = mbImpl->create_element( MBEDGE, &rad_dverts[rad_dverts.size() - 2], 2, new_edge );
            TRC tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &new_edge, 1, &is_dual );
            TRC

                // tag the new dual edge with the primal edge as it's dual entity; primal
                // edge IS NOT likewise tagged, since it's already tagged with the 2cell
                tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &new_edge, 1, &( *rit ) );
            TRC

                // add a graphics pt, position is center of primal edge
                tmp_result = add_graphics_point( dual_ent );
            TRC tmp_result = add_graphics_point( new_edge, avg_pos );
            TRC
        }

        else
        {
            // if inside, point goes on the 2cell, at primal edge mid-pt
            tmp_result = add_graphics_point( dual_ent, avg_pos );
            TRC
        }

        // check to see whether we have equiv entities; if we find any, save for later fixup
        Range dum_edges, dum_poly( dual_ent, dual_ent );
        tmp_result = mbImpl->get_adjacencies( dum_poly, 1, false, dum_edges );
        if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
        {
            // we do - need to add adjacencies to disambiguate; use the primal
            equiv_edges.merge( dum_edges );
        }
    }

    if( !equiv_edges.empty() ) result = check_dual_equiv_edges( equiv_edges );

    return result;
}

ErrorCode DualTool::check_dual_equiv_edges( Range& dual_edges )
{
    // fix equivalent dual edges (i.e. edges whose vertices define multiple edges)
    // by explicitly adding adjacencies to containing polygons; adjacent polygons
    // found by going through primal
    ErrorCode tmp_result, result = MB_SUCCESS;

    Range all_dedges( dual_edges );
    // first, go through all dual edges and find equivalent edges (by looking for
    // up-adjacent edges on the vertices of each edge)
    for( Range::iterator rit = dual_edges.begin(); rit != dual_edges.end(); ++rit )
    {
        Range connect, dum_range( *rit, *rit );
        tmp_result = mbImpl->get_adjacencies( dum_range, 0, false, connect );
        if( MB_SUCCESS != tmp_result ) continue;
        tmp_result = mbImpl->get_adjacencies( connect, 1, false, all_dedges, Interface::UNION );
        if( MB_SUCCESS != tmp_result ) continue;
    }

    // save a copy for checking later
    Range save_all_2cells;

    // go through each edge
    while( !all_dedges.empty() )
    {
        EntityHandle this_edge = *all_dedges.begin();
        all_dedges.erase( all_dedges.begin() );

        const EntityHandle* connect;
        int num_connect;
        result = mbImpl->get_connectivity( this_edge, connect, num_connect );
        if( MB_SUCCESS != result ) continue;

        Range dum_edges, verts;
        verts.insert( connect[0] );
        verts.insert( connect[1] );
        tmp_result = mbImpl->get_adjacencies( verts, 1, false, dum_edges );
        if( MB_SUCCESS != tmp_result )
        {
            result = tmp_result;
            continue;
        }
        if( dum_edges.size() == 1 )
        {
            // not an equiv edge - already removed from list, so just continue
            continue;
        }

        // ok, have an equiv entity - fix by looking through primal
        // pre-get the primal of these
        EntityHandle dedge_quad;
        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &this_edge, 1, &dedge_quad );
        if( MB_SUCCESS != tmp_result )
        {
            result = tmp_result;
            continue;
        }

        if( MBQUAD == mbImpl->type_from_handle( dedge_quad ) )
        {

            // get the primal edges adj to quad
            Range dum_quad_range( dedge_quad, dedge_quad ), adj_pedges;
            tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_pedges );
            if( MB_SUCCESS != tmp_result )
            {
                result = tmp_result;
                continue;
            }
            // get the dual 2cells corresponding to those pedges
            std::vector< EntityHandle > dcells;
            dcells.resize( adj_pedges.size() );
            tmp_result = mbImpl->tag_get_data( dualEntity_tag(), adj_pedges, &dcells[0] );
            if( MB_SUCCESS != tmp_result )
            {
                result = tmp_result;
                continue;
            }
            // now add explicit adjacencies from the dedge to those dcells
            std::vector< EntityHandle >::iterator vit;
            for( vit = dcells.begin(); vit != dcells.end(); ++vit )
            {
                save_all_2cells.insert( *vit );

                assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
                tmp_result = mbImpl->add_adjacencies( this_edge, &( *vit ), 1, false );
                if( MB_SUCCESS != tmp_result )
                {
                    result = tmp_result;
                    continue;
                }
                // check that there are really adjacencies and *vit is in them
                const EntityHandle* adjs;
                int num_adjs;
                tmp_result = reinterpret_cast< Core* >( mbImpl )->a_entity_factory()->get_adjacencies( this_edge, adjs,
                                                                                                       num_adjs );
                if( NULL == adjs || std::find( adjs, adjs + num_adjs, *vit ) == adjs + num_adjs )
                    std::cout << "Add_adjacencies failed in construct_dual_faces." << std::endl;
            }
        }
        else
        {
            // else, have a dual edge representing a bdy edge - tie directly to
            // dual entity if its dual entity
            EntityHandle bdy_dcell;
            tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &dedge_quad, 1, &bdy_dcell );
            TRC assert( MBPOLYGON == mbImpl->type_from_handle( bdy_dcell ) );

            tmp_result = mbImpl->add_adjacencies( this_edge, &bdy_dcell, 1, false );
            if( MB_SUCCESS != tmp_result )
            {
                result = tmp_result;
                continue;
            }
        }
    }

    // sanity check - look for adj edges again, and check for equiv entities
    for( Range::iterator vit = save_all_2cells.begin(); vit != save_all_2cells.end(); ++vit )
    {
        Range adj_edges, dum_quad_range;
        dum_quad_range.insert( *vit );
        assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
        tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_edges );
        if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
        {
            std::cout << "Multiple entities returned for polygon " << mbImpl->id_from_handle( *vit ) << "."
                      << std::endl;
            continue;
        }
    }
    // success!
    return result;
}

ErrorCode DualTool::construct_dual_cells( const Range& all_verts, Range& dual_ents )
{
    if( all_verts.empty() ) return MB_SUCCESS;

    // make sure they're all edges
    assert( 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.begin() ) ) &&
            0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.rbegin() ) ) );

    Range::const_iterator rit;
    EntityHandle dual_ent;
    unsigned int is_dual = 0x1;
    ErrorCode tmp_result = MB_SUCCESS;
    ErrorCode result     = MB_SUCCESS;
    std::vector< EntityHandle > edges, dfaces;

    for( rit = all_verts.begin(); rit != all_verts.end(); ++rit )
    {
        if( tmp_result != MB_SUCCESS ) result = tmp_result;

        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        if( MB_SUCCESS == tmp_result && 0 != dual_ent )
        {
            dual_ents.insert( dual_ent );
            continue;
        }

        // no dual entity; construct one; get the edges bounding the vertex
        edges.clear();
        dfaces.clear();
        tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 1, false, edges );
        if( MB_SUCCESS != tmp_result ) continue;

        // get the dual faces corresponding to the edges
        dfaces.resize( edges.size() );
        tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &edges[0], edges.size(), &dfaces[0] );
        if( MB_SUCCESS != tmp_result ) continue;

        // create the dual cell from those faces
        tmp_result = mbImpl->create_element( MBPOLYHEDRON, &dfaces[0], dfaces.size(), dual_ent );
        if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;

        // save it in the list of new dual ents
        dual_ents.insert( dual_ent );

        // tag it indicating it's a dual entity
        tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
        if( MB_SUCCESS != tmp_result ) continue;

        // tag the primal entity with its dual entity and vica versa
        tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
        if( MB_SUCCESS != tmp_result ) continue;
        tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
        if( MB_SUCCESS != tmp_result ) continue;
    }

    return result;
}

//! given an edge handle, return a list of dual vertices in radial order
//! around the edge
ErrorCode DualTool::get_radial_dverts( const EntityHandle edge,
                                       std::vector< EntityHandle >& rad_dverts,
                                       bool& bdy_edge )
{
    rad_dverts.clear();

    std::vector< EntityHandle > rad_faces, rad_ents;
    ErrorCode result = MeshTopoUtil( mbImpl ).star_entities( edge, rad_faces, bdy_edge, 0, &rad_ents );
    if( MB_SUCCESS != result ) return result;

    if( bdy_edge )
    {
        // if we're a bdy edge, change the order back to what DualTool expects
        rad_ents.push_back( *rad_faces.rbegin() );
        rad_ents.push_back( *rad_faces.begin() );
    }

    rad_dverts.resize( rad_ents.size() );
    for( unsigned int i = 0; i < rad_ents.size(); i++ )
    {
        EntityHandle dual_ent;
        result = mbImpl->tag_get_data( dualEntity_tag(), &rad_ents[i], 1, &dual_ent );
        if( !bdy_edge || i < rad_ents.size() - 2 )
            rad_dverts[i] = dual_ent;
        else
        {
            // fix up this entry
            assert( mbImpl->type_from_handle( dual_ent ) == MBEDGE );

            // get connectivity of that edge
            const EntityHandle* connect;
            int num_connect;
            result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
            if( MB_SUCCESS != result ) return result;

            // we want the one that's not already on the list; reuse last_face
            int last_hex           = ( i == rad_ents.size() - 1 ? 0 : i - 1 );
            EntityHandle last_face = ( connect[0] == rad_dverts[last_hex] ? connect[1] : connect[0] );
            rad_dverts[i]          = last_face;
        }
    }

    return result;
}

//! construct the dual entities for a hex mesh, including dual surfaces & curves
ErrorCode DualTool::construct_hex_dual( Range& entities )
{
    std::vector< EntityHandle > evec;
    std::copy( entities.begin(), entities.end(), std::back_inserter( evec ) );
    return construct_hex_dual( &evec[0], evec.size() );
}

//! construct the dual entities for a hex mesh, including dual surfaces & curves
ErrorCode DualTool::construct_hex_dual( EntityHandle* entities, const int num_entities )
{
    // really quite simple:

    // construct the dual...
    ErrorCode result = construct_dual( entities, num_entities );
    if( MB_SUCCESS != result )
    {
        std::cerr << "Error constructing dual entities for primal entities." << std::endl;
        return result;
    }

    // now traverse to build 1d and 2d hyperplanes
    result = construct_dual_hyperplanes( 1, entities, num_entities );
    if( MB_SUCCESS != result )
    {
        std::cerr << "Problem traversing 1d hyperplanes." << std::endl;
        return result;
    }

    result = construct_dual_hyperplanes( 2, entities, num_entities );
    if( MB_SUCCESS != result )
    {
        std::cerr << "Problem traversing 2d hyperplanes." << std::endl;
        return result;
    }

    result = construct_hp_parent_child();
    if( MB_SUCCESS != result )
    {
        std::cerr << "Problem constructing parent/child relations between hyperplanes." << std::endl;
        return result;
    }

    // see?  simple, just like I said
    return MB_SUCCESS;
}

//! get the cells of the dual
ErrorCode DualTool::get_dual_entities( const int dim, EntityHandle* entities, const int num_entities, Range& dual_ents )
{
    if( 0 == isDualCell_tag() ) return MB_SUCCESS;
    if( 0 > dim || 3 < dim ) return MB_INDEX_OUT_OF_RANGE;

    unsigned int dum              = 0x1;
    const void* dum_ptr           = &dum;
    static EntityType dual_type[] = { MBVERTEX, MBEDGE, MBPOLYGON, MBPOLYHEDRON };

    Range dim_ents;

    ErrorCode result;

    if( 0 == entities || 0 == num_entities )
    {
        // just get all the dual entities of this dimension
        result = mbImpl->get_entities_by_type_and_tag( 0, dual_type[dim], &isDualCellTag, &dum_ptr, 1, dual_ents );
    }
    else
    {
        // else look for specific dual entities
        result = mbImpl->get_adjacencies( entities, num_entities, 3 - dim, false, dim_ents, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
        std::vector< EntityHandle > dual_ents_vec( dim_ents.size() );
        result = mbImpl->tag_get_data( dualEntity_tag(), dim_ents, &dual_ents_vec[0] );
        if( MB_SUCCESS != result ) return result;
        std::copy( dual_ents_vec.begin(), dual_ents_vec.end(), range_inserter( dual_ents ) );
    }

    return result;
}

//! get the faces of the dual
ErrorCode DualTool::get_dual_entities( const int dim,
                                       EntityHandle* entities,
                                       const int num_entities,
                                       std::vector< EntityHandle >& dual_ents )
{
    Range tmp_range;
    ErrorCode result = get_dual_entities( dim, entities, num_entities, tmp_range );
    if( MB_SUCCESS != result ) return result;

    // dual_ents.insert(dual_ents.end(), tmp_range.begin(), tmp_range.end());
    dual_ents.reserve( dual_ents.size() + tmp_range.size() );
    for( Range::const_iterator it = tmp_range.begin(); it != tmp_range.end(); ++it )
    {
        dual_ents.push_back( *it );
    }
    return MB_SUCCESS;
}

ErrorCode DualTool::get_dual_hyperplanes( const Interface* impl, const int dim, Range& dual_ents )
{
    if( dim != 1 && dim != 2 ) return MB_INDEX_OUT_OF_RANGE;

    Tag dual_tag;
    ErrorCode result;

    if( dim == 1 )
        result = impl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
    else
        result = impl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );

    if( MB_SUCCESS == result )
        result = impl->get_entities_by_type_and_tag( 0, MBENTITYSET, &dual_tag, NULL, 1, dual_ents, Interface::UNION );

    return result;
}

ErrorCode DualTool::construct_dual_hyperplanes( const int dim, EntityHandle* entities, const int num_entities )
{
    // this function traverses dual faces of input dimension, constructing
    // dual hyperplanes of them in sets as it goes

    // check various inputs
    int num_quads, num_hexes;
    if(
        // this function only makes sense for dim == 1 or dim == 2
        ( dim != 1 && dim != 2 ) ||
        // should either be quads or hexes around
        mbImpl->get_number_entities_by_type( 0, MBQUAD, num_quads ) != MB_SUCCESS ||
        mbImpl->get_number_entities_by_type( 0, MBHEX, num_hexes ) != MB_SUCCESS ||
        // if we're asking for 1d dual ents, should be quads around
        ( num_quads == 0 && dim == 1 ) ||
        // if we're asking for 2d dual ents, should be hexes around
        ( num_hexes == 0 && dim == 2 ) )
        return MB_FAILURE;

    Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );

    // two stacks: one completely untreated entities, and the other untreated
    // entities on the current dual hyperplane
    std::vector< EntityHandle > tot_untreated;

    // put dual entities of this dimension on the untreated list
    ErrorCode result = get_dual_entities( dim, entities, num_entities, tot_untreated );
    if( MB_SUCCESS != result ) return result;

    // main part of traversal loop
    EntityHandle this_ent;
    EntityHandle this_hp;

    while( !tot_untreated.empty() )
    {
        if( debug && dim == 2 /*(tot_untreated.size()%report == 0)*/ )
            std::cout << "Untreated list size " << tot_untreated.size() << "." << std::endl;

        this_ent = tot_untreated.back();
        tot_untreated.pop_back();
        result = mbImpl->tag_get_data( hp_tag, &this_ent, 1, &this_hp );
        if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return result;

        // d for this entity having a hyperplane assignment already
        else if( this_hp != 0 )
            continue;

        if( 1 == dim && check_1d_loop_edge( this_ent ) ) continue;

        // inner loop: traverse the hyperplane 'till we don't have any more
        result = traverse_hyperplane( hp_tag, this_hp, this_ent );
        if( MB_SUCCESS != result )
        {
            std::cout << "Failed to traverse hyperplane ";
            if( this_hp )
                std::cout << mbImpl->id_from_handle( this_hp ) << "." << std::endl;
            else
                std::cout << "0." << std::endl;
            return result;
        }

        // ok, now order the edges if it's a chord
        if( 1 == dim ) order_chord( this_hp );
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::traverse_hyperplane( const Tag hp_tag, EntityHandle& this_hp, EntityHandle this_ent )
{
    Range tmp_star, star, tmp_range, new_hyperplane_ents;
    std::vector< EntityHandle > hp_untreated;
    int dim = mbImpl->dimension_from_handle( this_ent );
    MeshTopoUtil mtu( mbImpl );
    this_hp = 0;
    ErrorCode result;

    unsigned short mark_val = 0x0;
    Tag mark_tag;
    result = mbImpl->tag_get_handle( "__hyperplane_mark", 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT | MB_TAG_BIT );
    if( MB_SUCCESS != result ) return result;
    mark_val = 0x1;

    while( 0 != this_ent )
    {
        EntityHandle tmp_hp = get_dual_hyperplane( this_ent );
        if( 0 == this_hp && 0 != tmp_hp ) this_hp = tmp_hp;

        if( 0 == tmp_hp ) new_hyperplane_ents.insert( this_ent );

        if( debug && hp_untreated.size() % 10 == 0 )
            std::cout << "Dual surface " << this_hp << ", hp_untreated list size = " << hp_untreated.size() << "."
                      << std::endl;

        // get the 2nd order adjacencies through lower dimension
        tmp_range.clear();
        tmp_star.clear();
        star.clear();
        result = mtu.get_bridge_adjacencies( this_ent, dim - 1, dim, star );RR;

        // get the bridge adjacencies through higher dimension
        result = mtu.get_bridge_adjacencies( this_ent, dim + 1, dim, tmp_star );RR;
        tmp_range = subtract( star, tmp_star );

        for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
        {
            if( new_hyperplane_ents.find( *rit ) != new_hyperplane_ents.end() ) continue;

            // check for tag first, 'cuz it's probably faster than checking adjacencies
            // assign to avoid valgrind warning
            unsigned short tmp_mark = 0x0;
            result                  = mbImpl->tag_get_data( mark_tag, &( *rit ), 1, &tmp_mark );
            if( MB_SUCCESS == result && mark_val == tmp_mark ) continue;

            // if it's on the loop, it's not eligible
            if( 1 == dim && check_1d_loop_edge( *rit ) ) continue;

            // have one on this hp; just put it on the hp_untreated list for now,
            // will get tagged and put in the hp set later
            hp_untreated.push_back( *rit );
            result = mbImpl->tag_set_data( mark_tag, &( *rit ), 1, &mark_val );
            if( MB_SUCCESS != result ) return result;
        }

        // end of inner loop; get the next this_ent, or set to zero
        if( hp_untreated.empty() )
            this_ent = 0;
        else
        {
            this_ent = hp_untreated.back();
            hp_untreated.pop_back();
        }
    }

    if( debug_ap )
    {
        std::string hp_name;
        if( 2 == dim )
            hp_name = "sheet";
        else
            hp_name = "chord";

        if( 0 == this_hp )
            std::cout << "Constructed new " << hp_name << " with ";
        else
        {
            int this_id;
            result = mbImpl->tag_get_data( globalId_tag(), &this_hp, 1, &this_id );RR;
            std::cout << "Added to " << hp_name << " " << this_id << " ";
        }
        if( dim == 2 )
            std::cout << "edges:" << std::endl;
        else
            std::cout << "quads:" << std::endl;
        std::vector< EntityHandle > pents( new_hyperplane_ents.size() );
        result = mbImpl->tag_get_data( dualEntity_tag(), new_hyperplane_ents, &pents[0] );RR;
        for( std::vector< EntityHandle >::iterator vit = pents.begin(); vit != pents.end(); ++vit )
        {
            if( vit != pents.begin() ) std::cout << ", ";
            std::cout << mbImpl->id_from_handle( *vit );
        }
        std::cout << std::endl;
    }

    if( 0 == this_hp )
    {
        // ok, doesn't have one; make a new hyperplane
        int new_id = -1;
        result     = construct_new_hyperplane( dim, this_hp, new_id );
        if( MB_SUCCESS != result ) return result;

        if( debug_ap )
        {
            std::cout << "New ";
            if( 2 == dim )
                std::cout << " sheet ";
            else
                std::cout << " chord ";
            std::cout << new_id << " constructed." << std::endl;
        }
    }

    // set the hp_val for entities which didn't have one before
    std::vector< EntityHandle > hp_tags( new_hyperplane_ents.size() );
    std::fill( hp_tags.begin(), hp_tags.end(), this_hp );
    result = mbImpl->tag_set_data( hp_tag, new_hyperplane_ents, &hp_tags[0] );
    if( MB_SUCCESS != result ) return result;
    result = mbImpl->add_entities( this_hp, new_hyperplane_ents );
    if( MB_SUCCESS != result ) return result;

    // unmark the entities by removing the tag
    result = mbImpl->tag_delete( mark_tag );
    if( MB_SUCCESS != result ) return result;

    return MB_SUCCESS;
}

ErrorCode DualTool::order_chord( EntityHandle chord_set )
{
    // re-order the 1cells in the set so they are in order along the chord
    // start by finding the vertex dual to a quad
    Range verts, one_cells;
    ErrorCode result = mbImpl->get_entities_by_dimension( chord_set, 1, one_cells );
    if( MB_SUCCESS != result || one_cells.empty() ) return MB_FAILURE;

    result = mbImpl->get_adjacencies( one_cells, 0, false, verts, Interface::UNION );
    if( MB_SUCCESS != result || verts.empty() ) return MB_FAILURE;

    EntityHandle last_vert = 0;
    for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
    {
        if( TYPE_FROM_HANDLE( get_dual_entity( *rit ) ) == MBQUAD )
        {
            last_vert = *rit;
            break;
        }
    }
    // if there's no vertex owned by a quad, just start with 1st one
    if( 0 == last_vert ) last_vert = *verts.begin();

    // now, skip from vertex to vertex, building a list of 1cells
    std::vector< EntityHandle > ordered_1cells;
    EntityHandle last_1cell = 0;
    Range dum1, dum2;
    const EntityHandle* connect;
    int num_connect;
    ErrorCode tmp_result = MB_SUCCESS;
    while( ordered_1cells.size() != one_cells.size() )
    {
        dum1   = one_cells;
        result = mbImpl->get_adjacencies( &last_vert, 1, 1, false, dum1 );
        if( 0 != last_1cell ) dum1.erase( last_1cell );
        // assert(1 == dum1.size());
        if( 0 != last_1cell && 1 != dum1.size() )
        {
            std::cerr << "unexpected size traversing chord." << std::endl;
            tmp_result = MB_FAILURE;
        }

        last_1cell = *dum1.begin();
        ordered_1cells.push_back( last_1cell );
        result = mbImpl->get_connectivity( last_1cell, connect, num_connect );RR;
        if( last_vert == connect[0] )
            last_vert = connect[1];
        else
            last_vert = connect[0];
    }

    // now have the 1cells in order, replace them in the set
    if( MB_SUCCESS == tmp_result )
    {
        result = mbImpl->remove_entities( chord_set, one_cells );RR;
        result = mbImpl->add_entities( chord_set, &ordered_1cells[0], ordered_1cells.size() );RR;
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::construct_new_hyperplane( const int dim, EntityHandle& new_hyperplane, int& id )
{
    ErrorCode result;
    if( 1 == dim )
        result = mbImpl->create_meshset( ( MESHSET_ORDERED | MESHSET_TRACK_OWNER ), new_hyperplane );
    else
        result = mbImpl->create_meshset( ( MESHSET_SET | MESHSET_TRACK_OWNER ), new_hyperplane );
    if( MB_SUCCESS != result ) return result;

    if( -1 == id )
    {
        Range all_hyperplanes;
        result = get_dual_hyperplanes( mbImpl, dim, all_hyperplanes );RR;
        std::vector< int > gids( all_hyperplanes.size() );
        result = mbImpl->tag_get_data( globalIdTag, all_hyperplanes, ( gids.empty() ) ? NULL : &gids[0] );RR;
        for( unsigned int i = 0; i < gids.size(); i++ )
            if( gids[i] > id ) id = gids[i];
        id++;
        if( 0 == id ) id++;
    }

    result = mbImpl->tag_set_data( globalId_tag(), &new_hyperplane, 1, &id );RR;
    Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
    result     = mbImpl->tag_set_data( hp_tag, &new_hyperplane, 1, &new_hyperplane );RR;

    // assign a category name to these sets
    static const char dual_category_names[2][CATEGORY_TAG_SIZE] = { "Chord\0", "Sheet\0" };

    result = mbImpl->tag_set_data( categoryTag, &new_hyperplane, 1, dual_category_names[dim - 1] );

    return result;
}

bool DualTool::check_1d_loop_edge( EntityHandle this_ent )
{
    // make sure it's an edge
    if( MBEDGE != mbImpl->type_from_handle( this_ent ) ) return false;

    // also has to be a dual entity
    unsigned int dum;
    ErrorCode result = mbImpl->tag_get_data( isDualCell_tag(), &this_ent, 1, &dum );
    if( MB_SUCCESS != result || dum != 0x1 ) return false;

    const EntityHandle* verts;
    EntityHandle vert_tags[2];
    int num_verts;
    result = mbImpl->get_connectivity( this_ent, verts, num_verts );
    if( MB_SUCCESS != result ) return false;

    result = mbImpl->tag_get_data( dualEntity_tag(), verts, 2, vert_tags );
    if( MB_SUCCESS != result || mbImpl->type_from_handle( vert_tags[0] ) != MBQUAD ||
        mbImpl->type_from_handle( vert_tags[1] ) != MBQUAD )
        return false;

    else
        return true;
}

ErrorCode DualTool::construct_hp_parent_child()
{
    Range dual_surfs, dual_cells, dual_edges;
    ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );
    if( MB_SUCCESS != result || dual_surfs.empty() ) return result;
    std::vector< EntityHandle > dual_curve_sets;

    for( Range::iterator surf_it = dual_surfs.begin(); surf_it != dual_surfs.end(); ++surf_it )
    {
        // get all the cells, edges in those cells, and chords for those edges
        dual_cells.clear();
        result = mbImpl->get_entities_by_handle( *surf_it, dual_cells );
        if( MB_SUCCESS != result ) return result;
        dual_edges.clear();
        result = mbImpl->get_adjacencies( dual_cells, 1, false, dual_edges, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
        dual_curve_sets.resize( dual_edges.size() );
        result = mbImpl->tag_get_data( dualCurve_tag(), dual_edges, &dual_curve_sets[0] );
        if( MB_SUCCESS != result ) return result;

        // reuse dual_cells to get unique list of chord sets
        dual_cells.clear();
        for( unsigned int i = 0; i < dual_edges.size(); i++ )
            if( dual_curve_sets[i] != 0 ) dual_cells.insert( dual_curve_sets[i] );

        // now connect up this dual surf with all the 1d ones
        for( Range::iterator rit = dual_cells.begin(); rit != dual_cells.end(); ++rit )
        {
            result = mbImpl->add_parent_child( *surf_it, *rit );
            if( MB_SUCCESS != result ) return result;
        }
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::get_graphics_points( EntityHandle dual_ent,
                                         std::vector< int >& npts,
                                         std::vector< GraphicsPoint >& points )
{
    // shouldn't be a set
    assert( MBENTITYSET != mbImpl->type_from_handle( dual_ent ) );

    // get the graphics points comprising the given entity
    GraphicsPoint gp_array[DualTool::GP_SIZE];

    ErrorCode result = MB_SUCCESS;

    // switch based on topological dimension
    switch( mbImpl->dimension_from_handle( dual_ent ) )
    {
        case 0:
            // just return the vertex point
            result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
            if( MB_SUCCESS == result ) points.push_back( gp_array[0] );

            break;

        case 1:
            // get my graphics point then those of my vertices
            const EntityHandle* connect;
            int num_connect;
            result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
            if( MB_SUCCESS != result ) break;

            result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, gp_array );
            if( MB_SUCCESS == result )
            {
                points.push_back( gp_array[0] );
                points.push_back( gp_array[0] );
                points.push_back( gp_array[1] );
                result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
                if( MB_SUCCESS == result ) points[1] = gp_array[0];
            }

            npts.push_back( 3 );

            break;

        case 2:
            result = get_cell_points( dual_ent, npts, points );
            break;
    }

    return result;
}

ErrorCode DualTool::get_cell_points( EntityHandle dual_ent,
                                     std::vector< int >& npts,
                                     std::vector< GraphicsPoint >& points )
{
    assert( MBPOLYGON == mbImpl->type_from_handle( dual_ent ) );

    // get the 1cells in this 2cell
    Range one_cells;

    Range tc_range;
    tc_range.insert( dual_ent );
    ErrorCode result = mbImpl->get_adjacencies( tc_range, 1, false, one_cells, Interface::UNION );RR;

    int num_edges = one_cells.size();
    std::vector< GraphicsPoint > dum_gps( num_edges + 1 );

    // get graphics points for 0cells and for this cell
    result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), one_cells, &dum_gps[0] );RR;
    result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, &( dum_gps[num_edges] ) );RR;

    Range::iterator eit;
    const EntityHandle* connect;
    int num_connect;
    GraphicsPoint vert_gps[2];
    int i;
    for( i = 0, eit = one_cells.begin(); i < num_edges; i++, ++eit )
    {
        // get the vertices and the graphics points for them
        result = mbImpl->get_connectivity( *eit, connect, num_connect );RR;
        result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, vert_gps );RR;

        // make the 2 tris corresponding to this edge; don't worry about order
        // for now
        npts.push_back( 3 );
        points.push_back( dum_gps[num_edges] );
        points.push_back( vert_gps[0] );
        points.push_back( dum_gps[i] );

        npts.push_back( 3 );
        points.push_back( dum_gps[num_edges] );
        points.push_back( dum_gps[i] );
        points.push_back( vert_gps[1] );
    }

    return result;
}

ErrorCode DualTool::get_graphics_points( const Range& in_range,
                                         std::vector< GraphicsPoint >& points,
                                         const bool assign_ids,
                                         const int start_id )
{
    // return graphics points on dual entities in in_range or in entities
    // in sets in in_range
    ErrorCode result;

    // for each dual hyperplane set:
    Range::const_iterator rit;

    Range two_cells, all_cells;
    for( rit = in_range.begin(); rit != in_range.end(); ++rit )
    {
        // for each entity:
        two_cells.clear();
        EntityType this_type = mbImpl->type_from_handle( *rit );
        if( MBENTITYSET == this_type )
        {
            result = mbImpl->get_entities_by_handle( *rit, two_cells );RR;

            std::copy( two_cells.begin(), two_cells.end(), range_inserter( all_cells ) );
        }

        else
        {
            two_cells.insert( *rit );
            assert( this_type == MBVERTEX || this_type == MBEDGE || this_type == MBPOLYGON ||
                    this_type == MBPOLYHEDRON );
        }

        result = mbImpl->get_adjacencies( two_cells, 0, false, all_cells, Interface::UNION );RR;
        result = mbImpl->get_adjacencies( two_cells, 1, false, all_cells, Interface::UNION );RR;
    }

    // get graphics points
    points.resize( all_cells.size() );

    result = mbImpl->tag_get_data( dualGraphicsPointTag, all_cells, &points[0] );RR;

    if( assign_ids )
    {
        int i = start_id;

        for( std::vector< GraphicsPoint >::iterator vit = points.begin(); vit != points.end(); ++vit )
            vit->id = i++;

        result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), all_cells, &points[0] );RR;
    }

    return result;
}

EntityHandle DualTool::next_loop_vertex( const EntityHandle last_v,
                                         const EntityHandle this_v,
                                         const EntityHandle dual_surf )
{
    // given two vertices, find the next one on the loop; if one is a dual
    // surface, then just choose either one for that surface
    assert( ( 0 == last_v || mbImpl->type_from_handle( last_v ) == MBVERTEX ) &&
            mbImpl->type_from_handle( this_v ) == MBVERTEX && mbImpl->type_from_handle( dual_surf ) == MBENTITYSET );

    // get the connected vertices
    MeshTopoUtil tpu( mbImpl );
    Range other_verts;
    ErrorCode result = tpu.get_bridge_adjacencies( this_v, 1, 0, other_verts );
    if( MB_SUCCESS != result || other_verts.empty() ) return 0;

    // if (mbImpl->type_from_handle(last_v) == MBENTITYSET) {
    // dual surface, choose either; first get a 2cell on this surface
    Range tcells, tcells2, verts;
    result = mbImpl->get_entities_by_type( dual_surf, MBPOLYGON, tcells );
    if( MB_SUCCESS != result || tcells.empty() ) return 0;

    // ok, pay attention here: first get 2cells common to dual surface and this_v
    verts.insert( this_v );
    result = mbImpl->get_adjacencies( verts, 2, false, tcells );
    if( MB_SUCCESS != result || tcells.empty() ) return 0;

    // next get vertices common to both 2cells and subtract from other_verts; also
    // remove last_v if it's non-zero
    verts.clear();
    result = mbImpl->get_adjacencies( tcells, 0, false, verts );
    if( MB_SUCCESS != result || verts.empty() ) return 0;

    Range tmp_verts = subtract( other_verts, verts );
    other_verts.swap( tmp_verts );
    if( 0 != last_v ) other_verts.erase( last_v );

    // now get intersection of remaining vertices and 2 2cells vertices
    // look at each one successively; should find one, maybe not on both
    tmp_verts = other_verts;
    Range tmp_faces( *tcells.begin(), *tcells.begin() );
    result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );
    if( MB_SUCCESS == result && !tmp_verts.empty() ) return *tmp_verts.begin();
    tmp_faces.clear();
    tmp_faces.insert( *tcells.rbegin() );
    result = mbImpl->get_adjacencies( tmp_faces, 0, false, other_verts );
    if( MB_SUCCESS == result && !other_verts.empty() ) return *other_verts.begin();

    // if we got here, there isn't any
    return 0;
}

EntityHandle DualTool::get_dual_hyperplane( const EntityHandle ncell )
{
    // get the sheet or chord it's in
    std::vector< EntityHandle > adj_sets;
    ErrorCode result = mbImpl->get_adjacencies( &ncell, 1, 4, false, adj_sets );
    if( MB_SUCCESS != result ) return 0;

    EntityHandle dum_set;
    for( std::vector< EntityHandle >::iterator vit = adj_sets.begin(); vit != adj_sets.end(); ++vit )
    {
        if( mbImpl->tag_get_data( dualCurve_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND ||
            mbImpl->tag_get_data( dualSurface_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND )
            return *vit;
    }

    return 0;
}

//! set the dual surface or curve for an entity
ErrorCode DualTool::set_dual_surface_or_curve( EntityHandle entity,
                                               const EntityHandle dual_hyperplane,
                                               const int dual_entity_dimension )
{
    if( 1 == dual_entity_dimension )
        mbImpl->tag_set_data( dualCurve_tag(), &entity, 1, &dual_hyperplane );
    else if( 2 == dual_entity_dimension )
        mbImpl->tag_set_data( dualSurface_tag(), &entity, 1, &dual_hyperplane );
    else
        return MB_INDEX_OUT_OF_RANGE;

    return MB_SUCCESS;
}

//! return the corresponding dual entity
EntityHandle DualTool::get_dual_entity( const EntityHandle this_ent ) const
{
    EntityHandle dual_ent;
    ErrorCode result = mbImpl->tag_get_data( dualEntity_tag(), &this_ent, 1, &dual_ent );
    if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
        return 0;
    else
        return dual_ent;
}

//! return the corresponding dual entity
EntityHandle DualTool::get_extra_dual_entity( const EntityHandle this_ent )
{
    EntityHandle dual_ent;
    ErrorCode result = mbImpl->tag_get_data( extraDualEntity_tag(), &this_ent, 1, &dual_ent );
    if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
        return 0;
    else
        return dual_ent;
}

ErrorCode DualTool::atomic_pillow( EntityHandle odedge, EntityHandle& quad1, EntityHandle& quad2 )
{
    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    if( debug_ap )
    {
        Range sets;
        Tag ms_tag;

        ErrorCode result = mbImpl->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, ms_tag );
        if( MB_SUCCESS == result )
        {
            result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &ms_tag, NULL, 1, sets );
            if( MB_SUCCESS == result ) result = mbImpl->delete_entities( sets );
        }
    }

    std::cout << "-AP(";
    print_cell( odedge );
    std::cout << ")" << std::endl;

    // perform an atomic pillow operation around dedge

    // grab the quad before deleting the odedge
    quad1 = get_dual_entity( odedge );
    assert( 0 != quad1 );

    // 0. get star 2cells around odedge (before odedge changes) and 3cells around
    // those 2cells (going to delete the 2cells, therefore need to delete the 3cells
    // that depend on those too)
    MeshTopoUtil mtu( mbImpl );
    Range star_cells, tmp_cells;
    ErrorCode result = mbImpl->get_adjacencies( &odedge, 1, 2, false, star_cells );RR;
    result = mbImpl->get_adjacencies( star_cells, 3, false, tmp_cells, Interface::UNION );RR;
    star_cells.merge( tmp_cells );
    star_cells.insert( odedge );

    // tear down the dual entities which will be modified by the ap first
    result = delete_dual_entities( star_cells );RR;

    // now change the quad to an ap
    std::vector< EntityHandle > verts;
    result = mbImpl->get_connectivity( &quad1, 1, verts );RR;

    // get average position of vertices
    double coords[12], avg[3] = { 0.0, 0.0, 0.0 };
    result = mbImpl->get_coords( &verts[0], verts.size(), coords );RR;
    for( int i = 0; i < 4; i++ )
    {
        avg[0] += coords[3 * i];
        avg[1] += coords[3 * i + 1];
        avg[2] += coords[3 * i + 2];
    }
    for( int i = 0; i < 3; i++ )
        avg[i] *= 0.25;

    // for each position, get a corresponding position 1/2 way to avg
    double new_coords[12];
    for( int i = 0; i < 4; i++ )
    {
        new_coords[3 * i]     = avg[0] + .5 * ( coords[3 * i] - avg[0] );
        new_coords[3 * i + 1] = avg[1] + .5 * ( coords[3 * i + 1] - avg[1] );
        new_coords[3 * i + 2] = avg[2] + .5 * ( coords[3 * i + 2] - avg[2] );
    }

    // make the 4 new vertices; store in vector long enough for hex connectivity
    for( int i = 0; i < 4; i++ )
    {
        verts.push_back( 0 );
        result = mbImpl->create_vertex( &new_coords[3 * i], verts[4 + i] );RR;
    }

    // get the hexes connected to the quad
    Range hexes;
    result = mbImpl->get_adjacencies( &quad1, 1, 3, false, hexes );RR;
    assert( hexes.size() <= 2 );

    // remove any explicit adjacency from the first hex, since that'll get connected
    // to the new quad; add adjacency between quad and second hex, if there is a 2nd
    result = mbImpl->remove_adjacencies( quad1, &( *hexes.begin() ), 1 );RR;
    if( hexes.size() == 2 )
    {
        result = mbImpl->add_adjacencies( quad1, &( *hexes.rbegin() ), 1, false );RR;
    }

    // create the new, inner quad, and make it explicitly adjacent to 1st hex;
    // make the connectivity of this quad same as the original one
    std::vector< EntityHandle > tmp_verts;
    std::copy( verts.begin(), verts.end(), std::back_inserter( tmp_verts ) );

    result = mbImpl->create_element( MBQUAD, &tmp_verts[0], 4, quad2 );RR;
    result = mbImpl->add_adjacencies( quad2, &( *hexes.begin() ), 1, false );RR;

    // reverse the connectivity of the 1st hex
    std::reverse( verts.begin(), verts.begin() + 4 );
    std::reverse( verts.begin() + 4, verts.end() );

    // now make two inner hexes; note connectivity array is flipped for the two hexes
    EntityHandle new_hexes[2];
    result = mbImpl->create_element( MBHEX, &verts[0], 8, new_hexes[0] );RR;
    result = mbImpl->create_element( MBHEX, &tmp_verts[0], 8, new_hexes[1] );RR;

    // set the global id tag on the new hexes
    int new_hex_ids[2] = { maxHexId + 1, maxHexId + 2 };
    maxHexId += 2;
    result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 2, new_hex_ids );
    if( MB_SUCCESS != result ) return result;

    // by definition, quad1 is adj to new_hexes[0]
    result = mbImpl->add_adjacencies( quad1, &new_hexes[0], 1, false );RR;
    result = mbImpl->add_adjacencies( quad2, &new_hexes[1], 1, false );RR;

    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    // now update the dual
    result = construct_hex_dual( &new_hexes[0], 2 );RR;

    // get the new dual surface, by getting one of the edges between the center
    // and outer vertex rings
    Range new_edge;
    verts[1] = verts[4];
    result   = mbImpl->get_adjacencies( &verts[0], 2, 1, false, new_edge );
    if( MB_SUCCESS != result || new_edge.size() != 1 ) return result;

    return MB_SUCCESS;
}

//! effect reverse atomic pillow operation
ErrorCode DualTool::rev_atomic_pillow( EntityHandle pillow, Range& chords )
{
    // get the dual entities associated with elements in the pillow; go through
    // the elements instead of the pillow sheet so you get all of them, not just
    // the ones on the sheet
    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    std::cout << "-AP(";
    print_cell( pillow );
    std::cout << ")" << std::endl;

    Range dverts;
    ErrorCode result = get_dual_entities( pillow, NULL, NULL, &dverts, NULL, NULL );
    if( MB_SUCCESS != result ) return result;
    assert( 2 == dverts.size() );

    EntityHandle hexes[2];
    result = mbImpl->tag_get_data( dualEntity_tag(), dverts, hexes );RR;
    assert( hexes[0] != 0 && hexes[1] != 0 );

    std::vector< EntityHandle > dcells[4];
    Range pcells[4];
    std::copy( hexes, hexes + 2, range_inserter( pcells[3] ) );
    std::copy( dverts.begin(), dverts.end(), std::back_inserter( dcells[0] ) );
    for( int dim = 0; dim <= 2; dim++ )
    {
        result = mbImpl->get_adjacencies( hexes, 2, dim, false, pcells[dim], Interface::UNION );RR;
        dcells[3 - dim].resize( pcells[dim].size() );
        result = mbImpl->tag_get_data( dualEntity_tag(), pcells[dim], &dcells[3 - dim][0] );RR;
    }

    // delete the dual entities which are part of the original pillow
    result = mbImpl->delete_entities( &pillow, 1 );
    if( MB_SUCCESS != result ) return result;

    result = mbImpl->delete_entities( chords );
    if( MB_SUCCESS != result ) return result;

    for( int i = 3; i >= 0; i-- )
    {
        result = delete_dual_entities( &dcells[i][0], dcells[i].size() );RR;
    }

    // delete the primal entities inserted by the ap; be careful to get the right
    // faces, edges and vertices
    Range del_faces, del_edges, del_verts, tmp_faces, tmp_verts;
    // faces are the shared 5 and the 1 other one with greater handle (which
    // we know will be later in the range)
    result = mbImpl->get_adjacencies( hexes, 2, 2, false, del_faces );RR;
    assert( 5 == del_faces.size() );
    std::copy( pcells[2].begin(), pcells[2].end(), range_inserter( tmp_faces ) );
    tmp_faces = subtract( tmp_faces, del_faces );
    del_faces.insert( *tmp_faces.rbegin() );
    result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );RR;
    std::copy( pcells[0].begin(), pcells[0].end(), range_inserter( del_verts ) );
    del_verts = subtract( del_verts, tmp_verts );
    assert( 4 == del_verts.size() );
    result = mbImpl->get_adjacencies( del_verts, 1, false, del_edges, Interface::UNION );RR;
    assert( 8 == del_edges.size() );

    result = mbImpl->delete_entities( hexes, 2 );RR;
    result = mbImpl->delete_entities( del_faces );RR;
    result = mbImpl->delete_entities( del_edges );RR;
    result = mbImpl->delete_entities( del_verts );RR;

    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    // recompute the dual for the hexes on either side of the quad affected
    // by the ap removal
    Range tmp_hexes;
    result = mbImpl->get_adjacencies( tmp_verts, 3, false, tmp_hexes, Interface::UNION );RR;
    result = construct_hex_dual( tmp_hexes );RR;

    return MB_SUCCESS;
}

ErrorCode DualTool::delete_dual_entities( EntityHandle* entities, int num_entities )
{
    Range tmp_ents;
    std::copy( entities, entities + num_entities, range_inserter( tmp_ents ) );
    return delete_dual_entities( tmp_ents );
}

ErrorCode DualTool::delete_dual_entities( Range& entities )
{
    if( entities.empty() ) return delete_whole_dual();

    EntityHandle null_entity = 0;
    ErrorCode result;
    Range ents_to_delete;

    while( !entities.empty() )
    {
        EntityHandle this_entity = entities.pop_back();

        // reset the primal's dual entity
        EntityHandle primal = get_dual_entity( this_entity );
        if( get_dual_entity( primal ) == this_entity )
        {
            result = mbImpl->tag_set_data( dualEntity_tag(), &primal, 1, &null_entity );RR;
        }
        EntityHandle extra = get_extra_dual_entity( primal );
        if( 0 != extra )
        {
            result = mbImpl->tag_set_data( extraDualEntity_tag(), &primal, 1, &null_entity );RR;
        }

        ents_to_delete.insert( this_entity );

        // check for extra dual entities
        if( mbImpl->type_from_handle( this_entity ) == MBPOLYGON )
        {
            // for 2cell, might be a loop edge
            Range loop_edges;
            result = mbImpl->get_adjacencies( &this_entity, 1, 1, false, loop_edges );
            for( Range::iterator rit = loop_edges.begin(); rit != loop_edges.end(); ++rit )
                if( check_1d_loop_edge( *rit ) ) entities.insert( *rit );
        }
        else if( extra && extra != this_entity )
            // just put it on the list; primal for which we're extra has already been
            // reset to not point to extra entity
            ents_to_delete.insert( extra );
    }

    // now delete the entities (sheets and chords will be updated automatically)
    return mbImpl->delete_entities( ents_to_delete );
}

void DualTool::print_cell( EntityHandle cell )
{
    const EntityHandle* connect;
    int num_connect;
    ErrorCode result = mbImpl->get_connectivity( cell, connect, num_connect );
    if( MB_SUCCESS != result ) return;
    bool first = true;
    EntityHandle primals[20];
    std::vector< int > ids;

    assert( num_connect < 20 );
    result = mbImpl->tag_get_data( dualEntityTag, connect, num_connect, primals );
    if( MB_SUCCESS != result ) return;
    ids.resize( num_connect );
    result = mbImpl->tag_get_data( globalIdTag, primals, num_connect, &ids[0] );
    if( MB_SUCCESS != result ) return;
    for( int i = 0; i < num_connect; i++ )
    {
        if( !first ) std::cout << "-";
        EntityType this_type = mbImpl->type_from_handle( primals[i] );
        if( this_type == MBHEX )
            std::cout << "h";
        else if( this_type == MBQUAD )
            std::cout << "f";
        else
            std::cout << "u";

        if( ids[i] != 0 )
            std::cout << ids[i];
        else
            std::cout << mbImpl->id_from_handle( primals[i] );

        first = false;
    }
}

ErrorCode DualTool::face_open_collapse( EntityHandle ocl, EntityHandle ocr )
{
    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    MeshTopoUtil mtu( mbImpl );

    std::cout << "OC(";
    print_cell( ocl );
    std::cout << ")-(";
    print_cell( ocr );
    std::cout << ")" << std::endl;

    // get the primal entities we're dealing with
    EntityHandle split_quads[2] = { 0 }, split_edges[3] = { 0 }, split_nodes[2] = { 0 }, other_edges[6] = { 0 },
                 other_nodes[6] = { 0 };
    Range hexes;
    ErrorCode result = foc_get_ents( ocl, ocr, split_quads, split_edges, split_nodes, hexes, other_edges, other_nodes );RR;

    // get star entities around edges, separated into halves
    std::vector< EntityHandle > star_dp1[2], star_dp2[2];
    result = foc_get_stars( split_quads, split_edges, star_dp1, star_dp2 );RR;

    if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
        return MB_TYPE_OUT_OF_RANGE;

    result = foc_delete_dual( split_quads, split_edges, hexes );
    if( MB_SUCCESS != result ) return result;

    EntityHandle new_quads[2], new_edges[3], new_nodes[2];
    result = split_pair_nonmanifold( split_quads, split_edges, split_nodes, star_dp1, star_dp2, other_edges,
                                     other_nodes, new_quads, new_edges, new_nodes );
    if( MB_SUCCESS != result ) return result;

    // now merge entities, the C of foc
    EntityHandle keepit, deleteit;
#define MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
#define MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )
#define KEEP_DELETE( a, b, c, d ) \
    {                             \
        ( c ) = MIN( a, b );      \
        ( d ) = MAX( a, b );      \
    }

    // find how many shared edges there were
    int num_shared_edges = ( split_edges[2] ? 3 : ( split_edges[1] ? 2 : 1 ) );

    // first the node(s)
    for( int i = 0; i < 3 - num_shared_edges; i++ )
    {
        KEEP_DELETE( other_nodes[2 + 2 * i], other_nodes[3 + 2 * i], keepit, deleteit );
        result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
    }

    // now the edges
    for( int i = 0; i < 4 - num_shared_edges; i++ )
    {
        KEEP_DELETE( other_edges[2 * i], other_edges[2 * i + 1], keepit, deleteit );
        result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
    }

    // now the faces
    KEEP_DELETE( split_quads[0], split_quads[1], keepit, deleteit );
    result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;

    result = mbImpl->merge_entities( new_quads[0], new_quads[1], false, true );RR;

    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    // reconstruct dual
    result = construct_hex_dual( hexes );
    if( MB_SUCCESS != result ) return result;

    return check_dual_adjs();
}

ErrorCode DualTool::foc_get_ents( EntityHandle ocl,
                                  EntityHandle ocr,
                                  EntityHandle* split_quads,
                                  EntityHandle* split_edges,
                                  EntityHandle* split_nodes,
                                  Range& hexes,
                                  EntityHandle* other_edges,
                                  EntityHandle* other_nodes )
{
    // get the entities used for foc; ocl and ocr are dual 1-cells
    // representing quads to be split; returned from this function:
    // quads[2] - 2 quads to be split
    // split_edges[2] - edge(s) to be split (2nd is 0 if only one)
    // split_node - node to be split, if any (otherwise 0)
    // hexes - connected hexes to split_edges
    // other_edges[0], [1] - edges in quads[0] and [1] sharing node with
    //        one end of split_edges[0]
    // other_edges[2], [3] - other end of split_edges[0] (or [1] if 2
    //        split_edges)
    // other_edges[4], [5] - edges in quads[0], [1] opposite to split_edges[0]
    // other_nodes[0], [1] - nodes on other_edges[0], [1] not shared with
    //        split_edges[0]
    // other_nodes[2], [3] - nodes on other_edges[2], [3] not shared with
    //        split_edges[0] (if 2 split_edges, there's only 1 opposite node
    //        in each split quad)
    // (for diagram, see Tim's notes from 11/12/07)

    split_quads[0] = get_dual_entity( ocl );
    split_quads[1] = get_dual_entity( ocr );
    if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
        return MB_TYPE_OUT_OF_RANGE;

    Range common_edges;
    ErrorCode result = mbImpl->get_adjacencies( split_quads, 2, 1, false, common_edges );
    if( MB_SUCCESS != result ) return result;

    if( common_edges.empty() ) return MB_FAILURE;
    for( unsigned int i = 0; i < common_edges.size(); i++ )
        split_edges[i] = common_edges[i];

    MeshTopoUtil mtu( mbImpl );

    if( common_edges.size() == 3 )
    {
        // find other (non-shared) edges
        for( int i = 0; i < 2; i++ )
        {
            Range tmp_edges;
            result = mbImpl->get_adjacencies( &split_quads[i], 1, 1, false, tmp_edges );
            if( MB_SUCCESS != result ) return result;
            tmp_edges = subtract( tmp_edges, common_edges );
            assert( tmp_edges.size() == 1 );
            other_edges[i] = *tmp_edges.begin();
        }
        assert( other_edges[0] && other_edges[1] && other_edges[0] != other_edges[1] );

        // arrange common edges so middle is in middle
        result = mtu.opposite_entity( split_quads[0], other_edges[0], split_edges[1] );RR;
        common_edges.erase( split_edges[1] );
        split_edges[0] = *common_edges.begin();
        split_edges[2] = *common_edges.rbegin();
        common_edges.insert( split_edges[1] );

        // get split nodes and other nodes
        split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
        split_nodes[1] = mtu.common_entity( split_edges[2], split_edges[1], 0 );
        other_nodes[0] = mtu.common_entity( split_edges[0], other_edges[0], 0 );
        other_nodes[1] = mtu.common_entity( split_edges[2], other_edges[1], 0 );

        assert( other_nodes[0] && other_nodes[1] && split_nodes[0] && split_nodes[1] );
        assert( split_edges[0] && split_edges[1] && split_edges[2] && split_edges[0] != split_edges[1] &&
                split_edges[1] != split_edges[2] && split_edges[0] != split_edges[2] );
    }
    else if( common_edges.size() == 2 )
    {
        // split node is shared by split edges
        split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
        if( 0 == split_nodes[0] ) return MB_FAILURE;
        // first two other nodes are on split edges opposite split node
        result = mtu.opposite_entity( split_edges[0], split_nodes[0], other_nodes[0] );RR;
        result = mtu.opposite_entity( split_edges[1], split_nodes[0], other_nodes[1] );RR;
        // over split quads:
        for( int i = 0; i < 2; i++ )
        {
            // 1st other edge is opposite second split edge
            result = mtu.opposite_entity( split_quads[i], split_edges[1], other_edges[i] );RR;
            // 2nd other edge is opposite first split edge
            result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
            // last other node is opposite split node on split quad
            result = mtu.opposite_entity( split_quads[i], split_nodes[0], other_nodes[2 + i] );RR;
        }
    }
    else
    {
        const EntityHandle* connect;
        int num_connect;
        result = mbImpl->get_connectivity( split_edges[0], connect, num_connect );
        if( MB_SUCCESS != result ) return result;
        // other_nodes[0], [1] are on split edge
        other_nodes[0] = connect[0];
        other_nodes[1] = connect[1];

        // for each of the split quads
        for( int i = 0; i < 2; i++ )
        {
            // get the other edge on the split quad adj to node 0 on the split edge, by getting
            // edges adj to split quad and node and removing split edge; that's other_edge[i]
            Range tmp_range1, tmp_range2;
            tmp_range1.insert( connect[0] );
            tmp_range1.insert( split_quads[i] );
            result = mbImpl->get_adjacencies( tmp_range1, 1, false, tmp_range2 );
            if( MB_SUCCESS != result ) return result;
            tmp_range2.erase( split_edges[0] );
            assert( tmp_range2.size() == 1 );
            other_edges[i] = *tmp_range2.begin();
            // get edge connected to other node on split edge & split quad; that's
            // opposite prev other_edges on the split quad; that's other_edges[4+i]
            result = mtu.opposite_entity( split_quads[i], other_edges[i], other_edges[4 + i] );RR;
            // get the edge on the split quad opposite the split edge; that's other_edges[2+i]
            result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
            // get nodes on other side of split quad from split edge, by getting common
            // node between top/bottom edge and opposite edge
            other_nodes[2 + i] = mtu.common_entity( other_edges[i], other_edges[2 + i], 0 );
            other_nodes[4 + i] = mtu.common_entity( other_edges[4 + i], other_edges[2 + i], 0 );
            if( 0 == other_nodes[2 + i] || 0 == other_nodes[4 + i] ) return MB_FAILURE;
        }
    }

    result = mbImpl->get_adjacencies( split_edges, common_edges.size(), 3, false, hexes, Interface::UNION );
    if( MB_SUCCESS != result ) return result;

    assert( "split node not in other_nodes" && other_nodes[0] != split_nodes[0] && other_nodes[0] != split_nodes[1] &&
            other_nodes[1] != split_nodes[0] && other_nodes[1] != split_nodes[1] );
    assert( "each split node on an end of a split edge" && mtu.common_entity( other_nodes[0], split_edges[0], 0 ) &&
            ( ( ( split_edges[2] && mtu.common_entity( other_nodes[1], split_edges[2], 0 ) ) ||
                ( split_edges[1] && mtu.common_entity( other_nodes[1], split_edges[1], 0 ) ) ||
                mtu.common_entity( other_nodes[1], split_edges[0], 0 ) ) ) );
    assert( "opposite other edges meet at an other node" &&
            ( mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[0] ||
              ( split_edges[2] && mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[1] ) ) &&
            ( split_edges[2] ||
              ( split_edges[1] && mtu.common_entity( other_edges[2], other_edges[3], 0 ) == other_nodes[1] ) ||
              mtu.common_entity( other_edges[4], other_edges[5], 0 ) == other_nodes[1] ) );

    return MB_SUCCESS;
}

ErrorCode DualTool::split_pair_nonmanifold( EntityHandle* split_quads,
                                            EntityHandle* split_edges,
                                            EntityHandle* split_nodes,
                                            std::vector< EntityHandle >* star_dp1,
                                            std::vector< EntityHandle >* star_dp2,
                                            EntityHandle* /*other_edges*/,
                                            EntityHandle* /*other_nodes*/,
                                            EntityHandle* new_quads,
                                            EntityHandle* new_edges,
                                            EntityHandle* new_nodes )
{

    // if there's a bdy in the star around the shared edge(s), get the quads on that
    // bdy so we know which edges to merge after the split-nonmanifold
    MeshTopoUtil mtu( mbImpl );
    ErrorCode result;

    // get which star the split faces are in, and choose the other one
    int new_side = -1;
    if( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() )
        new_side = 1;
    else if( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() )
        new_side = 0;
    assert( -1 != new_side );
    if( -1 == new_side ) return MB_FAILURE;

    //=============== split faces

    for( int i = 0; i < 2; i++ )
    {
        // get a hex in star_dp2[new_side] that's adj to this split quad, to tell
        // mtu which one the new quad should go with; there should be at least one,
        // if we have any hexes connected to the split quad
        EntityHandle gowith_hex = 0;
        for( std::vector< EntityHandle >::iterator vit = star_dp2[new_side].begin(); vit != star_dp2[new_side].end();
             ++vit )
        {
            if( mtu.common_entity( *vit, split_quads[i], 2 ) )
            {
                gowith_hex = *vit;
                break;
            }
        }
        assert( 0 != gowith_hex );

        // split manifold each of the split_quads, and put the results on the merge list
        result =
            mtu.split_entities_manifold( split_quads + i, 1, new_quads + i, NULL, ( gowith_hex ? &gowith_hex : NULL ) );RR;
    }

    // make ranges of faces which need to be explicitly adj to old, new
    // edge; faces come from stars and new_quads (which weren't in the stars);
    // new_quads go with side j, which does not have split quads
    Range tmp_addl_faces[2], addl_faces[2];
    for( int i = 0; i < 2; i++ )
    {
        std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_faces[i] ) );
        tmp_addl_faces[new_side].insert( new_quads[i] );
    }
#ifndef NDEBUG
    bool cond1 = ( "split_quads on 1, new_quads on 0" &&
                   tmp_addl_faces[0].find( split_quads[0] ) == tmp_addl_faces[0].end() &&
                   tmp_addl_faces[0].find( split_quads[1] ) == tmp_addl_faces[0].end() &&
                   tmp_addl_faces[1].find( split_quads[0] ) != tmp_addl_faces[1].end() &&
                   tmp_addl_faces[1].find( split_quads[1] ) != tmp_addl_faces[1].end() &&
                   tmp_addl_faces[0].find( new_quads[0] ) != tmp_addl_faces[0].end() &&
                   tmp_addl_faces[0].find( new_quads[1] ) != tmp_addl_faces[0].end() &&
                   tmp_addl_faces[1].find( new_quads[0] ) == tmp_addl_faces[1].end() &&
                   tmp_addl_faces[1].find( new_quads[1] ) == tmp_addl_faces[1].end() ),
         cond2 = ( "split_quads on 0, new_quads on 1" &&
                   tmp_addl_faces[0].find( split_quads[0] ) != tmp_addl_faces[0].end() &&
                   tmp_addl_faces[0].find( split_quads[1] ) != tmp_addl_faces[0].end() &&
                   tmp_addl_faces[1].find( split_quads[0] ) == tmp_addl_faces[1].end() &&
                   tmp_addl_faces[1].find( split_quads[1] ) == tmp_addl_faces[1].end() &&
                   tmp_addl_faces[0].find( new_quads[0] ) == tmp_addl_faces[0].end() &&
                   tmp_addl_faces[0].find( new_quads[1] ) == tmp_addl_faces[0].end() &&
                   tmp_addl_faces[1].find( new_quads[0] ) != tmp_addl_faces[1].end() &&
                   tmp_addl_faces[1].find( new_quads[1] ) != tmp_addl_faces[1].end() );

    assert( cond1 || cond2 );
#endif

    //=============== split edge(s)
    for( int j = 0; j < 3; j++ )
    {
        if( !split_edges[j] ) break;

        // filter add'l faces to only those adj to split_edges[j]
        addl_faces[0] = tmp_addl_faces[0];
        addl_faces[1] = tmp_addl_faces[1];
        for( int i = 0; i < 2; i++ )
        {
            result = mbImpl->get_adjacencies( &split_edges[j], 1, 2, false, addl_faces[i] );RR;
        }

        // split edge
        result = mtu.split_entity_nonmanifold( split_edges[j], addl_faces[1 - new_side], addl_faces[new_side],
                                               new_edges[j] );RR;
    }

    //=============== split node(s)

    for( int j = 0; j < 2; j++ )
    {
        if( !split_nodes[j] ) break;

        // if we're splitting multiple edges, there might be other edges that have the split
        // node; also need to know which side they're on
        Range tmp_addl_edges[2];
        result = foc_get_addl_ents( star_dp1, star_dp2, split_edges, split_nodes[j], tmp_addl_edges );RR;

        // also, we need to know which of the split/new edges go
        // with the split/new node; new edges go with side 0, split with 1
        for( int i = 0; i < 3; i++ )
        {
            if( !split_edges[i] ) break;
            tmp_addl_edges[new_side].insert( new_edges[i] );
            tmp_addl_edges[1 - new_side].insert( split_edges[i] );
        }

        // same for star faces and hexes
        for( int i = 0; i < 2; i++ )
        {
            std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_edges[i] ) );
            std::copy( star_dp2[i].begin(), star_dp2[i].end(), range_inserter( tmp_addl_edges[i] ) );
        }

        // finally, new quads
        for( int i = 0; i < 2; i++ )
            tmp_addl_edges[new_side].insert( new_quads[i] );

        // filter the entities, keeping only the ones adjacent to this node
        Range addl_edges[2];
        for( int i = 0; i < 2; i++ )
        {
            for( Range::reverse_iterator rit = tmp_addl_edges[i].rbegin(); rit != tmp_addl_edges[i].rend(); ++rit )
            {
                if( mtu.common_entity( *rit, split_nodes[j], 0 ) ) addl_edges[i].insert( *rit );
            }
        }

        // now split the node too
        result = mtu.split_entity_nonmanifold( split_nodes[j], addl_edges[1 - new_side], addl_edges[new_side],
                                               new_nodes[j] );RR;
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::foc_get_addl_ents( std::vector< EntityHandle >* star_dp1,
                                       std::vector< EntityHandle >* /*star_dp2*/,
                                       EntityHandle* split_edges,
                                       EntityHandle split_node,
                                       Range* addl_ents )
{
    // if we're splitting 2 edges, there might be other edges that have the split
    // node; also need to know which side they're on

    // algorithm: for a given star_dp1 (faces) on a side:
    // - get all edges adj to all faces -> R1
    // - get all edges adj to split_node -> R2
    // - R3 = R1 & R2 (edges connected to split_node & adj to a star face)
    // - R3 -= split_edges (take split edges off addl_ents)

    Range R2;
    MeshTopoUtil mtu( mbImpl );
    ErrorCode result = mbImpl->get_adjacencies( &split_node, 1, 1, false, R2 );RR;
    Range::iterator rit;

    for( int i = 0; i < 2; i++ )
    {
        Range R1, R3;
        result = mbImpl->get_adjacencies( &star_dp1[i][0], star_dp1[i].size(), 1, false, R1, Interface::UNION );RR;
        R3 = intersect( R1, R2 );
        for( int j = 0; j < 3; j++ )
            if( split_edges[j] ) R3.erase( split_edges[j] );
        addl_ents[i].merge( R3 );
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::foc_get_stars( EntityHandle* split_quads,
                                   EntityHandle* split_edges,
                                   std::vector< EntityHandle >* star_dp1,
                                   std::vector< EntityHandle >* star_dp2 )
{
    bool on_bdy = false, on_bdy_tmp;
    ErrorCode result;
    MeshTopoUtil mtu( mbImpl );

    // get the star around the split_edge
    std::vector< EntityHandle > qstar, hstar;
    unsigned int qpos = 0;

    for( int i = 0; i < 3; i++ )
    {
        if( !split_edges[i] ) break;

        // get the star around this split edge
        unsigned int qpos_tmp = 0;
        std::vector< EntityHandle > qstar_tmp, hstar_tmp;
        result = mtu.star_entities( split_edges[i], qstar_tmp, on_bdy_tmp, 0, &hstar_tmp );RR;
        // if we're on the bdy, add a null to the hex star too
        if( on_bdy_tmp )
        {
            assert( hstar_tmp.size() == qstar_tmp.size() - 1 );
            hstar_tmp.push_back( 0 );
            on_bdy = true;
        }

        // get the position of first split quad in star
        while( qpos_tmp < qstar_tmp.size() && qstar_tmp[qpos_tmp] != split_quads[0] )
            qpos_tmp++;
        if( qpos_tmp == qstar_tmp.size() ) return MB_FAILURE;

        bool forward;
        // 1st iteration is forward by definition
        if( 0 == i ) forward = true;

        // need to be careful about direction on later iters
        else if( hstar[qpos] == hstar_tmp[qpos_tmp] )
            forward = true;
        else if( hstar[qpos] == hstar_tmp[( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size()] &&
                 hstar_tmp[qpos_tmp] == hstar[( qpos + qstar.size() - 1 ) % qstar.size()] )
            forward = false;
        else
            return MB_FAILURE;

        if( forward )
        {
            // 1st half of star
            // save hex right after split_quad[0] first
            star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
            qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
            while( qstar_tmp[qpos_tmp] != split_quads[1] )
            {
                star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
                star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
                qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
            }
            // 2nd half of star
            // save hex right after split_quad[1] first
            star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
            qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
            while( qstar_tmp[qpos_tmp] != split_quads[0] )
            {
                star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
                star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
                qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
            }
        }
        else
        {
            // go in reverse - take prev hex instead of current
            // one, and step in reverse

            // save hex right after split_quad[0] first
            qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
            star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
            while( qstar_tmp[qpos_tmp] != split_quads[1] )
            {
                star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
                qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
                star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
            }
            // 2nd half of star
            // save hex right after split_quad[1] first
            qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
            star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
            while( qstar_tmp[qpos_tmp] != split_quads[0] )
            {
                star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
                qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
                star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
            }
        }

        if( 0 == i )
        {
            // if we're on the first iteration, save results and continue, other iters
            // get compared to this one
            qstar.swap( qstar_tmp );
            hstar.swap( hstar_tmp );
            on_bdy = on_bdy_tmp;
            qpos   = qpos_tmp;
        }
    }

    // split quads go on list with NULLs, if any, otherwise on 2nd
    if( on_bdy )
    {
        if( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) != star_dp2[0].end() )
        {
            // remove *all* the zeros
            star_dp2[0].erase( std::remove( star_dp2[0].begin(), star_dp2[0].end(), 0 ), star_dp2[0].end() );
            // put the split quads on this half
            star_dp1[0].push_back( split_quads[0] );
            star_dp1[0].push_back( split_quads[1] );
        }
        else
        {
            star_dp2[1].erase( std::remove( star_dp2[1].begin(), star_dp2[1].end(), 0 ), star_dp2[1].end() );
            // put the split quads on this half
            star_dp1[1].push_back( split_quads[0] );
            star_dp1[1].push_back( split_quads[1] );
        }
    }
    else
    {
        star_dp1[1].push_back( split_quads[0] );
        star_dp1[1].push_back( split_quads[1] );
    }

    // some error checking
    if( !( ( ( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) == star_dp1[0].end() &&
               std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) == star_dp1[0].end() &&
               std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() &&
               std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) != star_dp1[1].end() ) ||

             ( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) == star_dp1[1].end() &&
               std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) == star_dp1[1].end() &&
               std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() &&
               std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) != star_dp1[0].end() ) ) ) )
    {
        std::cerr << "foc_get_stars: both split quads should be on the same star list half and not "
                  << "on the other, failed" << std::endl;
        return MB_FAILURE;
    }

    if( !( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) == star_dp2[0].end() &&
           std::find( star_dp2[1].begin(), star_dp2[1].end(), 0 ) == star_dp2[1].end() ) )
    {
        std::cerr << "foc_get_stars: no NULLs on the hstar lists, failed";
        return MB_FAILURE;
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::foc_delete_dual( EntityHandle* split_quads, EntityHandle* split_edges, Range& hexes )
{
    // special delete dual procedure, because in some cases we need to delete
    // a sheet too since it'll get merged into another

    // figure out whether we'll need to delete a sheet
    EntityHandle sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[0] ) );
    if( split_edges[1] ) sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[1] ) );
    EntityHandle chordl = get_dual_hyperplane( get_dual_entity( split_quads[0] ) );
    EntityHandle chordr = get_dual_hyperplane( get_dual_entity( split_quads[1] ) );
    assert( 0 != sheet1 && 0 != chordl && 0 != chordr );
    Range parentsl, parentsr;
    ErrorCode result = mbImpl->get_parent_meshsets( chordl, parentsl );
    if( MB_SUCCESS != result ) return result;
    result = mbImpl->get_parent_meshsets( chordr, parentsr );
    if( MB_SUCCESS != result ) return result;
    parentsl.erase( sheet1 );
    parentsr.erase( sheet1 );

    // before deciding which one to delete, collect the other cells which must
    // be deleted, and all the chords/sheets they're on
    Range adj_ents, dual_ents, cells1or2;
    for( int i = 0; i < 3; i++ )
    {
        result = mbImpl->get_adjacencies( hexes, i, false, adj_ents, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
    }

    // cache any adjacent hexes, for rebuilding the dual later
    result = mbImpl->get_adjacencies( adj_ents, 3, false, hexes, Interface::UNION );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator rit = adj_ents.begin(); rit != adj_ents.end(); ++rit )
    {
        EntityHandle this_ent = get_dual_entity( *rit );
        dual_ents.insert( this_ent );
        int dim = mbImpl->dimension_from_handle( this_ent );
        if( 1 == dim || 2 == dim ) cells1or2.insert( this_ent );
    }

    Range dual_hps;
    for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
        dual_hps.insert( get_dual_hyperplane( *rit ) );

    result = delete_dual_entities( dual_ents );
    if( MB_SUCCESS != result ) return result;

    // now decide which sheet to delete (to be merged into the other)
    EntityHandle sheet_delete = 0;
    if( is_blind( *parentsl.begin() ) )
        sheet_delete = *parentsl.begin();
    else if( is_blind( *parentsr.begin() ) )
        sheet_delete = *parentsr.begin();
    else
    {
        // neither is blind, take the one with fewer cells
        Range tmp_ents;
        int numl, numr;
        result = mbImpl->get_number_entities_by_handle( *parentsl.begin(), numl );
        if( MB_SUCCESS != result ) return result;
        result = mbImpl->get_number_entities_by_handle( *parentsr.begin(), numr );
        if( MB_SUCCESS != result ) return result;
        sheet_delete = ( numl > numr ? *parentsr.begin() : *parentsl.begin() );
    }
    assert( 0 != sheet_delete );

    // after deleting cells, check for empty chords & sheets, and delete those too
    for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
    {
        Range tmp_ents;
        result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
        if( MB_SUCCESS != result ) return result;
        if( tmp_ents.empty() )
        {
            result = mbImpl->delete_entities( &( *rit ), 1 );
            if( MB_SUCCESS != result ) return result;
        }
        else if( *rit == sheet_delete )
        {
            // delete the sheet
            result = mbImpl->delete_entities( &( *rit ), 1 );
            if( MB_SUCCESS != result ) return result;
        }
    }

    // now just to be safe, add the hexes bridge-adjacent across vertices
    // to the hexes we already have
    Range tmp_hexes;
    MeshTopoUtil mtu( mbImpl );
    for( Range::iterator rit = hexes.begin(); rit != hexes.end(); ++rit )
    {
        result = mtu.get_bridge_adjacencies( *rit, 0, 3, tmp_hexes );
        if( MB_SUCCESS != result ) return result;
    }
    hexes.merge( tmp_hexes );

    return MB_SUCCESS;
}

//! returns true if all vertices are dual to hexes (not faces)
bool DualTool::is_blind( const EntityHandle chord_or_sheet )
{
    // must be an entity set
    if( TYPE_FROM_HANDLE( chord_or_sheet ) != MBENTITYSET ) return false;

    // get the vertices
    Range verts, ents;
    ErrorCode result = mbImpl->get_entities_by_handle( chord_or_sheet, ents );
    if( MB_SUCCESS != result || ents.empty() ) return false;

    result = mbImpl->get_adjacencies( ents, 0, false, verts, Interface::UNION );
    if( MB_SUCCESS != result || verts.empty() ) return false;

    for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
    {
        // get dual entity for this vertex
        EntityHandle dual_ent = get_dual_entity( *rit );
        if( 0 == dual_ent ) continue;
        if( TYPE_FROM_HANDLE( dual_ent ) == MBQUAD ) return false;
    }

    // if none of the vertices' duals were quads, chord_or_sheet must be blind
    return true;
}

//! given a 1-cell and a chord, return the neighboring vertices on the
//! chord, in the same order as the 1-cell's vertices
ErrorCode DualTool::get_opposite_verts( const EntityHandle middle_edge, const EntityHandle chord, EntityHandle* verts )
{
    // get the edges on the chord, in order, and move to middle_edge
    std::vector< EntityHandle > chord_edges;
    const EntityHandle* connect;
    int num_connect;

    ErrorCode result = mbImpl->get_entities_by_handle( chord, chord_edges );RR;
    std::vector< EntityHandle >::iterator vit = std::find( chord_edges.begin(), chord_edges.end(), middle_edge );
    result                                    = mbImpl->get_connectivity( middle_edge, connect, num_connect );RR;

    if(
        // middle_edge isn't on this chord
        vit == chord_edges.end() ||
        // chord only has 1 edge
        chord_edges.size() == 1 ||
        // middle_edge is at beginning or end and chord isn't blind
        ( ( vit == chord_edges.begin() || vit == chord_edges.end() - 1 ) && !is_blind( chord ) ) )
        return MB_FAILURE;

    else if( chord_edges.size() == 2 )
    {
        // easier if it's a 2-edge blind chord, just get vertices in opposite order
        verts[0] = connect[1];
        verts[1] = connect[0];
        return MB_SUCCESS;
    }

    // get vertices with the prev edge & subtract vertices of 1-cell
    if( vit == chord_edges.begin() )
        vit = chord_edges.end() - 1;
    else
        --vit;
    Range dum_connect, middle_connect;
    result = mbImpl->get_connectivity( &middle_edge, 1, middle_connect );RR;
    result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
    dum_connect = subtract( dum_connect, middle_connect );
    if( dum_connect.size() != 1 )
    {
        std::cerr << "Trouble traversing chord." << std::endl;
        return MB_FAILURE;
    }

    // put in verts[0]
    verts[0] = *dum_connect.begin();

    // same with prev edge
    ++vit;
    if( vit == chord_edges.end() ) vit = chord_edges.begin();
    ++vit;
    dum_connect.clear();
    result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
    dum_connect = subtract( dum_connect, middle_connect );
    if( dum_connect.size() != 1 )
    {
        std::cerr << "Trouble traversing chord." << std::endl;
        return MB_FAILURE;
    }

    // put in verts[1]
    verts[1] = *dum_connect.begin();

    // if verts[0] and 1st vertex of 1cell don't have common edge, switch verts
    MeshTopoUtil mtu( mbImpl );
    if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
    {
        EntityHandle dum_h = verts[0];
        verts[0]           = verts[1];
        verts[1]           = dum_h;
    }

    if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
    {
        std::cerr << "Trouble traversing chord." << std::endl;
        return MB_FAILURE;
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::get_dual_entities( const EntityHandle dual_ent,
                                       Range* dcells,
                                       Range* dedges,
                                       Range* dverts,
                                       Range* dverts_loop,
                                       Range* dedges_loop )
{
    ErrorCode result = MB_SUCCESS;

    if( NULL != dcells )
    {
        result = mbImpl->get_entities_by_type( dual_ent, MBPOLYGON, *dcells );
        if( MB_SUCCESS != result ) return result;
    }

    if( NULL != dedges )
    {
        if( NULL != dcells )
            result = mbImpl->get_adjacencies( *dcells, 1, false, *dedges, Interface::UNION );
        else
            result = mbImpl->get_entities_by_type( dual_ent, MBEDGE, *dedges );

        if( MB_SUCCESS != result ) return result;
    }

    if( NULL != dverts )
    {
        if( NULL != dcells )
            result = mbImpl->get_adjacencies( *dcells, 0, false, *dverts, Interface::UNION );
        else if( NULL != dedges )
            result = mbImpl->get_adjacencies( *dedges, 0, false, *dverts, Interface::UNION );
        else
        {
            Range all_ents;
            result = mbImpl->get_entities_by_handle( dual_ent, all_ents );RR;
            result = mbImpl->get_adjacencies( all_ents, 0, false, *dverts, Interface::UNION );
        }

        if( MB_SUCCESS != result ) return result;
    }

    if( NULL != dverts_loop && NULL != dverts )
    {
        static std::vector< EntityHandle > dual_ents;
        dual_ents.resize( dverts->size() );
        result = mbImpl->tag_get_data( dualEntity_tag(), *dverts, &dual_ents[0] );
        if( MB_SUCCESS != result ) return result;
        Range::iterator rit;
        unsigned int i;
        for( rit = dverts->begin(), i = 0; rit != dverts->end(); ++rit, i++ )
            if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBQUAD ) dverts_loop->insert( *rit );
    }

    if( NULL != dedges_loop && NULL != dedges )
    {
        static std::vector< EntityHandle > dual_ents;
        dual_ents.resize( dedges->size() );
        result = mbImpl->tag_get_data( dualEntity_tag(), *dedges, &dual_ents[0] );
        if( MB_SUCCESS != result ) return result;
        Range::iterator rit;
        unsigned int i;
        for( rit = dedges->begin(), i = 0; rit != dedges->end(); ++rit, i++ )
            if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBEDGE ) dedges_loop->insert( *rit );
    }

    return result;
}

ErrorCode DualTool::list_entities( const EntityHandle* entities, const int num_entities ) const
{
    Range temp_range;
    ErrorCode result;
    if( NULL == entities && 0 == num_entities )
        return mbImpl->list_entities( entities, num_entities );

    else if( NULL == entities && 0 < num_entities )
    {

        // list all entities of all types
        std::cout << std::endl;
        for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
        {
            result = mbImpl->get_entities_by_type( 0, this_type, temp_range );
            if( MB_SUCCESS != result ) return result;
        }
    }

    else
    {
        std::copy( entities, entities + num_entities, range_inserter( temp_range ) );
    }

    return list_entities( temp_range );
}

ErrorCode DualTool::list_entities( const Range& entities ) const
{
    // now print each entity, listing the dual information first then calling Interface to do
    // the rest
    ErrorCode result = MB_SUCCESS, tmp_result;
    for( Range::const_iterator iter = entities.begin(); iter != entities.end(); ++iter )
    {
        EntityType this_type = TYPE_FROM_HANDLE( *iter );
        std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( *iter ) << ":" << std::endl;

        EntityHandle dual_ent = get_dual_entity( *iter );
        if( 0 != dual_ent )
        {
            std::cout << "Dual to " << CN::EntityTypeName( mbImpl->type_from_handle( dual_ent ) ) << " "
                      << mbImpl->id_from_handle( dual_ent ) << std::endl;
        }

        if( TYPE_FROM_HANDLE( *iter ) == MBENTITYSET )
        {
            EntityHandle chord = 0, sheet = 0;
            int id;
            result = mbImpl->tag_get_data( dualCurve_tag(), &( *iter ), 1, &chord );
            if( MB_SUCCESS != result ) return result;
            result = mbImpl->tag_get_data( dualSurface_tag(), &( *iter ), 1, &sheet );
            if( MB_SUCCESS != result ) return result;
            result = mbImpl->tag_get_data( globalId_tag(), &( *iter ), 1, &id );
            if( MB_SUCCESS != result ) return result;

            if( 0 != chord ) std::cout << "(Dual chord " << id << ")" << std::endl;
            if( 0 != sheet ) std::cout << "(Dual sheet " << id << ")" << std::endl;
        }

        tmp_result = mbImpl->list_entity( *iter );
        if( MB_SUCCESS != tmp_result ) result = tmp_result;
    }

    return result;
}

ErrorCode DualTool::face_shrink( EntityHandle odedge )
{
    // some preliminary checking
    if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;

    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    std::cout << "FS(";
    print_cell( odedge );
    std::cout << ")" << std::endl;

    EntityHandle quads[4], hexes[2];
    std::vector< EntityHandle > connects[4], side_quads[2];

    // get the quads along the chord through the 2 hexes, and the vertices
    // for those quads
    ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
    if( MB_SUCCESS != result ) return result;

    // flip/rotate connect arrays so they align & are same sense
    result = fs_check_quad_sense( hexes[0], quads[0], connects );
    if( MB_SUCCESS != result ) return result;

    // get the quad loops along the "side" surfaces
    result = fs_get_quad_loops( hexes, connects, side_quads );
    if( MB_SUCCESS != result ) return result;

    // ok, done with setup; now delete dual entities affected by this operation,
    // which is all the entities adjacent to vertices of dual edge
    Range adj_verts, adj_edges, dual_ents, cells1or2;
    MeshTopoUtil mtu( mbImpl );
    result = mtu.get_bridge_adjacencies( odedge, 0, 1, adj_edges );
    if( MB_SUCCESS != result ) return result;
    result = mbImpl->get_adjacencies( adj_edges, 0, false, adj_verts, Interface::UNION );
    if( MB_SUCCESS != result ) return result;
    for( int i = 1; i <= 3; i++ )
    {
        result = mbImpl->get_adjacencies( adj_verts, i, false, dual_ents, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
    }

    // before deleting dual, grab the 1- and 2-cells
    for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
    {
        int dim = mbImpl->dimension_from_handle( *rit );
        if( 1 == dim || 2 == dim ) cells1or2.insert( *rit );
    }
    Range dual_hps;
    for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
        dual_hps.insert( get_dual_hyperplane( *rit ) );

    dual_ents.insert( odedge );
    result = delete_dual_entities( dual_ents );
    if( MB_SUCCESS != result ) return result;

    // after deleting cells, check for empty chords & sheets, and delete those too
    for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
    {
        Range tmp_ents;
        result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
        if( MB_SUCCESS != result ) return result;
        if( tmp_ents.empty() )
        {
            result = mbImpl->delete_entities( &( *rit ), 1 );
            if( MB_SUCCESS != result ) return result;
        }
    }

    // remove any explicit adjacencies between side quads and hexes; don't check
    // for error, since there might not be adjacencies
    for( int i = 0; i < 4; i++ )
    {
        for( int j = 0; j < 2; j++ )
        {
            result = mbImpl->remove_adjacencies( side_quads[j][i], &hexes[j], 1 );
        }
    }

    // make inner ring of vertices
    // get centroid of quad2
    double q2coords[12], avg[3] = { 0.0, 0.0, 0.0 };
    result = mbImpl->get_coords( &connects[1][0], 4, q2coords );
    if( MB_SUCCESS != result ) return result;
    for( int i = 0; i < 4; i++ )
    {
        avg[0] += q2coords[3 * i];
        avg[1] += q2coords[3 * i + 1];
        avg[2] += q2coords[3 * i + 2];
    }
    avg[0] *= .25;
    avg[1] *= .25;
    avg[2] *= .25;
    // position new vertices
    connects[3].resize( 4 );
    for( int i = 0; i < 4; i++ )
    {
        q2coords[3 * i]     = avg[0] + .25 * ( q2coords[3 * i] - avg[0] );
        q2coords[3 * i + 1] = avg[1] + .25 * ( q2coords[3 * i + 1] - avg[1] );
        q2coords[3 * i + 2] = avg[2] + .25 * ( q2coords[3 * i + 2] - avg[2] );
        result              = mbImpl->create_vertex( &q2coords[3 * i], connects[3][i] );
        if( MB_SUCCESS != result ) return result;
    }

    // ok, now have the 4 connectivity arrays for 4 quads; construct hexes
    EntityHandle hconnect[8], new_hexes[4];
    int new_hex_ids[4];

    for( int i = 0; i < 4; i++ )
    {
        int i1 = i, i2 = ( i + 1 ) % 4;
        hconnect[0] = connects[0][i1];
        hconnect[1] = connects[0][i2];
        hconnect[2] = connects[3][i2];
        hconnect[3] = connects[3][i1];

        hconnect[4] = connects[1][i1];
        hconnect[5] = connects[1][i2];
        hconnect[6] = connects[2][i2];
        hconnect[7] = connects[2][i1];

        result = mbImpl->create_element( MBHEX, hconnect, 8, new_hexes[i] );
        if( MB_SUCCESS != result ) return result;

        // test for equiv entities from the side quads, and make explicit adjacencies
        // if there are any
        for( int j = 0; j < 2; j++ )
        {
            if( mtu.equivalent_entities( side_quads[j][i] ) )
            {
                result = mbImpl->add_adjacencies( side_quads[j][i], &new_hexes[i], 1, false );
                if( MB_SUCCESS != result ) return result;
            }
        }

        new_hex_ids[i] = ++maxHexId;
    }

    // set the global id tag on the new hexes
    result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 4, new_hex_ids );
    if( MB_SUCCESS != result ) return result;

    // now fixup other two hexes; start by getting hex through quads 0, 1
    // make this first hex switch to the other side, to make the dual look like
    // a hex push
    int tmp_ids[2];
    result = mbImpl->tag_get_data( globalId_tag(), hexes, 2, tmp_ids );
    if( MB_SUCCESS != result ) return result;

    result = mbImpl->delete_entities( hexes, 2 );
    if( MB_SUCCESS != result ) return result;
    result = mbImpl->delete_entities( &quads[1], 1 );
    if( MB_SUCCESS != result ) return result;
    for( int i = 0; i < 4; i++ )
    {
        hconnect[i]     = connects[3][i];
        hconnect[4 + i] = connects[2][i];
    }
    result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[0] );
    if( MB_SUCCESS != result ) return result;

    for( int i = 0; i < 4; i++ )
    {
        hconnect[i]     = connects[0][i];
        hconnect[4 + i] = connects[3][i];
    }
    result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[1] );
    if( MB_SUCCESS != result ) return result;

    // check for and fix any explicit adjacencies on either end quad
    if( mtu.equivalent_entities( quads[0] ) ) mbImpl->add_adjacencies( quads[0], &hexes[1], 1, false );
    if( mtu.equivalent_entities( quads[2] ) ) mbImpl->add_adjacencies( quads[2], &hexes[0], 1, false );

    // re-set the global ids for the hexes to what they were
    result = mbImpl->tag_set_data( globalId_tag(), hexes, 2, tmp_ids );
    if( MB_SUCCESS != result ) return result;

    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    // now update the dual
    Range tmph;
    result = mtu.get_bridge_adjacencies( hexes[0], 0, 3, tmph );
    if( MB_SUCCESS != result ) return result;
    result = mtu.get_bridge_adjacencies( hexes[1], 0, 3, tmph );
    if( MB_SUCCESS != result ) return result;
    tmph.insert( hexes[1] );
    result = construct_hex_dual( tmph );
    if( MB_SUCCESS != result ) return result;

    return result;
}

ErrorCode DualTool::fs_get_quad_loops( EntityHandle* hexes,
                                       std::vector< EntityHandle >* connects,
                                       std::vector< EntityHandle >* side_quads )
{
    for( int i = 0; i < 4; i++ )
    {
        for( int j = 0; j < 2; j++ )
        {
            Range adj_ents, dum_quads;
            adj_ents.insert( hexes[j] );
            adj_ents.insert( connects[j][i] );
            adj_ents.insert( connects[j][( i + 1 ) % 4] );
            adj_ents.insert( connects[j + 1][i] );
            adj_ents.insert( connects[j + 1][( i + 1 ) % 4] );

            ErrorCode result = mbImpl->get_adjacencies( adj_ents, 2, false, dum_quads );
            if( MB_SUCCESS != result ) return result;
            assert( 1 == dum_quads.size() );
            side_quads[j].push_back( *dum_quads.begin() );
        }
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::fs_check_quad_sense( EntityHandle hex0, EntityHandle quad0, std::vector< EntityHandle >* connects )
{
    // check sense of 0th quad wrt hex; since sense is out of element,
    // switch if quad is NOT reversed wrt hex
    int dum1, dum2, sense = 0;
    ErrorCode result = mbImpl->side_number( hex0, quad0, dum1, sense, dum2 );
    if( MB_SUCCESS != result ) return result;
    assert( 0 != sense );
    if( 1 == sense )
    {
        // just switch sense of this one; others will get switched next
        EntityHandle dum = connects[0][0];
        connects[0][0]   = connects[0][2];
        connects[0][2]   = dum;
    }

    // check sense of 1st, 2nd quads, rotate if necessary to align connect arrays
    int index0 = -1, index2 = -1, sense0 = 0, sense2 = 0;
    MeshTopoUtil mtu( mbImpl );
    for( int i = 0; i < 4; i++ )
    {
        if( 0 != mtu.common_entity( connects[0][0], connects[1][i], 1 ) )
        {
            index0 = i;
            if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 1 ) % 4], 1 ) )
                sense0 = 1;
            else if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 4 - 1 ) % 4], 1 ) )
                sense0 = -1;
            break;
        }
    }

    assert( index0 != -1 && sense0 != 0 );

    if( sense0 == -1 )
    {
        EntityHandle dumh = connects[1][0];
        connects[1][0]    = connects[1][2];
        connects[1][2]    = dumh;
        if( index0 % 2 == 0 ) index0 = ( index0 + 2 ) % 4;
    }

    if( index0 != 0 )
    {
        std::vector< EntityHandle > tmpc;
        for( int i = 0; i < 4; i++ )
            tmpc.push_back( connects[1][( index0 + i ) % 4] );
        connects[1].swap( tmpc );
    }

    for( int i = 0; i < 4; i++ )
    {
        if( 0 != mtu.common_entity( connects[1][0], connects[2][i], 1 ) )
        {
            index2 = i;
            if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 1 ) % 4], 1 ) )
                sense2 = 1;
            else if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 4 - 1 ) % 4], 1 ) )
                sense2 = -1;
            break;
        }
    }

    assert( index2 != -1 && sense2 != 0 );

    if( sense2 == -1 )
    {
        EntityHandle dumh = connects[2][0];
        connects[2][0]    = connects[2][2];
        connects[2][2]    = dumh;
        if( index2 % 2 == 0 ) index2 = ( index2 + 2 ) % 4;
    }

    if( index2 != 0 )
    {
        std::vector< EntityHandle > tmpc;
        for( int i = 0; i < 4; i++ )
            tmpc.push_back( connects[2][( index2 + i ) % 4] );
        connects[2].swap( tmpc );
    }

    return MB_SUCCESS;
}

//! effect reverse face shrink operation
ErrorCode DualTool::rev_face_shrink( EntityHandle odedge )
{
    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    // some preliminary checking
    if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;

    std::cout << "-FS(";
    print_cell( odedge );
    std::cout << ")" << std::endl;

    EntityHandle quads[4], hexes[2];
    std::vector< EntityHandle > connects[4], side_quads[2];

    // get three quads (shared quad & 2 end quads), hexes, and quad
    // connects
    ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
    if( MB_SUCCESS != result ) return result;

    // adjust sense & rotation so they're aligned, together & wrt first
    // hex
    result = fs_check_quad_sense( hexes[0], quads[0], connects );
    if( MB_SUCCESS != result ) return result;

    result = fsr_get_fourth_quad( connects, side_quads );
    if( MB_SUCCESS != result )
    {
        std::cout << "Can't do -FS here, two hexes must be adjacent to ring of 4 hexes." << std::endl;
        return result;
    }

    Range adj_ents, outer_hexes, all_adjs;

    // first get the entities connected to interior 4 verts
    for( int i = 1; i <= 3; i++ )
    {
        result = mbImpl->get_adjacencies( &connects[1][0], 4, i, false, adj_ents, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
    }

    // next get all entities adjacent to those; these will have their dual
    // entities deleted
    for( int i = 0; i < 3; i++ )
    {
        result = mbImpl->get_adjacencies( adj_ents, i, false, all_adjs, Interface::UNION );
        if( MB_SUCCESS != result ) return result;
    }

    // get the dual entities and delete them
    Range dual_ents, dual_hps;
    for( Range::iterator rit = all_adjs.begin(); rit != all_adjs.end(); ++rit )
    {
        EntityHandle this_ent = get_dual_entity( *rit );
        dual_ents.insert( this_ent );
    }

    // before deleting dual, grab the 1- and 2-cells
    for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
    {
        int dim = mbImpl->dimension_from_handle( *rit );
        if( 1 == dim || 2 == dim ) dual_hps.insert( get_dual_hyperplane( *rit ) );
    }

    result = delete_dual_entities( dual_ents );
    if( MB_SUCCESS != result ) return result;

    // after deleting cells, check for empty chords & sheets, and delete those too
    for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
    {
        Range tmp_ents;
        result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
        if( MB_SUCCESS != result ) return result;
        if( tmp_ents.empty() )
        {
            result = mbImpl->delete_entities( &( *rit ), 1 );
            if( MB_SUCCESS != result ) return result;
        }
    }

    // before re-connecting two hexes, check for existing quad on 4th quad vertices;
    // if there is a quad there, need to add explicit adjs to any adj hexes, since
    // by definition there'll be another quad on those vertices
    bool need_explicit = false;
    Range adj_quads;
    result = mbImpl->get_adjacencies( &connects[3][0], 4, 2, false, adj_quads );
    if( MB_MULTIPLE_ENTITIES_FOUND == result || !adj_quads.empty() )
    {
        // there's already a quad for these 4 vertices; by definition,
        // we'll be creating equivalent entities, so that original quad
        // needs explicit adj's to its bounding elements
        need_explicit = true;
        for( Range::iterator rit = adj_quads.begin(); rit != adj_quads.end(); ++rit )
        {
            Range adj_hexes;
            result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, adj_hexes );RR;
            result = mbImpl->add_adjacencies( *rit, adj_hexes, false );RR;
        }
    }

    // re-connect the two hexes
    std::vector< EntityHandle > new_connect;
    std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
    std::copy( connects[2].begin(), connects[2].end(), std::back_inserter( new_connect ) );
    result = mbImpl->set_connectivity( hexes[0], &new_connect[0], 8 );
    if( MB_SUCCESS != result ) return result;

    new_connect.clear();
    std::copy( connects[0].begin(), connects[0].end(), std::back_inserter( new_connect ) );
    std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
    result = mbImpl->set_connectivity( hexes[1], &new_connect[0], 8 );
    if( MB_SUCCESS != result ) return result;

    // test for equiv entities from the side quads, and make explicit
    // adjacencies if there are any
    MeshTopoUtil mtu( mbImpl );
    for( int j = 0; j < 2; j++ )
    {
        for( int i = 0; i < 4; i++ )
        {
            if( mtu.equivalent_entities( side_quads[j][i] ) )
            {
                result = mbImpl->add_adjacencies( side_quads[j][i], &hexes[j], 1, false );
                if( MB_SUCCESS != result ) return result;
            }
        }
    }

    // remove hexes we want to keep
    adj_ents.erase( hexes[0] );
    adj_ents.erase( hexes[1] );

    // delete the other interior entities
    result = mbImpl->delete_entities( adj_ents );
    if( MB_SUCCESS != result ) return result;

    EntityHandle new_quad;
    result = mbImpl->create_element( MBQUAD, &connects[3][0], 4, new_quad );RR;
    if( need_explicit )
    {
        result = mbImpl->add_adjacencies( new_quad, hexes, 2, false );RR;
    }

    if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();

    // now update the dual
    result = construct_hex_dual( hexes, 2 );
    if( MB_SUCCESS != result ) return result;

    return MB_SUCCESS;
}

ErrorCode DualTool::fsr_get_fourth_quad( std::vector< EntityHandle >* connects,
                                         std::vector< EntityHandle >* side_quads )
{
    // given the first three quad connectivities in ordered vectors, get the fourth,
    // where the fourth is really the 4 vertices originally shared by the 2 hexes
    // before the face shrink on them

    // vertex on 4th quad is in quad adj to other 3 verts
    for( int i = 0; i < 4; i++ )
    {
        Range start_verts, tmp_verts, quads;
        for( int j = 0; j < 3; j++ )
            start_verts.insert( connects[j][i] );
        ErrorCode result = mbImpl->get_adjacencies( start_verts, 2, false, quads );
        if( MB_SUCCESS != result ) return result;
        assert( quads.size() == 1 );
        result = mbImpl->get_adjacencies( &( *quads.begin() ), 1, 0, false, tmp_verts );RR;
        tmp_verts = subtract( tmp_verts, start_verts );
        assert( 1 == tmp_verts.size() );
        connects[3].push_back( *tmp_verts.begin() );
    }

    // now get the side quads
    for( int i = 0; i < 4; i++ )
    {
        Range dum_ents, hexes;
        dum_ents.insert( connects[1][i] );
        dum_ents.insert( connects[1][( i + 1 ) % 4] );
        dum_ents.insert( connects[3][i] );
        ErrorCode result = mbImpl->get_adjacencies( dum_ents, 3, false, hexes );
        if( MB_SUCCESS != result ) return result;
        assert( 1 == hexes.size() );

        hexes.insert( connects[0][i] );
        hexes.insert( connects[0][( i + 1 ) % 4] );
        hexes.insert( connects[3][i] );
        hexes.insert( connects[3][( i + 1 ) % 4] );
        dum_ents.clear();
        result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
        if( MB_SUCCESS != result ) return result;
        assert( dum_ents.size() == 1 );
        side_quads[0].push_back( *dum_ents.begin() );

        hexes.erase( connects[0][i] );
        hexes.erase( connects[0][( i + 1 ) % 4] );
        hexes.insert( connects[2][i] );
        hexes.insert( connects[2][( i + 1 ) % 4] );
        dum_ents.clear();
        result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
        if( MB_SUCCESS != result ) return result;
        side_quads[1].push_back( *dum_ents.begin() );
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::fs_get_quads( EntityHandle odedge,
                                  EntityHandle* quads,
                                  EntityHandle* hexes,
                                  std::vector< EntityHandle >* connects )
{
    // need to get the three quads along the chord
    EntityHandle chord = get_dual_hyperplane( odedge );
    if( 0 == chord ) return MB_FAILURE;

    std::vector< EntityHandle > edges;
    ErrorCode result = mbImpl->get_entities_by_handle( chord, edges );
    if( MB_FAILURE == result ) return result;

    std::vector< EntityHandle >::iterator vit = std::find( edges.begin(), edges.end(), odedge );
    // shouldn't be first or last edge on chord
    if( vit == edges.end() || *edges.begin() == *vit || *edges.rbegin() == *vit ) return MB_FAILURE;

    // get quads/connectivity for first 3 quads
    quads[0] = get_dual_entity( *( vit - 1 ) );
    quads[1] = get_dual_entity( *vit );
    quads[2] = get_dual_entity( *( vit + 1 ) );
    for( int i = 0; i < 3; i++ )
    {
        result = mbImpl->get_connectivity( &quads[i], 1, connects[i], true );
        if( MB_SUCCESS != result ) return result;
    }

    Range tmph;
    result = mbImpl->get_adjacencies( quads, 2, 3, false, tmph );
    if( MB_SUCCESS != result ) return result;
    assert( tmph.size() == 1 );
    hexes[0] = *tmph.begin();

    tmph.clear();
    result = mbImpl->get_adjacencies( &quads[1], 2, 3, false, tmph );
    if( MB_SUCCESS != result ) return result;
    assert( tmph.size() == 1 );
    hexes[1] = *tmph.begin();

    return MB_SUCCESS;
}

ErrorCode DualTool::delete_whole_dual()
{
    // delete dual hyperplanes
    Range dual_surfs, dual_curves;
    ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );RR;
    result = mbImpl->delete_entities( dual_surfs );RR;
    result = this->get_dual_hyperplanes( mbImpl, 1, dual_curves );RR;
    result = mbImpl->delete_entities( dual_curves );RR;

    // gather up all dual entities
    Range dual_ents;
    result = mbImpl->get_entities_by_type_and_tag( 0, MBVERTEX, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
    result = mbImpl->get_entities_by_type_and_tag( 0, MBEDGE, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
    result = mbImpl->get_entities_by_type_and_tag( 0, MBPOLYGON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
    result =
        mbImpl->get_entities_by_type_and_tag( 0, MBPOLYHEDRON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;

    // delete them, in reverse order of dimension
    ErrorCode tmp_result;
    for( Range::reverse_iterator rit = dual_ents.rbegin(); rit != dual_ents.rend(); ++rit )
    {
        tmp_result = mbImpl->delete_entities( &( *rit ), 1 );
        if( MB_SUCCESS != tmp_result ) result = tmp_result;
    }
    RR;

    // delete dual-related tags
    if( 0 != dualSurfaceTag )
    {
        tmp_result = mbImpl->tag_delete( dualSurfaceTag );
        if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
    }
    if( 0 != dualCurveTag )
    {
        tmp_result = mbImpl->tag_delete( dualCurveTag );
        if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
    }
    if( 0 != dualEntityTag )
    {
        tmp_result = mbImpl->tag_delete( dualEntityTag );
        if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
    }
    if( 0 != extraDualEntityTag )
    {
        tmp_result = mbImpl->tag_delete( extraDualEntityTag );
        if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
    }
    if( 0 != dualGraphicsPointTag )
    {
        tmp_result = mbImpl->tag_delete( dualGraphicsPointTag );
        if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
    }

    return MB_SUCCESS;
}

ErrorCode DualTool::check_dual_adjs()
{
    // check primal-dual correspondence

    // get the primal entities
    Range pents[4];
    ErrorCode result = mbImpl->get_entities_by_type( 0, MBHEX, pents[3] );RR;
    for( int i = 2; i >= 0; i-- )
    {
        result = mbImpl->get_adjacencies( pents[3], 2, false, pents[2], Interface::UNION );RR;
    }

    // for each primal entity of dimension pd
#define PRENT( ent ) CN::EntityTypeName( TYPE_FROM_HANDLE( ent ) ) << " " << ID_FROM_HANDLE( ent )
    ErrorCode overall_result = MB_SUCCESS;
    for( int pd = 1; pd <= 3; pd++ )
    {
        for( Range::iterator prit = pents[pd].begin(); prit != pents[pd].end(); ++prit )
        {
            // get corresponding dual entity of dimension dd = 3-pd
            EntityHandle dual_ent = get_dual_entity( *prit );
            if( 0 == dual_ent ) std::cerr << "Problem getting dual entity for " << PRENT( *prit ) << std::endl;

            // for each sub dimension sd = 0..pd-1
            for( int sd = 0; sd < pd; sd++ )
            {
                Range R1, R2, R3;
                //   R1 = entities bounding primal entity of dim sd
                result = mbImpl->get_adjacencies( &( *prit ), 1, sd, false, R1 );RR;

                //   R2 = entities bounded by dual entity, of dim 3-sd
                result = mbImpl->get_adjacencies( &dual_ent, 1, 3 - sd, false, R2 );RR;

                if( R1.size() != R2.size() )
                {
                    std::cerr << PRENT( *prit ) << ": number of adj ents in "
                              << "primal/dual don't agree for dimension " << sd << "." << std::endl;
                    overall_result = MB_FAILURE;
                }

                // for each entity in R1, get its dual and look for it in R2
                for( Range::iterator r1it = R1.begin(); r1it != R1.end(); ++r1it )
                {
                    EntityHandle tmp_dual = get_dual_entity( *r1it );
                    if( R2.find( tmp_dual ) == R2.end() )
                    {
                        std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r1it ) << " isn't adjacent in dual."
                                  << std::endl;
                        overall_result = MB_FAILURE;
                    }
                }
                // ditto for R2
                for( Range::iterator r2it = R2.begin(); r2it != R2.end(); ++r2it )
                {
                    EntityHandle tmp_prim = get_dual_entity( *r2it );
                    if( R1.find( tmp_prim ) == R1.end() )
                    {
                        std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r2it ) << " isn't adjacent in primal."
                                  << std::endl;
                        overall_result = MB_FAILURE;
                    }
                }
            }
        }
    }

    return overall_result;
}

}  // namespace moab