MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 <math.h> 00017 #include <limits> 00018 #include <algorithm> 00019 #include <stdio.h> // 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 }