cub_file_test.cpp

Go to the documentation of this file.

#include "TestUtil.hpp"
#include "moab/Core.hpp"
#include "MBTagConventions.hpp"
#include "moab/CN.hpp"
#include "moab/Range.hpp"
#include "moab/GeomTopoTool.hpp"
#include <cmath>
#include <algorithm>

using namespace moab;

/**\brief Input test file: test.cub
 * Cubit 10.2 file.
 * File contains:
 * Two merged 10x10x10 bricks sharing a single surface (surface 6).
 * - Each brick is meshed with 8 5x5x5 hexes.
 * - Each surface is meshed with 4 5x5 quads.
 * - Each curve is meshed with 2 5-unit edges.
 * A single block containing both bricks.
 * Two side sets
 * - sideset 1: surfaces 1 and 7
 * - sideset 2: surfaces 5 and 11
 * Two node sets:
 * - nodeset 1: surfaces 2 and 8
 * - nodeset 2: surfaces 3 and 9
 *
 * Surfaces:
 *                    2          8
 *                   /          /
 *          o----------o----------o
 *         /.      /  /.      /  /|
 *        / .  (5)/  / . (11)/  / |
 *       /  .    L  /  .    L  /  |
 *      o----------o----------o   |
 *      |   .      |   .      |(12)
 *   4--|-> o . . .|. .o. . . | . o
 *      |  . (1)   |  . (9)   |  /
 *      | .   ^    | .   ^    | /
 *      |.    |    |.    |    |/
 *      o----------o----------o
 *            |          |
 *            3          9
 *
 * Curves:
 *
 *     o----8-----o----20----o
 *    /.         /.         /|
 *  12 .       11 .       24 |
 *  /  7       /  5       /  17
 * o----2-----o----14----o   |
 * |   .      |   .      |   |
 * |   o . .6.|. .o. . 18| . o
 * 3  .       1  .      13  /
 * | 9        | 10       | 22
 * |.         |.         |/
 * o----4-----o----16----o
 *
 * Vertices:
 *
 *     8----------5----------13
 *    /.         /.         /|
 *   / .        / .        / |
 *  /  .       /  .       /  |
 * 3----------2----------10  |
 * |   .      |   .      |   |
 * |   7 . . .|. .6. . . | . 14
 * |  .       |  .       |  /
 * | .        | .        | /
 * |.         |.         |/
 * 4----------1----------9
 */

/* Input test file: ho_test.cub
 *
 * File is expected to contain at least one block for every
 * supported higher-order element type.  The coordinates of
 * every higher-order node are expected to be the mean of the
 * adjacent corner vertices of the element.
 */
static const std::string input_file_1 = std::string( TestDir + "unittest/io/test.cub" );
static const std::string ho_file      = std::string( TestDir + "unittest/io/ho_test.cub" );
static const std::string cubit12_file = std::string( TestDir + "unittest/io/cubtest12.cub" );
static const std::string cubit14_file = std::string( TestDir + "unittest/io/cubtest14.cub" );

void read_file( Interface& moab, const std::string& input_file );

// Check that adjacent lower-order entities have
// higher-order nodes consitent with input entity.
void check_adj_ho_nodes( Interface& moab, EntityHandle entity );

// Check that element has expected higher-order nodes
// and that each higher-order node is at the center
// of the sub-entity it is on.
void check_ho_element( Interface& moab, EntityHandle entity, int mid_nodes[4] );

// Validate elements of specified type.
// Looks for a block containing the specified entity type
// and with the specified mid-node flags set in its
// HAS_MID_NODES_TAG.
void test_ho_elements( EntityType type, int num_nodes );

void test_vertices();

void test_edges();

void test_quads();

void test_hexes();

void test_geometric_topology();

void test_geometric_sets();

void test_blocks();

void test_side_sets();

void test_node_sets();

void test_tri6()
{
    test_ho_elements( MBTRI, 6 );
}
void test_tri7()
{
    test_ho_elements( MBTRI, 7 );
}

void test_quad5()
{
    test_ho_elements( MBQUAD, 5 );
}
void test_quad8()
{
    test_ho_elements( MBQUAD, 8 );
}
void test_quad9()
{
    test_ho_elements( MBQUAD, 9 );
}

void test_tet8()
{
    test_ho_elements( MBTET, 8 );
}
void test_tet10()
{
    test_ho_elements( MBTET, 10 );
}
void test_tet14()
{
    test_ho_elements( MBTET, 14 );
}

void test_hex9()
{
    test_ho_elements( MBHEX, 9 );
}
void test_hex20()
{
    test_ho_elements( MBHEX, 20 );
}
void test_hex27()
{
    test_ho_elements( MBHEX, 27 );
}

void test_multiple_files();

void test_cubit12();
void test_cubit14();

int main()
{
    int result = 0;

    result += RUN_TEST( test_vertices );
    result += RUN_TEST( test_edges );
    result += RUN_TEST( test_quads );
    result += RUN_TEST( test_hexes );
    result += RUN_TEST( test_geometric_topology );
    result += RUN_TEST( test_geometric_sets );
    result += RUN_TEST( test_blocks );
    result += RUN_TEST( test_side_sets );
    result += RUN_TEST( test_node_sets );
    result += RUN_TEST( test_tri6 );
    result += RUN_TEST( test_tri7 );
    result += RUN_TEST( test_quad5 );
    result += RUN_TEST( test_quad8 );
    result += RUN_TEST( test_quad9 );
    result += RUN_TEST( test_tet8 );
    result += RUN_TEST( test_tet10 );
    result += RUN_TEST( test_tet14 );
    result += RUN_TEST( test_hex9 );
    result += RUN_TEST( test_hex20 );
    result += RUN_TEST( test_hex27 );
    result += RUN_TEST( test_multiple_files );
    result += RUN_TEST( test_cubit12 );
    result += RUN_TEST( test_cubit14 );
    return result;
}

void read_file( Interface& moab, const std::string& input_file )
{
    ErrorCode rval = moab.load_file( input_file.c_str() );CHECK_ERR( rval );
}

void test_vertices()
{
    /* Node coordinates, in order of node ID, beginning with 1. */
    const double node_coords[] = {
        5,  -5, 5,                                                   // 1
        5,  5,  5,  5,  0,  5,  -5, 5,  5,  0,  5,  5,               // 5
        -5, -5, 5,  -5, 0,  5,  0,  -5, 5,  0,  0,  5,  5,  5,  -5,  // 10
        5,  -5, -5, 5,  0,  -5, -5, -5, -5, 0,  -5, -5, -5, 5,  -5,  // 15
        -5, 0,  -5, 0,  5,  -5, 0,  0,  -5, -5, -5, 0,  5,  -5, 0,   // 20
        0,  -5, 0,  -5, 5,  0,  -5, 0,  0,  5,  5,  0,  0,  5,  0,   // 25
        5,  0,  0,  0,  0,  0,  15, -5, 5,  15, 5,  5,  15, 0,  5,   // 30
        10, 5,  5,  10, -5, 5,  10, 0,  5,  15, 5,  -5, 15, -5, -5,  // 35
        15, 0,  -5, 10, -5, -5, 10, 5,  -5, 10, 0,  -5, 15, -5, 0,   // 40
        10, -5, 0,  15, 5,  0,  10, 5,  0,  15, 0,  0,  10, 0,  0    // 45
    };

    ErrorCode rval;
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, input_file_1 );

    // get vertex handles and check correct number of vertices
    const size_t num_nodes = sizeof( node_coords ) / ( 3 * sizeof( double ) );
    Range verts;
    rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
    CHECK_EQUAL( num_nodes, (size_t)verts.size() );

    // check global ids (should be 1 to 45 for vertices.)
    Tag gid_tag;
    rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
    std::vector< int > ids( num_nodes );
    rval = mb.tag_get_data( gid_tag, verts, &ids[0] );CHECK_ERR( rval );
    std::vector< int > sorted( ids );
    std::sort( sorted.begin(), sorted.end() );
    for( size_t i = 0; i < num_nodes; ++i )
        CHECK_EQUAL( (int)( i + 1 ), sorted[i] );

    // check coordinates of each vertex
    std::vector< double > coords( 3 * num_nodes );
    rval = mb.get_coords( verts, &coords[0] );CHECK_ERR( rval );
    for( size_t i = 0; i < num_nodes; ++i )
    {
        const double* exp = node_coords + 3 * ( ids[i] - 1 );
        const double* act = &coords[3 * i];
        CHECK_REAL_EQUAL( exp[0], act[0], 1e-8 );
        CHECK_REAL_EQUAL( exp[1], act[1], 1e-8 );
        CHECK_REAL_EQUAL( exp[2], act[2], 1e-8 );
    }
}

void test_element( const std::string& filename, EntityType type, int num_elem, int node_per_elem, const int* conn_list )
{
    ErrorCode rval;
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, filename );

    Range elems;
    rval = mb.get_entities_by_type( 0, type, elems );CHECK_ERR( rval );
    CHECK_EQUAL( num_elem, (int)elems.size() );

    // get global ids
    Tag gid_tag;
    rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
    std::vector< int > ids( num_elem );
    rval = mb.tag_get_data( gid_tag, elems, &ids[0] );CHECK_ERR( rval );

    // check that global ids are consecutive, beginning with 1
    std::vector< int > sorted( ids );
    std::sort( sorted.begin(), sorted.end() );
    for( int i = 0; i < num_elem; ++i )
        CHECK_EQUAL( i + 1, sorted[i] );

    // check connectivity of each element
    std::vector< int > conn_ids( node_per_elem );
    std::vector< EntityHandle > conn_h;
    Range::iterator j = elems.begin();
    for( int i = 0; i < num_elem; ++i, ++j )
    {
        conn_h.clear();
        rval = mb.get_connectivity( &*j, 1, conn_h );CHECK_ERR( rval );
        CHECK_EQUAL( node_per_elem, (int)conn_h.size() );
        rval = mb.tag_get_data( gid_tag, &conn_h[0], node_per_elem, &conn_ids[0] );CHECK_ERR( rval );
        const int* exp = conn_list + node_per_elem * ( ids[i] - 1 );
        for( int k = 0; k < node_per_elem; ++k )
            CHECK_EQUAL( exp[k], conn_ids[k] );
    }
}

void test_edges()
{
    const int edge_conn[] = { 1,  3,                                                                           // 1
                              3,  2,  2,  5,  5,  4,  4,  7,  7,  6,  6,  8,  8,  1,  3,  9,  9,  8,           // 10
                              5,  9,  9,  7,  10, 12, 12, 11, 11, 14, 14, 13, 13, 16, 16, 15, 15, 17, 17, 10,  // 20
                              12, 18, 18, 17, 14, 18, 18, 16, 6,  19, 19, 13, 1,  20, 20, 11, 19, 21, 21, 8,   // 30
                              14, 21, 21, 20, 4,  22, 22, 15, 22, 23, 23, 7,  16, 23, 23, 19, 2,  24, 24, 10,  // 40
                              24, 25, 25, 5,  17, 25, 25, 22, 20, 26, 26, 3,  12, 26, 26, 24, 28, 30, 30, 29,  // 50
                              29, 31, 31, 2,  1,  32, 32, 28, 30, 33, 33, 32, 31, 33, 33, 3,  34, 36, 36, 35,  // 60
                              35, 37, 37, 11, 10, 38, 38, 34, 36, 39, 39, 38, 37, 39, 39, 12, 28, 40, 40, 35,  // 70
                              20, 41, 41, 32, 37, 41, 41, 40, 29, 42, 42, 34, 42, 43, 43, 31, 38, 43, 43, 24,  // 80
                              40, 44, 44, 30, 36, 44, 44, 42 };
    test_element( input_file_1, MBEDGE, 84, 2, edge_conn );
}

void test_quads()
{
    const int quad_conn[] = {
        1,  3,  9,  8,  // 1
        3,  2,  5,  9,  8,  9,  7,  6,  9,  5,  4,  7,  10, 12, 18, 17, 12, 11, 14, 18,
        17, 18, 16, 15, 18, 14, 13, 16, 6,  19, 21, 8,  19, 13, 14, 21,  // 10
        8,  21, 20, 1,  21, 14, 11, 20, 4,  22, 23, 7,  22, 15, 16, 23, 7,  23, 19, 6,
        23, 16, 13, 19, 2,  24, 25, 5,  24, 10, 17, 25, 5,  25, 22, 4,  25, 17, 15, 22,  // 20
        1,  20, 26, 3,  20, 11, 12, 26, 3,  26, 24, 2,  26, 12, 10, 24, 28, 30, 33, 32,
        30, 29, 31, 33, 32, 33, 3,  1,  33, 31, 2,  3,  34, 36, 39, 38, 36, 35, 37, 39,  // 30
        38, 39, 12, 10, 39, 37, 11, 12, 1,  20, 41, 32, 20, 11, 37, 41, 32, 41, 40, 28,
        41, 37, 35, 40, 29, 42, 43, 31, 42, 34, 38, 43, 31, 43, 24, 2,  43, 38, 10, 24,  // 40
        28, 40, 44, 30, 40, 35, 36, 44, 30, 44, 42, 29, 44, 36, 34, 42,
    };
    test_element( input_file_1, MBQUAD, 44, 4, quad_conn );
}

void test_hexes()
{
    const int hex_conn[] = { 6,  19, 23, 7,  8,  21, 27, 9,  19, 13, 16, 23, 21, 14, 18, 27, 7,  23, 22, 4,  9,  27,
                             25, 5,  23, 16, 15, 22, 27, 18, 17, 25, 8,  21, 27, 9,  1,  20, 26, 3,  21, 14, 18, 27,
                             20, 11, 12, 26, 9,  27, 25, 5,  3,  26, 24, 2,  27, 18, 17, 25, 26, 12, 10, 24, 1,  20,
                             26, 3,  32, 41, 45, 33, 20, 11, 12, 26, 41, 37, 39, 45, 3,  26, 24, 2,  33, 45, 43, 31,
                             26, 12, 10, 24, 45, 39, 38, 43, 32, 41, 45, 33, 28, 40, 44, 30, 41, 37, 39, 45, 40, 35,
                             36, 44, 33, 45, 43, 31, 30, 44, 42, 29, 45, 39, 38, 43, 44, 36, 34, 42 };
    test_element( input_file_1, MBHEX, 16, 8, hex_conn );
}

template < int L >
std::vector< int > find_parents( const int parent_conn[][L], int num_parent, int id )
{
    std::vector< int > results;
    for( int i = 0; i < num_parent; ++i )
    {
        for( int j = 0; j < L; ++j )
        {
            if( parent_conn[i][j] == id ) results.push_back( i + 1 );
        }
    }
    return results;
}

int check_geometric_set( Interface& moab,
                         int dim,
                         int id,
                         const int* children,
                         int num_children,
                         std::vector< int > parents )
{
    ErrorCode rval;
    Tag gid_tag, dim_tag;

    rval = moab.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
    rval = moab.tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, dim_tag );CHECK_ERR( rval );
    void* tag_vals[] = { &dim, &id };
    Tag tags[]       = { dim_tag, gid_tag };
    Range ents;
    rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, tag_vals, 2, ents );CHECK_ERR( rval );
    CHECK_EQUAL( 1u, (unsigned)ents.size() );

    const EntityHandle geom = ents.front();
    std::vector< int > exp_rel, act_rel;
    std::vector< EntityHandle > rel;

    if( num_children )
    {
        exp_rel.resize( num_children );
        std::copy( children, children + num_children, exp_rel.begin() );
        std::sort( exp_rel.begin(), exp_rel.end() );
        rel.clear();
        rval = moab.get_child_meshsets( geom, rel );CHECK_ERR( rval );
        CHECK_EQUAL( num_children, (int)rel.size() );
        act_rel.resize( rel.size() );
        rval = moab.tag_get_data( gid_tag, &rel[0], rel.size(), &act_rel[0] );CHECK_ERR( rval );
        std::sort( act_rel.begin(), act_rel.end() );
        CHECK( exp_rel == act_rel );
    }

    if( !parents.empty() )
    {
        exp_rel = parents;
        std::sort( exp_rel.begin(), exp_rel.end() );
        rel.clear();
        rval = moab.get_parent_meshsets( geom, rel );CHECK_ERR( rval );
        CHECK_EQUAL( parents.size(), rel.size() );
        act_rel.resize( rel.size() );
        rval = moab.tag_get_data( gid_tag, &rel[0], rel.size(), &act_rel[0] );CHECK_ERR( rval );
        std::sort( act_rel.begin(), act_rel.end() );
        CHECK( exp_rel == act_rel );
    }

    return 0;
}

void test_geometric_topology()
{
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, input_file_1 );
    // expected geometric vertices, specified by global ID
    const int vertex_ids[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 13, 14 };
    // List of global IDs of surfacs in geometric volumes, indexed by ID-1
    const int volume_surfs[2][6] = { { 1, 2, 3, 4, 5, 6 }, { 7, 8, 9, 6, 11, 12 } };
    // List of global IDs of curves in geometric surfaces, indexed by ID-1
    // Curve IDs of zero indicates that corresponding surface doesn't exist.
    const int surf_curves[12][4] = { { 1, 2, 3, 4 },     { 5, 6, 7, 8 },    { 9, 6, 10, 4 },   { 11, 7, 9, 3 },
                                     { 12, 8, 11, 2 },   { 10, 5, 12, 1 },  { 13, 14, 1, 16 }, { 17, 18, 5, 20 },
                                     { 10, 18, 22, 16 }, { 0, 0, 0, 0 },  // no surf 10
                                     { 24, 20, 12, 14 }, { 22, 17, 24, 13 } };
    // List of global IDs of vertices in geometric curves, indexed by ID-1
    // Vertex IDs of zero indicates that corresponding curve doesn't exist.
    const int curve_verts[24][2] = { { 1, 2 },  { 2, 3 },   { 3, 4 },  { 4, 1 }, { 5, 6 },  // 5
                                     { 6, 7 },  { 7, 8 },   { 8, 5 },  { 4, 7 }, { 1, 6 },  // 10
                                     { 3, 8 },  { 2, 5 },   { 9, 10 },                      // 13
                                     { 10, 2 }, { 0, 0 },                                   // no curve 15
                                     { 1, 9 },  { 13, 14 }, { 14, 6 }, { 0, 0 },            // no curve 19
                                     { 5, 13 }, { 0, 0 },                                   // no curve 21
                                     { 9, 14 }, { 0, 0 },                                   // no curve 23
                                     { 10, 13 } };

    // check all vertices
    for( unsigned i = 0; i < ( sizeof( vertex_ids ) / sizeof( vertex_ids[0] ) ); ++i )
        check_geometric_set( mb, 0, vertex_ids[i], 0, 0, find_parents< 2 >( curve_verts, 24, vertex_ids[i] ) );

    // check all curves
    for( int i = 1; i <= 24; ++i )
        if( curve_verts[i - 1][0] )
            check_geometric_set( mb, 1, i, curve_verts[i - 1], 2, find_parents< 4 >( surf_curves, 12, i ) );

    // check all surfs
    for( int i = 1; i <= 12; ++i )
        if( surf_curves[i - 1][0] )
            check_geometric_set( mb, 2, i, surf_curves[i - 1], 4, find_parents< 6 >( volume_surfs, 2, i ) );

    // check all volumes
    std::vector< int > empty;
    for( int i = 1; i <= 2; ++i )
        check_geometric_set( mb, 3, i, volume_surfs[i - 1], 6, empty );
}

void test_geometric_sets()
{
    ErrorCode rval;
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, input_file_1 );
    Tag gid_tag, dim_tag;
    rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
    rval = mb.tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, dim_tag );CHECK_ERR( rval );

    // verify mesh entity counts
    Range verts, curves, surfs, vols;
    int dim = 0;
    // Cppcheck warning (false positive): variable dim is assigned a value that is never used
    // Cppcheck warning (false positive): variable dim is reassigned a value before the old one has
    // been used
    const void* vals[] = { &dim };
    rval               = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, verts );CHECK_ERR( rval );
    dim  = 1;
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, curves );CHECK_ERR( rval );
    dim  = 2;
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, surfs );CHECK_ERR( rval );
    dim  = 3;
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, vols );CHECK_ERR( rval );

    CHECK_EQUAL( 12u, (unsigned)verts.size() );
    CHECK_EQUAL( 20u, (unsigned)curves.size() );
    CHECK_EQUAL( 11u, (unsigned)surfs.size() );
    CHECK_EQUAL( 2u, (unsigned)vols.size() );

    // check that each vertex has a single node, and that the
    // node is also contained in any parent curve
    Range ents;
    Range::iterator i;
    for( i = verts.begin(); i != verts.end(); ++i )
    {
        ents.clear();
        rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
        CHECK_EQUAL( 1u, (unsigned)ents.size() );
        CHECK( ents.all_of_type( MBVERTEX ) );
    }

    // check that each curve has one node and two edges
    for( i = curves.begin(); i != curves.end(); ++i )
    {
        ents.clear();
        rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
        CHECK_EQUAL( 1u, (unsigned)ents.num_of_type( MBVERTEX ) );
        CHECK_EQUAL( 2u, (unsigned)ents.num_of_type( MBEDGE ) );
        CHECK_EQUAL( 3u, (unsigned)ents.size() );
    }

    // check that each surface has 1 node, 4 edges, 4 quads
    for( i = surfs.begin(); i != surfs.end(); ++i )
    {
        ents.clear();
        rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
        CHECK_EQUAL( 1u, (unsigned)ents.num_of_type( MBVERTEX ) );
        CHECK_EQUAL( 4u, (unsigned)ents.num_of_type( MBEDGE ) );
        CHECK_EQUAL( 4u, (unsigned)ents.num_of_type( MBQUAD ) );
        CHECK_EQUAL( 9u, (unsigned)ents.size() );
    }

    // check that each volume has 1 node and 8 hexes.
    for( i = vols.begin(); i != vols.end(); ++i )
    {
        ents.clear();
        rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
        CHECK_EQUAL( 1u, (unsigned)ents.num_of_type( MBVERTEX ) );
        CHECK_EQUAL( 8u, (unsigned)ents.num_of_type( MBHEX ) );
        CHECK_EQUAL( 9u, (unsigned)ents.size() );
    }

    // Check that for each geometric entity, any contained vertices
    // are adjacent to some entity in one of its parents.
    Range parents, geom, nodes, tmp;
    for( int d = 0; d < 3; ++d )
    {
        const void* vals1[] = { &d };
        rval                = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals1, 1, geom );CHECK_ERR( rval );

        for( i = geom.begin(); i != geom.end(); ++i )
        {
            nodes.clear();
            ents.clear();
            parents.clear();
            rval = mb.get_entities_by_type( *i, MBVERTEX, nodes );CHECK_ERR( rval );
            rval = mb.get_parent_meshsets( *i, parents );CHECK_ERR( rval );
            for( Range::iterator j = parents.begin(); j != parents.end(); ++j )
            {
                tmp.clear();
                rval = mb.get_entities_by_dimension( *j, d + 1, tmp );CHECK_ERR( rval );
                ents.merge( tmp );
            }
            tmp.clear();
            rval = mb.get_adjacencies( ents, 0, false, tmp, Interface::UNION );CHECK_ERR( rval );
            nodes = subtract( nodes, tmp );
            CHECK( nodes.empty() );
        }
    }
}

// expect one block containing entire mesh, with id == 1
void test_blocks()
{
    ErrorCode rval;
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, input_file_1 );
    Tag mat_tag;
    rval = mb.tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mat_tag );CHECK_ERR( rval );

    Range blocks;
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &mat_tag, 0, 1, blocks );CHECK_ERR( rval );
    CHECK_EQUAL( 1u, (unsigned)blocks.size() );
    EntityHandle block = blocks.front();
    int id;
    rval = mb.tag_get_data( mat_tag, &block, 1, &id );CHECK_ERR( rval );
    CHECK_EQUAL( 1, id );

    Range block_hexes, mesh_hexes;
    rval = mb.get_entities_by_dimension( 0, 3, mesh_hexes );CHECK_ERR( rval );
    rval = mb.get_entities_by_dimension( block, 3, block_hexes, true );CHECK_ERR( rval );
    CHECK( mesh_hexes == block_hexes );
}

// Common code for test_side_sets and test_node sets
//\param tag_name  NEUMANN_SET_TAG_NAME or DIRICHLET_SET_TAG_NAME
//\param count     Number of expected sets
//\param ids       Expected IDs of sets
//\param set_surfs One list for each id in "ids" containing the
//                 ids of the geometric surfaces expected to be
//                 contained in the boundary condition set.
void test_bc_sets( const char* tag_name, unsigned count, const int* ids, const std::vector< int > set_surfs[] )
{
    ErrorCode rval;
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, input_file_1 );
    Tag ss_tag, gid_tag, dim_tag;
    rval = mb.tag_get_handle( tag_name, 1, MB_TYPE_INTEGER, ss_tag );CHECK_ERR( rval );
    rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
    rval = mb.tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, dim_tag );CHECK_ERR( rval );

    // check number of sidesets and IDs
    Range sidesets;
    rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &ss_tag, 0, 1, sidesets );CHECK_ERR( rval );
    CHECK_EQUAL( count, (unsigned)sidesets.size() );
    std::vector< EntityHandle > handles( count, 0 );
    for( Range::iterator i = sidesets.begin(); i != sidesets.end(); ++i )
    {
        int id;
        rval = mb.tag_get_data( ss_tag, &*i, 1, &id );CHECK_ERR( rval );
        unsigned idx;
        for( idx = 0; idx < count; ++idx )
            if( ids[idx] == id ) break;
        CHECK( idx != count );
        CHECK( handles[idx] == 0 );
        handles[idx] = *i;
    }

    // get surface faces
    std::vector< Range > exp( count );
    Range surfs, tmp;
    Tag tags[] = { dim_tag, gid_tag };
    for( unsigned i = 0; i < count; ++i )
    {
        exp[i].clear();
        surfs.clear();
        const int two = 2;
        for( unsigned j = 0; j < set_surfs[i].size(); ++j )
        {
            const void* vals[] = { &two, &set_surfs[i][j] };
            surfs.clear();
            rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, surfs );CHECK_ERR( rval );
            CHECK_EQUAL( 1u, (unsigned)surfs.size() );
            tmp.clear();
            rval = mb.get_entities_by_dimension( surfs.front(), 2, tmp, true );CHECK_ERR( rval );
            exp[i].merge( tmp );
        }
    }

    // check each bc set
    Range act;
    for( unsigned i = 0; i < count; ++i )
    {
        act.clear();
        rval = mb.get_entities_by_dimension( handles[i], 2, act, true );CHECK_ERR( rval );
        CHECK( exp[i] == act );
    }
}

// expect two sidesets containing two geometric surfaces each:
// sideset 1 : surfaces 1 and 7
// sideset 2 : surfaces 5 and 11
void test_side_sets()
{
    int ids[] = { 1, 2 };
    std::vector< int > surfs[2];
    surfs[0].push_back( 1 );
    surfs[0].push_back( 7 );
    surfs[1].push_back( 5 );
    surfs[1].push_back( 11 );
    test_bc_sets( NEUMANN_SET_TAG_NAME, 2, ids, surfs );
}

void test_node_sets()
{
    int ids[] = { 1, 2 };
    std::vector< int > surfs[2];
    surfs[0].push_back( 2 );
    surfs[0].push_back( 8 );
    surfs[1].push_back( 3 );
    surfs[1].push_back( 9 );
    test_bc_sets( DIRICHLET_SET_TAG_NAME, 2, ids, surfs );
}

static EntityHandle find_side( Interface& moab, EntityHandle entity, int side_dim, int side_num )
{
    ErrorCode rval;

    std::vector< EntityHandle > adj;
    rval = moab.get_adjacencies( &entity, 1, side_dim, false, adj );CHECK_ERR( rval );

    int sub_ent_indices[4];
    CN::SubEntityVertexIndices( TYPE_FROM_HANDLE( entity ), side_dim, side_num, sub_ent_indices );
    EntityType subtype  = CN::SubEntityType( TYPE_FROM_HANDLE( entity ), side_dim, side_num );
    int sub_ent_corners = CN::VerticesPerEntity( subtype );

    const EntityHandle* conn;
    int conn_len;
    rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );

    for( size_t i = 0; i < adj.size(); ++i )
    {
        if( TYPE_FROM_HANDLE( adj[i] ) != subtype ) continue;

        const EntityHandle* sub_conn;
        int sub_len;
        rval = moab.get_connectivity( adj[i], sub_conn, sub_len );CHECK_ERR( rval );

        int n = std::find( sub_conn, sub_conn + sub_len, conn[sub_ent_indices[0]] ) - sub_conn;
        if( n == sub_len )  // no vertex in common
            continue;

        // check forward direction
        int j;
        for( j = 1; j < sub_ent_corners; ++j )
            if( conn[sub_ent_indices[j]] != sub_conn[( j + n ) % sub_ent_corners] ) break;
        if( j == sub_ent_corners ) return adj[i];

        // check reverse direction
        for( j = 1; j < sub_ent_corners; ++j )
            if( conn[sub_ent_indices[j]] != sub_conn[( n + sub_ent_corners - j ) % sub_ent_corners] ) break;
        if( j == sub_ent_corners ) return adj[i];
    }

    // no match
    return 0;
}

void check_adj_ho_nodes( Interface& moab, EntityHandle entity )
{
    EntityType type = TYPE_FROM_HANDLE( entity );
    const EntityHandle* conn;
    int conn_len;
    ErrorCode rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );

    int ho[4];
    CN::HasMidNodes( type, conn_len, ho );
    for( int dim = CN::Dimension( type ) - 1; dim > 0; --dim )
    {
        if( !ho[dim] ) continue;

        for( int j = 0; j < CN::NumSubEntities( type, dim ); ++j )
        {
            EntityHandle side = find_side( moab, entity, dim, j );
            if( !side ) continue;

            const EntityHandle* side_conn;
            int side_len;
            rval = moab.get_connectivity( side, side_conn, side_len );CHECK_ERR( rval );

            int this_idx = CN::HONodeIndex( type, conn_len, dim, j );
            int side_idx = CN::HONodeIndex( TYPE_FROM_HANDLE( side ), side_len, dim, 0 );
            CHECK_EQUAL( side_conn[side_idx], conn[this_idx] );
        }
    }
}

// Check that element has expected higher-order nodes
// and that each higher-order node is at the center
// of the sub-entity it is on.
void check_ho_element( Interface& moab, EntityHandle entity, int mid_nodes[4] )
{
    // get element info
    const EntityType type = TYPE_FROM_HANDLE( entity );
    const EntityHandle* conn;
    int conn_len;
    ErrorCode rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );
    std::vector< double > coords( 3 * conn_len );
    rval = moab.get_coords( conn, conn_len, &coords[0] );CHECK_ERR( rval );

    // calculate and verify expected number of mid nodes
    int num_nodes = CN::VerticesPerEntity( type );
    for( int d = 1; d <= CN::Dimension( type ); ++d )
        if( mid_nodes[d] ) num_nodes += CN::NumSubEntities( type, d );
    CHECK_EQUAL( num_nodes, conn_len );

    // verify that each higher-order node is at the center
    // of its respective sub-entity.
    for( int i = CN::VerticesPerEntity( type ); i < num_nodes; ++i )
    {
        // get sub-entity owning ho-node
        int sub_dim, sub_num;
        CN::HONodeParent( type, num_nodes, i, sub_dim, sub_num );
        // get corner vertex indices
        int sub_conn[8], num_sub;
        if( sub_dim < CN::Dimension( type ) )
        {
            CN::SubEntityVertexIndices( type, sub_dim, sub_num, sub_conn );
            EntityType sub_type = CN::SubEntityType( type, sub_dim, sub_num );
            num_sub             = CN::VerticesPerEntity( sub_type );
        }
        else
        {
            num_sub = CN::VerticesPerEntity( type );
            for( int j = 0; j < num_sub; ++j )
                sub_conn[j] = j;
        }
        // calculate mean of corner vertices
        double mean[3] = { 0, 0, 0 };
        for( int j = 0; j < num_sub; ++j )
        {
            int co = 3 * sub_conn[j];
            mean[0] += coords[co];
            mean[1] += coords[co + 1];
            mean[2] += coords[co + 2];
        }
        mean[0] /= num_sub;
        mean[1] /= num_sub;
        mean[2] /= num_sub;
        // verify that higher-order node is at expected location
        CHECK_REAL_EQUAL( mean[0], coords[3 * i], 1e-6 );
        CHECK_REAL_EQUAL( mean[1], coords[3 * i + 1], 1e-6 );
        CHECK_REAL_EQUAL( mean[2], coords[3 * i + 2], 1e-6 );
    }
}

