fathom/moab-docs/ideas__test_8cpp_source.html

00001 #include "TestUtil.hpp"
00002 #include "moab/Core.hpp"
00003 #include "MBTagConventions.hpp"
00004 #define IS_BUILDING_MB
00005 #include "ReadIDEAS.hpp"
00006 #include "moab/Range.hpp"
00007 #include <cmath>
00008 #include <algorithm>
00009
00010 using namespace moab;
00011
00012 #ifdef MESHDIR
00013 static const char example[] = STRINGIFY( MESHDIR ) "/io/test.unv";
00014 #else
00015 static const char example[] = "test.unv";
00016 #endif
00017
00018 void read_file( Interface& moab, const char* input_file );
00019 void test_read_nodes();
00020 void test_read_tets();
00021 void test_read_hexes();
00022 void test_read_material_set();
00023 void test_read_physical_set();
00024
00025 int main()
00026 {
00027     int result = 0;
00028
00029     result += RUN_TEST( test_read_nodes );
00030     result += RUN_TEST( test_read_tets );
00031     result += RUN_TEST( test_read_hexes );
00032     result += RUN_TEST( test_read_material_set );
00033     result += RUN_TEST( test_read_physical_set );
00034
00035     return result;
00036 }
00037
00038 void read_file( Interface& moab, const char* input_file )
00039 {
00040     ErrorCode rval = moab.load_file( input_file );CHECK_ERR( rval );
00041 }
00042
00043 void test_read_nodes()
00044 {
00045     const double eps = 1e-6;
00046     ErrorCode rval;
00047     Core moab;
00048     Interface& mb = moab;
00049     read_file( moab, example );
00050
00051     std::vector< EntityHandle > nodes;
00052     rval = mb.get_entities_by_type( 0, MBVERTEX, nodes );CHECK_ERR( rval );
00053     CHECK_EQUAL( (size_t)17, nodes.size() );
00054
00055     Tag id_tag = mb.globalId_tag();
00056
00057     std::vector< int > ids( nodes.size() );
00058     rval = mb.tag_get_data( id_tag, &nodes[0], nodes.size(), &ids[0] );CHECK_ERR( rval );
00059
00060     std::vector< int > sorted_ids( ids );
00061     std::sort( sorted_ids.begin(), sorted_ids.end() );
00062
00063     std::vector< double > coords( 3 * nodes.size() );
00064     rval = mb.get_coords( &nodes[0], nodes.size(), &coords[0] );CHECK_ERR( rval );
00065
00066     int idx, pos = 0;
00067     // shared between 2 tets and 2 prisms
00068     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00069     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00070     CHECK_REAL_EQUAL( coords[3 * idx + 0], 0.0, eps );
00071     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0, eps );
00072     CHECK_REAL_EQUAL( coords[3 * idx + 2], 0.0, eps );
00073
00074     ++pos;
00075     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00076     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00077     CHECK_REAL_EQUAL( coords[3 * idx + 0], 0.0, eps );
00078     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0, eps );
00079     CHECK_REAL_EQUAL( coords[3 * idx + 2], 1.0, eps );
00080
00081     ++pos;
00082     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00083     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00084     CHECK_REAL_EQUAL( coords[3 * idx + 0], 0.0, eps );
00085     CHECK_REAL_EQUAL( coords[3 * idx + 1], 0.0, eps );
00086     CHECK_REAL_EQUAL( coords[3 * idx + 2], 1.0, eps );
00087
00088     // id=4
00089     ++pos;
00090     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00091     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00092     CHECK_REAL_EQUAL( coords[3 * idx + 0], 0.0, eps );
00093     CHECK_REAL_EQUAL( coords[3 * idx + 1], 0.0, eps );
00094     CHECK_REAL_EQUAL( coords[3 * idx + 2], 0.0, eps );
00095
00096     ++pos;
00097     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00098     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00099     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0, eps );
00100     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0, eps );
00101     CHECK_REAL_EQUAL( coords[3 * idx + 2], 0.0, eps );
00102
00103     ++pos;
00104     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00105     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00106     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0, eps );
00107     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0, eps );
00108     CHECK_REAL_EQUAL( coords[3 * idx + 2], 1.0, eps );
00109
00110     ++pos;
00111     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00112     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00113     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0, eps );
00114     CHECK_REAL_EQUAL( coords[3 * idx + 1], 0.0, eps );
00115     CHECK_REAL_EQUAL( coords[3 * idx + 2], 1.0, eps );
00116
00117     // id=8
00118     ++pos;
00119     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00120     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00121     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0, eps );
00122     CHECK_REAL_EQUAL( coords[3 * idx + 1], 0.0, eps );
00123     CHECK_REAL_EQUAL( coords[3 * idx + 2], 0.0, eps );
00124
00125     ++pos;
00126     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00127     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00128     CHECK_REAL_EQUAL( coords[3 * idx + 0], 0.0, eps );
00129     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0, eps );
00130     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00131
00132     ++pos;
00133     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00134     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00135     CHECK_REAL_EQUAL( coords[3 * idx + 0], 0.0, eps );
00136     CHECK_REAL_EQUAL( coords[3 * idx + 1], 0.0, eps );
00137     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00138
00139     ++pos;
00140     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00141     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00142     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0, eps );
00143     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0, eps );
00144     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00145
00146     // id=12
00147     ++pos;
00148     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00149     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00150     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0, eps );
00151     CHECK_REAL_EQUAL( coords[3 * idx + 1], 0.0, eps );
00152     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00153
00154     ++pos;
00155     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00156     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00157     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0000022760448197, eps );
00158     CHECK_REAL_EQUAL( coords[3 * idx + 1], -2.2760448196157412e-06, eps );
00159     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00160
00161     ++pos;
00162     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00163     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00164     CHECK_REAL_EQUAL( coords[3 * idx + 0], -2.2760448196157412e-6, eps );
00165     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0000022760448197, eps );
00166     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00167
00168     ++pos;
00169     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00170     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00171     CHECK_REAL_EQUAL( coords[3 * idx + 0], -2.2760448197267635e-6, eps );
00172     CHECK_REAL_EQUAL( coords[3 * idx + 1], -2.2760448196157412e-6, eps );
00173     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00174
00175     // id=16
00176     ++pos;
00177     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00178     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00179     CHECK_REAL_EQUAL( coords[3 * idx + 0], 5.0e-1, eps );
00180     CHECK_REAL_EQUAL( coords[3 * idx + 1], 5.0e-1, eps );
00181     CHECK_REAL_EQUAL( coords[3 * idx + 2], 3.0, eps );
00182
00183     ++pos;
00184     CHECK_EQUAL( pos + 1, sorted_ids[pos] );
00185     idx = std::find( ids.begin(), ids.end(), pos + 1 ) - ids.begin();
00186     CHECK_REAL_EQUAL( coords[3 * idx + 0], 1.0000022760448197, eps );
00187     CHECK_REAL_EQUAL( coords[3 * idx + 1], 1.0000022760448197, eps );
00188     CHECK_REAL_EQUAL( coords[3 * idx + 2], 2.0, eps );
00189 }
00190
00191 void test_read_tets()
00192 {
00193     ErrorCode rval;
00194     Core moab;
00195     Interface& mb = moab;
00196     read_file( moab, example );
00197
00198     std::vector< EntityHandle > tets;
00199     rval = mb.get_entities_by_type( 0, MBTET, tets );CHECK_ERR( rval );
00200     CHECK_EQUAL( (size_t)2, tets.size() );
00201
00202     Tag id_tag = mb.globalId_tag();
00203
00204     std::vector< int > ids( tets.size() );
00205     rval = mb.tag_get_data( id_tag, &tets[0], tets.size(), &ids[0] );CHECK_ERR( rval );
00206
00207     if( ids[0] != 3 )
00208     {
00209         std::swap( ids[0], ids[1] );
00210         std::swap( tets[0], tets[1] );
00211     }
00212
00213     int vtx_ids[4];
00214     const EntityHandle* conn;
00215     int len;
00216
00217     // The first tet has id=3
00218     const int conn1[] = { 13, 14, 15, 16 };
00219     int pos = 0, offset = 3;
00220     CHECK_EQUAL( pos + offset, ids[pos] );
00221     rval = mb.get_connectivity( tets[pos], conn, len );CHECK_ERR( rval );
00222     CHECK_EQUAL( 4, len );
00223     rval = mb.tag_get_data( id_tag, conn, len, vtx_ids );CHECK_ERR( rval );
00224     CHECK_ARRAYS_EQUAL( conn1, 4, vtx_ids, len );
00225
00226     // The second tet has id=4
00227     const int conn2[] = { 13, 17, 14, 16 };
00228     ++pos;
00229     CHECK_EQUAL( pos + offset, ids[pos] );
00230     rval = mb.get_connectivity( tets[pos], conn, len );CHECK_ERR( rval );
00231     CHECK_EQUAL( 4, len );
00232     rval = mb.tag_get_data( id_tag, conn, len, vtx_ids );CHECK_ERR( rval );
00233     CHECK_ARRAYS_EQUAL( conn2, 4, vtx_ids, len );
00234 }
00235
00236 void test_read_hexes()
00237 {
00238     ErrorCode rval;
00239     Core moab;
00240     Interface& mb = moab;
00241     read_file( moab, example );
00242
00243     std::vector< EntityHandle > hexes;
00244     rval = mb.get_entities_by_type( 0, MBHEX, hexes );CHECK_ERR( rval );
00245     CHECK_EQUAL( (size_t)2, hexes.size() );
00246
00247     Tag id_tag = mb.globalId_tag();
00248
00249     std::vector< int > ids( hexes.size() );
00250     rval = mb.tag_get_data( id_tag, &hexes[0], hexes.size(), &ids[0] );CHECK_ERR( rval );
00251
00252     int vtx_ids[8];
00253     const EntityHandle* conn;
00254     int len;
00255
00256     const int conn1[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
00257     int pos = 0, offset = 1;
00258     // Element id 1 is a hex
00259     CHECK_EQUAL( pos + offset, ids[pos] );
00260     rval = mb.get_connectivity( hexes[pos], conn, len );CHECK_ERR( rval );
00261     CHECK_EQUAL( 8, len );
00262     rval = mb.tag_get_data( id_tag, conn, len, vtx_ids );CHECK_ERR( rval );
00263     CHECK_ARRAYS_EQUAL( conn1, 8, vtx_ids, len );
00264
00265     const int conn2[] = { 2, 9, 10, 3, 6, 11, 12, 7 };
00266     ++pos;
00267     // Element id 2 is a hex
00268     CHECK_EQUAL( pos + offset, ids[pos] );
00269     rval = mb.get_connectivity( hexes[pos], conn, len );CHECK_ERR( rval );
00270     CHECK_EQUAL( 8, len );
00271     rval = mb.tag_get_data( id_tag, conn, len, vtx_ids );CHECK_ERR( rval );
00272     CHECK_ARRAYS_EQUAL( conn2, 8, vtx_ids, len );
00273 }
00274
00275 // Two tets and two hexes are in material set 100.
00276 void test_read_material_set()
00277 {
00278     ErrorCode rval;
00279     Core moab;
00280     Interface& mb = moab;
00281     read_file( moab, example );
00282
00283     Tag mat_tag;
00284     rval = mb.tag_get_handle( MAT_PROP_TABLE_TAG, 1, MB_TYPE_INTEGER, mat_tag );CHECK_ERR( rval );
00285
00286     Range mat_set;
00287     const int mat_set_id               = 100;
00288     const void* const mat_set_id_val[] = { &mat_set_id };
00289     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &mat_tag, mat_set_id_val, 1, mat_set );CHECK_ERR( rval );
00290     CHECK_EQUAL( 1, (int)mat_set.size() );
00291
00292     std::vector< EntityHandle > elements, contents;
00293     rval = mb.get_entities_by_type( 0, MBTET, elements );CHECK_ERR( rval );
00294     rval = mb.get_entities_by_type( 0, MBHEX, elements );CHECK_ERR( rval );
00295     rval = mb.get_entities_by_handle( mat_set.front(), contents );CHECK_ERR( rval );
00296     std::sort( elements.begin(), elements.end() );
00297     std::sort( contents.begin(), contents.end() );
00298     CHECK_EQUAL( elements, contents );
00299 }
00300
00301 // The tets are in physical set 4, which corresponds to volume 4 in Cubit.
00302 void test_read_physical_set()
00303 {
00304     ErrorCode rval;
00305     Core moab;
00306     Interface& mb = moab;
00307     read_file( moab, example );
00308
00309     Tag phys_tag;
00310     rval = mb.tag_get_handle( PHYS_PROP_TABLE_TAG, 1, MB_TYPE_INTEGER, phys_tag );CHECK_ERR( rval );
00311
00312     Range phys_set;
00313     const int phys_set_id               = 4;
00314     const void* const phys_set_id_val[] = { &phys_set_id };
00315     rval = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &phys_tag, phys_set_id_val, 1, phys_set );CHECK_ERR( rval );
00316     CHECK_EQUAL( 1, (int)phys_set.size() );
00317
00318     std::vector< EntityHandle > tets, contents;
00319     rval = mb.get_entities_by_type( 0, MBTET, tets );CHECK_ERR( rval );
00320     rval = mb.get_entities_by_handle( phys_set.front(), contents );CHECK_ERR( rval );
00321     std::sort( tets.begin(), tets.end() );
00322     std::sort( contents.begin(), contents.end() );
00323     CHECK_EQUAL( tets, contents );
00324 }