MOAB: Mesh Oriented datABase  (version 5.4.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 ) )            \
00229         {                                       \
00230             ++errors;                           \
00231         }                                       \
00232     } while( false )
00233 
00234 int main( int argc, char* argv[] )
00235 {
00236     ErrorCode rval;
00237     const char* filename = "test_geom.h5m";
00238 
00239 #ifdef MOAB_HAVE_MPI
00240     int fail = MPI_Init( &argc, &argv );
00241     if( fail ) return fail;
00242 #endif
00243 
00244     rval = write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file: " << filename );
00245 
00246     Interface* MBI = new Core();
00247 
00248     int errors = 0;
00249     rval       = MBI->load_file( filename );
00250     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file" );
00251 
00252     GeomTopoTool* gtt  = new GeomTopoTool( MBI );
00253     GeomQueryTool* gqt = new GeomQueryTool( gtt );
00254 
00255     errors += run_regular_tests( gqt );
00256 
00257     // clear out moab instance
00258     rval = MBI->delete_mesh();MB_CHK_SET_ERR( rval, "Failed to delete mesh" );
00259 
00260     delete gtt;
00261     delete gqt;
00262 
00263     // Now load a different geometry: two cubes that slightly overlap
00264     rval = overlap_write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file: " << filename );
00265 
00266     rval = MBI->load_file( filename );
00267     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file with overlaps" );
00268 
00269     gtt = new GeomTopoTool( MBI );
00270     gqt = new GeomQueryTool( gtt );
00271 
00272     errors += run_overlap_tests( gqt );
00273 
00274     // clear moab instance
00275     rval = MBI->delete_mesh();MB_CHK_SET_ERR( rval, "Failed to delete mesh" );
00276 
00277     delete gtt;
00278     delete gqt;
00279 
00280     // Re-run all tests with the alternate constructor
00281     std::cout << "-------------------------------------" << std::endl;
00282     std::cout << "Re-running tests with MBI constructor" << std::endl;
00283     std::cout << "-------------------------------------" << std::endl;
00284 
00285     rval = write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file" );
00286 
00287     rval = MBI->load_file( filename );
00288     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file" );
00289 
00290     gqt = new GeomQueryTool( MBI );
00291 
00292     errors += run_regular_tests( gqt );
00293 
00294     // clear moab and dagmc instance
00295     rval = MBI->delete_mesh();MB_CHK_SET_ERR( rval, "Failed to delete mesh" );
00296 
00297     delete gqt;
00298 
00299     // Now load a different geometry: two cubes that slightly overlap
00300     rval = overlap_write_geometry( filename );MB_CHK_SET_ERR( rval, "Failed to create input file: " );
00301 
00302     rval = MBI->load_file( filename );
00303     remove( filename );MB_CHK_SET_ERR( rval, "Failed to load file with overlaps." );
00304 
00305     gqt = new GeomQueryTool( MBI );
00306 
00307     errors += run_overlap_tests( gqt );
00308 
00309     // final cleanup
00310     delete gqt;
00311     delete MBI;
00312 
00313 #ifdef MOAB_HAVE_MPI
00314     fail = MPI_Finalize();
00315     if( fail ) return fail;
00316 #endif
00317 
00318     return errors;
00319 }
00320 
00321 int run_regular_tests( GeomQueryTool* gqt )
00322 {
00323     int errors = 0;
00324     ErrorCode rval;
00325 
00326     rval = gqt->initialize();MB_CHK_SET_ERR_CONT( rval, "Failed to initialize the GeomQueryTool" );
00327 
00328     RUN_TEST( test_ray_fire );
00329     RUN_TEST( test_closest_to_location );
00330     RUN_TEST( test_point_in_volume );
00331     RUN_TEST( test_measure_volume );
00332     RUN_TEST( test_measure_area );
00333     RUN_TEST( test_surface_sense );
00334 
00335     // change settings to use overlap-tolerant mode (arbitrary thickness)
00336     double overlap_thickness = 0.1;
00337     gqt->set_overlap_thickness( overlap_thickness );
00338     RUN_TEST( test_ray_fire );
00339     RUN_TEST( test_point_in_volume );
00340 
00341     return errors;
00342 }
00343 
00344 int run_overlap_tests( GeomQueryTool* gqt )
00345 {
00346     int errors = 0;
00347     ErrorCode rval;
00348 
00349     rval = gqt->initialize();MB_CHK_SET_ERR_CONT( rval, "Failed to initialize the GeomQueryTool" );
00350 
00351     // change settings to use overlap-tolerant mode (with a large enough thickness)
00352     double overlap_thickness = 3.0;
00353     gqt->set_overlap_thickness( overlap_thickness );
00354     RUN_TEST( overlap_test_ray_fire );
00355     RUN_TEST( overlap_test_point_in_volume );
00356     RUN_TEST( overlap_test_measure_volume );
00357     RUN_TEST( overlap_test_measure_area );
00358     RUN_TEST( overlap_test_surface_sense );
00359     RUN_TEST( overlap_test_tracking );
00360 
00361     return errors;
00362 }
00363 
00364 ErrorCode test_surface_sense( GeomQueryTool* gqt )
00365 {
00366     ErrorCode rval;
00367     Interface* moab = gqt->moab_instance();
00368 
00369     Tag dim_tag = gqt->gttool()->get_geom_tag();
00370     Range surfs, vols;
00371     const int two = 2, three = 3;
00372     const void* ptrs[] = { &two, &three };
00373     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00374     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00375 
00376     if( vols.size() != 2 )
00377     {
00378         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00379         return MB_FAILURE;
00380     }
00381     if( surfs.size() != 6 )
00382     {
00383         std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
00384         return MB_FAILURE;
00385     }
00386 
00387     for( Range::iterator i = surfs.begin(); i != surfs.end(); ++i )
00388     {
00389         int sense = 0;
00390         rval      = gqt->gttool()->get_sense( *i, vols.front(), sense );
00391         if( MB_SUCCESS != rval || sense != 1 )
00392         {
00393             std::cerr << "ERROR: Expected 1 for surface sense, got " << sense << std::endl;
00394             return MB_FAILURE;
00395         }
00396     }
00397 
00398     return MB_SUCCESS;
00399 }
00400 
00401 ErrorCode overlap_test_surface_sense( GeomQueryTool* gqt )
00402 {
00403     ErrorCode rval;
00404     Interface* moab = gqt->moab_instance();
00405 
00406     Tag dim_tag = gqt->gttool()->get_geom_tag();
00407     Range surfs, vols;
00408     const int two = 2, three = 3;
00409     const void* ptrs[] = { &two, &three };
00410     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00411     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00412 
00413     if( vols.size() != 2 )
00414     {
00415         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00416         return MB_FAILURE;
00417     }
00418     if( surfs.size() != 12 )
00419     {
00420         std::cerr << "ERROR: Expected 12 surfaces in input, found " << surfs.size() << std::endl;
00421         return MB_FAILURE;
00422     }
00423 
00424     for( Range::iterator i = surfs.begin(); i != surfs.end(); ++i )
00425     {
00426         int sense = 0;
00427         rval      = gqt->gttool()->get_sense( *i, vols.front(), sense );
00428         if( MB_SUCCESS != rval || sense != 1 )
00429         {
00430             std::cerr << "ERROR: Expected 1 for surface sense, got " << sense << std::endl;
00431             return MB_FAILURE;
00432         }
00433     }
00434 
00435     return MB_SUCCESS;
00436 }
00437 ErrorCode test_measure_volume( GeomQueryTool* gqt )
00438 {
00439     ErrorCode rval;
00440     Interface* moab = gqt->moab_instance();
00441 
00442     Tag dim_tag = gqt->gttool()->get_geom_tag();
00443     Range vols;
00444     const int three = 3;
00445     const void* ptr = &three;
00446     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00447 
00448     if( vols.size() != 2 )
00449     {
00450         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00451         return MB_FAILURE;
00452     }
00453 
00454     // input file is 2x2x2 cube with convexity in +Z face that touches the origin.
00455     // expected volume is 8 (2x2x2) less the volume of the pyrimid concavity
00456     double result;
00457     const double vol = 2 * 2 * 2 - 1 * 4. / 3;
00458 
00459     rval = gqt->measure_volume( vols.front(), result );CHKERR;
00460     if( fabs( result - vol ) > 10 * std::numeric_limits< double >::epsilon() )
00461     {
00462         std::cerr << "ERROR: Expected " << vol << " as measure of volume, got " << result << std::endl;
00463         return MB_FAILURE;
00464     }
00465 
00466     return MB_SUCCESS;
00467 }
00468 ErrorCode overlap_test_measure_volume( GeomQueryTool* gqt )
00469 {
00470     ErrorCode rval;
00471     Interface* moab = gqt->moab_instance();
00472 
00473     Tag dim_tag = gqt->gttool()->get_geom_tag();
00474     Range vols;
00475     const int three = 3;
00476     const void* ptr = &three;
00477     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00478 
00479     if( vols.size() != 2 )
00480     {
00481         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00482         return MB_FAILURE;
00483     }
00484 
00485     // The volume has two regions that overlap.
00486     double result;
00487     const double vol = ( 1 + 1.01 ) * 2 * 2;
00488 
00489     rval = gqt->measure_volume( vols.front(), result );CHKERR;
00490     if( fabs( result - vol ) > 2 * std::numeric_limits< double >::epsilon() )
00491     {
00492         std::cerr << "ERROR: Expected " << vol << " as measure of volume, got " << result << std::endl;
00493         return MB_FAILURE;
00494     }
00495 
00496     return MB_SUCCESS;
00497 }
00498 
00499 ErrorCode test_measure_area( GeomQueryTool* gqt )
00500 {
00501     ErrorCode rval;
00502     Interface* moab = gqt->moab_instance();
00503 
00504     Tag dim_tag = gqt->gttool()->get_geom_tag();
00505     Range surfs;
00506     const int two   = 2;
00507     const void* ptr = &two;
00508     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, surfs );CHKERR;
00509 
00510     if( surfs.size() != 6 )
00511     {
00512         std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
00513         return MB_FAILURE;
00514     }
00515 
00516     int ids[6];
00517     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00518 
00519     // expect area of 4 for all faces except face 6.
00520     // face 6 should have area == 4*sqrt(2)
00521     Range::iterator iter = surfs.begin();
00522     for( unsigned i = 0; i < 6; ++i, ++iter )
00523     {
00524         double expected = 4.0;
00525         if( ids[i] == 6 ) expected *= ROOT2;
00526 
00527         double result;
00528 
00529         rval = gqt->measure_area( *iter, result );CHKERR;
00530         if( fabs( result - expected ) > std::numeric_limits< double >::epsilon() )
00531         {
00532             std::cerr << "ERROR: Expected area of surface " << ids[i] << " to be " << expected << ".  Got " << result
00533                       << std::endl;
00534             return MB_FAILURE;
00535         }
00536     }
00537 
00538     return MB_SUCCESS;
00539 }
00540 
00541 ErrorCode overlap_test_measure_area( GeomQueryTool* gqt )
00542 {
00543     ErrorCode rval;
00544     Interface* moab = gqt->moab_instance();
00545 
00546     Tag dim_tag = gqt->gttool()->get_geom_tag();
00547     Range surfs;
00548     const int two   = 2;
00549     const void* ptr = &two;
00550     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, surfs );CHKERR;
00551 
00552     const unsigned num_surfs = 12;
00553     if( surfs.size() != num_surfs )
00554     {
00555         std::cerr << "ERROR: Expected " << num_surfs << " surfaces in input, found " << surfs.size() << std::endl;
00556         return MB_FAILURE;
00557     }
00558 
00559     int ids[num_surfs];
00560     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00561 
00562     const double x_area   = 2 * 2;
00563     const double yz_area0 = 2 * 1;
00564     const double yz_area1 = 2 * 1.01;
00565     Range::iterator iter  = surfs.begin();
00566     for( unsigned i = 0; i < num_surfs; ++i, ++iter )
00567     {
00568         double expected, result;
00569         if( 0 == i || 2 == i || 4 == i || 5 == i )
00570             expected = yz_area0;
00571         else if( 1 == i || 3 == i || 7 == i || 9 == i )
00572             expected = x_area;
00573         else if( 6 == i || 8 == i || 10 == i || 11 == i )
00574             expected = yz_area1;
00575 
00576         rval = gqt->measure_area( *iter, result );CHKERR;
00577         if( fabs( result - expected ) > std::numeric_limits< double >::epsilon() )
00578         {
00579             std::cerr << "ERROR: Expected area of surface " << ids[i] << " to be " << expected << ".  Got " << result
00580                       << std::endl;
00581             return MB_FAILURE;
00582         }
00583     }
00584 
00585     return MB_SUCCESS;
00586 }
00587 
00588 struct ray_fire
00589 {
00590     int prev_surf;
00591     double origin[3], direction[3];
00592     int hit_surf;
00593     double distance;
00594 };
00595 
00596 ErrorCode test_ray_fire( GeomQueryTool* gqt )
00597 {
00598     // Glancing ray-triangle intersections are not valid exit intersections.
00599     // Piercing ray-triangle intersections are valid exit intersections.
00600     // "0" destination surface implies that it is ambiguous.
00601     const struct ray_fire tests[] = { /* src_srf origin               direction                 dest dist */
00602                                       // piercing edge
00603                                       { 1, { 0.0, 0.0, -1. }, { -1.0 / ROOT2, 0.0, 1.0 / ROOT2 }, 4, ROOT2 },
00604                                       // piercing edge
00605                                       { 1, { 0.0, 0.0, -1. }, { 1.0 / ROOT2, 0.0, 1.0 / ROOT2 }, 2, ROOT2 },
00606                                       // piercing edge
00607                                       { 1, { 0.0, 0.0, -1. }, { 0.0, 1.0 / ROOT2, 1.0 / ROOT2 }, 3, ROOT2 },
00608                                       // piercing edge
00609                                       { 1, { 0.5, 0.5, -1. }, { 0.0, 0.0, 1.0 }, 6, 1.5 },
00610                                       // interior
00611                                       { 2, { 1.0, 0.0, 0.5 }, { -1.0, 0.0, 0.0 }, 6, 0.5 },
00612                                       // glancing node then piercing edge
00613                                       { 2, { 1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 2.0 },
00614                                       // piercing node
00615                                       { 1, { 0.0, 0.0, -1. }, { 0.0, 0.0, 1.0 }, 6, 1.0 },
00616                                       // glancing edge then interior
00617                                       { 2, { 1.0, 0.0, 0.5 }, { -1.0 / ROOT2, 1.0 / ROOT2, 0.0 }, 3, ROOT2 } };
00618 
00619     ErrorCode rval;
00620     Interface* moab = gqt->moab_instance();
00621 
00622     Tag dim_tag = gqt->gttool()->get_geom_tag();
00623     Range surfs, vols;
00624     const int two = 2, three = 3;
00625     const void* ptrs[] = { &two, &three };
00626     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00627     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00628 
00629     if( vols.size() != 2 )
00630     {
00631         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00632         return MB_FAILURE;
00633     }
00634     if( surfs.size() != 6 )
00635     {
00636         std::cerr << "ERROR: Expected 6 surfaces in input, found " << surfs.size() << std::endl;
00637         return MB_FAILURE;
00638     }
00639 
00640     int ids[6];
00641     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00642     EntityHandle surf[6];
00643     std::copy( surfs.begin(), surfs.end(), surf );
00644 
00645     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00646     for( int i = 0; i < num_test; ++i )
00647     {
00648         int* ptr = std::find( ids, ids + 6, tests[i].prev_surf );
00649         int idx  = ptr - ids;
00650         if( idx >= 6 )
00651         {
00652             std::cerr << "Surface " << tests[i].prev_surf << " not found." << std::endl;
00653             return MB_FAILURE;
00654         }
00655         // const EntityHandle src_surf = surf[idx];
00656 
00657         ptr = std::find( ids, ids + 6, tests[i].hit_surf );
00658         idx = ptr - ids;
00659         if( idx >= 6 )
00660         {
00661             std::cerr << "Surface " << tests[i].hit_surf << " not found." << std::endl;
00662             return MB_FAILURE;
00663         }
00664         const EntityHandle hit_surf = surf[idx];
00665 
00666         double dist;
00667         EntityHandle result;
00668         GeomQueryTool::RayHistory history;
00669         rval = gqt->ray_fire( vols.front(), tests[i].origin, tests[i].direction, result, dist, &history );
00670 
00671         if( result != hit_surf || fabs( dist - tests[i].distance ) > 1e-6 )
00672         {
00673             EntityHandle* p = std::find( surf, surf + 6, result );
00674             idx             = p - surf;
00675             int id          = idx > 5 ? 0 : ids[idx];
00676 
00677             std::cerr << "Rayfire test failed for " << std::endl
00678                       << "\t ray from (" << tests[i].origin[0] << ", " << tests[i].origin[1] << ", "
00679                       << tests[i].origin[2] << ") going [" << tests[i].direction[0] << ", " << tests[i].direction[1]
00680                       << ", " << tests[i].direction[2] << "]" << std::endl
00681                       << "\t Beginning on surface " << tests[i].prev_surf << std::endl
00682                       << "\t Expected to hit surface " << tests[i].hit_surf << " after " << tests[i].distance
00683                       << " units." << std::endl
00684                       << "\t Actually hit surface " << id << " after " << dist << " units." << std::endl;
00685             return MB_FAILURE;
00686         }
00687 
00688         CartVect loc = CartVect( tests[i].origin ) + ( dist * CartVect( tests[i].direction ) );
00689 
00690         std::vector< std::pair< int, GeomQueryTool::RayHistory* > > boundary_tests;
00691         boundary_tests.push_back( std::make_pair( 1, &history ) );
00692         boundary_tests.push_back( std::make_pair( 0, &history ) );
00693         boundary_tests.push_back( std::make_pair( 1, (GeomQueryTool::RayHistory*)NULL ) );
00694         boundary_tests.push_back( std::make_pair( 0, (GeomQueryTool::RayHistory*)NULL ) );
00695 
00696         for( unsigned int bt = 0; bt < boundary_tests.size(); ++bt )
00697         {
00698 
00699             int expected                 = boundary_tests[bt].first;
00700             GeomQueryTool::RayHistory* h = boundary_tests[bt].second;
00701 
00702             // pick the direction based on expected result of test. Either reuse the ray_fire
00703             // vector, or reverse it to check for a vector that enters the cell
00704             CartVect uvw( tests[i].direction );
00705             if( expected == 1 ) uvw = -uvw;
00706 
00707             int boundary_result = -1;
00708 
00709             rval = gqt->test_volume_boundary( vols.front(), result, loc.array(), uvw.array(), boundary_result, h );
00710 
00711             if( boundary_result != expected )
00712             {
00713                 std::cerr << "DagMC::test_volume_boundary failed (" << ( ( expected == 0 ) ? "+" : "-" ) << " dir,"
00714                           << ( ( h ) ? "+" : "-" ) << " history, i=" << i << ")" << std::endl;
00715                 return MB_FAILURE;
00716             }
00717         }
00718     }
00719 
00720     return MB_SUCCESS;
00721 }
00722 
00723 ErrorCode overlap_test_ray_fire( GeomQueryTool* gqt )
00724 {
00725     // Glancing ray-triangle intersections are not valid exit intersections.
00726     // Piercing ray-triangle intersections are valid exit intersections.
00727     // "0" destination surface implies that it is ambiguous.
00728     const struct ray_fire tests[] = { /*          ____________________
00729                                                   |                   |
00730                                          -x <---  | region1 | overlap | region0 |   ---> +x
00731                                                             |                   |
00732                                                             ---------------------
00733                                          -x <--- -1.0      0.0      0.01        1.0 ---> +x
00734                                        surf_id:   10        4         8         2
00735 
00736                                        The zero-distance advance would occur in the implicit volume between
00737                                        these two regions---not in this volume.
00738 
00739                                        src_srf origin               direction                 dest dist */
00740                                       // numerical location on surface 1 of region0
00741                                       { 4, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 2, 0.0 },
00742                                       { 2, { 1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 1.0 },
00743                                       // numerical location inside region0
00744                                       { 4, { 0.5, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 2, 0.5 },
00745                                       { 2, { 0.5, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.5 },
00746                                       // numerical location on surface 7 of region1
00747                                       { 4, { 0.01, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.0 },
00748                                       { 2, { 0.01, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.01 },
00749                                       // numerical location inside overlap
00750                                       { 10, { 0.005, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.005 },
00751                                       { 2, { 0.005, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.005 },
00752                                       // numerical location on surface 3 of region0
00753                                       { 10, { 0.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.01 },
00754                                       { 8, { 0.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 4, 0.0 },
00755                                       // numerical location inside region1
00756                                       { 10, { -0.5, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 0.51 },
00757                                       { 8, { -0.5, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 10, 0.5 },
00758                                       // numerical location on surface 9 of region1
00759                                       { 10, { -1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, 8, 1.01 },
00760                                       { 8, { -1.0, 0.0, 0.0 }, { -1.0, 0.0, 0.0 }, 10, 0.0 } };
00761 
00762     ErrorCode rval;
00763     Interface* moab = gqt->moab_instance();
00764 
00765     Tag dim_tag = gqt->gttool()->get_geom_tag();
00766     Range surfs, vols;
00767     const int two = 2, three = 3;
00768     const void* ptrs[] = { &two, &three };
00769     rval               = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
00770     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, vols );CHKERR;
00771 
00772     if( vols.size() != 2 )
00773     {
00774         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00775         return MB_FAILURE;
00776     }
00777     const unsigned num_surf = 12;
00778     if( surfs.size() != num_surf )
00779     {
00780         std::cerr << "ERROR: Expected " << num_surf << " surfaces in input, found " << surfs.size() << std::endl;
00781         return MB_FAILURE;
00782     }
00783 
00784     int ids[num_surf];
00785     rval = moab->tag_get_data( gqt->gttool()->get_gid_tag(), surfs, ids );CHKERR;
00786     EntityHandle surf[num_surf];
00787     std::copy( surfs.begin(), surfs.end(), surf );
00788 
00789     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00790     for( int i = 0; i < num_test; ++i )
00791     {
00792         int* ptr     = std::find( ids, ids + num_surf, tests[i].prev_surf );
00793         unsigned idx = ptr - ids;
00794         if( idx >= num_surf )
00795         {
00796             std::cerr << "Surface " << tests[i].prev_surf << " not found." << std::endl;
00797             return MB_FAILURE;
00798         }
00799         // const EntityHandle src_surf = surf[idx];
00800 
00801         ptr = std::find( ids, ids + num_surf, tests[i].hit_surf );
00802         idx = ptr - ids;
00803         if( idx >= num_surf )
00804         {
00805             std::cerr << "Surface " << tests[i].hit_surf << " not found." << std::endl;
00806             return MB_FAILURE;
00807         }
00808         const EntityHandle hit_surf = surf[idx];
00809 
00810         double dist;
00811         EntityHandle result;
00812         rval = gqt->ray_fire( vols.front(), tests[i].origin, tests[i].direction, result, dist );
00813 
00814         if( result != hit_surf || fabs( dist - tests[i].distance ) > 1e-6 )
00815         {
00816             EntityHandle* p = std::find( surf, surf + 6, result );
00817             idx             = p - surf;
00818             int id          = idx > num_surf - 1 ? 0 : ids[idx];
00819 
00820             std::cerr << "Rayfire test failed for " << std::endl
00821                       << "\t ray from (" << tests[i].origin[0] << ", " << tests[i].origin[1] << ", "
00822                       << tests[i].origin[2] << ") going [" << tests[i].direction[0] << ", " << tests[i].direction[1]
00823                       << ", " << tests[i].direction[2] << "]" << std::endl
00824                       << "\t Beginning on surface " << tests[i].prev_surf << std::endl
00825                       << "\t Expected to hit surface " << tests[i].hit_surf << " after " << tests[i].distance
00826                       << " units." << std::endl
00827                       << "\t Actually hit surface " << id << " after " << dist << " units." << std::endl;
00828             return MB_FAILURE;
00829         }
00830     }
00831 
00832     return MB_SUCCESS;
00833 }
00834 
00835 struct PointInVol
00836 {
00837     double coords[3];
00838     int result;
00839     double dir[3];
00840 };
00841 
00842 ErrorCode test_point_in_volume( GeomQueryTool* gqt )
00843 {
00844     const char* const NAME_ARR[] = { "Boundary", "Outside", "Inside" };
00845     const char* const* names     = NAME_ARR + 1;
00846     const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
00847     const struct PointInVol tests[] = { { { 0.0, 0.0, 0.5 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00848                                         { { 0.0, 0.0, -0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
00849                                         { { 0.7, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00850                                         { { -0.7, 0.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00851                                         { { 0.0, -0.7, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00852                                         { { 0.0, -0.7, 0.0 }, INSIDE, { 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                                         { { 1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00861                                         // Add some directions to test special cases of edge/node intersection
00862                                         { { 0.5, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
00863                                         { { 0.5, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
00864                                         { { 0.0, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
00865                                         { { 0.5, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
00866                                         { { 0.5, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } } };
00867 
00868     //    { { 1.0, 0.0, 0.0 }, BOUNDARY}, MCNP doesn't return on boundary
00869     //{ {-1.0, 0.0, 0.0 }, BOUNDARY},
00870     //{ { 0.0, 1.0, 0.0 }, BOUNDARY},
00871     //{ { 0.0,-1.0, 0.0 }, BOUNDARY},
00872     //{ { 0.0, 0.0, 0.0 }, BOUNDARY},
00873     //{ { 0.0, 0.0,-1.0 }, BOUNDARY} };
00874     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00875 
00876     ErrorCode rval;
00877     Interface* moab = gqt->moab_instance();
00878 
00879     Tag dim_tag = gqt->gttool()->get_geom_tag();
00880 
00881     Range vols;
00882     const int three = 3;
00883     const void* ptr = &three;
00884     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00885     if( vols.size() != 2 )
00886     {
00887         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00888         return MB_FAILURE;
00889     }
00890     const EntityHandle vol = vols.front();
00891 
00892     for( int i = 0; i < num_test; ++i )
00893     {
00894         int result;
00895         rval = gqt->point_in_volume( vol, tests[i].coords, result, tests[i].dir );CHKERR;
00896         if( result != tests[i].result )
00897         {
00898             std::cerr << "ERROR testing point_in_volume[" << i << "]:" << std::endl
00899                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
00900                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
00901             return MB_FAILURE;
00902         }
00903 
00904         // point_in_volume_slow doesn't to boundary.
00905         if( tests[i].result == BOUNDARY ) continue;
00906 
00907         rval = gqt->point_in_volume_slow( vol, tests[i].coords, result );CHKERR;
00908 
00909         if( result != tests[i].result )
00910         {
00911             std::cerr << "ERROR testing point_in_volume_slow[" << i << "]:" << std::endl
00912                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
00913                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
00914             return MB_FAILURE;
00915         }
00916     }
00917 
00918     return MB_SUCCESS;
00919 }
00920 
00921 ErrorCode test_closest_to_location( GeomQueryTool* gqt )
00922 {
00923     ErrorCode rval;
00924     Interface* moab = gqt->moab_instance();
00925 
00926     Range vols;
00927     const int three = 3;
00928     const void* ptr = &three;
00929     Tag dim_tag     = gqt->gttool()->get_geom_tag();
00930     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00931     if( vols.size() != 2 )
00932     {
00933         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00934         return MB_FAILURE;
00935     }
00936     const EntityHandle vol = vols.front();
00937 
00938     CartVect position( 10, 0, 0 );
00939     double result = 0.0;
00940     rval          = gqt->closest_to_location( vol, position.array(), result );
00941     if( result != 9.0 ) return MB_FAILURE;
00942 
00943     EntityHandle surface = 1;
00944     rval                 = gqt->closest_to_location( vol, position.array(), result, &surface );
00945 
00946     if( result != 9.0 ) return MB_FAILURE;
00947     if( surface == 1 ) return MB_FAILURE;
00948 
00949     return MB_SUCCESS;
00950 }
00951 
00952 ErrorCode overlap_test_point_in_volume( GeomQueryTool* gqt )
00953 {
00954     const char* const NAME_ARR[] = { "Boundary", "Outside", "Inside" };
00955     const char* const* names     = NAME_ARR + 1;
00956     const int INSIDE = 1, OUTSIDE = 0, BOUNDARY = -1;
00957     const struct PointInVol tests[] = { { { 0.5, 0.0, 0.5 }, INSIDE, { 0.0, 0.0, 0.0 } },
00958                                         { { 0.5, 0.0, -0.5 }, 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.0, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00961                                         { { 0.5, 0.5, 0.0 }, INSIDE, { 0.0, 0.0, 0.0 } },
00962                                         { { 0.5, -0.5, 0.0 }, INSIDE, { 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                                         { { 1.1, -1.1, -1.1 }, OUTSIDE, { 0.0, 0.0, 0.0 } },
00971                                         // Add some directions to test special cases of edge/node intersection
00972                                         { { 0.5, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
00973                                         { { 0.5, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
00974                                         { { 0.0, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
00975                                         { { 0.5, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
00976                                         { { 0.5, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } },
00977                                         // Test some points in the overlap
00978                                         { { 0.005, 0.0, 0.0 }, INSIDE, { -1.0, 0.0, 0.0 } },
00979                                         { { 0.005, 0.0, 0.0 }, INSIDE, { 1.0, 0.0, 0.0 } },
00980                                         { { 0.005, 0.0, 2.0 }, OUTSIDE, { 0.0, 0.0, -1.0 } },
00981                                         { { 0.005, 0.0, -0.5 }, INSIDE, { -1.0, 0.0, 0.0 } },
00982                                         { { 0.005, -0.5, -2.0 }, OUTSIDE, { 0.0, 0.0, 1.0 } } };
00983 
00984     const int num_test = sizeof( tests ) / sizeof( tests[0] );
00985 
00986     ErrorCode rval;
00987     Interface* moab = gqt->moab_instance();
00988 
00989     Tag dim_tag = gqt->gttool()->get_geom_tag();
00990     Range vols;
00991     const int three = 3;
00992     const void* ptr = &three;
00993     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, &ptr, 1, vols );CHKERR;
00994     if( vols.size() != 2 )
00995     {
00996         std::cerr << "ERROR: Expected 2 volumes in input, found " << vols.size() << std::endl;
00997         return MB_FAILURE;
00998     }
00999     const EntityHandle vol = vols.front();
01000     for( int i = 0; i < num_test; ++i )
01001     {
01002         int result;
01003         rval = gqt->point_in_volume( vol, tests[i].coords, result, tests[i].dir );CHKERR;
01004         if( result != tests[i].result )
01005         {
01006             std::cerr << "ERROR testing point_in_volume[" << i << "]:" << std::endl
01007                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
01008                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
01009             return MB_FAILURE;
01010         }
01011 
01012         // point_in_volume_slow doesn't to boundary.
01013         if( tests[i].result == BOUNDARY ) continue;
01014 
01015         rval = gqt->point_in_volume_slow( vol, tests[i].coords, result );CHKERR;
01016 
01017         if( result != tests[i].result )
01018         {
01019             std::cerr << "ERROR testing point_in_volume_slow[" << i << "]:" << std::endl
01020                       << "\tExpected " << names[tests[i].result] << " for (" << tests[i].coords[0] << ", "
01021                       << tests[i].coords[1] << ", " << tests[i].coords[2] << ").  Got " << names[result] << std::endl;
01022             return MB_FAILURE;
01023         }
01024     }
01025 
01026     return MB_SUCCESS;
01027 }
01028 
01029 ErrorCode overlap_test_tracking( GeomQueryTool* gqt )
01030 {
01031     /* Track a particle from left (-x) to right (+x) through an overlap.
01032                   ____________________
01033                   |                   |
01034          -x <---  | region1 | overlap | region0 |   ---> +x
01035                             |                   |
01036                             ---------------------
01037          -x <--- -1.0      0.0      0.01        1.0 ---> +x
01038        surf_id:   10        4         8         2                     */
01039 
01040     // get the surfaces and volumes
01041     Tag dim_tag = gqt->gttool()->get_geom_tag();
01042     Range surfs, explicit_vols;
01043     const int two = 2, three = 3;
01044     const void* ptrs[] = { &two, &three };
01045     ErrorCode rval;
01046     Interface* moab = gqt->moab_instance();
01047     rval            = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs, 1, surfs );CHKERR;
01048     rval = moab->get_entities_by_type_and_tag( 0, MBENTITYSET, &dim_tag, ptrs + 1, 1, explicit_vols );CHKERR;
01049 
01050     if( explicit_vols.size() != 2 )
01051     {
01052         std::cerr << "ERROR: Expected 2 explicit volumes in input, found " << explicit_vols.size() << std::endl;
01053         return MB_FAILURE;
01054     }
01055     const unsigned num_surf = 12;
01056     if( surfs.size() != num_surf )
01057     {
01058         std::cerr << "ERROR: Expected " << num_surf << " surfaces in input, found " << surfs.size() << std::endl;
01059         return MB_FAILURE;
01060     }
01061     const EntityHandle explicit_vol = explicit_vols.front();
01062 
01063     // start particle
01064     double point[]     = { -0.9, 0, 0 };
01065     const double dir[] = { 1, 0, 0 };
01066     EntityHandle vol   = explicit_vol;
01067     int result;
01068     const int INSIDE = 1;  // OUTSIDE = 0, BOUNDARY = -1;
01069     rval             = gqt->point_in_volume( explicit_vol, point, result, dir );CHKERR;
01070     if( result != INSIDE )
01071     {
01072         std::cerr << "ERROR: particle not inside explicit volume" << std::endl;
01073         return MB_FAILURE;
01074     }
01075 
01076     // get next surface
01077     double dist;
01078     EntityHandle next_surf;
01079     GeomQueryTool::RayHistory history;
01080     rval = gqt->ray_fire( vol, point, dir, next_surf, dist, &history );CHKERR;
01081     if( next_surf != surfs[7] || fabs( dist - 0.91 ) > 1e-6 )
01082     {
01083         std::cerr << "ERROR: failed on advance 1" << std::endl;
01084         return MB_FAILURE;
01085     }
01086     for( unsigned i = 0; i < 3; i++ )
01087         point[i] += dist * dir[i];
01088 
01089     // get the next volume (implicit complement)
01090     EntityHandle next_vol;
01091     rval = gqt->gttool()->next_vol( next_surf, vol, next_vol );CHKERR;
01092 
01093     // get the next surface (behind numerical location)
01094     vol  = next_vol;
01095     rval = gqt->ray_fire( vol, point, dir, next_surf, dist, &history );CHKERR;
01096     if( next_surf != surfs[3] || fabs( dist - 0.0 ) > 1e-6 )
01097     {
01098         std::cerr << "ERROR: failed on advance 2" << std::endl;
01099         return MB_FAILURE;
01100     }
01101     for( unsigned i = 0; i < 3; i++ )
01102         point[i] += dist * dir[i];
01103 
01104     // get the next volume (the explicit volume)
01105     rval = gqt->gttool()->next_vol( next_surf, vol, next_vol );CHKERR;
01106 
01107     // get the next surface
01108     vol  = next_vol;
01109     rval = gqt->ray_fire( vol, point, dir, next_surf, dist, &history );CHKERR;
01110     if( next_surf != surfs[1] || fabs( dist - 0.99 ) > 1e-6 )
01111     {
01112         std::cerr << "ERROR: failed on advance 3" << std::endl;
01113         return MB_FAILURE;
01114     }
01115     for( unsigned i = 0; i < 3; i++ )
01116         point[i] += dist * dir[i];
01117 
01118     return MB_SUCCESS;
01119 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines