exodus_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 "ReadNCDF.hpp"
#include "moab/FileOptions.hpp"
#define IS_BUILDING_MB
#include "ExoIIUtil.hpp"
#include <cmath>
#include <algorithm>
 
using namespace moab;
 
/* Input test file: ho_test.g
 *
 * 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.
 */
#ifdef MESHDIR
static const char ho_file[] = STRINGIFY( MESHDIR ) "/io/ho_test.g";
static const char file_one[] = STRINGIFY( MESHDIR ) "/mbtest1.g";
static const char alt_file[] = STRINGIFY( MESHDIR ) "/io/hex_2x2x2_ss.exo";
static const char polyg[] = STRINGIFY( MESHDIR ) "/io/poly8-10.vtk";
static const char polyh[] = STRINGIFY( MESHDIR ) "/io/polyhedra.vtk";
static const char rpolyg[] = STRINGIFY( MESHDIR ) "/io/polyg.exo";
static const char rpolyh[] = STRINGIFY( MESHDIR ) "/io/polyh.exo";
#else
static const char ho_file[] = "ho_test.g";
static const char file_one[] = "mbtest1.g";
static const char alt_file[] = "hex_2x2x2_ss.exo";
static const char polyg[] = "poly8-10.vtk";
static const char polyh[] = "polyhedra.vtk";
static const char rpolyg[] = "polyg.exo";
static const char rpolyh[] = "polyh.exo";
#endif
 
void read_file( Interface& moab, const char* input_file );
 
// 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] );
 
void test_read_side( int sideset_id, EntityType sideset_type, int sideset_nodes_per_elem, bool shell_side = false );
 
// 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 );
 
// Tests originally in MBTest.cpp
void mb_vertex_coordinate_test();
void mb_bar_connectivity_test();
void mb_tri_connectivity_test();
void mb_quad_connectivity_test();
void mb_hex_connectivity_test();
void mb_tet_connectivity_test();
void mb_write_mesh_test();
 
void test_types();
 
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_read_tri6_side()
{
 test_read_side( 1, MBEDGE, 3 );
} // sideset 1
void test_read_shell_side()
{
 test_read_side( 3, MBQUAD, 9, true );
} // sideset 3
void test_read_shell_edge()
{
 test_read_side( 4, MBEDGE, 3 );
} // sideset 4
void test_read_hex20_side()
{
 test_read_side( 2, MBQUAD, 8 );
} // sideset 2
 
void test_read_block_ids();
void test_read_sideset_ids();
void test_read_nodeset_ids();
void test_write_polygons();
void test_write_polyhedra();
void test_read_polygons();
void test_read_polyhedra();
 
void test_read_alternate_coord_format();
 
int main()
{
 int result = 0;
 
 result += RUN_TEST( mb_vertex_coordinate_test );
 result += RUN_TEST( mb_bar_connectivity_test );
 result += RUN_TEST( mb_tri_connectivity_test );
 result += RUN_TEST( mb_quad_connectivity_test );
 result += RUN_TEST( mb_hex_connectivity_test );
 result += RUN_TEST( mb_tet_connectivity_test );
 result += RUN_TEST( mb_write_mesh_test );
 
 result += RUN_TEST( test_types );
 
 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_read_tri6_side );
 result += RUN_TEST( test_read_shell_side );
 result += RUN_TEST( test_read_shell_edge );
 result += RUN_TEST( test_read_hex20_side );
 
 result += RUN_TEST( test_read_block_ids );
 result += RUN_TEST( test_read_sideset_ids );
 result += RUN_TEST( test_read_nodeset_ids );
 
 result += RUN_TEST( test_read_alternate_coord_format );
 
 result += RUN_TEST( test_write_polygons );
 result += RUN_TEST( test_write_polyhedra );
 result += RUN_TEST( test_read_polygons );
 result += RUN_TEST( test_read_polyhedra );
 
 return result;
}
 
void load_file_one( Interface* iface )
{
 ErrorCode error = iface->load_mesh( file_one );
 if( MB_SUCCESS != error )
 {
 std::cout << "Failed to load input file: " << file_one << std::endl;
 std::string error_reason;
 iface->get_last_error( error_reason );
 std::cout << error_reason << std::endl;
 }
 CHECK_ERR( error );
}
 
/*!
 @test
 Vertex Coordinates
 @li Get coordinates of vertex 1 correctly
 @li Get coordinates of vertex 8 correctly
 @li Get coordinates of vertex 6 correctly
*/
void mb_vertex_coordinate_test()
{
 double coords[3];
 EntityHandle handle;
 ErrorCode error;
 int err;
 
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 
 // coordinate 2 should be {1.5, -1.5, 3.5}
 
 handle = CREATE_HANDLE( MBVERTEX, 2, err );
 error = MB->get_coords( &handle, 1, coords );CHECK_ERR( error );
 const double exp2[] = { 1.5, -1.5, 3.5 };
 CHECK_ARRAYS_EQUAL( exp2, 3, coords, 3 );
 
 // coordinate 9 should be {1, -2, 3.5}
 handle = CREATE_HANDLE( MBVERTEX, 9, err );
 error = MB->get_coords( &handle, 1, coords );CHECK_ERR( error );
 const double exp9[] = { 1, -2, 3.5 };
 CHECK_ARRAYS_EQUAL( exp9, 3, coords, 3 );
 
 // coordinate 7 should be {0.5, -2, 3.5}
 handle = CREATE_HANDLE( MBVERTEX, 7, err );
 error = MB->get_coords( &handle, 1, coords );CHECK_ERR( error );
 const double exp7[] = { 0.5, -2, 3.5 };
 CHECK_ARRAYS_EQUAL( exp7, 3, coords, 3 );
 
 int node_count = 0;
 error = MB->get_number_entities_by_type( 0, MBVERTEX, node_count );CHECK_ERR( error );
 // Number of vertices (node_count) should be 83 assuming no gaps in the handle space
 CHECK_EQUAL( 47, node_count );
}
 
/*!
 @test
 MB Bar Element Connectivity Test
 @li Get coordinates for 2 node bar elements
*/
 
void mb_bar_connectivity_test()
{
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 
 std::vector< EntityHandle > conn;
 Range bars;
 
 ErrorCode error = MB->get_entities_by_type( 0, MBEDGE, bars );CHECK_ERR( error );
 
 // get the connectivity of the second bar
 EntityHandle handle = *( ++bars.begin() );
 
 error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
 
 CHECK_EQUAL( (size_t)2, conn.size() );
 
 // from ncdump the connectivity of bar 2 (0 based) is
 // 14, 13
 CHECK_EQUAL( (EntityHandle)14, conn[0] );
 CHECK_EQUAL( (EntityHandle)13, conn[1] );
 
 // Now try getting the connectivity of one of the vertices for fun.
 // just return the vertex in the connectivity
 handle = conn[0];
 error = MB->get_connectivity( &handle, 1, conn );
 CHECK_EQUAL( MB_FAILURE, error );
}
 
void mb_tri_connectivity_test()
{
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 
 std::vector< EntityHandle > conn;
 Range tris;
 ErrorCode error = MB->get_entities_by_type( 0, MBTRI, tris );CHECK_ERR( error );
 
 // get the connectivity of the second tri
 EntityHandle handle = *( ++tris.begin() );
 
 error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
 
 CHECK_EQUAL( (size_t)3, conn.size() );
 
 // from ncdump the connectivity of tri 2 (0 based) is
 // 45, 37, 38
 
 CHECK_EQUAL( (EntityHandle)45, conn[0] );
 CHECK_EQUAL( (EntityHandle)37, conn[1] );
 CHECK_EQUAL( (EntityHandle)38, conn[2] );
}
 
void mb_quad_connectivity_test()
{
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 
 std::vector< EntityHandle > conn;
 Range quads;
 
 ErrorCode error = MB->get_entities_by_type( 0, MBQUAD, quads );CHECK_ERR( error );
 
 // get the connectivity of the second quad
 EntityHandle handle = *( ++quads.begin() );
 
 error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
 
 CHECK_EQUAL( (size_t)4, conn.size() );
 
 // from ncdump the connectivity of quad 2 (0 based) is
 // 20, 11, 12, 26,
 
 CHECK_EQUAL( (EntityHandle)20, conn[0] );
 CHECK_EQUAL( (EntityHandle)11, conn[1] );
 CHECK_EQUAL( (EntityHandle)12, conn[2] );
 CHECK_EQUAL( (EntityHandle)26, conn[3] );
}
 
void mb_hex_connectivity_test()
{
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 
 std::vector< EntityHandle > conn;
 Range hexes;
 
 ErrorCode error = MB->get_entities_by_type( 0, MBHEX, hexes );CHECK_ERR( error );
 
 // get the connectivity of the second hex
 EntityHandle handle = *( ++hexes.begin() );
 
 error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
 
 CHECK_EQUAL( (size_t)8, conn.size() );
 
 // from ncdump the connectivity of hex 1 (0 based) is
 // 19, 13, 16, 23, 21, 14, 18, 27
 
 CHECK_EQUAL( (EntityHandle)19, conn[0] );
 CHECK_EQUAL( (EntityHandle)13, conn[1] );
 CHECK_EQUAL( (EntityHandle)16, conn[2] );
 CHECK_EQUAL( (EntityHandle)23, conn[3] );
 CHECK_EQUAL( (EntityHandle)21, conn[4] );
 CHECK_EQUAL( (EntityHandle)14, conn[5] );
 CHECK_EQUAL( (EntityHandle)18, conn[6] );
 CHECK_EQUAL( (EntityHandle)27, conn[7] );
}
 
void mb_tet_connectivity_test()
{
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 
 std::vector< EntityHandle > conn;
 Range tets;
 ErrorCode error = MB->get_entities_by_type( 0, MBTET, tets );CHECK_ERR( error );
 
 // get the connectivity of the second tet
 EntityHandle handle = *( ++tets.begin() );
 
 error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
 
 CHECK_EQUAL( (size_t)4, conn.size() );
 
 // from ncdump the connectivity of tet 2 (0 based) is:
 // 35, 34, 32, 43
 
 CHECK_EQUAL( (EntityHandle)35, conn[0] );
 CHECK_EQUAL( (EntityHandle)34, conn[1] );
 CHECK_EQUAL( (EntityHandle)32, conn[2] );
 CHECK_EQUAL( (EntityHandle)43, conn[3] );
}
 
void mb_write_mesh_test()
{
 Core moab;
 Interface* MB = &moab;
 load_file_one( MB );
 ErrorCode result;
 
 std::string file_name = "mb_write.g";
 
 // no need to get lists, write out the whole mesh
 result = MB->write_mesh( file_name.c_str() );CHECK_ERR( result );
 
 //---------The following tests outputting meshsets that are in meshsets of blocks ---/
 
 // lets create a block meshset and put some entities and meshsets into it
 EntityHandle block_ms;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, block_ms );CHECK_ERR( result );
 
 // make another meshset to put quads in, so SHELLs can be written out
 EntityHandle block_of_shells;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, block_of_shells );CHECK_ERR( result );
 
 // tag the meshset so it's a block, with id 100
 int id = 100;
 Tag tag_handle;
 const int negone = -1;
 result = MB->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
 result = MB->tag_set_data( tag_handle, &block_ms, 1, &id );CHECK_ERR( result );
 id = 101;
 result = MB->tag_set_data( tag_handle, &block_of_shells, 1, &id );CHECK_ERR( result );
 
 // set dimension tag on this to ensure shells get output; reuse id variable
 result = MB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
 id = 3;
 result = MB->tag_set_data( tag_handle, &block_of_shells, 1, &id );CHECK_ERR( result );
 
 // get some entities (tets)
 Range temp_range;
 result = MB->get_entities_by_type( 0, MBHEX, temp_range );CHECK_ERR( result );
 
 Range::iterator iter, end_iter;
 iter = temp_range.begin();
 end_iter = temp_range.end();
 
 // add evens to 'block_ms'
 std::vector< EntityHandle > temp_vec;
 for( ; iter != end_iter; ++iter )
 {
 if( ID_FROM_HANDLE( *iter ) % 2 == 0 ) temp_vec.push_back( *iter );
 }
 result = MB->add_entities( block_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
 
 // make another meshset
 EntityHandle ms_of_block_ms;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ms_of_block_ms );CHECK_ERR( result );
 
 // add some entities to it
 temp_vec.clear();
 iter = temp_range.begin();
 for( ; iter != end_iter; ++iter )
 {
 if( ID_FROM_HANDLE( *iter ) % 2 ) // add all odds
 temp_vec.push_back( *iter );
 }
 result = MB->add_entities( ms_of_block_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
 
 // add the other meshset to the block's meshset
 result = MB->add_entities( block_ms, &ms_of_block_ms, 1 );CHECK_ERR( result );
 
 //---------------testing sidesets----------------/
 
 // lets create a sideset meshset and put some entities and meshsets into it
 EntityHandle sideset_ms;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, sideset_ms );CHECK_ERR( result );
 
 // tag the meshset so it's a sideset, with id 104
 id = 104;
 result = MB->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
 
 result = MB->tag_set_data( tag_handle, &sideset_ms, 1, &id );CHECK_ERR( result );
 
 // get some entities (tris)
 temp_range.clear();
 result = MB->get_entities_by_type( 0, MBQUAD, temp_range );CHECK_ERR( result );
 
 iter = temp_range.begin();
 end_iter = temp_range.end();
 
 // add evens to 'sideset_ms'
 temp_vec.clear();
 for( ; iter != end_iter; ++iter )
 {
 if( ID_FROM_HANDLE( *iter ) % 2 == 0 ) temp_vec.push_back( *iter );
 }
 result = MB->add_entities( sideset_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
 
 // make another meshset
 EntityHandle ms_of_sideset_ms;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ms_of_sideset_ms );CHECK_ERR( result );
 
 // add some entities to it
 temp_vec.clear();
 iter = temp_range.begin();
 for( ; iter != end_iter; ++iter )
 {
 if( ID_FROM_HANDLE( *iter ) % 2 ) // add all odds
 temp_vec.push_back( *iter );
 }
 result = MB->add_entities( ms_of_sideset_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
 
 // add the other meshset to the sideset's meshset
 result = MB->add_entities( sideset_ms, &ms_of_sideset_ms, 1 );CHECK_ERR( result );
 
 //---------test sense on meshsets (reverse/foward)-------//
 
 // get all quads whose x-coord = 2.5 and put them into a meshset_a
 EntityHandle meshset_a;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_a );CHECK_ERR( result );
 
 temp_range.clear();
 result = MB->get_entities_by_type( 0, MBQUAD, temp_range );CHECK_ERR( result );
 
 std::vector< EntityHandle > nodes1, entity_vec;
 std::copy( temp_range.begin(), temp_range.end(), std::back_inserter( entity_vec ) );
 result = MB->get_connectivity( &entity_vec[0], entity_vec.size(), nodes1 );CHECK_ERR( result );
 assert( nodes1.size() == 4 * temp_range.size() );
 temp_vec.clear();
 std::vector< double > coords( 3 * nodes1.size() );
 result = MB->get_coords( &nodes1[0], nodes1.size(), &coords[0] );CHECK_ERR( result );
 
 unsigned int k = 0;
 for( Range::iterator it = temp_range.begin(); it != temp_range.end(); ++it )
 {
 if( coords[12 * k] == 2.5 && coords[12 * k + 3] == 2.5 && coords[12 * k + 6] == 2.5 &&
 coords[12 * k + 9] == 2.5 )
 temp_vec.push_back( *it );
 k++;
 }
 result = MB->add_entities( meshset_a, ( temp_vec.empty() ) ? NULL : &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
 result = MB->add_entities( block_of_shells, ( temp_vec.empty() ) ? NULL : &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
 
 // put these quads into a different meshset_b and tag them with a reverse sense tag
 EntityHandle meshset_b;
 result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_b );CHECK_ERR( result );
 
 result = MB->add_entities( meshset_b, &meshset_a, 1 );CHECK_ERR( result );
 
 result = MB->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, tag_handle, MB_TAG_SPARSE | MB_TAG_CREAT );CHECK_ERR( result );
 
 int reverse_value = -1;
 result = MB->tag_set_data( tag_handle, &meshset_b, 1, &reverse_value );CHECK_ERR( result );
 
 // get some random quad, whose x-coord != 2.5, and put it into a different meshset_c
 // and tag it with a reverse sense tag
 
 iter = temp_range.begin();
 end_iter = temp_range.end();
 
 temp_vec.clear();
 for( ; iter != end_iter; ++iter )
 {
 std::vector< EntityHandle > nodes;
 result = MB->get_connectivity( &( *iter ), 1, nodes );CHECK_ERR( result );
 
 bool not_equal_2_5 = true;
 for( unsigned int ku = 0; ku < nodes.size(); ku++ )
 {
 double coords2[3] = { 0 };
 
 result = MB->get_coords( &( nodes[ku] ), 1, coords2 );CHECK_ERR( result );
 
 if( coords2[0] == 2.5 )
 {
 not_equal_2_5 = false;
 break;
 }
 }
 
 if( not_equal_2_5 && nodes.size() > 0 )
 {
 temp_vec.push_back( *iter );
 break;
 }
 }
 
 EntityHandle meshset_c;
 MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_c );
 
 result = MB->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, tag_handle );CHECK_ERR( result );
 
 reverse_value = -1;
 result = MB->tag_set_data( tag_handle, &meshset_c, 1, &reverse_value );CHECK_ERR( result );
 
 MB->add_entities( meshset_c, &temp_vec[0], temp_vec.size() );
 MB->add_entities( block_of_shells, &temp_vec[0], temp_vec.size() );
 
 // create another meshset_abc, adding meshset_a, meshset_b, meshset_c
 EntityHandle meshset_abc;
 MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_abc );
 
 temp_vec.clear();
 temp_vec.push_back( meshset_a );
 temp_vec.push_back( meshset_b );
 temp_vec.push_back( meshset_c );
 
 MB->add_entities( meshset_abc, &temp_vec[0], temp_vec.size() );
 
 // tag it so it's a sideset
 id = 444;
 result = MB->tag_get_handle( "NEUMANN_SET", 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
 
 result = MB->tag_set_data( tag_handle, &meshset_abc, 1, &id );CHECK_ERR( result );
 
 //---------------do nodesets now -----------------//
 
 // lets create a nodeset meshset and put some entities and meshsets into it
 EntityHandle nodeset_ms;
 MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, nodeset_ms );
 
 // tag the meshset so it's a nodeset, with id 119
 id = 119;
 result = MB->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
 
 result = MB->tag_set_data( tag_handle, &nodeset_ms, 1, &id );CHECK_ERR( result );
 
 // get all Quads
 temp_range.clear();
 result = MB->get_entities_by_type( 0, MBQUAD, temp_range );CHECK_ERR( result );
 
 // get all the nodes of the tris
 Range nodes_of_quads;
 iter = temp_range.begin();
 end_iter = temp_range.end();
 
 for( ; iter != end_iter; ++iter )
 {
 std::vector< EntityHandle > nodes;
 result = MB->get_connectivity( &( *iter ), 1, nodes );CHECK_ERR( result );
 
 for( unsigned int ku = 0; ku < nodes.size(); ku++ )
 nodes_of_quads.insert( nodes[ku] );
 }
 
 iter = nodes_of_quads.begin();
 end_iter = nodes_of_quads.end();
 
 // add evens to 'nodeset_ms'
 temp_vec.clear();
 for( ; iter != end_iter; ++iter )
 {
 if( ID_FROM_HANDLE( *iter ) % 2 == 0 ) temp_vec.push_back( *iter );
 }
 MB->add_entities( nodeset_ms, &temp_vec[0], temp_vec.size() );
 
 // make another meshset
 EntityHandle ms_of_nodeset_ms;
 MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ms_of_nodeset_ms );
 
 // add some entities to it
 temp_vec.clear();
 iter = nodes_of_quads.begin();
 end_iter = nodes_of_quads.end();
 for( ; iter != end_iter; ++iter )
 {
 if( ID_FROM_HANDLE( *iter ) % 2 ) // add all odds
 temp_vec.push_back( *iter );
 }
 MB->add_entities( ms_of_nodeset_ms, &temp_vec[0], temp_vec.size() );
 
 // add the other meshset to the nodeset's meshset
 MB->add_entities( nodeset_ms, &ms_of_nodeset_ms, 1 );
 
 // no need to get lists, write out the whole mesh
 file_name = "mb_write2.g";
 std::vector< EntityHandle > output_list;
 output_list.push_back( block_ms );
 output_list.push_back( sideset_ms );
 output_list.push_back( meshset_abc );
 output_list.push_back( nodeset_ms );
 output_list.push_back( block_of_shells );
 result = MB->write_mesh( file_name.c_str(), &output_list[0], output_list.size() );CHECK_ERR( result );
}
 
struct TestType
{
 EntityType moab_type;
 ExoIIElementType exo_type;
 int num_nodes;
 std::string name;
};
 
void check_type( const TestType& type )
{
 int has_mid_nodes[4];
 CN::HasMidNodes( type.moab_type, type.num_nodes, has_mid_nodes );
 
 CHECK_EQUAL( type.moab_type, ExoIIUtil::ExoIIElementMBEntity[type.exo_type] );
 CHECK_EQUAL( type.name, std::string( ExoIIUtil::ElementTypeNames[type.exo_type] ) );
 CHECK_EQUAL( type.num_nodes, ExoIIUtil::VerticesPerElement[type.exo_type] );
 auto dim = CN::Dimension( type.moab_type );
 if( 3 == dim ) CHECK_EQUAL( has_mid_nodes[3], ExoIIUtil::HasMidNodes[type.exo_type][3] );
 if( 2 <= dim ) CHECK_EQUAL( has_mid_nodes[2], ExoIIUtil::HasMidNodes[type.exo_type][2] );
 if( 1 <= dim ) CHECK_EQUAL( has_mid_nodes[1], ExoIIUtil::HasMidNodes[type.exo_type][1] );
 
 Core moab;
 ExoIIUtil tool( &moab );
 CHECK_EQUAL( type.exo_type, tool.element_name_to_type( type.name.c_str() ) );
 CHECK_EQUAL( type.name, std::string( tool.element_type_name( type.exo_type ) ) );
}
 
void test_types()
{
 const TestType types[] = { { MBVERTEX, EXOII_SPHERE, 1, "SPHERE" },
 { MBEDGE, EXOII_SPRING, 1, "SPRING" },
 { MBEDGE, EXOII_BAR, 2, "BAR" },
 { MBEDGE, EXOII_BAR2, 2, "BAR2" },
 { MBEDGE, EXOII_BAR3, 3, "BAR3" },
 { MBEDGE, EXOII_BEAM, 2, "BEAM" },
 { MBEDGE, EXOII_BEAM2, 2, "BEAM2" },
 { MBEDGE, EXOII_BEAM3, 3, "BEAM3" },
 { MBEDGE, EXOII_TRUSS, 2, "TRUSS" },
 { MBEDGE, EXOII_TRUSS2, 2, "TRUSS2" },
 { MBEDGE, EXOII_TRUSS3, 3, "TRUSS3" },
 { MBTRI, EXOII_TRI, 3, "TRI" },
 { MBTRI, EXOII_TRI3, 3, "TRI3" },
 { MBTRI, EXOII_TRI6, 6, "TRI6" },
 { MBTRI, EXOII_TRI7, 7, "TRI7" },
 { MBQUAD, EXOII_QUAD, 4, "QUAD" },
 { MBQUAD, EXOII_QUAD4, 4, "QUAD4" },
 { MBQUAD, EXOII_QUAD5, 5, "QUAD5" },
 { MBQUAD, EXOII_QUAD8, 8, "QUAD8" },
 { MBQUAD, EXOII_QUAD9, 9, "QUAD9" },
 { MBQUAD, EXOII_SHELL, 4, "SHELL" },
 { MBQUAD, EXOII_SHELL4, 4, "SHELL4" },
 { MBQUAD, EXOII_SHELL5, 5, "SHELL5" },
 { MBQUAD, EXOII_SHELL8, 8, "SHELL8" },
 { MBQUAD, EXOII_SHELL9, 9, "SHELL9" },
 { MBTET, EXOII_TETRA, 4, "TETRA" },
 { MBTET, EXOII_TETRA4, 4, "TETRA4" },
 { MBTET, EXOII_TET4, 4, "TET4" },
 { MBTET, EXOII_TETRA8, 8, "TETRA8" },
 { MBTET, EXOII_TETRA10, 10, "TETRA10" },
 { MBTET, EXOII_TETRA14, 14, "TETRA14" },
 { MBPYRAMID, EXOII_PYRAMID, 5, "PYRAMID" },
 { MBPYRAMID, EXOII_PYRAMID5, 5, "PYRAMID5" },
 { MBPYRAMID, EXOII_PYRAMID10, 10, "PYRAMID10" },
 { MBPYRAMID, EXOII_PYRAMID13, 13, "PYRAMID13" },
 { MBPYRAMID, EXOII_PYRAMID18, 18, "PYRAMID18" },
 { MBPRISM, EXOII_WEDGE, 6, "WEDGE" },
 { MBKNIFE, EXOII_KNIFE, 7, "KNIFE" },
 { MBHEX, EXOII_HEX, 8, "HEX" },
 { MBHEX, EXOII_HEX8, 8, "HEX8" },
 { MBHEX, EXOII_HEX9, 9, "HEX9" },
 { MBHEX, EXOII_HEX20, 20, "HEX20" },
 { MBHEX, EXOII_HEX27, 27, "HEX27" },
 { MBHEX, EXOII_HEXSHELL, 12, "HEXSHELL" } };
 const int num_types = sizeof( types ) / sizeof( types[0] );
 for( int i = 0; i < num_types; ++i )
 check_type( types[i] );
}
 
void read_file( Interface& moab, const char* input_file )
{
 ErrorCode rval = moab.load_file( input_file );CHECK_ERR( rval );
}
 
void write_and_read( Interface& write_mb, Interface& read_mb, EntityHandle block = 0 )
{
 const char* tmp_file = "exodus_test_tmp.g";
 ErrorCode rval;
 
 EntityHandle* write_set_list = &block;
 int write_set_list_len = 0; //(block != 0);
 rval = write_mb.write_file( tmp_file, "EXODUS", 0, write_set_list, write_set_list_len );
 if( MB_SUCCESS != rval ) remove( tmp_file );CHECK_ERR( rval );
 
 rval = read_mb.load_file( tmp_file );
 remove( tmp_file );CHECK_ERR( rval );
}
 
void check_ho_elements( Interface& moab, EntityHandle block, EntityType type, int mid_nodes[4] )
{
 ErrorCode rval;
 Range elems;
 rval = moab.get_entities_by_handle( block, elems, ( type != MBENTITYSET ? true : false ) );CHECK_ERR( rval );
 CHECK( !elems.empty() );
 CHECK( elems.all_of_type( type ) );
 for( Range::const_iterator i = elems.begin(); i != elems.end(); ++i )
 check_ho_element( moab, *i, mid_nodes );
}
 
// 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 );
 }
}
 
EntityHandle find_block( Interface& mb, EntityType type, const int has_mid_nodes[4] )
{
 
 ErrorCode rval;
 Tag ho_tag, block_tag;
 const int negone = -1;
 rval = mb.tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag, 0, &negone );CHECK_ERR( rval );
 const int mids[] = { -1, -1, -1, -1 };
 rval = mb.tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, ho_tag, 0, mids );CHECK_ERR( rval );
 
 // get material sets with expected higher-order nodes
 Range blocks;
 Tag tags[2] = { ho_tag, block_tag };
 const void* vals[2] = { has_mid_nodes, NULL };
 rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, blocks );CHECK_ERR( rval );
 
 for( Range::iterator i = blocks.begin(); i != blocks.end(); ++i )
 {
 int n;
 rval = mb.get_number_entities_by_type( *i, type, n );CHECK_ERR( rval );
 if( n > 0 ) return *i;
 }
 
 CHECK( false ); // no block matching element type description
 return 0;
}
 
EntityHandle find_sideset( Interface& mb, int sideset_id, EntityType /*side_type*/ )
{
 ErrorCode rval;
 Tag ss_tag;
 const int negone = -1;
 rval = mb.tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ss_tag, 0, &negone );CHECK_ERR( rval );
 
 const void* tag_vals[] = { &sideset_id };
 Range side_sets;
 rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &ss_tag, tag_vals, 1, side_sets );CHECK_ERR( rval );
 CHECK_EQUAL( 1, (int)side_sets.size() );
 return side_sets.front();
}
 
// 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_impl1, mb_impl2;
 Interface &mb1 = mb_impl1, &mb2 = mb_impl2;
 int ho_flags[4];
 CN::HasMidNodes( type, num_nodes, ho_flags );
 
 // read file
 read_file( mb1, ho_file );
 // test element connectivity order
 EntityHandle block = find_block( mb1, type, ho_flags );
 CHECK( block != 0 );
 check_ho_elements( mb1, block, type, ho_flags );
 
 // write block and read it back in
 write_and_read( mb1, mb2, block );
 // test element connectivity order on re-read data
 block = find_block( mb2, type, ho_flags );
 CHECK( block != 0 );
 check_ho_elements( mb2, block, type, ho_flags );
}
 
void test_read_side( int id, EntityType sideset_type, int sideset_nodes_per_elem, bool shell_side )
{
 // read test file
 Core mb_impl;
 Interface& moab = mb_impl;
 read_file( moab, ho_file );
 
 // get side set
 EntityHandle set = find_sideset( moab, id, sideset_type );
 CHECK( set != 0 );
 
 // check expected element connectivity
 int ho_flags[4];
 CN::HasMidNodes( sideset_type, sideset_nodes_per_elem, ho_flags );
 check_ho_elements( moab, set, sideset_type, ho_flags );
 
 if( shell_side ) return;
 
 // check that each element is on the boundary of the mesh
 Range elems;
 ErrorCode rval = mb_impl.get_entities_by_handle( set, elems );CHECK_ERR( rval );
 
 int dim = CN::Dimension( sideset_type );
 for( Range::iterator i = elems.begin(); i != elems.end(); ++i )
 {
 Range adj;
 rval = mb_impl.get_adjacencies( &*i, 1, dim + 1, false, adj, Interface::UNION );CHECK_ERR( rval );
 CHECK_EQUAL( 1, (int)adj.size() );
 }
}
 
void test_read_ids_common( const char* file_name, const char* tag_name, const int* expected_vals, int num_expected )
{
 Core mb;
 ReadNCDF reader( &mb );
 
 FileOptions opts( "" );
 std::vector< int > values;
 ErrorCode rval = reader.read_tag_values( file_name, tag_name, opts, values );CHECK_ERR( rval );
 
 std::vector< int > expected( expected_vals, expected_vals + num_expected );
 std::sort( values.begin(), values.end() );
 std::sort( expected.begin(), expected.end() );
 CHECK_EQUAL( expected, values );
}
 
void test_read_block_ids()
{
 const int expected[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 };
 test_read_ids_common( ho_file, MATERIAL_SET_TAG_NAME, expected, sizeof( expected ) / sizeof( expected[0] ) );
}
 
void test_read_sideset_ids()
{
 const int expected[] = { 1, 2, 3, 4 };
 test_read_ids_common( ho_file, NEUMANN_SET_TAG_NAME, expected, sizeof( expected ) / sizeof( expected[0] ) );
}
 
void test_read_nodeset_ids()
{
 test_read_ids_common( ho_file, DIRICHLET_SET_TAG_NAME, 0, 0 );
}
 
void test_read_alternate_coord_format()
{
 Core moab;
 Interface& mb = moab;
 ErrorCode rval = mb.load_file( alt_file );CHECK_ERR( rval );
 
 const double exp[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1 };
 
 Range hexes;
 rval = mb.get_entities_by_type( 0, MBHEX, hexes );CHECK_ERR( rval );
 CHECK_EQUAL( (size_t)1, hexes.size() );
 EntityHandle hex = hexes.front();
 const EntityHandle* conn;
 int len;
 rval = mb.get_connectivity( hex, conn, len );CHECK_ERR( rval );
 CHECK_EQUAL( 8, len );
 double act[3 * 8];
 rval = mb.get_coords( conn, len, act );CHECK_ERR( rval );
 CHECK_ARRAYS_EQUAL( exp, 3 * 8, act, 3 * 8 );
}
void test_write_polygons()
{
 Core moab;
 Interface& mb = moab;
 ErrorCode rval = mb.load_file( polyg );CHECK_ERR( rval );
 rval = mb.write_file( "polyg.exo" );CHECK_ERR( rval );
}
 
void test_write_polyhedra()
{
 Core moab;
 Interface& mb = moab;
 ErrorCode rval = mb.load_file( polyh );CHECK_ERR( rval );
 rval = mb.write_file( "polyh.exo" );CHECK_ERR( rval );
}
void test_read_polygons()
{
 Core moab;
 Interface& mb = moab;
 ErrorCode rval = mb.load_file( rpolyg );CHECK_ERR( rval );
}
void test_read_polyhedra()
{
 Core moab;
 Interface& mb = moab;
 ErrorCode rval = mb.load_file( rpolyh );CHECK_ERR( rval );
}