MOAB: Mesh Oriented datABase  (version 5.4.0)
cub_file_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 "moab/GeomTopoTool.hpp"
00007 #include <cmath>
00008 #include <algorithm>
00009 
00010 using namespace moab;
00011 
00012 /**\brief Input test file: test.cub
00013  * Cubit 10.2 file.
00014  * File contains:
00015  * Two merged 10x10x10 bricks sharing a single surface (surface 6).
00016  * - Each brick is meshed with 8 5x5x5 hexes.
00017  * - Each surface is meshed with 4 5x5 quads.
00018  * - Each curve is meshed with 2 5-unit edges.
00019  * A single block containing both bricks.
00020  * Two side sets
00021  * - sideset 1: surfaces 1 and 7
00022  * - sideset 2: surfaces 5 and 11
00023  * Two node sets:
00024  * - nodeset 1: surfaces 2 and 8
00025  * - nodeset 2: surfaces 3 and 9
00026  *
00027  * Surfaces:
00028  *                    2          8
00029  *                   /          /
00030  *          o----------o----------o
00031  *         /.      /  /.      /  /|
00032  *        / .  (5)/  / . (11)/  / |
00033  *       /  .    L  /  .    L  /  |
00034  *      o----------o----------o   |
00035  *      |   .      |   .      |(12)
00036  *   4--|-> o . . .|. .o. . . | . o
00037  *      |  . (1)   |  . (9)   |  /
00038  *      | .   ^    | .   ^    | /
00039  *      |.    |    |.    |    |/
00040  *      o----------o----------o
00041  *            |          |
00042  *            3          9
00043  *
00044  * Curves:
00045  *
00046  *     o----8-----o----20----o
00047  *    /.         /.         /|
00048  *  12 .       11 .       24 |
00049  *  /  7       /  5       /  17
00050  * o----2-----o----14----o   |
00051  * |   .      |   .      |   |
00052  * |   o . .6.|. .o. . 18| . o
00053  * 3  .       1  .      13  /
00054  * | 9        | 10       | 22
00055  * |.         |.         |/
00056  * o----4-----o----16----o
00057  *
00058  * Vertices:
00059  *
00060  *     8----------5----------13
00061  *    /.         /.         /|
00062  *   / .        / .        / |
00063  *  /  .       /  .       /  |
00064  * 3----------2----------10  |
00065  * |   .      |   .      |   |
00066  * |   7 . . .|. .6. . . | . 14
00067  * |  .       |  .       |  /
00068  * | .        | .        | /
00069  * |.         |.         |/
00070  * 4----------1----------9
00071  */
00072 
00073 /* Input test file: ho_test.cub
00074  *
00075  * File is expected to contain at least one block for every
00076  * supported higher-order element type.  The coordinates of
00077  * every higher-order node are expected to be the mean of the
00078  * adjacent corner vertices of the element.
00079  */
00080 static const std::string input_file_1 = std::string( TestDir + "unittest/io/test.cub" );
00081 static const std::string ho_file      = std::string( TestDir + "unittest/io/ho_test.cub" );
00082 static const std::string cubit12_file = std::string( TestDir + "unittest/io/cubtest12.cub" );
00083 static const std::string cubit14_file = std::string( TestDir + "unittest/io/cubtest14.cub" );
00084 
00085 void read_file( Interface& moab, const std::string& input_file );
00086 
00087 // Check that adjacent lower-order entities have
00088 // higher-order nodes consitent with input entity.
00089 void check_adj_ho_nodes( Interface& moab, EntityHandle entity );
00090 
00091 // Check that element has expected higher-order nodes
00092 // and that each higher-order node is at the center
00093 // of the sub-entity it is on.
00094 void check_ho_element( Interface& moab, EntityHandle entity, int mid_nodes[4] );
00095 
00096 // Validate elements of specified type.
00097 // Looks for a block containing the specified entity type
00098 // and with the specified mid-node flags set in its
00099 // HAS_MID_NODES_TAG.
00100 void test_ho_elements( EntityType type, int num_nodes );
00101 
00102 void test_vertices();
00103 
00104 void test_edges();
00105 
00106 void test_quads();
00107 
00108 void test_hexes();
00109 
00110 void test_geometric_topology();
00111 
00112 void test_geometric_sets();
00113 
00114 void test_blocks();
00115 
00116 void test_side_sets();
00117 
00118 void test_node_sets();
00119 
00120 void test_tri6()
00121 {
00122     test_ho_elements( MBTRI, 6 );
00123 }
00124 void test_tri7()
00125 {
00126     test_ho_elements( MBTRI, 7 );
00127 }
00128 
00129 void test_quad5()
00130 {
00131     test_ho_elements( MBQUAD, 5 );
00132 }
00133 void test_quad8()
00134 {
00135     test_ho_elements( MBQUAD, 8 );
00136 }
00137 void test_quad9()
00138 {
00139     test_ho_elements( MBQUAD, 9 );
00140 }
00141 
00142 void test_tet8()
00143 {
00144     test_ho_elements( MBTET, 8 );
00145 }
00146 void test_tet10()
00147 {
00148     test_ho_elements( MBTET, 10 );
00149 }
00150 void test_tet14()
00151 {
00152     test_ho_elements( MBTET, 14 );
00153 }
00154 
00155 void test_hex9()
00156 {
00157     test_ho_elements( MBHEX, 9 );
00158 }
00159 void test_hex20()
00160 {
00161     test_ho_elements( MBHEX, 20 );
00162 }
00163 void test_hex27()
00164 {
00165     test_ho_elements( MBHEX, 27 );
00166 }
00167 
00168 void test_multiple_files();
00169 
00170 void test_cubit12();
00171 void test_cubit14();
00172 
00173 int main()
00174 {
00175     int result = 0;
00176 
00177     result += RUN_TEST( test_vertices );
00178     result += RUN_TEST( test_edges );
00179     result += RUN_TEST( test_quads );
00180     result += RUN_TEST( test_hexes );
00181     result += RUN_TEST( test_geometric_topology );
00182     result += RUN_TEST( test_geometric_sets );
00183     result += RUN_TEST( test_blocks );
00184     result += RUN_TEST( test_side_sets );
00185     result += RUN_TEST( test_node_sets );
00186     result += RUN_TEST( test_tri6 );
00187     result += RUN_TEST( test_tri7 );
00188     result += RUN_TEST( test_quad5 );
00189     result += RUN_TEST( test_quad8 );
00190     result += RUN_TEST( test_quad9 );
00191     result += RUN_TEST( test_tet8 );
00192     result += RUN_TEST( test_tet10 );
00193     result += RUN_TEST( test_tet14 );
00194     result += RUN_TEST( test_hex9 );
00195     result += RUN_TEST( test_hex20 );
00196     result += RUN_TEST( test_hex27 );
00197     result += RUN_TEST( test_multiple_files );
00198     result += RUN_TEST( test_cubit12 );
00199     result += RUN_TEST( test_cubit14 );
00200     return result;
00201 }
00202 
00203 void read_file( Interface& moab, const std::string& input_file )
00204 {
00205     ErrorCode rval = moab.load_file( input_file.c_str() );CHECK_ERR( rval );
00206 }
00207 
00208 void test_vertices()
00209 {
00210     /* Node coordinates, in order of node ID, beginning with 1. */
00211     const double node_coords[] = {
00212         5,  -5, 5,                                                   // 1
00213         5,  5,  5,  5,  0,  5,  -5, 5,  5,  0,  5,  5,               // 5
00214         -5, -5, 5,  -5, 0,  5,  0,  -5, 5,  0,  0,  5,  5,  5,  -5,  // 10
00215         5,  -5, -5, 5,  0,  -5, -5, -5, -5, 0,  -5, -5, -5, 5,  -5,  // 15
00216         -5, 0,  -5, 0,  5,  -5, 0,  0,  -5, -5, -5, 0,  5,  -5, 0,   // 20
00217         0,  -5, 0,  -5, 5,  0,  -5, 0,  0,  5,  5,  0,  0,  5,  0,   // 25
00218         5,  0,  0,  0,  0,  0,  15, -5, 5,  15, 5,  5,  15, 0,  5,   // 30
00219         10, 5,  5,  10, -5, 5,  10, 0,  5,  15, 5,  -5, 15, -5, -5,  // 35
00220         15, 0,  -5, 10, -5, -5, 10, 5,  -5, 10, 0,  -5, 15, -5, 0,   // 40
00221         10, -5, 0,  15, 5,  0,  10, 5,  0,  15, 0,  0,  10, 0,  0    // 45
00222     };
00223 
00224     ErrorCode rval;
00225     Core mb_impl;
00226     Interface& mb = mb_impl;
00227     read_file( mb, input_file_1 );
00228 
00229     // get vertex handles and check correct number of vertices
00230     const size_t num_nodes = sizeof( node_coords ) / ( 3 * sizeof( double ) );
00231     Range verts;
00232     rval = mb.get_entities_by_type( 0, MBVERTEX, verts );CHECK_ERR( rval );
00233     CHECK_EQUAL( num_nodes, (size_t)verts.size() );
00234 
00235     // check global ids (should be 1 to 45 for vertices.)
00236     Tag gid_tag;
00237     rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
00238     std::vector< int > ids( num_nodes );
00239     rval = mb.tag_get_data( gid_tag, verts, &ids[0] );CHECK_ERR( rval );
00240     std::vector< int > sorted( ids );
00241     std::sort( sorted.begin(), sorted.end() );
00242     for( size_t i = 0; i < num_nodes; ++i )
00243         CHECK_EQUAL( (int)( i + 1 ), sorted[i] );
00244 
00245     // check coordinates of each vertex
00246     std::vector< double > coords( 3 * num_nodes );
00247     rval = mb.get_coords( verts, &coords[0] );CHECK_ERR( rval );
00248     for( size_t i = 0; i < num_nodes; ++i )
00249     {
00250         const double* exp = node_coords + 3 * ( ids[i] - 1 );
00251         const double* act = &coords[3 * i];
00252         CHECK_REAL_EQUAL( exp[0], act[0], 1e-8 );
00253         CHECK_REAL_EQUAL( exp[1], act[1], 1e-8 );
00254         CHECK_REAL_EQUAL( exp[2], act[2], 1e-8 );
00255     }
00256 }
00257 
00258 void test_element( const std::string& filename, EntityType type, int num_elem, int node_per_elem, const int* conn_list )
00259 {
00260     ErrorCode rval;
00261     Core mb_impl;
00262     Interface& mb = mb_impl;
00263     read_file( mb, filename );
00264 
00265     Range elems;
00266     rval = mb.get_entities_by_type( 0, type, elems );CHECK_ERR( rval );
00267     CHECK_EQUAL( num_elem, (int)elems.size() );
00268 
00269     // get global ids
00270     Tag gid_tag;
00271     rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
00272     std::vector< int > ids( num_elem );
00273     rval = mb.tag_get_data( gid_tag, elems, &ids[0] );CHECK_ERR( rval );
00274 
00275     // check that global ids are consecutive, beginning with 1
00276     std::vector< int > sorted( ids );
00277     std::sort( sorted.begin(), sorted.end() );
00278     for( int i = 0; i < num_elem; ++i )
00279         CHECK_EQUAL( i + 1, sorted[i] );
00280 
00281     // check connectivity of each element
00282     std::vector< int > conn_ids( node_per_elem );
00283     std::vector< EntityHandle > conn_h;
00284     Range::iterator j = elems.begin();
00285     for( int i = 0; i < num_elem; ++i, ++j )
00286     {
00287         conn_h.clear();
00288         rval = mb.get_connectivity( &*j, 1, conn_h );CHECK_ERR( rval );
00289         CHECK_EQUAL( node_per_elem, (int)conn_h.size() );
00290         rval = mb.tag_get_data( gid_tag, &conn_h[0], node_per_elem, &conn_ids[0] );CHECK_ERR( rval );
00291         const int* exp = conn_list + node_per_elem * ( ids[i] - 1 );
00292         for( int k = 0; k < node_per_elem; ++k )
00293             CHECK_EQUAL( exp[k], conn_ids[k] );
00294     }
00295 }
00296 
00297 void test_edges()
00298 {
00299     const int edge_conn[] = { 1,  3,                                                                           // 1
00300                               3,  2,  2,  5,  5,  4,  4,  7,  7,  6,  6,  8,  8,  1,  3,  9,  9,  8,           // 10
00301                               5,  9,  9,  7,  10, 12, 12, 11, 11, 14, 14, 13, 13, 16, 16, 15, 15, 17, 17, 10,  // 20
00302                               12, 18, 18, 17, 14, 18, 18, 16, 6,  19, 19, 13, 1,  20, 20, 11, 19, 21, 21, 8,   // 30
00303                               14, 21, 21, 20, 4,  22, 22, 15, 22, 23, 23, 7,  16, 23, 23, 19, 2,  24, 24, 10,  // 40
00304                               24, 25, 25, 5,  17, 25, 25, 22, 20, 26, 26, 3,  12, 26, 26, 24, 28, 30, 30, 29,  // 50
00305                               29, 31, 31, 2,  1,  32, 32, 28, 30, 33, 33, 32, 31, 33, 33, 3,  34, 36, 36, 35,  // 60
00306                               35, 37, 37, 11, 10, 38, 38, 34, 36, 39, 39, 38, 37, 39, 39, 12, 28, 40, 40, 35,  // 70
00307                               20, 41, 41, 32, 37, 41, 41, 40, 29, 42, 42, 34, 42, 43, 43, 31, 38, 43, 43, 24,  // 80
00308                               40, 44, 44, 30, 36, 44, 44, 42 };
00309     test_element( input_file_1, MBEDGE, 84, 2, edge_conn );
00310 }
00311 
00312 void test_quads()
00313 {
00314     const int quad_conn[] = {
00315         1,  3,  9,  8,  // 1
00316         3,  2,  5,  9,  8,  9,  7,  6,  9,  5,  4,  7,  10, 12, 18, 17, 12, 11, 14, 18,
00317         17, 18, 16, 15, 18, 14, 13, 16, 6,  19, 21, 8,  19, 13, 14, 21,  // 10
00318         8,  21, 20, 1,  21, 14, 11, 20, 4,  22, 23, 7,  22, 15, 16, 23, 7,  23, 19, 6,
00319         23, 16, 13, 19, 2,  24, 25, 5,  24, 10, 17, 25, 5,  25, 22, 4,  25, 17, 15, 22,  // 20
00320         1,  20, 26, 3,  20, 11, 12, 26, 3,  26, 24, 2,  26, 12, 10, 24, 28, 30, 33, 32,
00321         30, 29, 31, 33, 32, 33, 3,  1,  33, 31, 2,  3,  34, 36, 39, 38, 36, 35, 37, 39,  // 30
00322         38, 39, 12, 10, 39, 37, 11, 12, 1,  20, 41, 32, 20, 11, 37, 41, 32, 41, 40, 28,
00323         41, 37, 35, 40, 29, 42, 43, 31, 42, 34, 38, 43, 31, 43, 24, 2,  43, 38, 10, 24,  // 40
00324         28, 40, 44, 30, 40, 35, 36, 44, 30, 44, 42, 29, 44, 36, 34, 42,
00325     };
00326     test_element( input_file_1, MBQUAD, 44, 4, quad_conn );
00327 }
00328 
00329 void test_hexes()
00330 {
00331     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,
00332                              25, 5,  23, 16, 15, 22, 27, 18, 17, 25, 8,  21, 27, 9,  1,  20, 26, 3,  21, 14, 18, 27,
00333                              20, 11, 12, 26, 9,  27, 25, 5,  3,  26, 24, 2,  27, 18, 17, 25, 26, 12, 10, 24, 1,  20,
00334                              26, 3,  32, 41, 45, 33, 20, 11, 12, 26, 41, 37, 39, 45, 3,  26, 24, 2,  33, 45, 43, 31,
00335                              26, 12, 10, 24, 45, 39, 38, 43, 32, 41, 45, 33, 28, 40, 44, 30, 41, 37, 39, 45, 40, 35,
00336                              36, 44, 33, 45, 43, 31, 30, 44, 42, 29, 45, 39, 38, 43, 44, 36, 34, 42 };
00337     test_element( input_file_1, MBHEX, 16, 8, hex_conn );
00338 }
00339 
00340 template < int L >
00341 std::vector< int > find_parents( const int parent_conn[][L], int num_parent, int id )
00342 {
00343     std::vector< int > results;
00344     for( int i = 0; i < num_parent; ++i )
00345     {
00346         for( int j = 0; j < L; ++j )
00347         {
00348             if( parent_conn[i][j] == id ) results.push_back( i + 1 );
00349         }
00350     }
00351     return results;
00352 }
00353 
00354 int check_geometric_set( Interface& moab,
00355                          int dim,
00356                          int id,
00357                          const int* children,
00358                          int num_children,
00359                          std::vector< int > parents )
00360 {
00361     ErrorCode rval;
00362     Tag gid_tag, dim_tag;
00363 
00364     rval = moab.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
00365     rval = moab.tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, dim_tag );CHECK_ERR( rval );
00366     void* tag_vals[] = { &dim, &id };
00367     Tag tags[]       = { dim_tag, gid_tag };
00368     Range ents;
00369     rval = moab.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, tag_vals, 2, ents );CHECK_ERR( rval );
00370     CHECK_EQUAL( 1u, (unsigned)ents.size() );
00371 
00372     const EntityHandle geom = ents.front();
00373     std::vector< int > exp_rel, act_rel;
00374     std::vector< EntityHandle > rel;
00375 
00376     if( num_children )
00377     {
00378         exp_rel.resize( num_children );
00379         std::copy( children, children + num_children, exp_rel.begin() );
00380         std::sort( exp_rel.begin(), exp_rel.end() );
00381         rel.clear();
00382         rval = moab.get_child_meshsets( geom, rel );CHECK_ERR( rval );
00383         CHECK_EQUAL( num_children, (int)rel.size() );
00384         act_rel.resize( rel.size() );
00385         rval = moab.tag_get_data( gid_tag, &rel[0], rel.size(), &act_rel[0] );CHECK_ERR( rval );
00386         std::sort( act_rel.begin(), act_rel.end() );
00387         CHECK( exp_rel == act_rel );
00388     }
00389 
00390     if( !parents.empty() )
00391     {
00392         exp_rel = parents;
00393         std::sort( exp_rel.begin(), exp_rel.end() );
00394         rel.clear();
00395         rval = moab.get_parent_meshsets( geom, rel );CHECK_ERR( rval );
00396         CHECK_EQUAL( parents.size(), rel.size() );
00397         act_rel.resize( rel.size() );
00398         rval = moab.tag_get_data( gid_tag, &rel[0], rel.size(), &act_rel[0] );CHECK_ERR( rval );
00399         std::sort( act_rel.begin(), act_rel.end() );
00400         CHECK( exp_rel == act_rel );
00401     }
00402 
00403     return 0;
00404 }
00405 
00406 void test_geometric_topology()
00407 {
00408     Core mb_impl;
00409     Interface& mb = mb_impl;
00410     read_file( mb, input_file_1 );
00411     // expected geometric vertices, specified by global ID
00412     const int vertex_ids[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 13, 14 };
00413     // List of global IDs of surfacs in geometric volumes, indexed by ID-1
00414     const int volume_surfs[2][6] = { { 1, 2, 3, 4, 5, 6 }, { 7, 8, 9, 6, 11, 12 } };
00415     // List of global IDs of curves in geometric surfaces, indexed by ID-1
00416     // Curve IDs of zero indicates that corresponding surface doesn't exist.
00417     const int surf_curves[12][4] = { { 1, 2, 3, 4 },     { 5, 6, 7, 8 },    { 9, 6, 10, 4 },   { 11, 7, 9, 3 },
00418                                      { 12, 8, 11, 2 },   { 10, 5, 12, 1 },  { 13, 14, 1, 16 }, { 17, 18, 5, 20 },
00419                                      { 10, 18, 22, 16 }, { 0, 0, 0, 0 },  // no surf 10
00420                                      { 24, 20, 12, 14 }, { 22, 17, 24, 13 } };
00421     // List of global IDs of vertices in geometric curves, indexed by ID-1
00422     // Vertex IDs of zero indicates that corresponding curve doesn't exist.
00423     const int curve_verts[24][2] = { { 1, 2 },  { 2, 3 },   { 3, 4 },  { 4, 1 }, { 5, 6 },  // 5
00424                                      { 6, 7 },  { 7, 8 },   { 8, 5 },  { 4, 7 }, { 1, 6 },  // 10
00425                                      { 3, 8 },  { 2, 5 },   { 9, 10 },                      // 13
00426                                      { 10, 2 }, { 0, 0 },                                   // no curve 15
00427                                      { 1, 9 },  { 13, 14 }, { 14, 6 }, { 0, 0 },            // no curve 19
00428                                      { 5, 13 }, { 0, 0 },                                   // no curve 21
00429                                      { 9, 14 }, { 0, 0 },                                   // no curve 23
00430                                      { 10, 13 } };
00431 
00432     // check all vertices
00433     for( unsigned i = 0; i < ( sizeof( vertex_ids ) / sizeof( vertex_ids[0] ) ); ++i )
00434         check_geometric_set( mb, 0, vertex_ids[i], 0, 0, find_parents< 2 >( curve_verts, 24, vertex_ids[i] ) );
00435 
00436     // check all curves
00437     for( int i = 1; i <= 24; ++i )
00438         if( curve_verts[i - 1][0] )
00439             check_geometric_set( mb, 1, i, curve_verts[i - 1], 2, find_parents< 4 >( surf_curves, 12, i ) );
00440 
00441     // check all surfs
00442     for( int i = 1; i <= 12; ++i )
00443         if( surf_curves[i - 1][0] )
00444             check_geometric_set( mb, 2, i, surf_curves[i - 1], 4, find_parents< 6 >( volume_surfs, 2, i ) );
00445 
00446     // check all volumes
00447     std::vector< int > empty;
00448     for( int i = 1; i <= 2; ++i )
00449         check_geometric_set( mb, 3, i, volume_surfs[i - 1], 6, empty );
00450 }
00451 
00452 void test_geometric_sets()
00453 {
00454     ErrorCode rval;
00455     Core mb_impl;
00456     Interface& mb = mb_impl;
00457     read_file( mb, input_file_1 );
00458     Tag gid_tag, dim_tag;
00459     rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
00460     rval = mb.tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, dim_tag );CHECK_ERR( rval );
00461 
00462     // verify mesh entity counts
00463     Range verts, curves, surfs, vols;
00464     int dim = 0;
00465     // Cppcheck warning (false positive): variable dim is assigned a value that is never used
00466     // Cppcheck warning (false positive): variable dim is reassigned a value before the old one has
00467     // been used
00468     const void* vals[] = { &dim };
00469     rval               = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, verts );CHECK_ERR( rval );
00470     dim  = 1;
00471     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, curves );CHECK_ERR( rval );
00472     dim  = 2;
00473     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, surfs );CHECK_ERR( rval );
00474     dim  = 3;
00475     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals, 1, vols );CHECK_ERR( rval );
00476 
00477     CHECK_EQUAL( 12u, (unsigned)verts.size() );
00478     CHECK_EQUAL( 20u, (unsigned)curves.size() );
00479     CHECK_EQUAL( 11u, (unsigned)surfs.size() );
00480     CHECK_EQUAL( 2u, (unsigned)vols.size() );
00481 
00482     // check that each vertex has a single node, and that the
00483     // node is also contained in any parent curve
00484     Range ents;
00485     Range::iterator i;
00486     for( i = verts.begin(); i != verts.end(); ++i )
00487     {
00488         ents.clear();
00489         rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
00490         CHECK_EQUAL( 1u, (unsigned)ents.size() );
00491         CHECK( ents.all_of_type( MBVERTEX ) );
00492     }
00493 
00494     // check that each curve has one node and two edges
00495     for( i = curves.begin(); i != curves.end(); ++i )
00496     {
00497         ents.clear();
00498         rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
00499         CHECK_EQUAL( 1u, (unsigned)ents.num_of_type( MBVERTEX ) );
00500         CHECK_EQUAL( 2u, (unsigned)ents.num_of_type( MBEDGE ) );
00501         CHECK_EQUAL( 3u, (unsigned)ents.size() );
00502     }
00503 
00504     // check that each surface has 1 node, 4 edges, 4 quads
00505     for( i = surfs.begin(); i != surfs.end(); ++i )
00506     {
00507         ents.clear();
00508         rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
00509         CHECK_EQUAL( 1u, (unsigned)ents.num_of_type( MBVERTEX ) );
00510         CHECK_EQUAL( 4u, (unsigned)ents.num_of_type( MBEDGE ) );
00511         CHECK_EQUAL( 4u, (unsigned)ents.num_of_type( MBQUAD ) );
00512         CHECK_EQUAL( 9u, (unsigned)ents.size() );
00513     }
00514 
00515     // check that each volume has 1 node and 8 hexes.
00516     for( i = vols.begin(); i != vols.end(); ++i )
00517     {
00518         ents.clear();
00519         rval = mb.get_entities_by_handle( *i, ents );CHECK_ERR( rval );
00520         CHECK_EQUAL( 1u, (unsigned)ents.num_of_type( MBVERTEX ) );
00521         CHECK_EQUAL( 8u, (unsigned)ents.num_of_type( MBHEX ) );
00522         CHECK_EQUAL( 9u, (unsigned)ents.size() );
00523     }
00524 
00525     // Check that for each geometric entity, any contained vertices
00526     // are adjacent to some entity in one of its parents.
00527     Range parents, geom, nodes, tmp;
00528     for( int d = 0; d < 3; ++d )
00529     {
00530         const void* vals1[] = { &d };
00531         rval                = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, vals1, 1, geom );CHECK_ERR( rval );
00532 
00533         for( i = geom.begin(); i != geom.end(); ++i )
00534         {
00535             nodes.clear();
00536             ents.clear();
00537             parents.clear();
00538             rval = mb.get_entities_by_type( *i, MBVERTEX, nodes );CHECK_ERR( rval );
00539             rval = mb.get_parent_meshsets( *i, parents );CHECK_ERR( rval );
00540             for( Range::iterator j = parents.begin(); j != parents.end(); ++j )
00541             {
00542                 tmp.clear();
00543                 rval = mb.get_entities_by_dimension( *j, d + 1, tmp );CHECK_ERR( rval );
00544                 ents.merge( tmp );
00545             }
00546             tmp.clear();
00547             rval = mb.get_adjacencies( ents, 0, false, tmp, Interface::UNION );CHECK_ERR( rval );
00548             nodes = subtract( nodes, tmp );
00549             CHECK( nodes.empty() );
00550         }
00551     }
00552 }
00553 
00554 // expect one block containing entire mesh, with id == 1
00555 void test_blocks()
00556 {
00557     ErrorCode rval;
00558     Core mb_impl;
00559     Interface& mb = mb_impl;
00560     read_file( mb, input_file_1 );
00561     Tag mat_tag;
00562     rval = mb.tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mat_tag );CHECK_ERR( rval );
00563 
00564     Range blocks;
00565     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &mat_tag, 0, 1, blocks );CHECK_ERR( rval );
00566     CHECK_EQUAL( 1u, (unsigned)blocks.size() );
00567     EntityHandle block = blocks.front();
00568     int id;
00569     rval = mb.tag_get_data( mat_tag, &block, 1, &id );CHECK_ERR( rval );
00570     CHECK_EQUAL( 1, id );
00571 
00572     Range block_hexes, mesh_hexes;
00573     rval = mb.get_entities_by_dimension( 0, 3, mesh_hexes );CHECK_ERR( rval );
00574     rval = mb.get_entities_by_dimension( block, 3, block_hexes, true );CHECK_ERR( rval );
00575     CHECK( mesh_hexes == block_hexes );
00576 }
00577 
00578 // Common code for test_side_sets and test_node sets
00579 //\param tag_name  NEUMANN_SET_TAG_NAME or DIRICHLET_SET_TAG_NAME
00580 //\param count     Number of expected sets
00581 //\param ids       Expected IDs of sets
00582 //\param set_surfs One list for each id in "ids" containing the
00583 //                 ids of the geometric surfaces expected to be
00584 //                 contained in the boundary condition set.
00585 void test_bc_sets( const char* tag_name, unsigned count, const int* ids, const std::vector< int > set_surfs[] )
00586 {
00587     ErrorCode rval;
00588     Core mb_impl;
00589     Interface& mb = mb_impl;
00590     read_file( mb, input_file_1 );
00591     Tag ss_tag, gid_tag, dim_tag;
00592     rval = mb.tag_get_handle( tag_name, 1, MB_TYPE_INTEGER, ss_tag );CHECK_ERR( rval );
00593     rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
00594     rval = mb.tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, dim_tag );CHECK_ERR( rval );
00595 
00596     // check number of sidesets and IDs
00597     Range sidesets;
00598     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &ss_tag, 0, 1, sidesets );CHECK_ERR( rval );
00599     CHECK_EQUAL( count, (unsigned)sidesets.size() );
00600     std::vector< EntityHandle > handles( count, 0 );
00601     for( Range::iterator i = sidesets.begin(); i != sidesets.end(); ++i )
00602     {
00603         int id;
00604         rval = mb.tag_get_data( ss_tag, &*i, 1, &id );CHECK_ERR( rval );
00605         unsigned idx;
00606         for( idx = 0; idx < count; ++idx )
00607             if( ids[idx] == id ) break;
00608         CHECK( idx != count );
00609         CHECK( handles[idx] == 0 );
00610         handles[idx] = *i;
00611     }
00612 
00613     // get surface faces
00614     std::vector< Range > exp( count );
00615     Range surfs, tmp;
00616     Tag tags[] = { dim_tag, gid_tag };
00617     for( unsigned i = 0; i < count; ++i )
00618     {
00619         exp[i].clear();
00620         surfs.clear();
00621         const int two = 2;
00622         for( unsigned j = 0; j < set_surfs[i].size(); ++j )
00623         {
00624             const void* vals[] = { &two, &set_surfs[i][j] };
00625             surfs.clear();
00626             rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, surfs );CHECK_ERR( rval );
00627             CHECK_EQUAL( 1u, (unsigned)surfs.size() );
00628             tmp.clear();
00629             rval = mb.get_entities_by_dimension( surfs.front(), 2, tmp, true );CHECK_ERR( rval );
00630             exp[i].merge( tmp );
00631         }
00632     }
00633 
00634     // check each bc set
00635     Range act;
00636     for( unsigned i = 0; i < count; ++i )
00637     {
00638         act.clear();
00639         rval = mb.get_entities_by_dimension( handles[i], 2, act, true );CHECK_ERR( rval );
00640         CHECK( exp[i] == act );
00641     }
00642 }
00643 
00644 // expect two sidesets containing two geometric surfaces each:
00645 // sideset 1 : surfaces 1 and 7
00646 // sideset 2 : surfaces 5 and 11
00647 void test_side_sets()
00648 {
00649     int ids[] = { 1, 2 };
00650     std::vector< int > surfs[2];
00651     surfs[0].push_back( 1 );
00652     surfs[0].push_back( 7 );
00653     surfs[1].push_back( 5 );
00654     surfs[1].push_back( 11 );
00655     test_bc_sets( NEUMANN_SET_TAG_NAME, 2, ids, surfs );
00656 }
00657 
00658 void test_node_sets()
00659 {
00660     int ids[] = { 1, 2 };
00661     std::vector< int > surfs[2];
00662     surfs[0].push_back( 2 );
00663     surfs[0].push_back( 8 );
00664     surfs[1].push_back( 3 );
00665     surfs[1].push_back( 9 );
00666     test_bc_sets( DIRICHLET_SET_TAG_NAME, 2, ids, surfs );
00667 }
00668 
00669 static EntityHandle find_side( Interface& moab, EntityHandle entity, int side_dim, int side_num )
00670 {
00671     ErrorCode rval;
00672 
00673     std::vector< EntityHandle > adj;
00674     rval = moab.get_adjacencies( &entity, 1, side_dim, false, adj );CHECK_ERR( rval );
00675 
00676     int sub_ent_indices[4];
00677     CN::SubEntityVertexIndices( TYPE_FROM_HANDLE( entity ), side_dim, side_num, sub_ent_indices );
00678     EntityType subtype  = CN::SubEntityType( TYPE_FROM_HANDLE( entity ), side_dim, side_num );
00679     int sub_ent_corners = CN::VerticesPerEntity( subtype );
00680 
00681     const EntityHandle* conn;
00682     int conn_len;
00683     rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );
00684 
00685     for( size_t i = 0; i < adj.size(); ++i )
00686     {
00687         if( TYPE_FROM_HANDLE( adj[i] ) != subtype ) continue;
00688 
00689         const EntityHandle* sub_conn;
00690         int sub_len;
00691         rval = moab.get_connectivity( adj[i], sub_conn, sub_len );CHECK_ERR( rval );
00692 
00693         int n = std::find( sub_conn, sub_conn + sub_len, conn[sub_ent_indices[0]] ) - sub_conn;
00694         if( n == sub_len )  // no vertex in common
00695             continue;
00696 
00697         // check forward direction
00698         int j;
00699         for( j = 1; j < sub_ent_corners; ++j )
00700             if( conn[sub_ent_indices[j]] != sub_conn[( j + n ) % sub_ent_corners] ) break;
00701         if( j == sub_ent_corners ) return adj[i];
00702 
00703         // check reverse direction
00704         for( j = 1; j < sub_ent_corners; ++j )
00705             if( conn[sub_ent_indices[j]] != sub_conn[( n + sub_ent_corners - j ) % sub_ent_corners] ) break;
00706         if( j == sub_ent_corners ) return adj[i];
00707     }
00708 
00709     // no match
00710     return 0;
00711 }
00712 
00713 void check_adj_ho_nodes( Interface& moab, EntityHandle entity )
00714 {
00715     EntityType type = TYPE_FROM_HANDLE( entity );
00716     const EntityHandle* conn;
00717     int conn_len;
00718     ErrorCode rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );
00719 
00720     int ho[4];
00721     CN::HasMidNodes( type, conn_len, ho );
00722     for( int dim = CN::Dimension( type ) - 1; dim > 0; --dim )
00723     {
00724         if( !ho[dim] ) continue;
00725 
00726         for( int j = 0; j < CN::NumSubEntities( type, dim ); ++j )
00727         {
00728             EntityHandle side = find_side( moab, entity, dim, j );
00729             if( !side ) continue;
00730 
00731             const EntityHandle* side_conn;
00732             int side_len;
00733             rval = moab.get_connectivity( side, side_conn, side_len );CHECK_ERR( rval );
00734 
00735             int this_idx = CN::HONodeIndex( type, conn_len, dim, j );
00736             int side_idx = CN::HONodeIndex( TYPE_FROM_HANDLE( side ), side_len, dim, 0 );
00737             CHECK_EQUAL( side_conn[side_idx], conn[this_idx] );
00738         }
00739     }
00740 }
00741 
00742 // Check that element has expected higher-order nodes
00743 // and that each higher-order node is at the center
00744 // of the sub-entity it is on.
00745 void check_ho_element( Interface& moab, EntityHandle entity, int mid_nodes[4] )
00746 {
00747     // get element info
00748     const EntityType type = TYPE_FROM_HANDLE( entity );
00749     const EntityHandle* conn;
00750     int conn_len;
00751     ErrorCode rval = moab.get_connectivity( entity, conn, conn_len );CHECK_ERR( rval );
00752     std::vector< double > coords( 3 * conn_len );
00753     rval = moab.get_coords( conn, conn_len, &coords[0] );CHECK_ERR( rval );
00754 
00755     // calculate and verify expected number of mid nodes
00756     int num_nodes = CN::VerticesPerEntity( type );
00757     for( int d = 1; d <= CN::Dimension( type ); ++d )
00758         if( mid_nodes[d] ) num_nodes += CN::NumSubEntities( type, d );
00759     CHECK_EQUAL( num_nodes, conn_len );
00760 
00761     // verify that each higher-order node is at the center
00762     // of its respective sub-entity.
00763     for( int i = CN::VerticesPerEntity( type ); i < num_nodes; ++i )
00764     {
00765         // get sub-entity owning ho-node
00766         int sub_dim, sub_num;
00767         CN::HONodeParent( type, num_nodes, i, sub_dim, sub_num );
00768         // get corner vertex indices
00769         int sub_conn[8], num_sub;
00770         if( sub_dim < CN::Dimension( type ) )
00771         {
00772             CN::SubEntityVertexIndices( type, sub_dim, sub_num, sub_conn );
00773             EntityType sub_type = CN::SubEntityType( type, sub_dim, sub_num );
00774             num_sub             = CN::VerticesPerEntity( sub_type );
00775         }
00776         else
00777         {
00778             num_sub = CN::VerticesPerEntity( type );
00779             for( int j = 0; j < num_sub; ++j )
00780                 sub_conn[j] = j;
00781         }
00782         // calculate mean of corner vertices
00783         double mean[3] = { 0, 0, 0 };
00784         for( int j = 0; j < num_sub; ++j )
00785         {
00786             int co = 3 * sub_conn[j];
00787             mean[0] += coords[co];
00788             mean[1] += coords[co + 1];
00789             mean[2] += coords[co + 2];
00790         }
00791         mean[0] /= num_sub;
00792         mean[1] /= num_sub;
00793         mean[2] /= num_sub;
00794         // verify that higher-order node is at expected location
00795         CHECK_REAL_EQUAL( mean[0], coords[3 * i], 1e-6 );
00796         CHECK_REAL_EQUAL( mean[1], coords[3 * i + 1], 1e-6 );
00797         CHECK_REAL_EQUAL( mean[2], coords[3 * i + 2], 1e-6 );
00798     }
00799 }
00800 
00801 // Validate elements of specified type.
00802 // Looks for a block containing the specified entity type
00803 // and with the specified mid-node flags set in its
00804 // HAS_MID_NODES_TAG.
00805 void test_ho_elements( EntityType type, int num_nodes )
00806 {
00807     Core mb_impl;
00808     Interface& mb = mb_impl;
00809     read_file( mb, ho_file );
00810 
00811     ErrorCode rval;
00812     Tag ho_tag, block_tag;
00813     rval = mb.tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, block_tag );CHECK_ERR( rval );
00814     rval = mb.tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, ho_tag );CHECK_ERR( rval );
00815 
00816     // get material sets with expected higher-order nodes
00817     Range blocks;
00818     int ho_flags[4];
00819     CN::HasMidNodes( type, num_nodes, ho_flags );
00820     Tag tags[2]   = { ho_tag, block_tag };
00821     void* vals[2] = { ho_flags, NULL };
00822     rval          = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, blocks );CHECK_ERR( rval );
00823 
00824     Range::iterator i;
00825     Range entities;
00826     for( i = blocks.begin(); i != blocks.end(); ++i )
00827     {
00828         rval = mb.get_entities_by_type( *i, type, entities, true );CHECK_ERR( rval );
00829     }
00830     // test file should contain all types of HO entities
00831     CHECK( !entities.empty() );
00832     // test ho node positions -- expect to be at center of adj corners
00833     for( i = entities.begin(); i != entities.end(); ++i )
00834         check_ho_element( mb, *i, ho_flags );
00835     // test ho node handles consistent with adjacent entities
00836     for( i = entities.begin(); i != entities.end(); ++i )
00837         check_adj_ho_nodes( mb, *i );
00838 }
00839 
00840 void test_multiple_files()
00841 {
00842     // load two surface meshes, one at z=+5 at z=-5.
00843     ErrorCode rval;
00844     Core mb_impl;
00845     Interface& mb = mb_impl;
00846     Range file1_ents, file2_ents;
00847     read_file( mb, input_file_1 );
00848     mb.get_entities_by_handle( 0, file1_ents );
00849     read_file( mb, input_file_1 );
00850     mb.get_entities_by_handle( 0, file2_ents );
00851     file2_ents = subtract( file2_ents, file1_ents );
00852     EntityHandle file1, file2;
00853     mb.create_meshset( MESHSET_SET, file1 );
00854     mb.create_meshset( MESHSET_SET, file2 );
00855     mb.add_entities( file1, file1_ents );
00856     mb.add_entities( file2, file2_ents );
00857 
00858     // first check that we get the same number of verts from
00859     // each file and that they are distinct vertices
00860     Range file1_verts, file2_verts;
00861     rval = mb.get_entities_by_type( file1, MBVERTEX, file1_verts );CHECK_ERR( rval );
00862     CHECK( !file1_verts.empty() );
00863     rval = mb.get_entities_by_type( file2, MBVERTEX, file2_verts );CHECK_ERR( rval );
00864     CHECK( !file2_verts.empty() );
00865     CHECK_EQUAL( file1_verts.size(), file2_verts.size() );
00866     CHECK( intersect( file1_verts, file2_verts ).empty() );
00867 
00868     // now check that we get the same number of elements from
00869     // each file and that they are distinct
00870     Range file1_elems, file2_elems;
00871     rval = mb.get_entities_by_dimension( file1, 3, file1_elems );CHECK_ERR( rval );
00872     CHECK( !file1_elems.empty() );
00873     rval = mb.get_entities_by_dimension( file2, 3, file2_elems );CHECK_ERR( rval );
00874     CHECK( !file2_elems.empty() );
00875     CHECK_EQUAL( file1_elems.size(), file2_elems.size() );
00876     CHECK( intersect( file1_elems, file2_elems ).empty() );
00877 
00878     // now check that the connectivity for each element is
00879     // defined using the appropriate vertex instances
00880     Range file1_elem_verts, file2_elem_verts;
00881     rval = mb.get_adjacencies( file1_elems, 0, false, file1_elem_verts, Interface::UNION );CHECK_ERR( rval );
00882     rval = mb.get_adjacencies( file2_elems, 0, false, file2_elem_verts, Interface::UNION );CHECK_ERR( rval );
00883     CHECK_EQUAL( file1_elem_verts.size(), file2_elem_verts.size() );
00884     CHECK( intersect( file1_elem_verts, file1_verts ) == file1_elem_verts );
00885     CHECK( intersect( file2_elem_verts, file2_verts ) == file2_elem_verts );
00886 }
00887 
00888 void test_cubit12()
00889 {
00890     Core mb_impl;
00891     Interface& mb = mb_impl;
00892     read_file( mb, cubit12_file );
00893 }
00894 
00895 void test_cubit14()
00896 {
00897     Core mb_impl;
00898     Interface& mb = mb_impl;
00899     read_file( mb, cubit14_file );
00900     // check the global id for some geometry sets
00901     GeomTopoTool gtt( &mb_impl );
00902     Range ranges[5];
00903     ErrorCode rval = gtt.find_geomsets( ranges );CHECK_ERR( rval );
00904     EntityHandle set0 = ranges[0][0];  // does it have a global id > 0?
00905     Tag gid_tag;
00906     rval = mb.tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gid_tag );CHECK_ERR( rval );
00907 
00908     int val;
00909     rval = mb.tag_get_data( gid_tag, &set0, 1, &val );
00910     CHECK( val != 0 );
00911 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines