h5partial.cpp

Go to the documentation of this file.

#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "TestRunner.hpp"
#include "ReadHDF5.hpp"
#include "MBTagConventions.hpp"
#include "moab/FileOptions.hpp"

#ifdef MOAB_HAVE_MPI
#include "moab_mpi.h"
#endif

#include <vector>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <limits>

using namespace moab;

const char TEST_FILE[] = "partial.h5m";
#define READ_OPTS "BUFFER_SIZE=256"
const char ID_TAG_NAME[] = "test_id_tag";

static void test_read_nothing_common( bool non_existant );
static void test_read_nodes_common( int num_read_sets, bool blocked_coordinate_io );
static void test_read_handle_tag_common( bool var_len );

const int MBQUAD_INT       = 20;
const int NUM_SETS         = 10;
const int SET_WIDTH        = ( MBQUAD_INT + NUM_SETS - 1 ) / NUM_SETS;  // ceil(MBQUAD_INT/NUM_SETS)
const char LOGICAL_NAME[]  = "logical";                                 // tag storing logical (i,j) coordinates
const char CENTROID_NAME[] = "centroid";                                // tag storing position of centroid (x,y,0)
//! Create a regular MBQUAD_INT^2 element quad mesh with regularly
//! spaced coordinates in the range [1,100].  Group elements
//! into 10 vertical strips MBQUAD_INT/10 elements wide.  Tag elements,
//! vertices and/or sets with ID in [1,10] stored in ID_TAG_NAME
//! tag.  Write new mesh to TEST_FILE.
void create_mesh( bool create_element_sets,
                  bool create_vertex_sets,
                  bool tag_elements_with_id,
                  bool tag_vertices_with_id,
                  const char* adj_elem_tag_name = 0,
                  bool var_len_adj_elems        = false );
// Given a list of vertices adjacent to a quad strip, identify it as one of the
// NUM_SETS strips of quads written by create_mesh.
int identify_set( Interface& mb, const Range& verts );
int identify_set( Interface& mb, EntityHandle set );

static Tag check_tag( Interface& mb, const char* name, TagType storage, DataType type, int size );

enum GatherTestMode
{
    GATHER_SETS,
    GATHER_CONTENTS,
    GATHER_NONE
};
void test_gather_sets_common( bool contained_sets, GatherTestMode mode, bool no_parent_containing_sets = false );
void test_gather_sets_ranged( bool contained_sets, GatherTestMode mode, bool no_parent_containing_sets = false );

//! Read a set containing no entities
void test_read_empty_set()
{
    test_read_nothing_common( false );
}

//! Specify ID that doesn't exist in file
void test_read_non_existant_set()
{
    test_read_nothing_common( true );
}

//! Read in the nodes contained in a set.
void test_read_one_set_nodes()
{
    test_read_nodes_common( 1, false );
}

//! Read in the nodes contained in a set.
void test_read_one_set_nodes_blocked()
{
    test_read_nodes_common( 1, true );
}

//! Read in the elems contained in a set
void test_read_one_set_elems();

//! Read in the polyhedra contained in a set
void test_read_one_set_polyhedra();

//! Read in the sets contained in a set.
//! Should read all sets containing read elements or nodes
//! and all sets that are contained the the specified "read"
//! set.  Test the later here.
void test_read_set_sets();

//! Read in the nodes contained in a sets.
void test_read_two_sets_nodes()
{
    test_read_nodes_common( 2, false );
}

//! Read in the elems contained in a sets
void test_read_two_sets_elems();

//! For any set selected to be read by either explicit designation,
//! containing read entities, or contained in an explcitly designated
//! set, any child sets are also read.  Check that here.
void test_read_child_sets_only()
{
    test_gather_sets_common( false, GATHER_SETS );
    test_gather_sets_ranged( false, GATHER_SETS );
}
void test_read_child_set_contents()
{
    test_gather_sets_common( false, GATHER_CONTENTS );
    test_gather_sets_ranged( false, GATHER_CONTENTS );
}
void test_read_no_child_sets()
{
    test_gather_sets_common( false, GATHER_NONE );
    test_gather_sets_ranged( false, GATHER_NONE );
}

//! For any set selected to be read by either explicit designation,
//! containing read entities, or contained in an explcitly designated
//! set, any contained sets are also read.  Check that here.
void test_read_contained_sets_only()
{
    test_gather_sets_common( true, GATHER_SETS, true );
    test_gather_sets_ranged( true, GATHER_SETS );
}
void test_read_contained_set_contents()
{
    test_gather_sets_common( true, GATHER_CONTENTS, true );
    test_gather_sets_ranged( true, GATHER_CONTENTS );
}
void test_read_no_contained_sets()
{
    test_gather_sets_common( true, GATHER_NONE, true );
    test_gather_sets_ranged( true, GATHER_NONE );
}

//! Read in the sets contained in a set.
//! Should read all sets containing read elements or nodes
//! and all sets that are contained the the specified "read"
//! set.  Test the former here.
void test_read_containing_sets();

//! Test reading of explicit adjacencies
void test_read_adjacencies();

//! Test reading of sparse double tag data
void test_read_double_tag();

//! Test reading of sparse opaque tag data
void test_read_opaque_tag();

//! Test reading of sparse handle tag data
void test_read_handle_tag()
{
    test_read_handle_tag_common( false );
}

//! Test reading of variable-length tag data
void test_var_len_tag()
{
    test_read_handle_tag_common( true );
}

void test_read_tagged_elems();

void test_read_tagged_nodes();

void test_read_sides();

void test_read_ids();

void test_read_partial_ids();

int main( int argc, char* argv[] )
{
#ifdef MOAB_HAVE_MPI
    int fail = MPI_Init( &argc, &argv );
    if( fail ) return fail;
#endif

    REGISTER_TEST( test_read_empty_set );
    REGISTER_TEST( test_read_non_existant_set );
    REGISTER_TEST( test_read_one_set_nodes );
    REGISTER_TEST( test_read_one_set_nodes_blocked );
    REGISTER_TEST( test_read_one_set_elems );
    REGISTER_TEST( test_read_one_set_polyhedra );
    REGISTER_TEST( test_read_set_sets );
    REGISTER_TEST( test_read_two_sets_nodes );
    REGISTER_TEST( test_read_two_sets_elems );
    REGISTER_TEST( test_read_child_sets_only );
    REGISTER_TEST( test_read_child_set_contents );
    REGISTER_TEST( test_read_no_child_sets );
    REGISTER_TEST( test_read_contained_sets_only );
    REGISTER_TEST( test_read_contained_set_contents );
    REGISTER_TEST( test_read_no_contained_sets );
    REGISTER_TEST( test_read_containing_sets );
    REGISTER_TEST( test_read_double_tag );
    REGISTER_TEST( test_read_opaque_tag );
    REGISTER_TEST( test_read_handle_tag );
    REGISTER_TEST( test_var_len_tag );
    REGISTER_TEST( test_read_adjacencies );
    REGISTER_TEST( test_read_tagged_elems );
    REGISTER_TEST( test_read_tagged_nodes );
    REGISTER_TEST( test_read_sides );
    REGISTER_TEST( test_read_ids );
    REGISTER_TEST( test_read_partial_ids );
    int result = RUN_TESTS( argc, argv );

#ifdef MOAB_HAVE_MPI
    fail = MPI_Finalize();
    if( fail ) return fail;
#endif

    return result;
}

void test_read_nothing_common( bool non_existant )
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    // create a few nodes to write to file
    std::vector< double > coords( 3000 );
    Range verts;
    rval = mb.create_vertices( &coords[0], coords.size() / 3, verts );CHECK_ERR( rval );

    // create three entity sets
    EntityHandle sets[3];
    rval = mb.create_meshset( MESHSET_SET, sets[0] );CHECK_ERR( rval );
    rval = mb.create_meshset( MESHSET_SET, sets[1] );CHECK_ERR( rval );
    rval = mb.create_meshset( MESHSET_ORDERED, sets[2] );CHECK_ERR( rval );

    // put all vertices into two of the sets
    rval = mb.add_entities( sets[0], verts );CHECK_ERR( rval );
    rval = mb.add_entities( sets[2], verts );CHECK_ERR( rval );

    // tag all three sets
    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );
    int ids[3] = { 5, 7, 9 };
    rval       = mb.tag_set_data( id_tag, sets, 3, ids );CHECK_ERR( rval );

    // write mesh
    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );
    rval = mb.delete_mesh();CHECK_ERR( rval );

    // now read back in only the empty set
    EntityHandle file_set;
    int id = non_existant ? 8 : 7;
    rval   = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
    rval = mb.load_file( TEST_FILE, &file_set, READ_OPTS, ID_TAG_NAME, &id, 1 );
    if( non_existant )
    {
        CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
        return;
    }
    else
        CHECK_ERR( rval );

    // the file should contain exactly two sets (the specified one and the new
    // file set, and nothing else.)
    for( EntityType t = MBVERTEX; t < MBENTITYSET; ++t )
    {
        int count = -1;
        rval      = mb.get_number_entities_by_type( 0, t, count );CHECK_ERR( rval );
        CHECK_EQUAL( 0, count );
    }
    Range setrange;
    rval = mb.get_entities_by_type( 0, MBENTITYSET, setrange );CHECK_ERR( rval );
    CHECK_EQUAL( ( non_existant ? 1 : 2 ), (int)setrange.size() );
    CHECK( setrange.find( file_set ) != setrange.end() );
}

static void vtx_coords( int set_id, int j, int num_sets, double coords[3] )
{
    int i     = num_sets * j + set_id;
    coords[0] = i;
    coords[1] = i + 0.25;
    coords[2] = i + 0.5;
}

void test_read_nodes_common( int num_read_sets, bool blocked )
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    // create 1000 nodes
    const int num_sets = 2 * num_read_sets;
    std::vector< EntityHandle > verts( 1000 );
    std::vector< std::vector< EntityHandle > > set_verts( num_sets );
    for( size_t i = 0; i < verts.size(); ++i )
    {
        double coords[3];
        int j = i % num_sets;
        vtx_coords( j + 1, set_verts[j].size(), num_sets, coords );
        rval = mb.create_vertex( coords, verts[i] );
        set_verts[j].push_back( verts[i] );CHECK_ERR( rval );
    }

    // create two sets, each containing half of the nodes
    std::vector< EntityHandle > sets( num_sets );
    for( int i = 0; i < num_sets; ++i )
    {
        rval = mb.create_meshset( MESHSET_ORDERED, sets[i] );CHECK_ERR( rval );
        rval = mb.add_entities( sets[i], &set_verts[i][0], set_verts[i].size() );CHECK_ERR( rval );
    }

    // tag both sets
    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );
    std::vector< int > values( num_sets );
    for( int i = 0; i < num_sets; ++i )
        values[i] = i + 1;
    rval = mb.tag_set_data( id_tag, &sets[0], num_sets, &values[0] );CHECK_ERR( rval );

    // write file
    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );
    rval = mb.delete_mesh();CHECK_ERR( rval );

    // now read back in only the specified number of sets
    std::string opts( READ_OPTS );
    if( !opts.empty() ) opts += ';';
    if( blocked )
        opts += "BLOCKED_COORDINATE_IO=yes";
    else
        opts += "BLOCKED_COORDINATE_IO=no";

    values.resize( num_read_sets );
    for( int i = 0; i < num_read_sets; ++i )
        values[i] = 2 * ( i + 1 );
    EntityHandle file_set;
    rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
    rval = mb.load_file( TEST_FILE, &file_set, opts.c_str(), ID_TAG_NAME, &values[0], num_read_sets );CHECK_ERR( rval );

    int count, expected = 0;
    rval = mb.get_number_entities_by_dimension( 0, 0, count );CHECK_ERR( rval );
    for( int i = 0; i < num_sets; ++i )
        if( i % 2 ) expected += set_verts[i].size();
    CHECK_EQUAL( expected, count );

    Range sets2;
    rval = mb.get_entities_by_type( 0, MBENTITYSET, sets2 );CHECK_ERR( rval );
    CHECK_EQUAL( 1 + num_read_sets, (int)sets2.size() );
    Range::iterator it = sets2.find( file_set );
    CHECK( it != sets2.end() );
    sets2.erase( it );

    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag );CHECK_ERR( rval );
    while( !sets2.empty() )
    {
        EntityHandle set = sets2.pop_front();
        int id;
        rval = mb.tag_get_data( id_tag, &set, 1, &id );CHECK_ERR( rval );
        CHECK( std::find( values.begin(), values.end(), id ) != values.end() );
        CHECK( id > 0 );
        CHECK( (unsigned)id <= set_verts.size() );

        std::vector< EntityHandle > verts2;
        rval = mb.get_entities_by_handle( set, verts2 );CHECK_ERR( rval );
        CHECK_EQUAL( set_verts[id - 1].size(), verts2.size() );

        for( size_t i = 0; i < verts2.size(); ++i )
        {
            double exp_coords[3], coords[3];
            vtx_coords( id, i, num_sets, exp_coords );
            rval = mb.get_coords( &verts2[i], 1, coords );CHECK_ERR( rval );
            CHECK_REAL_EQUAL( exp_coords[0], coords[0], 1e-12 );
            CHECK_REAL_EQUAL( exp_coords[1], coords[1], 1e-12 );
            CHECK_REAL_EQUAL( exp_coords[2], coords[2], 1e-12 );
        }
    }
}

//! Create a regular MBQUAD_INT^2 element quad mesh with regularly
//! spaced coordinates in the range [1,100].  Group elements
//! into 10 vertical strips MBQUAD_INT/10 elements wide.  Tag elements,
//! vertices and/or sets with ID in [1,10] stored in ID_TAG_NAME
//! tag.  Write new mesh to TEST_FILE.
void create_mesh( bool create_element_sets,
                  bool create_vertex_sets,
                  bool tag_elements_with_id,
                  bool tag_vertices_with_id,
                  const char* adj_elem_tag_name,
                  bool var_len_adj_elems )
{
    Core moab;
    Interface& mb = moab;
    ErrorCode rval;

    // create tags
    Tag logical_tag, centroid_tag, id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );
    rval =
        mb.tag_get_handle( LOGICAL_NAME, 2 * sizeof( int ), MB_TYPE_OPAQUE, logical_tag, MB_TAG_DENSE | MB_TAG_EXCL );CHECK_ERR( rval );
    rval = mb.tag_get_handle( CENTROID_NAME, 3, MB_TYPE_DOUBLE, centroid_tag, MB_TAG_DENSE | MB_TAG_EXCL );CHECK_ERR( rval );

    EntityHandle sets[NUM_SETS];
    if( create_element_sets || create_vertex_sets )
    {
        for( int i = 0; i < NUM_SETS; ++i )
        {
            rval = mb.create_meshset( MESHSET_ORDERED, sets[i] );CHECK_ERR( rval );
            int id = i + 1;
            rval   = mb.tag_set_data( id_tag, &sets[i], 1, &id );CHECK_ERR( rval );
        }
    }

    // create elements
    EntityHandle verts[MBQUAD_INT + 1][MBQUAD_INT + 1], quads[MBQUAD_INT][MBQUAD_INT];
    for( int i = 0; i <= MBQUAD_INT; ++i )
        for( int j = 0; j <= MBQUAD_INT; ++j )
        {
            double coords[3] = { static_cast< double >( i ), static_cast< double >( j ), 0 };
            rval             = mb.create_vertex( coords, verts[j][i] );CHECK_ERR( rval );
            int logical[2] = { i, j };
            rval           = mb.tag_set_data( logical_tag, &verts[j][i], 1, logical );CHECK_ERR( rval );
            rval = mb.tag_set_data( centroid_tag, &verts[j][i], 1, coords );CHECK_ERR( rval );
            int id = ( i - 1 ) / SET_WIDTH + 1;  // Note: assumes SET_WIDTH > 1
            if( tag_vertices_with_id )
            {
                rval = mb.tag_set_data( id_tag, &verts[j][i], 1, &id );CHECK_ERR( rval );
            }
            if( create_vertex_sets )
            {
                rval = mb.add_entities( sets[id - 1], &verts[j][i], 1 );CHECK_ERR( rval );
                // Some vertices are shared by quads in different sets.
                // put such vertices in both sets.
                int id2 = i / SET_WIDTH + 1;
                if( id2 != id && id2 <= NUM_SETS )
                {
                    rval = mb.add_entities( sets[id2 - 1], &verts[j][i], 1 );CHECK_ERR( rval );
                }
            }
        }
    for( int i = 0; i < MBQUAD_INT; ++i )
        for( int j = 0; j < MBQUAD_INT; ++j )
        {
            EntityHandle conn[4] = { verts[j][i], verts[j][i + 1], verts[j + 1][i + 1], verts[j + 1][i] };
            rval                 = mb.create_element( MBQUAD, conn, 4, quads[j][i] );CHECK_ERR( rval );
            int logical[2] = { i, j };
            rval           = mb.tag_set_data( logical_tag, &quads[j][i], 1, logical );CHECK_ERR( rval );
            double centroid[3] = { i + 0.5, j + 0.5, 0 };
            rval               = mb.tag_set_data( centroid_tag, &quads[j][i], 1, centroid );CHECK_ERR( rval );
            int id = i / SET_WIDTH + 1;
            if( tag_elements_with_id )
            {
                rval = mb.tag_set_data( id_tag, &quads[j][i], 1, &id );CHECK_ERR( rval );
            }
            if( create_element_sets )
            {
                rval = mb.add_entities( sets[id - 1], &quads[j][i], 1 );CHECK_ERR( rval );
            }
        }

    if( adj_elem_tag_name && !var_len_adj_elems )
    {
        Tag handle_tag;
        rval = mb.tag_get_handle( adj_elem_tag_name, 4, MB_TYPE_HANDLE, handle_tag, MB_TAG_DENSE | MB_TAG_EXCL );CHECK_ERR( rval );
        for( int i = 0; i <= MBQUAD_INT; ++i )
            for( int j = 0; j <= MBQUAD_INT; ++j )
            {
                EntityHandle val[4] = { ( i > 0 && j > 0 ) ? quads[j - 1][i - 1] : 0,
                                        ( i > 0 && j < MBQUAD_INT ) ? quads[j][i - 1] : 0,
                                        ( i < MBQUAD_INT && j < MBQUAD_INT ) ? quads[j][i] : 0,
                                        ( i < MBQUAD_INT && j > 0 ) ? quads[j - 1][i] : 0 };
                rval = mb.tag_set_data( handle_tag, &verts[j][i], 1, val );CHECK_ERR( rval );
            }
    }
    else if( adj_elem_tag_name && var_len_adj_elems )
    {
        Tag handle_tag;
        rval = mb.tag_get_handle( adj_elem_tag_name, 0, MB_TYPE_HANDLE, handle_tag,
                                  MB_TAG_DENSE | MB_TAG_VARLEN | MB_TAG_EXCL );CHECK_ERR( rval );
        for( int i = 0; i <= MBQUAD_INT; ++i )
            for( int j = 0; j <= MBQUAD_INT; ++j )
            {
                EntityHandle val[4];
                int num = 0;
                if( i > 0 && j > 0 ) val[num++] = quads[j - 1][i - 1];
                if( i > 0 && j < MBQUAD_INT ) val[num++] = quads[j][i - 1];
                if( i < MBQUAD_INT && j < MBQUAD_INT ) val[num++] = quads[j][i];
                if( i < MBQUAD_INT && j > 0 ) val[num++] = quads[j - 1][i];
                const void* ptr = val;
                rval            = mb.tag_set_by_ptr( handle_tag, &verts[j][i], 1, &ptr, &num );CHECK_ERR( rval );
            }
    }

    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );
}

// Given a list of vertices adjacent to a quad strip, identify it as one of the
// NUM_SETS strips of quads written by create_mesh.
int identify_set( Interface& mb, const Range& verts )
{
    const int COL = SET_WIDTH + 1;
    CHECK_EQUAL( ( 1 + MBQUAD_INT ) * COL, (int)verts.size() );

    // Get X range of vertices
    int min_x = std::numeric_limits< int >::max();
    int max_x = std::numeric_limits< int >::min();
    for( Range::const_iterator i = verts.begin(); i != verts.end(); ++i )
    {
        double coords[3];
        ErrorCode rval = mb.get_coords( &*i, 1, coords );CHECK_ERR( rval );
        // Expect whole-valued coorindates
        int int_x = (int)coords[0];
        CHECK( fabs( coords[0] - (double)int_x ) < 1e-12 );

        if( int_x < min_x ) min_x = int_x;
        if( int_x > max_x ) max_x = int_x;
    }
    CHECK( max_x - min_x < COL );

    // Calculate ID (return value) from coordinate range)
    const int ID = min_x / SET_WIDTH + 1;

    // Now verify that all vertices correctly form a grid
    EntityHandle grid[MBQUAD_INT + 1][COL];
    memset( grid, 0, sizeof( grid ) );
    for( Range::const_iterator i = verts.begin(); i != verts.end(); ++i )
    {
        double coords[3];
        ErrorCode rval = mb.get_coords( &*i, 1, coords );CHECK_ERR( rval );
        // Expect whole-valued coorindates
        int x = (int)coords[0] - ( ID - 1 ) * SET_WIDTH, y = (int)coords[1];
        CHECK( fabs( coords[1] - (double)y ) < 1e-12 );
        CHECK( fabs( coords[2] ) < 1e-12 );
        CHECK( y >= 0 && y <= MBQUAD_INT );
        CHECK_EQUAL( (EntityHandle)0, grid[y][x] );
        grid[y][x] = *i;
    }

    return ID;
}
int identify_set( Interface& mb, EntityHandle set )
{
    ErrorCode rval;
    Range verts, elems;
    rval = mb.get_entities_by_handle( set, elems );CHECK_ERR( rval );
    Range::iterator it = elems.upper_bound( MBVERTEX );
    verts.merge( elems.begin(), it );
    elems.erase( elems.begin(), it );
    it = elems.lower_bound( MBENTITYSET );
    elems.erase( it, elems.end() );
    rval = mb.get_adjacencies( elems, 0, false, verts, Interface::UNION );CHECK_ERR( rval );
    return identify_set( mb, verts );
}

//! Read in the elems contained in a set
void test_read_one_set_elems()
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    create_mesh( true, false, false, false );

    for( int id = 1; id <= NUM_SETS; ++id )
    {
        rval = mb.delete_mesh();CHECK_ERR( rval );
        rval = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, &id, 1 );CHECK_ERR( rval );
        Range verts;
        rval       = mb.get_entities_by_type( 0, MBVERTEX, verts );
        int act_id = identify_set( mb, verts );
        CHECK_EQUAL( id, act_id );
    }
}

//! Read in the elems contained in a sets
void test_read_two_sets_elems()
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    create_mesh( true, false, false, false );
    int ids[2] = { 2, 8 };
    EntityHandle file_set;
    rval = mb.create_meshset( MESHSET_SET, file_set );CHECK_ERR( rval );
    rval = mb.load_file( TEST_FILE, &file_set, READ_OPTS, ID_TAG_NAME, ids, 2 );CHECK_ERR( rval );

    Range sets;
    rval = mb.get_entities_by_type( 0, MBENTITYSET, sets );CHECK_ERR( rval );
    CHECK_EQUAL( 3, (int)sets.size() );
    Range::iterator it = sets.find( file_set );
    CHECK( it != sets.end() );
    sets.erase( it );

    int id1 = identify_set( mb, sets.front() );
    int id2 = identify_set( mb, sets.back() );
    if( id1 == ids[0] )
    {
        CHECK_EQUAL( ids[1], id2 );
    }
    else
    {
        CHECK_EQUAL( ids[1], id1 );
        CHECK_EQUAL( ids[0], id2 );
    }
}

Tag check_tag( Interface& mb, const char* name, TagType storage, DataType type, int size )
{

    Tag tag;
    ErrorCode rval = mb.tag_get_handle( name, size, type, tag );CHECK_ERR( rval );

    TagType storage1;
    rval = mb.tag_get_type( tag, storage1 );CHECK_ERR( rval );
    CHECK_EQUAL( storage, storage1 );

    DataType type1;
    rval = mb.tag_get_data_type( tag, type1 );CHECK_ERR( rval );
    CHECK_EQUAL( type, type1 );

    int size1;
    rval = mb.tag_get_length( tag, size1 );
    if( size <= 0 )
    {  // variable-length tag
        CHECK_EQUAL( MB_VARIABLE_DATA_LENGTH, rval );
    }
    else
    {
        CHECK_ERR( rval );
        CHECK_EQUAL( size, size1 );
    }

    return tag;
}

//! Test reading of sparse double tag data
void test_read_double_tag()
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    create_mesh( true, false, false, false );
    int ids[2] = { 1, 4 };
    rval       = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, 2 );CHECK_ERR( rval );

    Tag tag = check_tag( mb, CENTROID_NAME, MB_TAG_DENSE, MB_TYPE_DOUBLE, 3 );
    Range verts;
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    CHECK( !verts.empty() );
    for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
    {
        double coords[3], data[3];
        rval = mb.get_coords( &*i, 1, coords );CHECK_ERR( rval );
        rval = mb.tag_get_data( tag, &*i, 1, data );CHECK_ERR( rval );
        CHECK_REAL_EQUAL( coords[0], data[0], 1e-12 );
        CHECK_REAL_EQUAL( coords[1], data[1], 1e-12 );
        CHECK_REAL_EQUAL( coords[2], data[2], 1e-12 );
    }
}

//! Test reading of sparse opaque tag data
void test_read_opaque_tag()
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    create_mesh( true, false, false, false );
    int ids[2] = { 1, 4 };
    rval       = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, 2 );CHECK_ERR( rval );

    Tag tag = check_tag( mb, LOGICAL_NAME, MB_TAG_DENSE, MB_TYPE_OPAQUE, 2 * sizeof( int ) );
    Range verts;
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    CHECK( !verts.empty() );
    for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
    {
        double coords[3];
        int data[2];
        rval = mb.get_coords( &*i, 1, coords );CHECK_ERR( rval );
        rval = mb.tag_get_data( tag, &*i, 1, data );CHECK_ERR( rval );
        CHECK_REAL_EQUAL( coords[0], (double)data[0], 1e-12 );
        CHECK_REAL_EQUAL( coords[1], (double)data[1], 1e-12 );
    }
}

static void test_read_handle_tag_common( bool var_len )
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    const char tag_name[] = "VTX_ADJ";
    create_mesh( true, false, false, false, tag_name, var_len );
    int ids[2] = { 7, 10 };
    rval       = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, 2 );CHECK_ERR( rval );

    Tag tag = check_tag( mb, tag_name, MB_TAG_DENSE, MB_TYPE_HANDLE, var_len ? 0 : 4 );
    Range verts;
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    CHECK( !verts.empty() );
    for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
    {
        std::vector< EntityHandle > adj, val;
        rval = mb.get_adjacencies( &*i, 1, 2, false, adj, Interface::UNION );CHECK_ERR( rval );
        CHECK( !adj.empty() );

        int num;
        const void* ptr;
        rval = mb.tag_get_by_ptr( tag, &*i, 1, &ptr, &num );CHECK_ERR( rval );

        if( var_len )
        {
            CHECK( num > 0 );
            CHECK( num < 5 );
        }
        else
        {
            CHECK_EQUAL( 4, num );
        }

        val.clear();
        const EntityHandle* dat = (const EntityHandle*)ptr;
        for( const EntityHandle* end = dat + num; dat != end; ++dat )
            if( *dat ) val.push_back( *dat );

        CHECK_EQUAL( adj.size(), val.size() );
        std::sort( adj.begin(), adj.end() );
        std::sort( val.begin(), val.end() );
        CHECK( adj == val );
    }
}

void test_read_tagged_elems()
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    create_mesh( false, false, true, false );
    int id = 5;
    rval   = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, &id, 1 );CHECK_ERR( rval );

    Range verts;
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    int id2 = identify_set( mb, verts );
    CHECK_EQUAL( id, id2 );

    int elems;
    rval = mb.get_number_entities_by_type( 0, MBQUAD, elems );CHECK_ERR( rval );
    CHECK_EQUAL( MBQUAD_INT * MBQUAD_INT / NUM_SETS, elems );
}

void test_read_tagged_nodes()
{
    ErrorCode rval;
    Core moab;
    Interface& mb = moab;

    create_mesh( false, false, false, true );
    int id = 1;  // NOTE: this test will only succeed for ID == 1
    rval   = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, &id, 1 );CHECK_ERR( rval );

    Range verts;
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    int id2 = identify_set( mb, verts );
    CHECK_EQUAL( id, id2 );

    int elems;
    rval = mb.get_number_entities_by_type( 0, MBQUAD, elems );CHECK_ERR( rval );
    CHECK_EQUAL( MBQUAD_INT * MBQUAD_INT / NUM_SETS, elems );
}

//! Read in the polyhedra contained in a set
void test_read_one_set_polyhedra()
{
    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    // create a 2x2x1 block of hexes, splitting each hex face
    // into two triangles to form an 12-sided polyhedron
    EntityHandle verts[18], hexes[4];
    double coords[18][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 2, 0, 0 }, { 0, 1, 0 }, { 1, 1, 0 }, { 2, 1, 0 },
                             { 0, 0, 1 }, { 1, 0, 1 }, { 2, 0, 1 }, { 0, 1, 1 }, { 1, 1, 1 }, { 2, 1, 1 },
                             { 0, 0, 2 }, { 1, 0, 2 }, { 2, 0, 2 }, { 0, 1, 2 }, { 1, 1, 2 }, { 2, 1, 2 } };
    int hexconn[4][8]    = { { 0, 1, 4, 3, 6, 7, 10, 9 },
                          { 1, 2, 5, 4, 7, 8, 11, 10 },
                          { 6, 7, 10, 9, 12, 13, 16, 15 },
                          { 7, 8, 11, 10, 13, 14, 17, 16 } };
    for( int i = 0; i < 18; ++i )
    {
        rval = mb.create_vertex( coords[i], verts[i] );CHECK_ERR( rval );
    }
    for( int i = 0; i < 4; ++i )
    {
        EntityHandle conn[8];
        for( int j = 0; j < 8; ++j )
            conn[j] = verts[hexconn[i][j]];
        rval = mb.create_element( MBHEX, conn, 8, hexes[i] );CHECK_ERR( rval );
    }

    Tag tri_tag;
    rval = mb.tag_get_handle( "tris", 2, MB_TYPE_HANDLE, tri_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );

    std::vector< EntityHandle > quads;
    EntityHandle tris[12], poly[4];
    for( int i = 0; i < 4; ++i )
    {
        quads.clear();
        rval = mb.get_adjacencies( &hexes[i], 1, 2, true, quads );CHECK_ERR( rval );
        CHECK_EQUAL( 6, (int)quads.size() );

        for( int j = 0; j < 6; ++j )
        {
            rval = mb.tag_get_data( tri_tag, &quads[j], 1, tris + 2 * j );
            if( MB_SUCCESS == rval ) continue;
            CHECK_EQUAL( MB_TAG_NOT_FOUND, rval );
            const EntityHandle* conn;
            int len;
            rval = mb.get_connectivity( quads[j], conn, len );CHECK_ERR( rval );
            CHECK_EQUAL( 4, len );
            EntityHandle tri_conn[2][3] = { { conn[0], conn[1], conn[2] }, { conn[2], conn[3], conn[0] } };
            rval                        = mb.create_element( MBTRI, tri_conn[0], 3, tris[2 * j] );CHECK_ERR( rval );
            rval = mb.create_element( MBTRI, tri_conn[1], 3, tris[2 * j + 1] );CHECK_ERR( rval );
            rval = mb.tag_set_data( tri_tag, &quads[j], 1, tris + 2 * j );CHECK_ERR( rval );
        }

        rval = mb.create_element( MBPOLYHEDRON, tris, 12, poly[i] );CHECK_ERR( rval );
    }

    Range all_tri;
    rval = mb.get_entities_by_type( 0, MBTRI, all_tri );CHECK_ERR( rval );
    CHECK_EQUAL( 40, (int)all_tri.size() );

    rval = mb.delete_entities( hexes, 4 );CHECK_ERR( rval );
    rval = mb.delete_entities( &quads[0], quads.size() );CHECK_ERR( rval );

    EntityHandle sets[2];
    rval = mb.create_meshset( 0, sets[0] );CHECK_ERR( rval );
    rval = mb.add_entities( sets[0], poly, 2 );CHECK_ERR( rval );
    rval = mb.create_meshset( 0, sets[1] );CHECK_ERR( rval );
    rval = mb.add_entities( sets[1], poly + 2, 2 );CHECK_ERR( rval );

    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );
    int ids[2] = { 2, 3 };
    rval       = mb.tag_set_data( id_tag, sets, 2, ids );CHECK_ERR( rval );

    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );
    rval = mb.delete_mesh();CHECK_ERR( rval );

    rval = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, 1 );CHECK_ERR( rval );

    Range rpoly;
    rval = mb.get_entities_by_type( 0, MBPOLYHEDRON, rpoly );CHECK_ERR( rval );
    CHECK_EQUAL( 2, (int)rpoly.size() );

    Range polyverts;
    rval = mb.get_adjacencies( rpoly, 0, false, polyverts, Interface::UNION );CHECK_ERR( rval );
    CHECK_EQUAL( 12, (int)polyverts.size() );

    for( Range::iterator it = polyverts.begin(); it != polyverts.end(); ++it )
    {
        double coords2[3];
        rval = mb.get_coords( &*it, 1, coords2 );CHECK_ERR( rval );
        CHECK( coords2[0] > -1e-12 && coords2[0] - 2 < 1e-12 );
        CHECK( coords2[1] > -1e-12 && coords2[1] - 1 < 1e-12 );
        CHECK( coords2[2] > -1e-12 && coords2[2] - 1 < 1e-12 );
    }
}

//! Read in the sets contained in a set.
//! Should read all sets containing read elements or nodes
//! and all sets that are contained the the specified "read"
//! set.  Test the later here.
void test_read_set_sets()
{
    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );

    // create sets and assign an ID to each
    const int len = 5;
    EntityHandle set[2 * len + 2];
    for( int i = 0; i < 2 * len + 2; ++i )
    {
        rval = mb.create_meshset( MESHSET_SET, set[i] );CHECK_ERR( rval );
        int id = i + 1;
        rval   = mb.tag_set_data( id_tag, set + i, 1, &id );CHECK_ERR( rval );
    }

    // make set containment as follows (values are assigned IDs):
    int cont_ids[2][len] = { { 3, 4, 5, 9, 10 }, { 6, 7, 8, 11, 12 } };
    for( int i = 0; i < 2; ++i )
    {
        EntityHandle contents[len] = { set[cont_ids[i][0] - 1], set[cont_ids[i][1] - 1], set[cont_ids[i][2] - 1],
                                       set[cont_ids[i][3] - 1], set[cont_ids[i][4] - 1] };
        rval                       = mb.add_entities( set[i], contents, len );CHECK_ERR( rval );
    }

    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );

    for( int i = 0; i < 2; ++i )
    {
        rval = mb.delete_mesh();CHECK_ERR( rval );

        EntityHandle file;
        rval = mb.create_meshset( MESHSET_SET, file );CHECK_ERR( rval );
        int id = i + 1;
        rval   = mb.load_file( TEST_FILE, &file, READ_OPTS ";SETS=NONE", ID_TAG_NAME, &id, 1 );CHECK_ERR( rval );

        // check that the total number of sets read is as expected
        Range sets;
        rval = mb.get_entities_by_type( 0, MBENTITYSET, sets );CHECK_ERR( rval );
        Range::iterator it = sets.find( file );
        if( it != sets.end() ) sets.erase( it );
        CHECK_EQUAL( len + 1, (int)sets.size() );

        // check that we read in the set specified by ID to the reader
        rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag );CHECK_ERR( rval );
        sets.clear();
        const void* data[] = { &id };
        rval               = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &id_tag, data, 1, sets );CHECK_ERR( rval );
        CHECK_EQUAL( 1, (int)sets.size() );

        // check that it contains the expected sets
        EntityHandle owner = sets.front();
        sets.clear();
        rval = mb.get_entities_by_type( owner, MBENTITYSET, sets );CHECK_ERR( rval );
        CHECK_EQUAL( len, (int)sets.size() );

        std::vector< int > expected( cont_ids[i], cont_ids[i] + len );
        std::vector< int > actual( len );
        rval = mb.tag_get_data( id_tag, sets, &actual[0] );CHECK_ERR( rval );
        std::sort( expected.begin(), expected.end() );
        std::sort( actual.begin(), actual.end() );
        CHECK( expected == actual );
    }
}

static void check_children( bool contents, GatherTestMode mode, Interface& mb, int id, Tag id_tag, EntityHandle file )
{
    // Increase number of expected sets by one if contents is true because
    // we always read immediately contained (depth 1) sets.
    const int exp_num_sets  = ( mode == GATHER_NONE ) ? 1 + contents : id;
    const int exp_num_edges = ( mode == GATHER_CONTENTS ) ? id : 1;

    ErrorCode rval;
    Range range;
    rval = mb.get_entities_by_type( 0, MBEDGE, range );CHECK_ERR( rval );
    CHECK_EQUAL( exp_num_edges, (int)range.size() );
    range.clear();
    rval = mb.get_entities_by_type( 0, MBENTITYSET, range );CHECK_ERR( rval );
    Range::iterator it = range.find( file );
    CHECK( it != range.end() );
    range.erase( it );
    CHECK_EQUAL( exp_num_sets, (int)range.size() );

    EntityHandle set;
    const void* val[] = { &id };
    range.clear();
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &id_tag, val, 1, range );CHECK_ERR( rval );
    CHECK_EQUAL( 1, (int)range.size() );
    set = range.front();

    if( mode == GATHER_NONE )
    {
        range.clear();
        rval = mb.get_entities_by_type( set, MBEDGE, range );CHECK_ERR( rval );
        CHECK_EQUAL( 1, (int)range.size() );
        return;
    }

    for( int i = id; i > 0; --i )
    {
        int act_id;
        rval = mb.tag_get_data( id_tag, &set, 1, &act_id );CHECK_ERR( rval );
        CHECK_EQUAL( i, act_id );

        range.clear();
        rval = mb.get_entities_by_type( set, MBEDGE, range );CHECK_ERR( rval );
        if( mode == GATHER_CONTENTS || i == id )
        {
            CHECK_EQUAL( 1, (int)range.size() );
            const EntityHandle* conn;
            int len;
            rval = mb.get_connectivity( range.front(), conn, len );CHECK_ERR( rval );
            CHECK_EQUAL( 2, len );
            double coords[3];
            rval = mb.get_coords( conn + 1, 1, coords );CHECK_ERR( rval );
            CHECK_EQUAL( i, (int)coords[0] );
        }
        else
        {
            CHECK( range.empty() );
        }

        std::vector< EntityHandle > children;
        if( contents )
            rval = mb.get_entities_by_type( set, MBENTITYSET, children );
        else
            rval = mb.get_child_meshsets( set, children );CHECK_ERR( rval );
        if( i == 1 )
        {
            CHECK( children.empty() );
        }
        else
        {
            CHECK_EQUAL( 1, (int)children.size() );
            set = children[0];
        }
    }
}

const char* set_read_opts[] = { "SETS", "CONTENTS", "NONE" };
void test_gather_sets_common( bool contents, GatherTestMode mode, bool no_parent_containing_sets )
{
    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );

    // Create a string of edges from [0,INT] along the X axis each 1 unit in length.
    // Create a set for edge edge containing the edge and make it the parent of the
    // set containing the previous (closer to origin) edge.  Assign each set an
    // ID that is the X coordinate of the larger of the two vertices of the edge
    // contained in the set.
    const int INT = 64;
    EntityHandle verts[INT + 1], edges[INT], sets[INT];
    double coords[] = { 0, 0, 0 };
    rval            = mb.create_vertex( coords, verts[0] );CHECK_ERR( rval );
    for( int i = 0; i < INT; ++i )
    {
        const int id = i + 1;
        coords[0]    = id;
        rval         = mb.create_vertex( coords, verts[id] );CHECK_ERR( rval );
        rval = mb.create_element( MBEDGE, verts + i, 2, edges[i] );CHECK_ERR( rval );
        rval = mb.create_meshset( MESHSET_SET, sets[i] );CHECK_ERR( rval );
        rval = mb.add_entities( sets[i], edges + i, 1 );CHECK_ERR( rval );
        rval = mb.tag_set_data( id_tag, sets + i, 1, &id );CHECK_ERR( rval );
        if( i > 0 )
        {
            if( contents )
                rval = mb.add_entities( sets[i], sets + ( i - 1 ), 1 );
            else
                rval = mb.add_child_meshset( sets[i], sets[i - 1] );CHECK_ERR( rval );
        }
    }

    // Write the data
    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );

    EntityHandle file;
    std::string opt( READ_OPTS );
    if( contents )
        opt += ";CHILDREN=NONE;SETS=";
    else
        opt += ";SETS=NONE;CHILDREN=";
    opt += set_read_opts[mode];

    if( no_parent_containing_sets ) opt += ";NO_SET_CONTAINING_PARENTS";

    const int test_ids[]   = { 2, 7, INT / 3 - 1, INT / 2 + 1, INT - 3 };
    const int num_test_ids = sizeof( test_ids ) / sizeof( int );
    for( int i = 0; i < num_test_ids; ++i )
    {
        CHECK( test_ids[i] <= INT );

        rval = mb.delete_mesh();CHECK_ERR( rval );

        rval = mb.create_meshset( MESHSET_SET, file );CHECK_ERR( rval );

        rval = mb.load_file( TEST_FILE, 0, opt.c_str(), ID_TAG_NAME, test_ids + i, 1 );CHECK_ERR( rval );
        rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag );CHECK_ERR( rval );

        check_children( contents, mode, mb, test_ids[i], id_tag, file );
    }
}

void test_gather_sets_ranged( bool contents, GatherTestMode mode, bool no_parent_containing_sets )
{
    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    Range verts;
    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );

    // create four groups of vertices, where all vertices in the same group
    // have the same x-coordinate
    const int NUM_GRP_VTX = 20;
    const int NUM_GRP     = 4;
    EntityHandle sets[NUM_GRP];
    for( int i = 0; i < NUM_GRP; ++i )
    {
        double coords[3 * NUM_GRP_VTX];
        for( int j = 0; j < NUM_GRP_VTX; ++j )
        {
            coords[3 * j]     = i;
            coords[3 * j + 1] = j;
            coords[3 * j + 2] = 0;
        }
        rval = mb.create_vertices( coords, NUM_GRP_VTX, verts );CHECK_ERR( rval );

        rval = mb.create_meshset( MESHSET_SET, sets[i] );CHECK_ERR( rval );
        rval = mb.add_entities( sets[i], verts );CHECK_ERR( rval );
        int id = i + 1;
        rval   = mb.tag_set_data( id_tag, sets + i, 1, &id );CHECK_ERR( rval );
    }

    // place two of the sets inside the others
    if( contents )
    {
        rval = mb.add_entities( sets[0], &sets[1], 1 );CHECK_ERR( rval );
        rval = mb.add_entities( sets[2], &sets[3], 1 );CHECK_ERR( rval );
    }
    else
    {
        rval = mb.add_child_meshset( sets[0], sets[1] );CHECK_ERR( rval );
        rval = mb.add_child_meshset( sets[2], sets[3] );CHECK_ERR( rval );
    }

    // Write the data
    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );

    // Read the data
    std::string opt( READ_OPTS );
    if( contents )
        opt += ";CHILDREN=NONE;SETS=";
    else
        opt += ";SETS=NONE;CHILDREN=";
    opt += set_read_opts[mode];

    if( no_parent_containing_sets ) opt += ";NO_PARENT_CONTAINING_SETS";

    EntityHandle file;
    const int read_id = 3;
    rval              = mb.delete_mesh();CHECK_ERR( rval );
    rval = mb.create_meshset( MESHSET_SET, file );CHECK_ERR( rval );
    rval = mb.load_file( TEST_FILE, &file, opt.c_str(), ID_TAG_NAME, &read_id, 1 );CHECK_ERR( rval );

    // get any sets that were read it
    Range read_sets;
    rval = mb.get_entities_by_type( file, MBENTITYSET, read_sets );CHECK_ERR( rval );

    // count number of vertices in each group
    int counts[NUM_GRP];
    memset( counts, 0, sizeof( counts ) );
    verts.clear();
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    for( Range::iterator it = verts.begin(); it != verts.end(); ++it )
    {
        double coords[3];
        rval = mb.get_coords( &*it, 1, coords );CHECK_ERR( rval );
        int i = (int)( coords[0] + 1e-12 );
        CHECK( i >= 0 && i < NUM_GRP );
        counts[i]++;
    }

    // check expected counts
    CHECK_EQUAL( 0, counts[0] );
    CHECK_EQUAL( 0, counts[1] );
    CHECK_EQUAL( NUM_GRP_VTX, counts[2] );
    switch( mode )
    {
        case GATHER_NONE:
            CHECK_EQUAL( 0, counts[3] );
            CHECK_EQUAL( 1 + contents, (int)read_sets.size() );
            break;
        case GATHER_SETS:
            CHECK_EQUAL( 0, counts[3] );
            CHECK_EQUAL( 2, (int)read_sets.size() );
            break;
        case GATHER_CONTENTS:
            CHECK_EQUAL( NUM_GRP_VTX, counts[3] );
            CHECK_EQUAL( 2, (int)read_sets.size() );
            break;
    }
}

static void check_num_verts( Interface& mb, Tag tag, int id, int num_vtx )
{
    ErrorCode rval;
    const void* val[] = { &id };
    Range range;
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, val, 1, range );CHECK_ERR( rval );
    CHECK_EQUAL( 1, (int)range.size() );

    EntityHandle set = range.front();
    range.clear();
    rval = mb.get_entities_by_type( set, MBVERTEX, range );CHECK_ERR( rval );
    CHECK_EQUAL( num_vtx, (int)range.size() );
}

//! Read in the sets contained in a set.
//! Should read all sets containing read elements or nodes
//! and all sets that are contained the the specified "read"
//! set.  Test the former here.
void test_read_containing_sets()
{
    // create mesh decomposed by elements but create
    // sets containing all vertices of decomposed elements
    // such that adjacent sets share vertices.
    create_mesh( false, true, false, false );

    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    // read some sets
    const int ids[]    = { 1, 5, 9 };
    const int num_sets = sizeof( ids ) / sizeof( int );
    rval               = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, num_sets );CHECK_ERR( rval );

    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag );CHECK_ERR( rval );

    // expect all sets adjacent to the specified sets because
    // they share vertices.
    Range verts;
    for( int i = 0; i < num_sets; ++i )
    {
        if( ids[i] > 1 ) check_num_verts( mb, id_tag, ids[i] - 1, MBQUAD_INT + 1 );
        check_num_verts( mb, id_tag, ids[i], ( MBQUAD_INT + 1 ) * ( SET_WIDTH + 1 ) );
        if( ids[i] < NUM_SETS ) check_num_verts( mb, id_tag, ids[i] + 1, MBQUAD_INT + 1 );
    }
}

//! Test reading of explicit adjacencies
void test_read_adjacencies()
{
    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    // create four hexes sharing an edge
    EntityHandle verts[3][3][2], hexes[2][2];
    for( int k = 0; k < 2; ++k )
    {
        for( int j = 0; j < 3; ++j )
        {
            for( int i = 0; i < 3; ++i )
            {
                double coords[] = { static_cast< double >( i ), static_cast< double >( j ),
                                    static_cast< double >( k ) };
                rval            = mb.create_vertex( coords, verts[i][j][k] );CHECK_ERR( rval );
            }
        }
    }
    for( int j = 0; j < 2; ++j )
    {
        for( int i = 0; i < 2; ++i )
        {
            EntityHandle conn[] = { verts[i][j][0], verts[i + 1][j][0], verts[i + 1][j + 1][0], verts[i][j + 1][0],
                                    verts[i][j][1], verts[i + 1][j][1], verts[i + 1][j + 1][1], verts[i][j + 1][1] };
            rval                = mb.create_element( MBHEX, conn, 8, hexes[i][j] );CHECK_ERR( rval );
        }
    }

    // create two duplicate edges that connect the vertices common to all four hexes
    EntityHandle edge_conn[2] = { verts[1][1][0], verts[1][1][1] };
    EntityHandle edges[2];
    rval = mb.create_element( MBEDGE, edge_conn, 2, edges[0] );CHECK_ERR( rval );
    rval = mb.create_element( MBEDGE, edge_conn, 2, edges[1] );CHECK_ERR( rval );
    // mark one edge as adjacent to the left two hexes and the
    // other as adjacent to the right two
    rval = mb.add_adjacencies( edges[0], hexes[0], 2, true );CHECK_ERR( rval );
    rval = mb.add_adjacencies( edges[1], hexes[1], 2, true );CHECK_ERR( rval );
    // create two sets containing the front two and the rear two
    // hexes, respectively.
    EntityHandle sets[2];
    rval = mb.create_meshset( MESHSET_SET, sets[0] );CHECK_ERR( rval );
    rval = mb.create_meshset( MESHSET_SET, sets[1] );CHECK_ERR( rval );
    EntityHandle set1[4] = { hexes[0][0], hexes[1][0], edges[0], edges[1] };
    EntityHandle set2[4] = { hexes[0][1], hexes[1][1], edges[0], edges[1] };
    rval                 = mb.add_entities( sets[0], set1, 4 );CHECK_ERR( rval );
    rval = mb.add_entities( sets[1], set2, 4 );CHECK_ERR( rval );

    // assign IDs to sets
    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );
    int ids[2] = { 1, 2 };
    rval       = mb.tag_set_data( id_tag, sets, 2, ids );CHECK_ERR( rval );

    // write mesh
    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );

    // read mesh
    rval = mb.delete_mesh();CHECK_ERR( rval );
    rval = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, 1 );CHECK_ERR( rval );

    // expect two hexes and two edges
    Range range;
    rval = mb.get_entities_by_type( 0, MBHEX, range );CHECK_ERR( rval );
    CHECK_EQUAL( 2, (int)range.size() );
    EntityHandle h1 = range.front(), h2 = range.back();
    range.clear();
    rval = mb.get_entities_by_type( 0, MBEDGE, range );CHECK_ERR( rval );
    CHECK_EQUAL( 2, (int)range.size() );

    // expecte each edge to have one of the hexes
    range.clear();
    rval = mb.get_adjacencies( &h1, 1, 1, false, range );CHECK_ERR( rval );
    CHECK_EQUAL( 1, (int)range.size() );
    EntityHandle e1 = range.front();
    range.clear();
    rval = mb.get_adjacencies( &h2, 1, 1, false, range );CHECK_ERR( rval );
    CHECK_EQUAL( 1, (int)range.size() );
    EntityHandle e2 = range.front();

    CHECK( e1 != e2 );
}

void test_read_sides()
{
    ErrorCode rval;
    Core instance;
    Interface& mb = instance;

    // create 4x4 grid of quads with edges
    const int INT = 4;
    EntityHandle verts[INT + 1][INT + 1];
    for( int j = 0; j <= INT; ++j )
    {
        for( int i = 0; i <= INT; ++i )
        {
            double coords[3] = { static_cast< double >( i ), static_cast< double >( j ), 0 };
            rval             = mb.create_vertex( coords, verts[INT - j][i] );CHECK_ERR( rval );
        }
    }
    EntityHandle quads[INT][INT];
    for( int j = 0; j < INT; ++j )
    {
        for( int i = 0; i < INT; ++i )
        {
            EntityHandle conn[4] = { verts[INT - j][i], verts[INT - j][i + 1], verts[INT - j - 1][i + 1],
                                     verts[INT - j - 1][i] };
            rval                 = mb.create_element( MBQUAD, conn, 4, quads[INT - j - 1][i] );CHECK_ERR( rval );
        }
    }
    Range edges;
    rval = mb.get_adjacencies( &quads[0][0], INT * INT, 1, true, edges, Interface::UNION );CHECK_ERR( rval );
    CHECK_EQUAL( 40, (int)edges.size() );

    // group quads into two sets
    EntityHandle sets[2];
    rval = mb.create_meshset( MESHSET_SET, sets[0] );CHECK_ERR( rval );
    rval = mb.create_meshset( MESHSET_SET, sets[1] );CHECK_ERR( rval );
    rval = mb.add_entities( sets[0], quads[0], INT );CHECK_ERR( rval );
    rval = mb.add_entities( sets[1], quads[1], INT );CHECK_ERR( rval );
    rval = mb.add_entities( sets[0], quads[2], INT );CHECK_ERR( rval );
    rval = mb.add_entities( sets[1], quads[3], INT );CHECK_ERR( rval );

    // assign IDS
    Tag id_tag;
    rval = mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag, MB_TAG_SPARSE | MB_TAG_EXCL );CHECK_ERR( rval );
    int ids[2] = { 4, 5 };
    rval       = mb.tag_set_data( id_tag, sets, 2, ids );CHECK_ERR( rval );

    // write mesh
    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );

    // read first set back in
    rval = mb.delete_mesh();CHECK_ERR( rval );
    rval = mb.load_file( TEST_FILE, 0, READ_OPTS, ID_TAG_NAME, ids, 1 );CHECK_ERR( rval );

    // check expected counts
    int count;
    rval = mb.get_number_entities_by_type( 0, MBVERTEX, count );CHECK_ERR( rval );
    CHECK_EQUAL( ( INT + 1 ) * INT, count );
    rval = mb.get_number_entities_by_type( 0, MBQUAD, count );CHECK_ERR( rval );
    CHECK_EQUAL( INT * INT / 2, count );
    rval = mb.get_number_entities_by_type( 0, MBEDGE, count );CHECK_ERR( rval );
    CHECK_EQUAL( 2 * ( INT + 1 ) + INT * INT, count );

    // check edges adjacent to each quad
    Range elems;
    rval = mb.get_entities_by_type( 0, MBQUAD, elems );CHECK_ERR( rval );
    for( Range::iterator it = elems.begin(); it != elems.end(); ++it )
    {
        edges.clear();
        rval = mb.get_adjacencies( &*it, 1, 1, false, edges );CHECK_ERR( rval );
        CHECK_EQUAL( 4, (int)edges.size() );
    }
}

const int expected_ids[]  = { 2, 4, 6, 8, 10, 12, 14, 16, 18 };
const int expected_vols[] = { 3, 7, 10 };

void write_id_test_file()
{
    Core moab;
    Interface& mb = moab;
    ErrorCode rval;

    // create 12 entity sets
    EntityHandle sets[12];
    for( int i = 0; i < 12; ++i )
    {
        rval = mb.create_meshset( MESHSET_SET, sets[i] );CHECK_ERR( rval );
    }

    // create tag handles
    Tag id = 0, gid = 0, dim = 0;
    mb.tag_get_handle( ID_TAG_NAME, 1, MB_TYPE_INTEGER, id, MB_TAG_SPARSE | MB_TAG_EXCL );
    mb.tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim, MB_TAG_SPARSE | MB_TAG_EXCL );
    gid = mb.globalId_tag();

    // set ID tag on first 10 sets
    rval = mb.tag_set_data( id, sets, sizeof( expected_ids ) / sizeof( int ), expected_ids );CHECK_ERR( rval );
    // set geom dim on all sets, only three of them have dim == 3
    int num_vol = sizeof( expected_vols ) / sizeof( int );
    int dims[12], ids[12];
    int v = 0;
    for( int i = 0; i < 12; ++i )
    {
        dims[i] = i % 3 + 1;
        if( dims[i] == 3 )
        {
            if( v < num_vol )
                ids[i] = expected_vols[v++];
            else
                ids[i] = expected_vols[0];
        }
        else
            ids[i] = 100;
    }
    rval = mb.tag_set_data( gid, sets, 12, ids );CHECK_ERR( rval );
    rval = mb.tag_set_data( dim, sets, 12, dims );CHECK_ERR( rval );

    rval = mb.write_file( TEST_FILE, "MOAB" );CHECK_ERR( rval );
}

void test_read_ids()
{
    write_id_test_file();

    Core moab;
    ReadHDF5 reader( &moab );
    FileOptions opts( "" );
    ErrorCode rval;
    std::vector< int > values;
    rval = reader.read_tag_values( TEST_FILE, ID_TAG_NAME, opts, values );
    remove( TEST_FILE );CHECK_ERR( rval );

    std::sort( values.begin(), values.end() );
    std::vector< int > expected( expected_ids, expected_ids + sizeof( expected_ids ) / sizeof( int ) );
    CHECK_EQUAL( expected, values );
}

void test_read_partial_ids()
{
    write_id_test_file();

    const int three                = 3;
    ReaderIface::IDTag vols        = { GEOM_DIMENSION_TAG_NAME, &three, 1 };
    ReaderIface::SubsetList subset = { &vols, 1, 0, 0 };

    Core moab;
    ReadHDF5 reader( &moab );
    FileOptions opts( "" );
    ErrorCode rval;
    std::vector< int > values;
    rval = reader.read_tag_values( TEST_FILE, GLOBAL_ID_TAG_NAME, opts, values, &subset );
    remove( TEST_FILE );CHECK_ERR( rval );

    std::sort( values.begin(), values.end() );
    std::vector< int > expected( expected_ids, expected_ids + sizeof( expected_ids ) / sizeof( int ) );
}