// Validate elements of specified type.
// Looks for a block containing the specified entity type
// and with the specified mid-node flags set in its
// HAS_MID_NODES_TAG.
void test_ho_elements( EntityType type, int num_nodes )
{
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, ho_file );

    ErrorCode rval;
    Tag ho_tag, block_tag;
    rval = mb.tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag );CHECK_ERR( rval );
    rval = mb.tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, ho_tag );CHECK_ERR( rval );

    // get material sets with expected higher-order nodes
    Range blocks;
    int ho_flags[4];
    CN::HasMidNodes( type, num_nodes, ho_flags );
    Tag tags[2]   = { ho_tag, block_tag };
    void* vals[2] = { ho_flags, NULL };
    rval          = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, blocks );CHECK_ERR( rval );

    Range::iterator i;
    Range entities;
    for( i = blocks.begin(); i != blocks.end(); ++i )
    {
        rval = mb.get_entities_by_type( *i, type, entities, true );CHECK_ERR( rval );
    }
    // test file should contain all types of HO entities
    CHECK( !entities.empty() );
    // test ho node positions -- expect to be at center of adj corners
    for( i = entities.begin(); i != entities.end(); ++i )
        check_ho_element( mb, *i, ho_flags );
    // test ho node handles consistent with adjacent entities
    for( i = entities.begin(); i != entities.end(); ++i )
        check_adj_ho_nodes( mb, *i );
}

void test_multiple_files()
{
    // load two surface meshes, one at z=+5 at z=-5.
    ErrorCode rval;
    Core mb_impl;
    Interface& mb = mb_impl;
    Range file1_ents, file2_ents;
    read_file( mb, input_file_1 );
    mb.get_entities_by_handle( 0, file1_ents );
    read_file( mb, input_file_1 );
    mb.get_entities_by_handle( 0, file2_ents );
    file2_ents = subtract( file2_ents, file1_ents );
    EntityHandle file1, file2;
    mb.create_meshset( MESHSET_SET, file1 );
    mb.create_meshset( MESHSET_SET, file2 );
    mb.add_entities( file1, file1_ents );
    mb.add_entities( file2, file2_ents );

    // first check that we get the same number of verts from
    // each file and that they are distinct vertices
    Range file1_verts, file2_verts;
    rval = mb.get_entities_by_type( file1, MBVERTEX, file1_verts );CHECK_ERR( rval );
    CHECK( !file1_verts.empty() );
    rval = mb.get_entities_by_type( file2, MBVERTEX, file2_verts );CHECK_ERR( rval );
    CHECK( !file2_verts.empty() );
    CHECK_EQUAL( file1_verts.size(), file2_verts.size() );
    CHECK( intersect( file1_verts, file2_verts ).empty() );

    // now check that we get the same number of elements from
    // each file and that they are distinct
    Range file1_elems, file2_elems;
    rval = mb.get_entities_by_dimension( file1, 3, file1_elems );CHECK_ERR( rval );
    CHECK( !file1_elems.empty() );
    rval = mb.get_entities_by_dimension( file2, 3, file2_elems );CHECK_ERR( rval );
    CHECK( !file2_elems.empty() );
    CHECK_EQUAL( file1_elems.size(), file2_elems.size() );
    CHECK( intersect( file1_elems, file2_elems ).empty() );

    // now check that the connectivity for each element is
    // defined using the appropriate vertex instances
    Range file1_elem_verts, file2_elem_verts;
    rval = mb.get_adjacencies( file1_elems, 0, false, file1_elem_verts, Interface::UNION );CHECK_ERR( rval );
    rval = mb.get_adjacencies( file2_elems, 0, false, file2_elem_verts, Interface::UNION );CHECK_ERR( rval );
    CHECK_EQUAL( file1_elem_verts.size(), file2_elem_verts.size() );
    CHECK( intersect( file1_elem_verts, file1_verts ) == file1_elem_verts );
    CHECK( intersect( file2_elem_verts, file2_verts ) == file2_elem_verts );
}

void test_cubit12()
{
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, cubit12_file );
}

void test_cubit14()
{
    Core mb_impl;
    Interface& mb = mb_impl;
    read_file( mb, cubit14_file );
    // check the global id for some geometry sets
    GeomTopoTool gtt( &mb_impl );
    Range ranges[5];
    ErrorCode rval = gtt.find_geomsets( ranges );CHECK_ERR( rval );
    EntityHandle set0 = ranges[0][0];  // does it have a global id > 0?
    Tag gid_tag;
    rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );

    int val;
    rval = mb.tag_get_data( gid_tag, &set0, 1, &val );
    CHECK( val != 0 );
}