MOAB: Mesh Oriented datABase  (version 5.3.0)
test_geom_gqt.cpp
Go to the documentation of this file.
00001 #include "MBTagConventions.hpp"
00002 
00003 #include "moab/Interface.hpp"
00004 #include "moab/Core.hpp"
00005 #include "moab/Range.hpp"
00006 #include "moab/CartVect.hpp"
00007 #include "moab/GeomQueryTool.hpp"
00008 #include "moab/GeomTopoTool.hpp"
00009 
00010 #ifdef MOAB_HAVE_MPI
00011 #include "moab_mpi.h"
00012 #endif
00013 
00014 #include <vector>
00015 #include <iostream>
00016 #include <cmath>
00017 #include <limits>
00018 #include <algorithm>
00019 #include <cstdio>  // for remove()
00020 
00021 #define CHKERR \
00022     if( MB_SUCCESS != rval ) return rval
00023 
00024 const double ROOT2 = 1.4142135623730951;
00025 
00026 using namespace moab;
00027 
00028 // Create file containing geometry for 2x2x2 cube
00029 // centered at origin and having a convex +Z face
00030 // (center of face at origin).
00031 ErrorCode write_geometry( const char* output_file_name );
00032 
00033 ErrorCode test_ray_fire( GeomQueryTool* );
00034 
00035 ErrorCode test_point_in_volume( GeomQueryTool* );
00036 
00037 ErrorCode test_closest_to_location( GeomQueryTool* );
00038 
00039 ErrorCode test_measure_volume( GeomQueryTool* );
00040 
00041 ErrorCode test_measure_area( GeomQueryTool* );
00042 
00043 ErrorCode test_surface_sense( GeomQueryTool* );
00044 
00045 ErrorCode overlap_write_geometry( const char* output_file_name );
00046 ErrorCode overlap_test_ray_fire( GeomQueryTool* );
00047 ErrorCode overlap_test_point_in_volume( GeomQueryTool* );
00048 ErrorCode overlap_test_measure_volume( GeomQueryTool* );
00049 ErrorCode overlap_test_measure_area( GeomQueryTool* );
00050 ErrorCode overlap_test_surface_sense( GeomQueryTool* );
00051 ErrorCode overlap_test_tracking( GeomQueryTool* );
00052 
00053 int run_regular_tests( GeomQueryTool* gqt );
00054 int run_overlap_tests( GeomQueryTool* gqt );
00055 
00056 ErrorCode write_geometry( const char* output_file_name )
00057 {
00058     ErrorCode rval;
00059 
00060     Interface* moab = new Core();
00061 
00062     // Define a 2x2x2 cube centered at orgin
00063     // with concavity in +Z face.
00064     const double coords[]    = { 1, -1, -1, 1, 1,  -1, -1, 1,  -1, -1, -1, -1, 1, -1,
00065                               1, 1,  1,  1, -1, 1,  1,  -1, -1, 1,  0,  0,  0 };
00066     const int connectivity[] = {
00067         0, 3, 1, 3, 2, 1,  // -Z
00068         0, 1, 4, 5, 4, 1,  // +X
00069         1, 2, 6, 6, 5, 1,  // +Y
00070         6, 2, 3, 7, 6, 3,  // -X
00071         0, 4, 3, 7, 3, 4,  // -Y
00072         4, 5, 8, 5, 6, 8,  // +Z
00073         6, 7, 8, 7, 4, 8   // +Z
00074     };
00075     const unsigned tris_per_surf[] = { 2, 2, 2, 2, 2, 4 };
00076 
00077     // Create the geometry
00078     const unsigned num_verts = sizeof( coords ) / ( 3 * sizeof( double ) );
00079     const unsigned num_tris  = sizeof( connectivity ) / ( 3 * sizeof( int ) );
00080     const unsigned num_surfs = sizeof( tris_per_surf ) / sizeof( unsigned );
00081     EntityHandle verts[num_verts], tris[num_tris], surfs[num_surfs];
00082     for( unsigned i = 0; i < num_verts; ++i )
00083     {
00084         rval = moab->create_vertex( coords + 3 * i, verts[i] );CHKERR;
00085     }
00086     for( unsigned i = 0; i < num_tris; ++i )
00087     {
00088         const EntityHandle conn[] = { verts[connectivity[3 * i]], verts[connectivity[3 * i + 1]],
00089                                       verts[connectivity[3 * i + 2]] };
00090         rval                      = moab->create_element( MBTRI, conn, 3, tris[i] );CHKERR;
00091     }
00092 
00093     // create CAD topology
00094     EntityHandle* tri_iter = tris;
00095     for( unsigned i = 0; i < num_surfs; ++i )
00096     {
00097         rval = moab->create_meshset( MESHSET_SET, surfs[i] );CHKERR;
00098         rval = moab->add_entities( surfs[i], tri_iter, tris_per_surf[i] );CHKERR;
00099         tri_iter += tris_per_surf[i];
00100     }
00101 
00102     Tag dim_tag, id_tag, sense_tag;
00103     rval = moab->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag, MB_TAG_SPARSE | MB_TAG_CREAT );CHKERR;
00104     id_tag = moab->globalId_tag();
00105     rval   = moab->tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, sense_tag, MB_TAG_SPARSE | MB_TAG_CREAT );CHKERR;
00106 
00107     std::vector< int > dims( num_surfs, 2 );
00108     rval = moab->tag_set_data( dim_tag, surfs, num_surfs, &dims[0] );CHKERR;
00109     std::vector< int > ids( num_surfs );
00110     for( size_t i = 0; i < ids.size(); ++i )
00111         ids[i] = i + 1;
00112     rval = moab->tag_set_data( id_tag, surfs, num_surfs, &ids[0] );CHKERR;
00113 
00114     EntityHandle volume;
00115     rval = moab->create_meshset( MESHSET_SET, volume );CHKERR;
00116     for( unsigned i = 0; i < num_surfs; ++i )
00117     {
00118         rval = moab->add_parent_child( volume, surfs[i] );CHKERR;
00119     }
00120 
00121     std::vector< EntityHandle > senses( 2 * num_surfs, 0 );
00122     for( size_t i = 0; i < senses.size(); i += 2 )
00123         senses[i] = volume;
00124     rval = moab->tag_set_data( sense_tag, surfs, num_surfs, &senses[0] );CHKERR;
00125 
00126     const int three = 3;
00127     const int one   = 1;
00128     rval            = moab->tag_set_data( dim_tag, &volume, 1, &three );CHKERR;
00129     rval = moab->tag_set_data( id_tag, &volume, 1, &one );CHKERR;
00130 
00131     rval = moab->write_mesh( output_file_name );CHKERR;
00132     delete moab;
00133     return MB_SUCCESS;
00134 }
00135 
00136 ErrorCode overlap_write_geometry( const char* output_file_name )
00137 {
00138     ErrorCode rval;
00139     Interface* moab = new Core();
00140 
00141     // Define two 1x2x2 cubes that overlap from 0 <= x <= 0.01
00142     // cube 0 centered at (0.5,0,0)
00143     const double coords[]    = { 1, -1, -1, 1, 1, -1, 0, 1, -1, 0, -1, -1, 1, -1, 1, 1, 1, 1, 0, 1, 1, 0, -1, 1,
00144                               // cube 1 centered near (-0.5,0,0)
00145                               0.01, -1, -1, 0.01, 1, -1, -1, 1, -1, -1, -1, -1, 0.01, -1, 1, 0.01, 1, 1, -1, 1, 1, -1,
00146                               -1, 1 };
00147     const int connectivity[] = { 0, 3, 1, 3, 2, 1,    // -Z
00148                                  0, 1, 4, 5, 4, 1,    // +X
00149                                  1, 2, 6, 6, 5, 1,    // +Y
00150                                  6, 2, 3, 7, 6, 3,    // -X
00151                                  0, 4, 3, 7, 3, 4,    // -Y
00152                                  4, 5, 6, 6, 7, 4 };  // +Z
00153 
00154     // Create the geometry
00155     const unsigned tris_per_surf = 2;
00156     const unsigned num_cubes     = 2;
00157     const unsigned num_verts     = sizeof( coords ) / ( 3 * sizeof( double ) );
00158     const unsigned num_tris      = num_cubes * sizeof( connectivity ) / ( 3 * sizeof( int ) );
00159     const unsigned num_surfs     = num_tris / tris_per_surf;
00160     EntityHandle verts[num_verts], tris[num_tris], surfs[num_surfs];
00161 
00162     for( unsigned i = 0; i < num_verts; ++i )
00163     {
00164         rval = moab->create_vertex( coords + 3 * i, verts[i] );CHKERR;
00165     }
00166     // cube0
00167     for( unsigned i = 0; i < num_tris / 2; ++i )
00168     {
00169         const EntityHandle conn[] = { verts[connectivity[3 * i]], verts[connectivity[3 * i + 1]],
00170                                       verts[connectivity[3 * i + 2]] };
00171         rval                      = moab->create_element( MBTRI, conn, 3, tris[i] );CHKERR;
00172     }
00173     // cube1
00174     for( unsigned i = 0; i < num_tris / 2; ++i )
00175     {
00176         const EntityHandle conn[] = { verts[8 + connectivity[3 * i]], verts[8 + connectivity[3 * i + 1]],
00177                                       verts[8 + connectivity[3 * i + 2]] };
00178         rval                      = moab->create_element( MBTRI, conn, 3, tris[num_tris / 2 + i] );CHKERR;
00179     }
00180 
00181     // create CAD topology
00182     EntityHandle* tri_iter = tris;
00183     for( unsigned i = 0; i < num_surfs; ++i )
00184     {
00185         rval = moab->create_meshset( MESHSET_SET, surfs[i] );CHKERR;
00186         rval = moab->add_entities( surfs[i], tri_iter, tris_per_surf );CHKERR;
00187         tri_iter += tris_per_surf;
00188     }
00189 
00190     Tag dim_tag, id_tag, sense_tag;
00191     rval = moab->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dim_tag, MB_TAG_SPARSE | MB_TAG_CREAT );CHKERR;
00192     id_tag = moab->globalId_tag();
00193     rval   = moab->tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, sense_tag, MB_TAG_SPARSE | MB_TAG_CREAT );CHKERR;
00194 
00195     std::vector< int > dims( num_surfs, 2 );
00196     rval = moab->tag_set_data( dim_tag, surfs, num_surfs, &dims[0] );CHKERR;
00197     std::vector< int > ids( num_surfs );
00198     for( size_t i = 0; i < ids.size(); ++i )
00199         ids[i] = i + 1;
00200     rval = moab->tag_set_data( id_tag, surfs, num_surfs, &ids[0] );CHKERR;
00201 
00202     EntityHandle volume;
00203     rval = moab->create_meshset( MESHSET_SET, volume );CHKERR;
00204     for( unsigned i = 0; i < num_surfs; ++i )
00205     {
00206         rval = moab->add_parent_child( volume, surfs[i] );CHKERR;
00207     }
00208 
00209     std::vector< EntityHandle > senses( 2 * num_surfs, 0 );
00210     for( size_t i = 0; i < senses.size(); i += 2 )
00211         senses[i] = volume;
00212     rval = moab->tag_set_data( sense_tag, surfs, num_surfs, &senses[0] );CHKERR;
00213 
00214     const int three = 3;
00215     const int one   = 1;
00216     rval            = moab->tag_set_data( dim_tag, &volume, 1, &three );CHKERR;
00217     rval = moab->tag_set_data( id_tag, &volume, 1, &one );CHKERR;
00218 
00219     rval = moab->write_mesh( output_file_name );CHKERR;
00220     delete moab;
00221     return MB_SUCCESS;
00222 }
00223 
00224 #define RUN_TEST( A )                              \
00225     do                                             \
00226     {                                              \
00227         std::cout << #A << "... " << std::endl;    \
00228         if( MB_SUCCESS != A( gqt ) ) { ++errors; } \
00229     } while( false )
00230 
00231 int main( int argc, char* argv[] )
00232 {
00233     ErrorCode rval;
00234     const char* filename = "test_geom.h5m";
00235 
00236 #ifdef MOAB_HAVE_MPI
00237     int fail = MPI_Init( &argc, &argv );
00238     if( fail ) return fail;
00239 #endif
00240 
00241     rval = write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file: " << filename );
00242 
00243     Interface* MBI = new Core();
00244 
00245     int errors = 0;
00246     rval       = MBI->load_file( filename );
00247     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file" );
00248 
00249     GeomTopoTool* gtt  = new GeomTopoTool( MBI );
00250     GeomQueryTool* gqt = new GeomQueryTool( gtt );
00251 
00252     errors += run_regular_tests( gqt );
00253 
00254     // clear out moab instance
00255     rval = MBI->delete_mesh();MB_CHK_SET_ERR( rval, "Failed to delete mesh" );
00256 
00257     delete gtt;
00258     delete gqt;
00259 
00260     // Now load a different geometry: two cubes that slightly overlap
00261     rval = overlap_write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file: " << filename );
00262 
00263     rval = MBI->load_file( filename );
00264     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file with overlaps" );
00265 
00266     gtt = new GeomTopoTool( MBI );
00267     gqt = new GeomQueryTool( gtt );
00268 
00269     errors += run_overlap_tests( gqt );
00270 
00271     // clear moab instance
00272     rval = MBI->delete_mesh();MB_CHK_SET_ERR( rval, "Failed to delete mesh" );
00273 
00274     delete gtt;
00275     delete gqt;
00276 
00277     // Re-run all tests with the alternate constructor
00278     std::cout << "-------------------------------------" << std::endl;
00279     std::cout << "Re-running tests with MBI constructor" << std::endl;
00280     std::cout << "-------------------------------------" << std::endl;
00281 
00282     rval = write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file" );
00283 
00284     rval = MBI->load_file( filename );
00285     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file" );
00286 
00287     gqt = new GeomQueryTool( MBI );
00288 
00289     errors += run_regular_tests( gqt );
00290 
00291     // clear moab and dagmc instance
00292     rval = MBI->delete_mesh();MB_CHK_SET_ERR( rval, "Failed to delete mesh" );
00293 
00294     delete gqt;
00295 
00296     // Now load a different geometry: two cubes that slightly overlap
00297     rval = overlap_write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file: " );
00298 
00299     rval = MBI->load_file( filename );
00300     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file with overlaps." );
00301 
00302     gqt = new GeomQueryTool( MBI );
00303 
00304     errors += run_overlap_tests( gqt );
00305 
00306     // final cleanup
00307     delete gqt;
00308     delete MBI;
00309 
00310 #ifdef MOAB_HAVE_MPI
00311     fail = MPI_Finalize();
00312     if( fail ) return fail;
00313 #endif
00314 
00315     return errors;
00316 }
00317 
00318 int run_regular_tests( GeomQueryTool* gqt )
00319 {
00320     int errors = 0;
00321     ErrorCode rval;
00322 
00323     rval = gqt->initialize();MB_CHK_SET_ERR_CONT( rval, "Failed to initialize the GeomQueryTool" );
00324 
00325     RUN_TEST( test_ray_fire );
00326     RUN_TEST( test_closest_to_location );
00327     RUN_TEST( test_point_in_volume );
00328     RUN_TEST( test_measure_volume );
00329     RUN_TEST( test_measure_area );
00330     RUN_TEST( test_surface_sense );
00331 
00332     // change settings to use overlap-tolerant mode (arbitrary thickness)
00333     double overlap_thickness = 0.1;
00334     gqt->set_overlap_thickness( overlap_thickness );
00335     RUN_TEST( test_ray_fire );
00336     RUN_TEST( test_point_in_volume );
00337 
00338     return errors;
00339 }
00340 
00341 int run_overlap_tests( GeomQueryTool* gqt )
00342 {
00343     int errors = 0;
00344     ErrorCode rval;
00345 
00346     rval = gqt->initialize();MB_CHK_SET_ERR_CONT( rval, "Failed to initialize the GeomQueryTool" );
00347 
00348     // change settings to use overlap-tolerant mode (with a large enough thickness)
00349     double overlap_thickness = 3.0;
00350     gqt->set_overlap_thickness( overlap_thickness );
00351     RUN_TEST( overlap_test_ray_fire );
00352     RUN_TEST( overlap_test_point_in_volume );
00353     RUN_TEST( overlap_test_measure_volume );
00354     RUN_TEST( overlap_test_measure_area );
00355     RUN_TEST( overlap_test_surface_sense );
00356     RUN_TEST( overlap_test_tracking );
00357 
00358     return errors;
00359 }
00360 
00361 ErrorCode test_surface_sense( GeomQueryTool* gqt )
00362 {
00363     ErrorCode rval;
00364     Interface* moab = gqt->moab_instance();
00365 
00366     Tag dim_tag = gqt->gttool()->get_geom_tag();
00367     Range surfs, vols;
00368     const int two = 2, three = 3;
00369     const void* ptrs[] = { &two, &three };
00370     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00371     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00372 
00373     if( vols.size() != 2 )
00374     {
00375         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00376         return MB_FAILURE;
00377     }
00378     if( surfs.size() != 6 )
00379     {
00380         std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
00381         return MB_FAILURE;
00382     }
00383 
00384     for( Range::iterator i = surfs.begin(); i != surfs.end(); ++i )
00385     {
00386         int sense = 0;
00387         rval      = gqt->gttool()->get_sense( *i, vols.front(), sense );
00388         if( MB_SUCCESS != rval || sense != 1 )
00389         {
00390             std::cerr << "ERROR: Expected 1 for surface sense, got " << sense << std::endl;
00391             return MB_FAILURE;
00392         }
00393     }
00394 
00395     return MB_SUCCESS;
00396 }
00397 
00398 ErrorCode overlap_test_surface_sense( GeomQueryTool* gqt )
00399 {
00400     ErrorCode rval;
00401     Interface* moab = gqt->moab_instance();
00402 
00403     Tag dim_tag = gqt->gttool()->get_geom_tag();
00404     Range surfs, vols;
00405     const int two = 2, three = 3;
00406     const void* ptrs[] = { &two, &three };
00407     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00408     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00409 
00410     if( vols.size() != 2 )
00411     {
00412         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00413         return MB_FAILURE;
00414     }
00415     if( surfs.size() != 12 )
00416     {
00417         std::cerr << "ERROR: Expected 12 surfaces in input, found " << surfs.size() << std::endl;
00418         return MB_FAILURE;
00419     }
00420 
00421     for( Range::iterator i = surfs.begin(); i != surfs.end(); ++i )
00422     {
00423         int sense = 0;
00424         rval      = gqt->gttool()->get_sense( *i, vols.front(), sense );
00425         if( MB_SUCCESS != rval || sense != 1 )
00426         {
00427             std::cerr << "ERROR: Expected 1 for surface sense, got " << sense << std::endl;
00428             return MB_FAILURE;
00429         }
00430     }
00431 
00432     return MB_SUCCESS;
00433 }
00434 ErrorCode test_measure_volume( GeomQueryTool* gqt )
00435 {
00436     ErrorCode rval;
00437     Interface* moab = gqt->moab_instance();
00438 
00439     Tag dim_tag = gqt->gttool()->get_geom_tag();
00440     Range vols;
00441     const int three = 3;
00442     const void* ptr = &three;
00443     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00444 
00445     if( vols.size() != 2 )
00446     {
00447         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00448         return MB_FAILURE;
00449     }
00450 
00451     // input file is 2x2x2 cube with convexity in +Z face that touches the origin.
00452     // expected volume is 8 (2x2x2) less the volume of the pyrimid concavity
00453     double result;
00454     const double vol = 2 * 2 * 2 - 1 * 4. / 3;
00455 
00456     rval = gqt->measure_volume( vols.front(), result );CHKERR;
00457     if( fabs( result - vol ) > 10 * std::numeric_limits< double >::epsilon() )
00458     {
00459         std::cerr << "ERROR: Expected " << vol << " as measure of volume, got " << result << std::endl;
00460         return MB_FAILURE;
00461     }
00462 
00463     return MB_SUCCESS;
00464 }
00465 ErrorCode overlap_test_measure_volume( GeomQueryTool* gqt )
00466 {
00467     ErrorCode rval;
00468     Interface* moab = gqt->moab_instance();
00469 
00470     Tag dim_tag = gqt->gttool()->get_geom_tag();
00471     Range vols;
00472     const int three = 3;
00473     const void* ptr = &three;
00474     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00475 
00476     if( vols.size() != 2 )
00477     {
00478         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00479         return MB_FAILURE;
00480     }
00481 
00482     // The volume has two regions that overlap.
00483     double result;
00484     const double vol = ( 1 + 1.01 ) * 2 * 2;
00485 
00486     rval = gqt->measure_volume( vols.front(), result );CHKERR;
00487     if( fabs( result - vol ) > 2 * std::numeric_limits< double >::epsilon() )
00488     {
00489         std::cerr << "ERROR: Expected " << vol << " as measure of volume, got " << result << std::endl;
00490         return MB_FAILURE;
00491     }
00492 
00493     return MB_SUCCESS;
00494 }
00495 
00496 ErrorCode test_measure_area( GeomQueryTool* gqt )
00497 {
00498     ErrorCode rval;
00499     Interface* moab = gqt->moab_instance();
00500 
00501     Tag dim_tag = gqt->gttool()->get_geom_tag();
00502     Range surfs;
00503     const int two   = 2;
00504     const void* ptr = &two;
00505     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, surfs );CHKERR;
00506 
00507     if( surfs.size() != 6 )
00508     {
00509         std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
00510         return MB_FAILURE;
00511     }
00512 
00513     int ids[6];
00514     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00515 
00516     // expect area of 4 for all faces except face 6.
00517     // face 6 should have area == 4*sqrt(2)
00518     Range::iterator iter = surfs.begin();
00519     for( unsigned i = 0; i < 6; ++i, ++iter )
00520     {
00521         double expected = 4.0;
00522         if( ids[i] == 6 ) expected *= ROOT2;
00523 
00524         double result;
00525 
00526         rval = gqt->measure_area( *iter, result );CHKERR;
00527         if( fabs( result - expected ) > std::numeric_limits< double >::epsilon() )
00528         {
00529             std::cerr << "ERROR: Expected area of surface " << ids[i] << " to be " << expected << ".  Got " << result
00530                       << std::endl;
00531             return MB_FAILURE;
00532         }
00533     }
00534 
00535     return MB_SUCCESS;
00536 }
00537 
00538 ErrorCode overlap_test_measure_area( GeomQueryTool* gqt )
00539 {
00540     ErrorCode rval;
00541     Interface* moab = gqt->moab_instance();
00542 
00543     Tag dim_tag = gqt->gttool()->get_geom_tag();
00544     Range surfs;
00545     const int two   = 2;
00546     const void* ptr = &two;
00547     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, surfs );CHKERR;
00548 
00549     const unsigned num_surfs = 12;
00550     if( surfs.size() != num_surfs )
00551     {
00552         std::cerr << "ERROR: Expected " << num_surfs << " surfaces in input, found " << surfs.size() << std::endl;
00553         return MB_FAILURE;
00554     }
00555 
00556     int ids[num_surfs];
00557     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00558 
00559     const double x_area   = 2 * 2;
00560     const double yz_area0 = 2 * 1;
00561     const double yz_area1 = 2 * 1.01;
00562     Range::iterator iter  = surfs.begin();
00563     for( unsigned i = 0; i < num_surfs; ++i, ++iter )
00564     {
00565         double expected, result;
00566         if( 0 == i || 2 == i || 4 == i || 5 == i )
00567             expected = yz_area0;
00568         else if( 1 == i || 3 == i || 7 == i || 9 == i )
00569             expected = x_area;
00570         else if( 6 == i || 8 == i || 10 == i || 11 == i )
00571             expected = yz_area1;
00572 
00573         rval = gqt->measure_area( *iter, result );CHKERR;
00574         if( fabs( result - expected ) > std::numeric_limits< double >::epsilon() )
00575         {
00576             std::cerr << "ERROR: Expected area of surface " << ids[i] << " to be " << expected << ".  Got " << result
00577                       << std::endl;
00578             return MB_FAILURE;
00579         }
00580     }
00581 
00582     return MB_SUCCESS;
00583 }
00584 
00585 struct ray_fire
00586 {
00587     int prev_surf;
00588     double origin[3], direction[3];
00589     int hit_surf;
00590     double distance;
00591 };
00592 
00593 ErrorCode test_ray_fire( GeomQueryTool* gqt )
00594 {
00595     // Glancing ray-triangle intersections are not valid exit intersections.
00596     // Piercing ray-triangle intersections are valid exit intersections.
00597     // "0" destination surface implies that it is ambiguous.
00598     const struct ray_fire tests[] = { /* src_srf origin               direction                 dest dist */
00599                                       // piercing edge
00600                                       { 1, { 0.0, 0.0, -1. }, { -1.0 / ROOT2, 0.0, 1.0 / ROOT2 }, 4, ROOT2 },
00601                                       // piercing edge
00602                                       { 1, { 0.0, 0.0, -1. }, { 1.0 / ROOT2, 0.0, 1.0 / ROOT2 }, 2, ROOT2 },
00603                                       // piercing edge
00604                                       { 1, { 0.0, 0.0, -1. }, { 0.0, 1.0 / ROOT2, 1.0 / ROOT2 }, 3, ROOT2 },
00605                                       // piercing edge
00606                                       { 1, { 0.5, 0.5, -1. }, { 0.0, 0.0, 1.0 }, 6, 1.5 },
00607                                       // interior
00608                                       { 2, { 1.0, 0.0, 0.5 }, { -1.0, 0.0, 0.0 }, 6, 0.5 },
00609                                       // glancing node then piercing edge
00610                                       { 2, { 1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 2.0 },
00611                                       // piercing node
00612                                       { 1, { 0.0, 0.0, -1. }, { 0.0, 0.0, 1.0 }, 6, 1.0 },
00613                                       // glancing edge then interior
00614                                       { 2, { 1.0, 0.0, 0.5 }, { -1.0 / ROOT2, 1.0 / ROOT2, 0.0 }, 3, ROOT2 }
00615     };
00616 
00617     ErrorCode rval;
00618     Interface* moab = gqt->moab_instance();
00619 
00620     Tag dim_tag = gqt->gttool()->get_geom_tag();
00621     Range surfs, vols;
00622     const int two = 2, three = 3;
00623     const void* ptrs[] = { &two, &three };
00624     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00625     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00626 
00627     if( vols.size() != 2 )
00628     {
00629         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00630         return MB_FAILURE;
00631     }
00632     if( surfs.size() != 6 )
00633     {
00634         std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
00635         return MB_FAILURE;
00636     }
00637 
00638     int ids[6];
00639     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00640     EntityHandle surf[6];
00641     std::copy( surfs.begin(), surfs.end(), surf );
00642 
00643     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00644     for( int i = 0; i < num_test; ++i )
00645     {
00646         int* ptr = std::find( ids, ids + 6, tests[i].prev_surf );
00647         int idx  = ptr - ids;
00648         if( idx >= 6 )
00649         {
00650             std::cerr << "Surface " << tests[i].prev_surf << " not found." << std::endl;
00651             return MB_FAILURE;
00652         }
00653         // const EntityHandle src_surf = surf[idx];
00654 
00655         ptr = std::find( ids, ids + 6, tests[i].hit_surf );
00656         idx = ptr - ids;
00657         if( idx >= 6 )
00658         {
00659             std::cerr << "Surface " << tests[i].hit_surf << " not found." << std::endl;
00660             return MB_FAILURE;
00661         }
00662         const EntityHandle hit_surf = surf[idx];
00663 
00664         double dist;
00665         EntityHandle result;
00666         GeomQueryTool::RayHistory history;
00667         rval = gqt->ray_fire( vols.front(), tests[i].origin, tests[i].direction, result, dist, &history );
00668 
00669         if( result != hit_surf || fabs( dist - tests[i].distance ) > 1e-6 )
00670         {
00671             EntityHandle* p = std::find( surf, surf + 6, result );
00672             idx             = p - surf;
00673             int id          = idx > 5 ? 0 : ids[idx];
00674 
00675             std::cerr << "Rayfire test failed for " << std::endl
00676                       << "\t ray from (" << tests[i].origin[0] << ", " << tests[i].origin[1] << ", "
00677                       << tests[i].origin[2] << ") going [" << tests[i].direction[0] << ", " << tests[i].direction[1]
00678                       << ", " << tests[i].direction[2] << "]" << std::endl
00679                       << "\t Beginning on surface " << tests[i].prev_surf << std::endl
00680                       << "\t Expected to hit surface " << tests[i].hit_surf << " after " << tests[i].distance
00681                       << " units." << std::endl
00682                       << "\t Actually hit surface " << id << " after " << dist << " units." << std::endl;
00683             return MB_FAILURE;
00684         }
00685 
00686         CartVect loc = CartVect( tests[i].origin ) + ( dist * CartVect( tests[i].direction ) );
00687 
00688         std::vector< std::pair< int, GeomQueryTool::RayHistory* > > boundary_tests;
00689         boundary_tests.push_back( std::make_pair( 1, &history ) );
00690         boundary_tests.push_back( std::make_pair( 0, &history ) );
00691         boundary_tests.push_back( std::make_pair( 1, (GeomQueryTool::RayHistory*)NULL ) );
00692         boundary_tests.push_back( std::make_pair( 0, (GeomQueryTool::RayHistory*)NULL ) );
00693 
00694         for( unsigned int bt = 0; bt < boundary_tests.size(); ++bt )
00695         {
00696 
00697             int expected                 = boundary_tests[bt].first;
00698             GeomQueryTool::RayHistory* h = boundary_tests[bt].second;
00699 
00700             // pick the direction based on expected result of test. Either reuse the ray_fire
00701             // vector, or reverse it to check for a vector that enters the cell
00702             CartVect uvw( tests[i].direction );
00703             if( expected == 1 ) uvw = -uvw;
00704 
00705             int boundary_result = -1;
00706 
00707             rval = gqt->test_volume_boundary( vols.front(), result, loc.array(), uvw.array(), boundary_result, h );
00708 
00709             if( boundary_result != expected )
00710             {
00711                 std::cerr << "DagMC::test_volume_boundary failed (" << ( ( expected == 0 ) ? "+" : "-" ) << " dir,"
00712                           << ( ( h ) ? "+" : "-" ) << " history, i=" << i << ")" << std::endl;
00713                 return MB_FAILURE;
00714             }
00715         }
00716     }
00717 
00718     return MB_SUCCESS;
00719 }
00720 
00721 ErrorCode overlap_test_ray_fire( GeomQueryTool* gqt )
00722 {
00723     // Glancing ray-triangle intersections are not valid exit intersections.
00724     // Piercing ray-triangle intersections are valid exit intersections.
00725     // "0" destination surface implies that it is ambiguous.
00726     const struct ray_fire tests[] = { /*          ____________________
00727                                                   |                   |
00728                                          -x <---  | region1 | overlap | region0 |   ---> +x
00729                                                             |                   |
00730                                                             ---------------------
00731                                          -x <--- -1.0      0.0      0.01        1.0 ---> +x
00732                                        surf_id:   10        4         8         2
00733 
00734                                        The zero-distance advance would occur in the implicit volume between
00735                                        these two regions---not in this volume.
00736 
00737                                        src_srf origin               direction                 dest dist */
00738                                       // numerical location on surface 1 of region0
00739                                       { 4, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 2, 0.0 },
00740                                       { 2, { 1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 1.0 },
00741                                       // numerical location inside region0
00742                                       { 4, { 0.5, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 2, 0.5 },
00743                                       { 2, { 0.5, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.5 },
00744                                       // numerical location on surface 7 of region1
00745                                       { 4, { 0.01, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.0 },
00746                                       { 2, { 0.01, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.01 },
00747                                       // numerical location inside overlap
00748                                       { 10, { 0.005, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.005 },
00749                                       { 2, { 0.005, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.005 },
00750                                       // numerical location on surface 3 of region0
00751                                       { 10, { 0.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.01 },
00752                                       { 8, { 0.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.0 },
00753                                       // numerical location inside region1
00754                                       { 10, { -0.5, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.51 },
00755                                       { 8, { -0.5, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 10, 0.5 },
00756                                       // numerical location on surface 9 of region1
00757                                       { 10, { -1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 1.01 },
00758                                       { 8, { -1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 10, 0.0 }
00759     };
00760 
00761     ErrorCode rval;
00762     Interface* moab = gqt->moab_instance();
00763 
00764     Tag dim_tag = gqt->gttool()->get_geom_tag();
00765     Range surfs, vols;
00766     const int two = 2, three = 3;
00767     const void* ptrs[] = { &two, &three };
00768     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00769     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00770 
00771     if( vols.size() != 2 )
00772     {
00773         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00774         return MB_FAILURE;
00775     }
00776     const unsigned num_surf = 12;
00777     if( surfs.size() != num_surf )
00778     {
00779         std::cerr << "ERROR: Expected " << num_surf << " surfaces in input, found " << surfs.size() << std::endl;
00780         return MB_FAILURE;
00781     }
00782 
00783     int ids[num_surf];
00784     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00785     EntityHandle surf[num_surf];
00786     std::copy( surfs.begin(), surfs.end(), surf );
00787 
00788     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00789     for( int i = 0; i < num_test; ++i )
00790     {
00791         int* ptr     = std::find( ids, ids + num_surf, tests[i].prev_surf );
00792         unsigned idx = ptr - ids;
00793         if( idx >= num_surf )
00794         {
00795             std::cerr << "Surface " << tests[i].prev_surf << " not found." << std::endl;
00796             return MB_FAILURE;
00797         }
00798         // const EntityHandle src_surf = surf[idx];
00799 
00800         ptr = std::find( ids, ids + num_surf, tests[i].hit_surf );
00801         idx = ptr - ids;
00802         if( idx >= num_surf )
00803         {
00804             std::cerr << "Surface " << tests[i].hit_surf << " not found." << std::endl;
00805             return MB_FAILURE;
00806         }
00807         const EntityHandle hit_surf = surf[idx];
00808 
00809         double dist;
00810         EntityHandle result;
00811         rval = gqt->ray_fire( vols.front(), tests[i].origin, tests[i].direction, result, dist );
00812 
00813         if( result != hit_surf || fabs( dist - tests[i].distance ) > 1e-6 )
00814         {
00815             EntityHandle* p = std::find( surf, surf + 6, result );
00816             idx             = p - surf;
00817             int id          = idx > num_surf - 1 ? 0 : ids[idx];
00818 
00819             std::cerr << "Rayfire test failed for " << std::endl
00820                       << "\t ray from (" << tests[i].origin[0] << ", " << tests[i].origin[1] << ", "
00821                       << tests[i].origin[2] << ") going [" << tests[i].direction[0] << ", " << tests[i].direction[1]
00822                       << ", " << tests[i].direction[2] << "]" << std::endl
00823                       << "\t Beginning on surface " << tests[i].prev_surf << std::endl
00824                       << "\t Expected to hit surface " << tests[i].hit_surf << " after " << tests[i].distance
00825                       << " units." << std::endl
00826                       << "\t Actually hit surface " << id << " after " << dist << " units." << std::endl;
00827             return MB_FAILURE;
00828         }
00829     }
00830 
00831     return MB_SUCCESS;
00832 }
00833 
00834 struct PointInVol
00835 {
00836     double coords[3];
00837     int result;
00838     double dir[3];
00839 };
00840 
00841 ErrorCode test_point_in_volume( GeomQueryTool* gqt )
00842 {
00843     const char* const NAME_ARR[] = { "Boundary", "Outside", "Inside" };
00844     const char* const* names     = NAME_ARR + 1;
00845     const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
00846     const struct PointInVol tests[] = { { { 0.0, 0.0, 0.5 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00847                                         { { 0.0, 0.0, -0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
00848                                         { { 0.7, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00849                                         { { -0.7, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00850                                         { { 0.0, -0.7, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00851                                         { { 0.0, -0.7, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00852                                         { { 1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00853                                         { { -1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00854                                         { { -1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00855                                         { { 1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00856                                         { { 1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00857                                         { { -1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00858                                         { { -1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00859                                         { { 1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00860                                         // Add some directions to test special cases of edge/node intersection
00861                                         { { 0.5, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
00862                                         { { 0.5, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
00863                                         { { 0.0, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
00864                                         { { 0.5, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
00865                                         { { 0.5, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } } };
00866 
00867     //    { { 1.0, 0.0, 0.0 }, BOUNDARY}, MCNP doesn't return on boundary
00868     //{ {-1.0, 0.0, 0.0 }, BOUNDARY},
00869     //{ { 0.0, 1.0, 0.0 }, BOUNDARY},
00870     //{ { 0.0,-1.0, 0.0 }, BOUNDARY},
00871     //{ { 0.0, 0.0, 0.0 }, BOUNDARY},
00872     //{ { 0.0, 0.0,-1.0 }, BOUNDARY} };
00873     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00874 
00875     ErrorCode rval;
00876     Interface* moab = gqt->moab_instance();
00877 
00878     Tag dim_tag = gqt->gttool()->get_geom_tag();
00879 
00880     Range vols;
00881     const int three = 3;
00882     const void* ptr = &three;
00883     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00884     if( vols.size() != 2 )
00885     {
00886         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00887         return MB_FAILURE;
00888     }
00889     const EntityHandle vol = vols.front();
00890 
00891     for( int i = 0; i < num_test; ++i )
00892     {
00893         int result;
00894         rval = gqt->point_in_volume( vol, tests[i].coords, result, tests[i].dir );CHKERR;
00895         if( result != tests[i].result )
00896         {
00897             std::cerr << "ERROR testing point_in_volume[" << i << "]:" << std::endl
00898                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
00899                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
00900             return MB_FAILURE;
00901         }
00902 
00903         // point_in_volume_slow doesn't to boundary.
00904         if( tests[i].result == BOUNDARY ) continue;
00905 
00906         rval = gqt->point_in_volume_slow( vol, tests[i].coords, result );CHKERR;
00907 
00908         if( result != tests[i].result )
00909         {
00910             std::cerr << "ERROR testing point_in_volume_slow[" << i << "]:" << std::endl
00911                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
00912                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
00913             return MB_FAILURE;
00914         }
00915     }
00916 
00917     return MB_SUCCESS;
00918 }
00919 
00920 ErrorCode test_closest_to_location( GeomQueryTool* gqt )
00921 {
00922     ErrorCode rval;
00923     Interface* moab = gqt->moab_instance();
00924 
00925     Range vols;
00926     const int three = 3;
00927     const void* ptr = &three;
00928     Tag dim_tag     = gqt->gttool()->get_geom_tag();
00929     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00930     if( vols.size() != 2 )
00931     {
00932         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00933         return MB_FAILURE;
00934     }
00935     const EntityHandle vol = vols.front();
00936 
00937     CartVect position( 10, 0, 0 );
00938     double result = 0.0;
00939     rval          = gqt->closest_to_location( vol, position.array(), result );
00940     if( result != 9.0 ) return MB_FAILURE;
00941 
00942     EntityHandle surface = 1;
00943     rval                 = gqt->closest_to_location( vol, position.array(), result, &surface );
00944 
00945     if( result != 9.0 ) return MB_FAILURE;
00946     if( surface == 1 ) return MB_FAILURE;
00947 
00948     return MB_SUCCESS;
00949 }
00950 
00951 ErrorCode overlap_test_point_in_volume( GeomQueryTool* gqt )
00952 {
00953     const char* const NAME_ARR[] = { "Boundary", "Outside", "Inside" };
00954     const char* const* names     = NAME_ARR + 1;
00955     const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
00956     const struct PointInVol tests[] = { { { 0.5, 0.0, 0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
00957                                         { { 0.5, 0.0, -0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
00958                                         { { 0.5, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00959                                         { { -0.5, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00960                                         { { 0.5, 0.5, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00961                                         { { 0.5, -0.5, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00962                                         { { 1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00963                                         { { -1.1, 1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00964                                         { { -1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00965                                         { { 1.1, -1.1, 1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00966                                         { { 1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00967                                         { { -1.1, 1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00968                                         { { -1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00969                                         { { 1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00970                                         // Add some directions to test special cases of edge/node intersection
00971                                         { { 0.5, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
00972                                         { { 0.5, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
00973                                         { { 0.0, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
00974                                         { { 0.5, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
00975                                         { { 0.5, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } },
00976                                         // Test some points in the overlap
00977                                         { { 0.005, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
00978                                         { { 0.005, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
00979                                         { { 0.005, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
00980                                         { { 0.005, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
00981                                         { { 0.005, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } } };
00982 
00983     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00984 
00985     ErrorCode rval;
00986     Interface* moab = gqt->moab_instance();
00987 
00988     Tag dim_tag = gqt->gttool()->get_geom_tag();
00989     Range vols;
00990     const int three = 3;
00991     const void* ptr = &three;
00992     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00993     if( vols.size() != 2 )
00994     {
00995         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00996         return MB_FAILURE;
00997     }
00998     const EntityHandle vol = vols.front();
00999     for( int i = 0; i < num_test; ++i )
01000     {
01001         int result;
01002         rval = gqt->point_in_volume( vol, tests[i].coords, result, tests[i].dir );CHKERR;
01003         if( result != tests[i].result )
01004         {
01005             std::cerr << "ERROR testing point_in_volume[" << i << "]:" << std::endl
01006                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
01007                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
01008             return MB_FAILURE;
01009         }
01010 
01011         // point_in_volume_slow doesn't to boundary.
01012         if( tests[i].result == BOUNDARY ) continue;
01013 
01014         rval = gqt->point_in_volume_slow( vol, tests[i].coords, result );CHKERR;
01015 
01016         if( result != tests[i].result )
01017         {
01018             std::cerr << "ERROR testing point_in_volume_slow[" << i << "]:" << std::endl
01019                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
01020                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
01021             return MB_FAILURE;
01022         }
01023     }
01024 
01025     return MB_SUCCESS;
01026 }
01027 
01028 ErrorCode overlap_test_tracking( GeomQueryTool* gqt )
01029 {
01030     /* Track a particle from left (-x) to right (+x) through an overlap.
01031                   ____________________
01032                   |                   |
01033          -x <---  | region1 | overlap | region0 |   ---> +x
01034                             |                   |
01035                             ---------------------
01036          -x <--- -1.0      0.0      0.01        1.0 ---> +x
01037        surf_id:   10        4         8         2                     */
01038 
01039     // get the surfaces and volumes
01040     Tag dim_tag = gqt->gttool()->get_geom_tag();
01041     Range surfs, explicit_vols;
01042     const int two = 2, three = 3;
01043     const void* ptrs[] = { &two, &three };
01044     ErrorCode rval;
01045     Interface* moab = gqt->moab_instance();
01046     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
01047     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, explicit_vols );CHKERR;
01048 
01049     if( explicit_vols.size() != 2 )
01050     {
01051         std::cerr << "ERROR: Expected 2 explicit volumes in input, found " << explicit_vols.size() << std::endl;
01052         return MB_FAILURE;
01053     }
01054     const unsigned num_surf = 12;
01055     if( surfs.size() != num_surf )
01056     {
01057         std::cerr << "ERROR: Expected " << num_surf << " surfaces in input, found " << surfs.size() << std::endl;
01058         return MB_FAILURE;
01059     }
01060     const EntityHandle explicit_vol = explicit_vols.front();
01061 
01062     // start particle
01063     double point[]     = { -0.9, 0, 0 };
01064     const double dir[] = { 1, 0, 0 };
01065     EntityHandle vol   = explicit_vol;
01066     int result;
01067     const int INSIDE = 1;  // OUTSIDE = 0, BOUNDARY = -1;
01068     rval             = gqt->point_in_volume( explicit_vol, point, result, dir );CHKERR;
01069     if( result != INSIDE )
01070     {
01071         std::cerr << "ERROR: particle not inside explicit volume" << std::endl;
01072         return MB_FAILURE;
01073     }
01074 
01075     // get next surface
01076     double dist;
01077     EntityHandle next_surf;
01078     GeomQueryTool::RayHistory history;
01079     rval = gqt->ray_fire( vol, point, dir, next_surf, dist, &history );CHKERR;
01080     if( next_surf != surfs[7] || fabs( dist - 0.91 ) > 1e-6 )
01081     {
01082         std::cerr << "ERROR: failed on advance 1" << std::endl;
01083         return MB_FAILURE;
01084     }
01085     for( unsigned i = 0; i < 3; i++ )
01086         point[i] += dist * dir[i];
01087 
01088     // get the next volume (implicit complement)
01089     EntityHandle next_vol;
01090     rval = gqt->gttool()->next_vol( next_surf, vol, next_vol );CHKERR;
01091 
01092     // get the next surface (behind numerical location)
01093     vol  = next_vol;
01094     rval = gqt->ray_fire( vol, point, dir, next_surf, dist, &history );CHKERR;
01095     if( next_surf != surfs[3] || fabs( dist - 0.0 ) > 1e-6 )
01096     {
01097         std::cerr << "ERROR: failed on advance 2" << std::endl;
01098         return MB_FAILURE;
01099     }
01100     for( unsigned i = 0; i < 3; i++ )
01101         point[i] += dist * dir[i];
01102 
01103     // get the next volume (the explicit volume)
01104     rval = gqt->gttool()->next_vol( next_surf, vol, next_vol );CHKERR;
01105 
01106     // get the next surface
01107     vol  = next_vol;
01108     rval = gqt->ray_fire( vol, point, dir, next_surf, dist, &history );CHKERR;
01109     if( next_surf != surfs[1] || fabs( dist - 0.99 ) > 1e-6 )
01110     {
01111         std::cerr << "ERROR: failed on advance 3" << std::endl;
01112         return MB_FAILURE;
01113     }
01114     for( unsigned i = 0; i < 3; i++ )
01115         point[i] += dist * dir[i];
01116 
01117     return MB_SUCCESS;
01118 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines