MOAB: Mesh Oriented datABase  (version 5.3.0)
exodus_test.cpp
Go to the documentation of this file.
00001 #include "TestUtil.hpp"
00002 #include "moab/Core.hpp"
00003 #include "MBTagConventions.hpp"
00004 #include "moab/CN.hpp"
00005 #include "moab/Range.hpp"
00006 #include "ReadNCDF.hpp"
00007 #include "moab/FileOptions.hpp"
00008 #define IS_BUILDING_MB
00009 #include "ExoIIUtil.hpp"
00010 #include <cmath>
00011 #include <algorithm>
00012 
00013 using namespace moab;
00014 
00015 /* Input test file: ho_test.g
00016  *
00017  * File is expected to contain at least one block for every
00018  * supported higher-order element type.  The coordinates of
00019  * every higher-order node are expected to be the mean of the
00020  * adjacent corner vertices of the element.
00021  */
00022 #ifdef MESHDIR
00023 static const char ho_file[]  = STRINGIFY( MESHDIR ) "/io/ho_test.g";
00024 static const char file_one[] = STRINGIFY( MESHDIR ) "/mbtest1.g";
00025 static const char alt_file[] = STRINGIFY( MESHDIR ) "/io/hex_2x2x2_ss.exo";
00026 static const char polyg[]    = STRINGIFY( MESHDIR ) "/io/poly8-10.vtk";
00027 static const char polyh[]    = STRINGIFY( MESHDIR ) "/io/polyhedra.vtk";
00028 static const char rpolyg[]   = STRINGIFY( MESHDIR ) "/io/polyg.exo";
00029 static const char rpolyh[]   = STRINGIFY( MESHDIR ) "/io/polyh.exo";
00030 #else
00031 static const char ho_file[]  = "ho_test.g";
00032 static const char file_one[] = "mbtest1.g";
00033 static const char alt_file[] = "hex_2x2x2_ss.exo";
00034 static const char polyg[]    = "poly8-10.vtk";
00035 static const char polyh[]    = "polyhedra.vtk";
00036 static const char rpolyg[]   = "polyg.exo";
00037 static const char rpolyh[]   = "polyh.exo";
00038 #endif
00039 
00040 void read_file( Interface& moab, const char* input_file );
00041 
00042 // Check that element has expected higher-order nodes
00043 // and that each higher-order node is at the center
00044 // of the sub-entity it is on.
00045 void check_ho_element( Interface& moab, EntityHandle entity, int mid_nodes[4] );
00046 
00047 void test_read_side( int sideset_id, EntityType sideset_type, int sideset_nodes_per_elem, bool shell_side = false );
00048 
00049 // Validate elements of specified type.
00050 // Looks for a block containing the specified entity type
00051 // and with the specified mid-node flags set in its
00052 // HAS_MID_NODES_TAG.
00053 void test_ho_elements( EntityType type, int num_nodes );
00054 
00055 // Tests originally in MBTest.cpp
00056 void mb_vertex_coordinate_test();
00057 void mb_bar_connectivity_test();
00058 void mb_tri_connectivity_test();
00059 void mb_quad_connectivity_test();
00060 void mb_hex_connectivity_test();
00061 void mb_tet_connectivity_test();
00062 void mb_write_mesh_test();
00063 
00064 void test_types();
00065 
00066 void test_tri6()
00067 {
00068     test_ho_elements( MBTRI, 6 );
00069 }
00070 void test_tri7()
00071 {
00072     test_ho_elements( MBTRI, 7 );
00073 }
00074 
00075 void test_quad5()
00076 {
00077     test_ho_elements( MBQUAD, 5 );
00078 }
00079 void test_quad8()
00080 {
00081     test_ho_elements( MBQUAD, 8 );
00082 }
00083 void test_quad9()
00084 {
00085     test_ho_elements( MBQUAD, 9 );
00086 }
00087 
00088 void test_tet8()
00089 {
00090     test_ho_elements( MBTET, 8 );
00091 }
00092 void test_tet10()
00093 {
00094     test_ho_elements( MBTET, 10 );
00095 }
00096 void test_tet14()
00097 {
00098     test_ho_elements( MBTET, 14 );
00099 }
00100 
00101 void test_hex9()
00102 {
00103     test_ho_elements( MBHEX, 9 );
00104 }
00105 void test_hex20()
00106 {
00107     test_ho_elements( MBHEX, 20 );
00108 }
00109 void test_hex27()
00110 {
00111     test_ho_elements( MBHEX, 27 );
00112 }
00113 
00114 void test_read_tri6_side()
00115 {
00116     test_read_side( 1, MBEDGE, 3 );
00117 }  // sideset 1
00118 void test_read_shell_side()
00119 {
00120     test_read_side( 3, MBQUAD, 9, true );
00121 }  // sideset 3
00122 void test_read_shell_edge()
00123 {
00124     test_read_side( 4, MBEDGE, 3 );
00125 }  // sideset 4
00126 void test_read_hex20_side()
00127 {
00128     test_read_side( 2, MBQUAD, 8 );
00129 }  // sideset 2
00130 
00131 void test_read_block_ids();
00132 void test_read_sideset_ids();
00133 void test_read_nodeset_ids();
00134 void test_write_polygons();
00135 void test_write_polyhedra();
00136 void test_read_polygons();
00137 void test_read_polyhedra();
00138 
00139 void test_read_alternate_coord_format();
00140 
00141 int main()
00142 {
00143     int result = 0;
00144 
00145     result += RUN_TEST( mb_vertex_coordinate_test );
00146     result += RUN_TEST( mb_bar_connectivity_test );
00147     result += RUN_TEST( mb_tri_connectivity_test );
00148     result += RUN_TEST( mb_quad_connectivity_test );
00149     result += RUN_TEST( mb_hex_connectivity_test );
00150     result += RUN_TEST( mb_tet_connectivity_test );
00151     result += RUN_TEST( mb_write_mesh_test );
00152 
00153     result += RUN_TEST( test_types );
00154 
00155     result += RUN_TEST( test_tri6 );
00156     result += RUN_TEST( test_tri7 );
00157     result += RUN_TEST( test_quad5 );
00158     result += RUN_TEST( test_quad8 );
00159     result += RUN_TEST( test_quad9 );
00160     result += RUN_TEST( test_tet8 );
00161     result += RUN_TEST( test_tet10 );
00162     result += RUN_TEST( test_tet14 );
00163     result += RUN_TEST( test_hex9 );
00164     result += RUN_TEST( test_hex20 );
00165     result += RUN_TEST( test_hex27 );
00166 
00167     result += RUN_TEST( test_read_tri6_side );
00168     result += RUN_TEST( test_read_shell_side );
00169     result += RUN_TEST( test_read_shell_edge );
00170     result += RUN_TEST( test_read_hex20_side );
00171 
00172     result += RUN_TEST( test_read_block_ids );
00173     result += RUN_TEST( test_read_sideset_ids );
00174     result += RUN_TEST( test_read_nodeset_ids );
00175 
00176     result += RUN_TEST( test_read_alternate_coord_format );
00177 
00178     result += RUN_TEST( test_write_polygons );
00179     result += RUN_TEST( test_write_polyhedra );
00180     result += RUN_TEST( test_read_polygons );
00181     result += RUN_TEST( test_read_polyhedra );
00182 
00183     return result;
00184 }
00185 
00186 void load_file_one( Interface* iface )
00187 {
00188     ErrorCode error = iface->load_mesh( file_one );
00189     if( MB_SUCCESS != error )
00190     {
00191         std::cout << "Failed to load input file: " << file_one << std::endl;
00192         std::string error_reason;
00193         iface->get_last_error( error_reason );
00194         std::cout << error_reason << std::endl;
00195     }
00196     CHECK_ERR( error );
00197 }
00198 
00199 /*!
00200   @test
00201   Vertex Coordinates
00202   @li Get coordinates of vertex 1 correctly
00203   @li Get coordinates of vertex 8 correctly
00204   @li Get coordinates of vertex 6 correctly
00205 */
00206 void mb_vertex_coordinate_test()
00207 {
00208     double coords[3];
00209     EntityHandle handle;
00210     ErrorCode error;
00211     int err;
00212 
00213     Core moab;
00214     Interface* MB = &moab;
00215     load_file_one( MB );
00216 
00217     // coordinate 2 should be {1.5, -1.5, 3.5}
00218 
00219     handle = CREATE_HANDLE( MBVERTEX, 2, err );
00220     error  = MB->get_coords( &handle, 1, coords );CHECK_ERR( error );
00221     const double exp2[] = { 1.5, -1.5, 3.5 };
00222     CHECK_ARRAYS_EQUAL( exp2, 3, coords, 3 );
00223 
00224     // coordinate 9 should be {1, -2, 3.5}
00225     handle = CREATE_HANDLE( MBVERTEX, 9, err );
00226     error  = MB->get_coords( &handle, 1, coords );CHECK_ERR( error );
00227     const double exp9[] = { 1, -2, 3.5 };
00228     CHECK_ARRAYS_EQUAL( exp9, 3, coords, 3 );
00229 
00230     // coordinate 7 should be {0.5, -2, 3.5}
00231     handle = CREATE_HANDLE( MBVERTEX, 7, err );
00232     error  = MB->get_coords( &handle, 1, coords );CHECK_ERR( error );
00233     const double exp7[] = { 0.5, -2, 3.5 };
00234     CHECK_ARRAYS_EQUAL( exp7, 3, coords, 3 );
00235 
00236     int node_count = 0;
00237     error          = MB->get_number_entities_by_type( 0, MBVERTEX, node_count );CHECK_ERR( error );
00238     // Number of vertices (node_count) should be 83 assuming no gaps in the handle space
00239     CHECK_EQUAL( 47, node_count );
00240 }
00241 
00242 /*!
00243   @test
00244   MB Bar Element Connectivity Test
00245   @li Get coordinates for 2 node bar elements
00246 */
00247 
00248 void mb_bar_connectivity_test()
00249 {
00250     Core moab;
00251     Interface* MB = &moab;
00252     load_file_one( MB );
00253 
00254     std::vector< EntityHandle > conn;
00255     Range bars;
00256 
00257     ErrorCode error = MB->get_entities_by_type( 0, MBEDGE, bars );CHECK_ERR( error );
00258 
00259     // get the connectivity of the second bar
00260     EntityHandle handle = *( ++bars.begin() );
00261 
00262     error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
00263 
00264     CHECK_EQUAL( (size_t)2, conn.size() );
00265 
00266     // from ncdump the connectivity of bar 2 (0 based) is
00267     //  14, 13
00268     CHECK_EQUAL( (EntityHandle)14, conn[0] );
00269     CHECK_EQUAL( (EntityHandle)13, conn[1] );
00270 
00271     // Now try getting the connectivity of one of the vertices for fun.
00272     // just return the vertex in the connectivity
00273     handle = conn[0];
00274     error  = MB->get_connectivity( &handle, 1, conn );
00275     CHECK_EQUAL( MB_FAILURE, error );
00276 }
00277 
00278 void mb_tri_connectivity_test()
00279 {
00280     Core moab;
00281     Interface* MB = &moab;
00282     load_file_one( MB );
00283 
00284     std::vector< EntityHandle > conn;
00285     Range tris;
00286     ErrorCode error = MB->get_entities_by_type( 0, MBTRI, tris );CHECK_ERR( error );
00287 
00288     // get the connectivity of the second tri
00289     EntityHandle handle = *( ++tris.begin() );
00290 
00291     error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
00292 
00293     CHECK_EQUAL( (size_t)3, conn.size() );
00294 
00295     // from ncdump the connectivity of tri 2 (0 based) is
00296     //  45, 37, 38
00297 
00298     CHECK_EQUAL( (EntityHandle)45, conn[0] );
00299     CHECK_EQUAL( (EntityHandle)37, conn[1] );
00300     CHECK_EQUAL( (EntityHandle)38, conn[2] );
00301 }
00302 
00303 void mb_quad_connectivity_test()
00304 {
00305     Core moab;
00306     Interface* MB = &moab;
00307     load_file_one( MB );
00308 
00309     std::vector< EntityHandle > conn;
00310     Range quads;
00311 
00312     ErrorCode error = MB->get_entities_by_type( 0, MBQUAD, quads );CHECK_ERR( error );
00313 
00314     // get the connectivity of the second quad
00315     EntityHandle handle = *( ++quads.begin() );
00316 
00317     error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
00318 
00319     CHECK_EQUAL( (size_t)4, conn.size() );
00320 
00321     // from ncdump the connectivity of quad 2 (0 based) is
00322     // 20, 11, 12, 26,
00323 
00324     CHECK_EQUAL( (EntityHandle)20, conn[0] );
00325     CHECK_EQUAL( (EntityHandle)11, conn[1] );
00326     CHECK_EQUAL( (EntityHandle)12, conn[2] );
00327     CHECK_EQUAL( (EntityHandle)26, conn[3] );
00328 }
00329 
00330 void mb_hex_connectivity_test()
00331 {
00332     Core moab;
00333     Interface* MB = &moab;
00334     load_file_one( MB );
00335 
00336     std::vector< EntityHandle > conn;
00337     Range hexes;
00338 
00339     ErrorCode error = MB->get_entities_by_type( 0, MBHEX, hexes );CHECK_ERR( error );
00340 
00341     // get the connectivity of the second hex
00342     EntityHandle handle = *( ++hexes.begin() );
00343 
00344     error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
00345 
00346     CHECK_EQUAL( (size_t)8, conn.size() );
00347 
00348     // from ncdump the connectivity of hex 1 (0 based) is
00349     // 19, 13, 16, 23, 21, 14, 18, 27
00350 
00351     CHECK_EQUAL( (EntityHandle)19, conn[0] );
00352     CHECK_EQUAL( (EntityHandle)13, conn[1] );
00353     CHECK_EQUAL( (EntityHandle)16, conn[2] );
00354     CHECK_EQUAL( (EntityHandle)23, conn[3] );
00355     CHECK_EQUAL( (EntityHandle)21, conn[4] );
00356     CHECK_EQUAL( (EntityHandle)14, conn[5] );
00357     CHECK_EQUAL( (EntityHandle)18, conn[6] );
00358     CHECK_EQUAL( (EntityHandle)27, conn[7] );
00359 }
00360 
00361 void mb_tet_connectivity_test()
00362 {
00363     Core moab;
00364     Interface* MB = &moab;
00365     load_file_one( MB );
00366 
00367     std::vector< EntityHandle > conn;
00368     Range tets;
00369     ErrorCode error = MB->get_entities_by_type( 0, MBTET, tets );CHECK_ERR( error );
00370 
00371     // get the connectivity of the second tet
00372     EntityHandle handle = *( ++tets.begin() );
00373 
00374     error = MB->get_connectivity( &handle, 1, conn );CHECK_ERR( error );
00375 
00376     CHECK_EQUAL( (size_t)4, conn.size() );
00377 
00378     // from ncdump the connectivity of tet 2 (0 based) is:
00379     // 35, 34, 32, 43
00380 
00381     CHECK_EQUAL( (EntityHandle)35, conn[0] );
00382     CHECK_EQUAL( (EntityHandle)34, conn[1] );
00383     CHECK_EQUAL( (EntityHandle)32, conn[2] );
00384     CHECK_EQUAL( (EntityHandle)43, conn[3] );
00385 }
00386 
00387 void mb_write_mesh_test()
00388 {
00389     Core moab;
00390     Interface* MB = &moab;
00391     load_file_one( MB );
00392     ErrorCode result;
00393 
00394     std::string file_name = "mb_write.g";
00395 
00396     // no need to get lists, write out the whole mesh
00397     result = MB->write_mesh( file_name.c_str() );CHECK_ERR( result );
00398 
00399     //---------The following tests outputting meshsets that are in meshsets of blocks ---/
00400 
00401     // lets create a block meshset and put some entities and meshsets into it
00402     EntityHandle block_ms;
00403     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, block_ms );CHECK_ERR( result );
00404 
00405     // make another meshset to put quads in, so SHELLs can be written out
00406     EntityHandle block_of_shells;
00407     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, block_of_shells );CHECK_ERR( result );
00408 
00409     // tag the meshset so it's a block, with id 100
00410     int id = 100;
00411     Tag tag_handle;
00412     const int negone = -1;
00413     result           = MB->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
00414     result = MB->tag_set_data( tag_handle, &block_ms, 1, &id );CHECK_ERR( result );
00415     id     = 101;
00416     result = MB->tag_set_data( tag_handle, &block_of_shells, 1, &id );CHECK_ERR( result );
00417 
00418     // set dimension tag on this to ensure shells get output; reuse id variable
00419     result = MB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
00420     id     = 3;
00421     result = MB->tag_set_data( tag_handle, &block_of_shells, 1, &id );CHECK_ERR( result );
00422 
00423     // get some entities (tets)
00424     Range temp_range;
00425     result = MB->get_entities_by_type( 0, MBHEX, temp_range );CHECK_ERR( result );
00426 
00427     Range::iterator iter, end_iter;
00428     iter     = temp_range.begin();
00429     end_iter = temp_range.end();
00430 
00431     // add evens to 'block_ms'
00432     std::vector< EntityHandle > temp_vec;
00433     for( ; iter != end_iter; ++iter )
00434     {
00435         if( ID_FROM_HANDLE( *iter ) % 2 == 0 ) temp_vec.push_back( *iter );
00436     }
00437     result = MB->add_entities( block_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
00438 
00439     // make another meshset
00440     EntityHandle ms_of_block_ms;
00441     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ms_of_block_ms );CHECK_ERR( result );
00442 
00443     // add some entities to it
00444     temp_vec.clear();
00445     iter = temp_range.begin();
00446     for( ; iter != end_iter; ++iter )
00447     {
00448         if( ID_FROM_HANDLE( *iter ) % 2 )  // add all odds
00449             temp_vec.push_back( *iter );
00450     }
00451     result = MB->add_entities( ms_of_block_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
00452 
00453     // add the other meshset to the block's meshset
00454     result = MB->add_entities( block_ms, &ms_of_block_ms, 1 );CHECK_ERR( result );
00455 
00456     //---------------testing sidesets----------------/
00457 
00458     // lets create a sideset meshset and put some entities and meshsets into it
00459     EntityHandle sideset_ms;
00460     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, sideset_ms );CHECK_ERR( result );
00461 
00462     // tag the meshset so it's a sideset, with id 104
00463     id     = 104;
00464     result = MB->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
00465 
00466     result = MB->tag_set_data( tag_handle, &sideset_ms, 1, &id );CHECK_ERR( result );
00467 
00468     // get some entities (tris)
00469     temp_range.clear();
00470     result = MB->get_entities_by_type( 0, MBQUAD, temp_range );CHECK_ERR( result );
00471 
00472     iter     = temp_range.begin();
00473     end_iter = temp_range.end();
00474 
00475     // add evens to 'sideset_ms'
00476     temp_vec.clear();
00477     for( ; iter != end_iter; ++iter )
00478     {
00479         if( ID_FROM_HANDLE( *iter ) % 2 == 0 ) temp_vec.push_back( *iter );
00480     }
00481     result = MB->add_entities( sideset_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
00482 
00483     // make another meshset
00484     EntityHandle ms_of_sideset_ms;
00485     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ms_of_sideset_ms );CHECK_ERR( result );
00486 
00487     // add some entities to it
00488     temp_vec.clear();
00489     iter = temp_range.begin();
00490     for( ; iter != end_iter; ++iter )
00491     {
00492         if( ID_FROM_HANDLE( *iter ) % 2 )  // add all odds
00493             temp_vec.push_back( *iter );
00494     }
00495     result = MB->add_entities( ms_of_sideset_ms, &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
00496 
00497     // add the other meshset to the sideset's meshset
00498     result = MB->add_entities( sideset_ms, &ms_of_sideset_ms, 1 );CHECK_ERR( result );
00499 
00500     //---------test sense on meshsets (reverse/foward)-------//
00501 
00502     // get all quads whose x-coord = 2.5 and put them into a meshset_a
00503     EntityHandle meshset_a;
00504     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_a );CHECK_ERR( result );
00505 
00506     temp_range.clear();
00507     result = MB->get_entities_by_type( 0, MBQUAD, temp_range );CHECK_ERR( result );
00508 
00509     std::vector< EntityHandle > nodes1, entity_vec;
00510     std::copy( temp_range.begin(), temp_range.end(), std::back_inserter( entity_vec ) );
00511     result = MB->get_connectivity( &entity_vec[0], entity_vec.size(), nodes1 );CHECK_ERR( result );
00512     assert( nodes1.size() == 4 * temp_range.size() );
00513     temp_vec.clear();
00514     std::vector< double > coords( 3 * nodes1.size() );
00515     result = MB->get_coords( &nodes1[0], nodes1.size(), &coords[0] );CHECK_ERR( result );
00516 
00517     unsigned int k = 0;
00518     for( Range::iterator it = temp_range.begin(); it != temp_range.end(); ++it )
00519     {
00520         if( coords[12 * k] == 2.5 && coords[12 * k + 3] == 2.5 && coords[12 * k + 6] == 2.5 &&
00521             coords[12 * k + 9] == 2.5 )
00522             temp_vec.push_back( *it );
00523         k++;
00524     }
00525     result = MB->add_entities( meshset_a, ( temp_vec.empty() ) ? NULL : &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
00526     result = MB->add_entities( block_of_shells, ( temp_vec.empty() ) ? NULL : &temp_vec[0], temp_vec.size() );CHECK_ERR( result );
00527 
00528     // put these quads into a different meshset_b and tag them with a reverse sense tag
00529     EntityHandle meshset_b;
00530     result = MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_b );CHECK_ERR( result );
00531 
00532     result = MB->add_entities( meshset_b, &meshset_a, 1 );CHECK_ERR( result );
00533 
00534     result = MB->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, tag_handle, MB_TAG_SPARSE | MB_TAG_CREAT );CHECK_ERR( result );
00535 
00536     int reverse_value = -1;
00537     result            = MB->tag_set_data( tag_handle, &meshset_b, 1, &reverse_value );CHECK_ERR( result );
00538 
00539     // get some random quad, whose x-coord != 2.5, and put it into a different meshset_c
00540     // and tag it with a reverse sense tag
00541 
00542     iter     = temp_range.begin();
00543     end_iter = temp_range.end();
00544 
00545     temp_vec.clear();
00546     for( ; iter != end_iter; ++iter )
00547     {
00548         std::vector< EntityHandle > nodes;
00549         result = MB->get_connectivity( &( *iter ), 1, nodes );CHECK_ERR( result );
00550 
00551         bool not_equal_2_5 = true;
00552         for( unsigned int ku = 0; ku < nodes.size(); ku++ )
00553         {
00554             double coords2[3] = { 0 };
00555 
00556             result = MB->get_coords( &( nodes[ku] ), 1, coords2 );CHECK_ERR( result );
00557 
00558             if( coords2[0] == 2.5 )
00559             {
00560                 not_equal_2_5 = false;
00561                 break;
00562             }
00563         }
00564 
00565         if( not_equal_2_5 && nodes.size() > 0 )
00566         {
00567             temp_vec.push_back( *iter );
00568             break;
00569         }
00570     }
00571 
00572     EntityHandle meshset_c;
00573     MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_c );
00574 
00575     result = MB->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, tag_handle );CHECK_ERR( result );
00576 
00577     reverse_value = -1;
00578     result        = MB->tag_set_data( tag_handle, &meshset_c, 1, &reverse_value );CHECK_ERR( result );
00579 
00580     MB->add_entities( meshset_c, &temp_vec[0], temp_vec.size() );
00581     MB->add_entities( block_of_shells, &temp_vec[0], temp_vec.size() );
00582 
00583     // create another meshset_abc, adding meshset_a, meshset_b, meshset_c
00584     EntityHandle meshset_abc;
00585     MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, meshset_abc );
00586 
00587     temp_vec.clear();
00588     temp_vec.push_back( meshset_a );
00589     temp_vec.push_back( meshset_b );
00590     temp_vec.push_back( meshset_c );
00591 
00592     MB->add_entities( meshset_abc, &temp_vec[0], temp_vec.size() );
00593 
00594     // tag it so it's a sideset
00595     id     = 444;
00596     result = MB->tag_get_handle( "NEUMANN_SET", 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
00597 
00598     result = MB->tag_set_data( tag_handle, &meshset_abc, 1, &id );CHECK_ERR( result );
00599 
00600     //---------------do nodesets now -----------------//
00601 
00602     // lets create a nodeset meshset and put some entities and meshsets into it
00603     EntityHandle nodeset_ms;
00604     MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, nodeset_ms );
00605 
00606     // tag the meshset so it's a nodeset, with id 119
00607     id     = 119;
00608     result = MB->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, tag_handle, 0, &negone );CHECK_ERR( result );
00609 
00610     result = MB->tag_set_data( tag_handle, &nodeset_ms, 1, &id );CHECK_ERR( result );
00611 
00612     // get all Quads
00613     temp_range.clear();
00614     result = MB->get_entities_by_type( 0, MBQUAD, temp_range );CHECK_ERR( result );
00615 
00616     // get all the nodes of the tris
00617     Range nodes_of_quads;
00618     iter     = temp_range.begin();
00619     end_iter = temp_range.end();
00620 
00621     for( ; iter != end_iter; ++iter )
00622     {
00623         std::vector< EntityHandle > nodes;
00624         result = MB->get_connectivity( &( *iter ), 1, nodes );CHECK_ERR( result );
00625 
00626         for( unsigned int ku = 0; ku < nodes.size(); ku++ )
00627             nodes_of_quads.insert( nodes[ku] );
00628     }
00629 
00630     iter     = nodes_of_quads.begin();
00631     end_iter = nodes_of_quads.end();
00632 
00633     // add evens to 'nodeset_ms'
00634     temp_vec.clear();
00635     for( ; iter != end_iter; ++iter )
00636     {
00637         if( ID_FROM_HANDLE( *iter ) % 2 == 0 ) temp_vec.push_back( *iter );
00638     }
00639     MB->add_entities( nodeset_ms, &temp_vec[0], temp_vec.size() );
00640 
00641     // make another meshset
00642     EntityHandle ms_of_nodeset_ms;
00643     MB->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ms_of_nodeset_ms );
00644 
00645     // add some entities to it
00646     temp_vec.clear();
00647     iter     = nodes_of_quads.begin();
00648     end_iter = nodes_of_quads.end();
00649     for( ; iter != end_iter; ++iter )
00650     {
00651         if( ID_FROM_HANDLE( *iter ) % 2 )  // add all odds
00652             temp_vec.push_back( *iter );
00653     }
00654     MB->add_entities( ms_of_nodeset_ms, &temp_vec[0], temp_vec.size() );
00655 
00656     // add the other meshset to the nodeset's meshset
00657     MB->add_entities( nodeset_ms, &ms_of_nodeset_ms, 1 );
00658 
00659     // no need to get lists, write out the whole mesh
00660     file_name = "mb_write2.g";
00661     std::vector< EntityHandle > output_list;
00662     output_list.push_back( block_ms );
00663     output_list.push_back( sideset_ms );
00664     output_list.push_back( meshset_abc );
00665     output_list.push_back( nodeset_ms );
00666     output_list.push_back( block_of_shells );
00667     result = MB->write_mesh( file_name.c_str(), &output_list[0], output_list.size() );CHECK_ERR( result );
00668 }
00669 
00670 struct TestType
00671 {
00672     EntityType moab_type;
00673     ExoIIElementType exo_type;
00674     int num_nodes;
00675     std::string name;
00676 };
00677 
00678 void check_type( const TestType& type )
00679 {
00680     int has_mid_nodes[4];
00681     CN::HasMidNodes( type.moab_type, type.num_nodes, has_mid_nodes );
00682 
00683     CHECK_EQUAL( type.moab_type, ExoIIUtil::ExoIIElementMBEntity[type.exo_type] );
00684     CHECK_EQUAL( type.name, std::string( ExoIIUtil::ElementTypeNames[type.exo_type] ) );
00685     CHECK_EQUAL( type.num_nodes, ExoIIUtil::VerticesPerElement[type.exo_type] );
00686     switch( CN::Dimension( type.moab_type ) )
00687     {
00688         case 3:
00689             CHECK_EQUAL( has_mid_nodes[3], ExoIIUtil::HasMidNodes[type.exo_type][3] );
00690         case 2:
00691             CHECK_EQUAL( has_mid_nodes[2], ExoIIUtil::HasMidNodes[type.exo_type][2] );
00692         case 1:
00693             CHECK_EQUAL( has_mid_nodes[1], ExoIIUtil::HasMidNodes[type.exo_type][1] );
00694     }
00695 
00696     Core moab;
00697     ExoIIUtil tool( &moab );
00698     CHECK_EQUAL( type.exo_type, tool.element_name_to_type( type.name.c_str() ) );
00699     CHECK_EQUAL( type.name, std::string( tool.element_type_name( type.exo_type ) ) );
00700 }
00701 
00702 void test_types()
00703 {
00704     const TestType types[] = { { MBVERTEX, EXOII_SPHERE, 1, "SPHERE" },
00705                                { MBEDGE, EXOII_SPRING, 1, "SPRING" },
00706                                { MBEDGE, EXOII_BAR, 2, "BAR" },
00707                                { MBEDGE, EXOII_BAR2, 2, "BAR2" },
00708                                { MBEDGE, EXOII_BAR3, 3, "BAR3" },
00709                                { MBEDGE, EXOII_BEAM, 2, "BEAM" },
00710                                { MBEDGE, EXOII_BEAM2, 2, "BEAM2" },
00711                                { MBEDGE, EXOII_BEAM3, 3, "BEAM3" },
00712                                { MBEDGE, EXOII_TRUSS, 2, "TRUSS" },
00713                                { MBEDGE, EXOII_TRUSS2, 2, "TRUSS2" },
00714                                { MBEDGE, EXOII_TRUSS3, 3, "TRUSS3" },
00715                                { MBTRI, EXOII_TRI, 3, "TRI" },
00716                                { MBTRI, EXOII_TRI3, 3, "TRI3" },
00717                                { MBTRI, EXOII_TRI6, 6, "TRI6" },
00718                                { MBTRI, EXOII_TRI7, 7, "TRI7" },
00719                                { MBQUAD, EXOII_QUAD, 4, "QUAD" },
00720                                { MBQUAD, EXOII_QUAD4, 4, "QUAD4" },
00721                                { MBQUAD, EXOII_QUAD5, 5, "QUAD5" },
00722                                { MBQUAD, EXOII_QUAD8, 8, "QUAD8" },
00723                                { MBQUAD, EXOII_QUAD9, 9, "QUAD9" },
00724                                { MBQUAD, EXOII_SHELL, 4, "SHELL" },
00725                                { MBQUAD, EXOII_SHELL4, 4, "SHELL4" },
00726                                { MBQUAD, EXOII_SHELL5, 5, "SHELL5" },
00727                                { MBQUAD, EXOII_SHELL8, 8, "SHELL8" },
00728                                { MBQUAD, EXOII_SHELL9, 9, "SHELL9" },
00729                                { MBTET, EXOII_TETRA, 4, "TETRA" },
00730                                { MBTET, EXOII_TETRA4, 4, "TETRA4" },
00731                                { MBTET, EXOII_TET4, 4, "TET4" },
00732                                { MBTET, EXOII_TETRA8, 8, "TETRA8" },
00733                                { MBTET, EXOII_TETRA10, 10, "TETRA10" },
00734                                { MBTET, EXOII_TETRA14, 14, "TETRA14" },
00735                                { MBPYRAMID, EXOII_PYRAMID, 5, "PYRAMID" },
00736                                { MBPYRAMID, EXOII_PYRAMID5, 5, "PYRAMID5" },
00737                                { MBPYRAMID, EXOII_PYRAMID10, 10, "PYRAMID10" },
00738                                { MBPYRAMID, EXOII_PYRAMID13, 13, "PYRAMID13" },
00739                                { MBPYRAMID, EXOII_PYRAMID18, 18, "PYRAMID18" },
00740                                { MBPRISM, EXOII_WEDGE, 6, "WEDGE" },
00741                                { MBKNIFE, EXOII_KNIFE, 7, "KNIFE" },
00742                                { MBHEX, EXOII_HEX, 8, "HEX" },
00743                                { MBHEX, EXOII_HEX8, 8, "HEX8" },
00744                                { MBHEX, EXOII_HEX9, 9, "HEX9" },
00745                                { MBHEX, EXOII_HEX20, 20, "HEX20" },
00746                                { MBHEX, EXOII_HEX27, 27, "HEX27" },
00747                                { MBHEX, EXOII_HEXSHELL, 12, "HEXSHELL" } };
00748     const int num_types    = sizeof( types ) / sizeof( types[0] );
00749     for( int i = 0; i < num_types; ++i )
00750         check_type( types[i] );
00751 }
00752 
00753 void read_file( Interface& moab, const char* input_file )
00754 {
00755     ErrorCode rval = moab.load_file( input_file );CHECK_ERR( rval );
00756 }
00757 
00758 void write_and_read( Interface& write_mb, Interface& read_mb, EntityHandle block = 0 )
00759 {
00760     const char* tmp_file = "exodus_test_tmp.g";
00761     ErrorCode rval;
00762 
00763     EntityHandle* write_set_list = &block;
00764     int write_set_list_len       = 0;  //(block != 0);
00765     rval                         = write_mb.write_file( tmp_file, "EXODUS", 0, write_set_list, write_set_list_len );
00766     if( MB_SUCCESS != rval ) remove( tmp_file );CHECK_ERR( rval );
00767 
00768     rval = read_mb.load_file( tmp_file );
00769     remove( tmp_file );CHECK_ERR( rval );
00770 }
00771 
00772 void check_ho_elements( Interface& moab, EntityHandle block, EntityType type, int mid_nodes[4] )
00773 {
00774     ErrorCode rval;
00775     Range elems;
00776     rval = moab.get_entities_by_handle( block, elems, ( type != MBENTITYSET ? true : false ) );CHECK_ERR( rval );
00777     CHECK( !elems.empty() );
00778     CHECK( elems.all_of_type( type ) );
00779     for( Range::const_iterator i = elems.begin(); i != elems.end(); ++i )
00780         check_ho_element( moab, *i, mid_nodes );
00781 }
00782 
00783 // Check that element has expected higher-order nodes
00784 // and that each higher-order node is at the center
00785 // of the sub-entity it is on.
00786 void check_ho_element( Interface& moab, EntityHandle entity, int mid_nodes[4] )
00787 {
00788     // get element info
00789     const EntityType type = TYPE_FROM_HANDLE( entity );
00790     const EntityHandle* conn;
00791     int conn_len;
00792     ErrorCode rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );
00793     std::vector< double > coords( 3 * conn_len );
00794     rval = moab.get_coords( conn, conn_len, &coords[0] );CHECK_ERR( rval );
00795 
00796     // calculate and verify expected number of mid nodes
00797     int num_nodes = CN::VerticesPerEntity( type );
00798     for( int d = 1; d <= CN::Dimension( type ); ++d )
00799         if( mid_nodes[d] ) num_nodes += CN::NumSubEntities( type, d );
00800     CHECK_EQUAL( num_nodes, conn_len );
00801 
00802     // verify that each higher-order node is at the center
00803     // of its respective sub-entity.
00804     for( int i = CN::VerticesPerEntity( type ); i < num_nodes; ++i )
00805     {
00806         // get sub-entity owning ho-node
00807         int sub_dim, sub_num;
00808         CN::HONodeParent( type, num_nodes, i, sub_dim, sub_num );
00809         // get corner vertex indices
00810         int sub_conn[8], num_sub;
00811         if( sub_dim < CN::Dimension( type ) )
00812         {
00813             CN::SubEntityVertexIndices( type, sub_dim, sub_num, sub_conn );
00814             EntityType sub_type = CN::SubEntityType( type, sub_dim, sub_num );
00815             num_sub             = CN::VerticesPerEntity( sub_type );
00816         }
00817         else
00818         {
00819             num_sub = CN::VerticesPerEntity( type );
00820             for( int j = 0; j < num_sub; ++j )
00821                 sub_conn[j] = j;
00822         }
00823         // calculate mean of corner vertices
00824         double mean[3] = { 0, 0, 0 };
00825         for( int j = 0; j < num_sub; ++j )
00826         {
00827             int co = 3 * sub_conn[j];
00828             mean[0] += coords[co];
00829             mean[1] += coords[co + 1];
00830             mean[2] += coords[co + 2];
00831         }
00832         mean[0] /= num_sub;
00833         mean[1] /= num_sub;
00834         mean[2] /= num_sub;
00835         // verify that higher-order node is at expected location
00836         CHECK_REAL_EQUAL( mean[0], coords[3 * i], 1e-6 );
00837         CHECK_REAL_EQUAL( mean[1], coords[3 * i + 1], 1e-6 );
00838         CHECK_REAL_EQUAL( mean[2], coords[3 * i + 2], 1e-6 );
00839     }
00840 }
00841 
00842 EntityHandle find_block( Interface& mb, EntityType type, const int has_mid_nodes[4] )
00843 {
00844 
00845     ErrorCode rval;
00846     Tag ho_tag, block_tag;
00847     const int negone = -1;
00848     rval             = mb.tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag, 0, &negone );CHECK_ERR( rval );
00849     const int mids[] = { -1, -1, -1, -1 };
00850     rval             = mb.tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, ho_tag, 0, mids );CHECK_ERR( rval );
00851 
00852     // get material sets with expected higher-order nodes
00853     Range blocks;
00854     Tag tags[2]         = { ho_tag, block_tag };
00855     const void* vals[2] = { has_mid_nodes, NULL };
00856     rval                = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, blocks );CHECK_ERR( rval );
00857 
00858     for( Range::iterator i = blocks.begin(); i != blocks.end(); ++i )
00859     {
00860         int n;
00861         rval = mb.get_number_entities_by_type( *i, type, n );CHECK_ERR( rval );
00862         if( n > 0 ) return *i;
00863     }
00864 
00865     CHECK( false );  // no block matching element type description
00866     return 0;
00867 }
00868 
00869 EntityHandle find_sideset( Interface& mb, int sideset_id, EntityType /*side_type*/ )
00870 {
00871     ErrorCode rval;
00872     Tag ss_tag;
00873     const int negone = -1;
00874     rval             = mb.tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ss_tag, 0, &negone );CHECK_ERR( rval );
00875 
00876     const void* tag_vals[] = { &sideset_id };
00877     Range side_sets;
00878     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &ss_tag, tag_vals, 1, side_sets );CHECK_ERR( rval );
00879     CHECK_EQUAL( 1, (int)side_sets.size() );
00880     return side_sets.front();
00881 }
00882 
00883 // Validate elements of specified type.
00884 // Looks for a block containing the specified entity type
00885 // and with the specified mid-node flags set in its
00886 // HAS_MID_NODES_TAG.
00887 void test_ho_elements( EntityType type, int num_nodes )
00888 {
00889     Core mb_impl1, mb_impl2;
00890     Interface &mb1 = mb_impl1, &mb2 = mb_impl2;
00891     int ho_flags[4];
00892     CN::HasMidNodes( type, num_nodes, ho_flags );
00893 
00894     // read file
00895     read_file( mb1, ho_file );
00896     // test element connectivity order
00897     EntityHandle block = find_block( mb1, type, ho_flags );
00898     CHECK( block != 0 );
00899     check_ho_elements( mb1, block, type, ho_flags );
00900 
00901     // write block and read it back in
00902     write_and_read( mb1, mb2, block );
00903     // test element connectivity order on re-read data
00904     block = find_block( mb2, type, ho_flags );
00905     CHECK( block != 0 );
00906     check_ho_elements( mb2, block, type, ho_flags );
00907 }
00908 
00909 void test_read_side( int id, EntityType sideset_type, int sideset_nodes_per_elem, bool shell_side )
00910 {
00911     // read test file
00912     Core mb_impl;
00913     Interface& moab = mb_impl;
00914     read_file( moab, ho_file );
00915 
00916     // get side set
00917     EntityHandle set = find_sideset( moab, id, sideset_type );
00918     CHECK( set != 0 );
00919 
00920     // check expected element connectivity
00921     int ho_flags[4];
00922     CN::HasMidNodes( sideset_type, sideset_nodes_per_elem, ho_flags );
00923     check_ho_elements( moab, set, sideset_type, ho_flags );
00924 
00925     if( shell_side ) return;
00926 
00927     // check that each element is on the boundary of the mesh
00928     Range elems;
00929     ErrorCode rval = mb_impl.get_entities_by_handle( set, elems );CHECK_ERR( rval );
00930 
00931     int dim = CN::Dimension( sideset_type );
00932     for( Range::iterator i = elems.begin(); i != elems.end(); ++i )
00933     {
00934         Range adj;
00935         rval = mb_impl.get_adjacencies( &*i, 1, dim + 1, false, adj, Interface::UNION );CHECK_ERR( rval );
00936         CHECK_EQUAL( 1, (int)adj.size() );
00937     }
00938 }
00939 
00940 void test_read_ids_common( const char* file_name, const char* tag_name, const int* expected_vals, int num_expected )
00941 {
00942     Core mb;
00943     ReadNCDF reader( &mb );
00944 
00945     FileOptions opts( "" );
00946     std::vector< int > values;
00947     ErrorCode rval = reader.read_tag_values( file_name, tag_name, opts, values );CHECK_ERR( rval );
00948 
00949     std::vector< int > expected( expected_vals, expected_vals + num_expected );
00950     std::sort( values.begin(), values.end() );
00951     std::sort( expected.begin(), expected.end() );
00952     CHECK_EQUAL( expected, values );
00953 }
00954 
00955 void test_read_block_ids()
00956 {
00957     const int expected[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 };
00958     test_read_ids_common( ho_file, MATERIAL_SET_TAG_NAME, expected, sizeof( expected ) / sizeof( expected[0] ) );
00959 }
00960 
00961 void test_read_sideset_ids()
00962 {
00963     const int expected[] = { 1, 2, 3, 4 };
00964     test_read_ids_common( ho_file, NEUMANN_SET_TAG_NAME, expected, sizeof( expected ) / sizeof( expected[0] ) );
00965 }
00966 
00967 void test_read_nodeset_ids()
00968 {
00969     test_read_ids_common( ho_file, DIRICHLET_SET_TAG_NAME, 0, 0 );
00970 }
00971 
00972 void test_read_alternate_coord_format()
00973 {
00974     Core moab;
00975     Interface& mb  = moab;
00976     ErrorCode rval = mb.load_file( alt_file );CHECK_ERR( rval );
00977 
00978     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 };
00979 
00980     Range hexes;
00981     rval = mb.get_entities_by_type( 0, MBHEX, hexes );CHECK_ERR( rval );
00982     CHECK_EQUAL( (size_t)1, hexes.size() );
00983     EntityHandle hex = hexes.front();
00984     const EntityHandle* conn;
00985     int len;
00986     rval = mb.get_connectivity( hex, conn, len );CHECK_ERR( rval );
00987     CHECK_EQUAL( 8, len );
00988     double act[3 * 8];
00989     rval = mb.get_coords( conn, len, act );CHECK_ERR( rval );
00990     CHECK_ARRAYS_EQUAL( exp, 3 * 8, act, 3 * 8 );
00991 }
00992 void test_write_polygons()
00993 {
00994     Core moab;
00995     Interface& mb  = moab;
00996     ErrorCode rval = mb.load_file( polyg );CHECK_ERR( rval );
00997     rval = mb.write_file( "polyg.exo" );CHECK_ERR( rval );
00998 }
00999 
01000 void test_write_polyhedra()
01001 {
01002     Core moab;
01003     Interface& mb  = moab;
01004     ErrorCode rval = mb.load_file( polyh );CHECK_ERR( rval );
01005     rval = mb.write_file( "polyh.exo" );CHECK_ERR( rval );
01006 }
01007 void test_read_polygons()
01008 {
01009     Core moab;
01010     Interface& mb  = moab;
01011     ErrorCode rval = mb.load_file( rpolyg );CHECK_ERR( rval );
01012 }
01013 void test_read_polyhedra()
01014 {
01015     Core moab;
01016     Interface& mb  = moab;
01017     ErrorCode rval = mb.load_file( rpolyh );CHECK_ERR( rval );
01018 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines