MOAB: Mesh Oriented datABase  (version 5.4.0)
obb_test.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/Range.hpp"
00003 #include "moab/OrientedBoxTreeTool.hpp"
00004 #include "moab/OrientedBox.hpp"
00005 #include "MBTagConventions.hpp"
00006 #include "moab/GeomUtil.hpp"
00007 #include "moab/CN.hpp"
00008 
00009 #ifdef MOAB_HAVE_MPI
00010 #include "moab_mpi.h"
00011 #endif
00012 
00013 #include "../TestUtil.hpp"
00014 #include <iostream>
00015 #include <sstream>
00016 #include <cstdlib>
00017 #include <limits>
00018 #include <cstdio>
00019 #include <set>
00020 #include <map>
00021 
00022 using namespace moab;
00023 
00024 struct RayTest
00025 {
00026     const char* description;
00027     unsigned expected_hits;
00028     CartVect point, direction;
00029 };
00030 
00031 static const char* NAME = "obb_test";
00032 
00033 std::map< std::string, std::vector< RayTest > > default_files_tests;
00034 typedef std::map< std::string, std::vector< RayTest > >::iterator ray_test_itr;
00035 
00036 void initialize_default_files()
00037 {
00038     size_t num_tests;
00039     std::vector< RayTest > tests;
00040     std::string file = STRINGIFY( MESHDIR ) "/3k-tri-sphere.vtk";
00041     RayTest set1[]   = { { "triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 99.8792, -5, 0.121729 ) },
00042                        { "triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( 4.99167, 0, 99.7502 ) },
00043                        { "triangle node ", 6, CartVect( 0, 0, 0 ), CartVect( 0, 0, 100 ) } };
00044 
00045     num_tests = sizeof( set1 ) / sizeof( set1[0] );
00046     tests.insert( tests.begin(), &set1[0], &set1[num_tests] );
00047     default_files_tests[file] = tests;
00048     tests.clear();
00049 
00050 #ifdef MOAB_HAVE_HDF5
00051     file = STRINGIFY( MESHDIR ) "/3k-tri-cube.h5m";
00052 #else
00053     file = STRINGIFY( MESHDIR ) "/3k-tri-cube.vtk";
00054 #endif
00055 
00056     RayTest set2[] = { { "interior triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 0, 0, 5 ) },
00057                        { "interior triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( -0.25, -0.05, 5 ) },
00058                        { "interior triangle node ", 5, CartVect( 0, 0, 0 ), CartVect( -0.3, -0.3, 5 ) },
00059                        { "edge triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 5, 2.9, 2.9 ) },
00060                        { "edge triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( 5, 5, 2.9 ) },
00061                        { "edge triangle node ", 6, CartVect( 0, 0, 0 ), CartVect( 5, 5, 3 ) },
00062                        { "corner triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 5, 4.9, 4.9 ) },
00063                        { "corner triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( 5, 5, 4.9 ) },
00064                        { "corner triangle node ", 3, CartVect( 0, 0, 0 ), CartVect( 5, 5, 5 ) } };
00065 
00066     num_tests = sizeof( set2 ) / sizeof( set2[0] );
00067     tests.insert( tests.begin(), &set2[0], &set2[num_tests] );
00068     default_files_tests[file] = tests;
00069     tests.clear();
00070 }
00071 
00072 // Another possible file
00073 // STRINGIFY(MESHDIR) "/4k-tri-plane.vtk",
00074 
00075 static void usage( const char* error, const char* opt )
00076 {
00077     const char* default_message = "Invalid option";
00078     if( opt && !error ) error = default_message;
00079 
00080     std::ostream& str = error ? std::cerr : std::cout;
00081     if( error )
00082     {
00083         str << error;
00084         if( opt ) str << ": " << opt;
00085         str << std::endl;
00086     }
00087 
00088     str << "Usage: " << NAME << " [output opts.] [settings] [file ...]" << std::endl;
00089     str << "       " << NAME << " -h" << std::endl;
00090     if( !error )
00091         str << " If no file is specified the defautl test files will be used" << std::endl
00092             << " -h  print help text. " << std::endl
00093             << " -v  verbose output (may be specified multiple times) " << std::endl
00094             << " -q  quiet (minimal output) " << std::endl
00095             << " -f <x>:<y>:<z>:<i>:<j>:<k> Do ray fire" << std::endl
00096             << " -c  write box geometry to Cubit journal file." << std::endl
00097             << " -k  write leaf contents to vtk files." << std::endl
00098             << " -K  write contents of leaf boxes intersected by rays to vtk file." << std::endl
00099             << " -o <name> Save file containing tree and triangles.  Mesh tag \"OBB_ROOT\"." << std::endl
00100             << " -t <real> specify tolerance" << std::endl
00101             << " -n <int>  specify max entities per leaf node " << std::endl
00102             << " -l <int>  specify max tree levels" << std::endl
00103             << " -r <real> specify worst cell split ratio" << std::endl
00104             << " -R <real> specify best cell split ratio" << std::endl
00105             << " -s force construction of surface tree" << std::endl
00106             << " -S do not build surface tree." << std::endl
00107             << "    (Default: surface tree if file contains multiple surfaces" << std::endl
00108             << " -u  use unordered (Range) meshsets for tree nodes" << std::endl
00109             << " -U  use ordered (vector) meshsets for tree nodes" << std::endl
00110             << " Verbosity (-q sets to 0, each -v increments, default is 1):" << std::endl
00111             << "  0 - no output" << std::endl
00112             << "  1 - status messages and error summary" << std::endl
00113             << "  2 - print tree statistics " << std::endl
00114             << "  3 - output errors for each node" << std::endl
00115             << "  4 - print tree" << std::endl
00116             << "  5 - print tree w/ contents of each node" << std::endl
00117             << " See documentation for OrientedBoxTreeTool::Settings for " << std::endl
00118             << " a description of tree generation settings." << std::endl;
00119     exit( !!error );
00120 }
00121 
00122 static const char* get_option( int& i, int argc, char* argv[] )
00123 {
00124     ++i;
00125     if( i == argc ) usage( "Expected argument following option", argv[i - 1] );
00126     return argv[i];
00127 }
00128 
00129 static int get_int_option( int& i, int argc, char* argv[] )
00130 {
00131     const char* str = get_option( i, argc, argv );
00132     char* end_ptr;
00133     long val = strtol( str, &end_ptr, 0 );
00134     if( !*str || *end_ptr ) usage( "Expected integer following option", argv[i - 1] );
00135     return val;
00136 }
00137 
00138 static double get_double_option( int& i, int argc, char* argv[] )
00139 {
00140     const char* str = get_option( i, argc, argv );
00141     char* end_ptr;
00142     double val = strtod( str, &end_ptr );
00143     if( !*str || *end_ptr ) usage( "Expected real number following option", argv[i - 1] );
00144     return val;
00145 }
00146 
00147 static bool do_file( const char* filename );
00148 
00149 enum TriOption
00150 {
00151     DISABLE,
00152     ENABLE,
00153     AUTO
00154 };
00155 
00156 static int verbosity = 1;
00157 static OrientedBoxTreeTool::Settings settings;
00158 static double tolerance   = 1e-6;
00159 static bool write_cubit   = false;
00160 static bool write_vtk     = false;
00161 static bool write_ray_vtk = false;
00162 static std::vector< CartVect > rays;
00163 static const char* save_file_name = 0;
00164 static TriOption surfTree         = AUTO;
00165 
00166 static void parse_ray( int& i, int argc, char* argv[] );
00167 
00168 int main( int argc, char* argv[] )
00169 {
00170 #ifdef MOAB_HAVE_MPI
00171     int fail = MPI_Init( &argc, &argv );
00172     if( fail ) return fail;
00173 #endif
00174 
00175     initialize_default_files();
00176 
00177     std::vector< const char* > file_names;
00178     bool flags = true;
00179     for( int i = 1; i < argc; ++i )
00180     {
00181         if( flags && argv[i][0] == '-' )
00182         {
00183             if( !argv[i][1] || argv[i][2] ) usage( 0, argv[i] );
00184             switch( argv[i][1] )
00185             {
00186                 default:
00187                     usage( 0, argv[i] );
00188                     break;
00189                 case '-':
00190                     flags = false;
00191                     break;
00192                 case 'v':
00193                     ++verbosity;
00194                     break;
00195                 case 'q':
00196                     verbosity = 0;
00197                     break;
00198                 case 'h':
00199                     usage( 0, 0 );
00200                     break;
00201                 case 'c':
00202                     write_cubit = true;
00203                     break;
00204                 case 'k':
00205                     write_vtk = true;
00206                     break;
00207                 case 'K':
00208                     write_ray_vtk = true;
00209                     break;
00210                 case 's':
00211                     surfTree = ENABLE;
00212                     break;
00213                 case 'S':
00214                     surfTree = DISABLE;
00215                     break;
00216                 case 'u':
00217                     settings.set_options = MESHSET_SET;
00218                     break;
00219                 case 'U':
00220                     settings.set_options = MESHSET_ORDERED;
00221                     break;
00222                 case 'o':
00223                     save_file_name = get_option( i, argc, argv );
00224                     break;
00225                 case 'n':
00226                     settings.max_leaf_entities = get_int_option( i, argc, argv );
00227                     break;
00228                 case 'l':
00229                     settings.max_depth = get_int_option( i, argc, argv );
00230                     break;
00231                 case 'r':
00232                     settings.worst_split_ratio = get_double_option( i, argc, argv );
00233                     break;
00234                 case 'R':
00235                     settings.best_split_ratio = get_double_option( i, argc, argv );
00236                     break;
00237                 case 't':
00238                     tolerance = get_double_option( i, argc, argv );
00239                     break;
00240                 case 'f':
00241                     parse_ray( i, argc, argv );
00242                     break;
00243             }
00244         }
00245         else
00246         {
00247             file_names.push_back( argv[i] );
00248         }
00249     }
00250 
00251     if( verbosity )
00252     {
00253         Core core;
00254         std::string version;
00255         core.impl_version( &version );
00256         std::cout << version << std::endl;
00257         if( verbosity > 1 )
00258             std::cout << "max_leaf_entities:      " << settings.max_leaf_entities << std::endl
00259                       << "max_depth:              " << settings.max_depth << std::endl
00260                       << "worst_split_ratio:      " << settings.worst_split_ratio << std::endl
00261                       << "best_split_ratio:       " << settings.best_split_ratio << std::endl
00262                       << "tolerance:              " << tolerance << std::endl
00263                       << "set type:               "
00264                       << ( ( settings.set_options & MESHSET_ORDERED ) ? "ordered" : "set" ) << std::endl
00265                       << std::endl;
00266     }
00267 
00268     if( !settings.valid() || tolerance < 0.0 )
00269     {
00270         std::cerr << "Invalid settings specified." << std::endl;
00271         return 2;
00272     }
00273 
00274     if( file_names.empty() )
00275     {
00276         std::cerr << "No file(s) specified." << std::endl;
00277         for( ray_test_itr file = default_files_tests.begin(); file != default_files_tests.end(); file++ )
00278         {
00279             std::cerr << "Using default file \"" << file->first << '"' << std::endl;
00280             file_names.push_back( file->first.c_str() );
00281         }
00282     }
00283 
00284     if( save_file_name && file_names.size() != 1 )
00285     {
00286         std::cout << "Only one input file allowed if \"-o\" option is specified." << std::endl;
00287         std::cout << "Only testing with single file " << file_names[0] << std::endl;
00288         file_names.erase( ++file_names.begin(), file_names.end() );
00289     }
00290 
00291     int exit_val = 0;
00292     for( unsigned j = 0; j < file_names.size(); ++j )
00293         if( !do_file( file_names[j] ) ) ++exit_val;
00294 
00295 #ifdef MOAB_HAVE_MPI
00296     fail = MPI_Finalize();
00297     if( fail ) return fail;
00298 #endif
00299 
00300     return exit_val ? exit_val + 2 : 0;
00301 }
00302 
00303 static void parse_ray( int& i, int argc, char* argv[] )
00304 {
00305     CartVect point, direction;
00306     if( 6 != sscanf( get_option( i, argc, argv ), "%lf:%lf:%lf:%lf:%lf:%lf", &point[0], &point[1], &point[2],
00307                      &direction[0], &direction[1], &direction[2] ) )
00308         usage( "Expected ray specified as <x>:<y>:<z>:<i>:<j>:<k>", 0 );
00309     direction.normalize();
00310     rays.push_back( point );
00311     rays.push_back( direction );
00312 }
00313 
00314 class TreeValidator : public OrientedBoxTreeTool::Op
00315 {
00316   private:
00317     Interface* const instance;
00318     OrientedBoxTreeTool* const tool;
00319     const bool printing;
00320     const double epsilon;
00321     bool surfaces;
00322     std::ostream& stream;
00323     OrientedBoxTreeTool::Settings settings;
00324 
00325     void print( EntityHandle handle, const char* string )
00326     {
00327         if( printing ) stream << instance->id_from_handle( handle ) << ": " << string << std::endl;
00328     }
00329 
00330     ErrorCode error( EntityHandle handle, const char* string )
00331     {
00332         ++error_count;
00333         print( handle, string );
00334         return MB_SUCCESS;
00335     }
00336 
00337   public:
00338     unsigned entity_count;
00339 
00340     unsigned loose_box_count;
00341     unsigned child_outside_count;
00342     unsigned entity_outside_count;
00343     unsigned num_entities_outside;
00344     unsigned non_ortho_count;
00345     unsigned error_count;
00346     unsigned empty_leaf_count;
00347     unsigned non_empty_non_leaf_count;
00348     unsigned entity_invalid_count;
00349     unsigned unsorted_axis_count;
00350     unsigned non_unit_count;
00351     unsigned duplicate_entity_count;
00352     unsigned bad_outer_radius_count;
00353     unsigned missing_surface_count;
00354     unsigned multiple_surface_count;
00355     std::set< EntityHandle > seen;
00356     int surface_depth;
00357     EntityHandle surface_handle;
00358 
00359     TreeValidator( Interface* instance_ptr,
00360                    OrientedBoxTreeTool* tool_ptr,
00361                    bool print_errors,
00362                    std::ostream& str,
00363                    double tol,
00364                    bool surfs,
00365                    OrientedBoxTreeTool::Settings s )
00366         : instance( instance_ptr ), tool( tool_ptr ), printing( print_errors ), epsilon( tol ), surfaces( surfs ),
00367           stream( str ), settings( s ), entity_count( 0 ), loose_box_count( 0 ), child_outside_count( 0 ),
00368           entity_outside_count( 0 ), num_entities_outside( 0 ), non_ortho_count( 0 ), error_count( 0 ),
00369           empty_leaf_count( 0 ), non_empty_non_leaf_count( 0 ), entity_invalid_count( 0 ), unsorted_axis_count( 0 ),
00370           non_unit_count( 0 ), duplicate_entity_count( 0 ), bad_outer_radius_count( 0 ), missing_surface_count( 0 ),
00371           multiple_surface_count( 0 ), surface_depth( -1 ), surface_handle( 0 )
00372     {
00373     }
00374 
00375     bool is_valid() const
00376     {
00377         return 0 == loose_box_count + child_outside_count + entity_outside_count + num_entities_outside +
00378                         non_ortho_count + error_count + empty_leaf_count + non_empty_non_leaf_count +
00379                         entity_invalid_count + unsorted_axis_count + non_unit_count + duplicate_entity_count +
00380                         bad_outer_radius_count + missing_surface_count + multiple_surface_count;
00381     }
00382 
00383     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
00384 
00385     virtual ErrorCode leaf( EntityHandle )
00386     {
00387         return MB_SUCCESS;
00388     }
00389 };
00390 
00391 ErrorCode TreeValidator::visit( EntityHandle node, int depth, bool& descend )
00392 {
00393     ErrorCode rval;
00394     descend = true;
00395 
00396     Range contents;
00397     rval = instance->get_entities_by_handle( node, contents );
00398     if( MB_SUCCESS != rval ) return error( node, "Error getting contents of tree node.  Corrupt tree?" );
00399     entity_count += contents.size();
00400 
00401     if( surfaces )
00402     {
00403         // if no longer in subtree for previous surface, clear
00404         if( depth <= surface_depth ) surface_depth = -1;
00405 
00406         EntityHandle surface      = 0;
00407         Range::iterator surf_iter = contents.lower_bound( MBENTITYSET );
00408         if( surf_iter != contents.end() )
00409         {
00410             surface = *surf_iter;
00411             contents.erase( surf_iter );
00412         }
00413 
00414         if( surface )
00415         {
00416             if( surface_depth >= 0 )
00417             {
00418                 ++multiple_surface_count;
00419                 print( node, "Multiple surfaces in encountered in same subtree." );
00420             }
00421             else
00422             {
00423                 surface_depth  = depth;
00424                 surface_handle = surface;
00425             }
00426         }
00427     }
00428 
00429     std::vector< EntityHandle > children;
00430     rval = tool->get_moab_instance()->get_child_meshsets( node, children );
00431     if( MB_SUCCESS != rval || ( !children.empty() && children.size() != 2 ) )
00432         return error( node, "Error getting children.  Corrupt tree?" );
00433 
00434     OrientedBox box;
00435     rval = tool->box( node, box );
00436     if( MB_SUCCESS != rval ) return error( node, "Error getting oriented box from tree node.  Corrupt tree?" );
00437 
00438     if( children.empty() && contents.empty() )
00439     {
00440         ++empty_leaf_count;
00441         print( node, "Empty leaf node.\n" );
00442     }
00443     else if( !children.empty() && !contents.empty() )
00444     {
00445         ++non_empty_non_leaf_count;
00446         print( node, "Non-leaf node is not empty." );
00447     }
00448 
00449     if( surfaces && children.empty() && surface_depth < 0 )
00450     {
00451         ++missing_surface_count;
00452         print( node, "Reached leaf node w/out encountering any surface set." );
00453     }
00454 
00455     double dot_epsilon = epsilon * ( box.axis( 0 ) + box.axis( 1 ) + box.axis( 2 ) ).length();
00456     if( box.axis( 0 ) % box.axis( 1 ) > dot_epsilon || box.axis( 0 ) % box.axis( 2 ) > dot_epsilon ||
00457         box.axis( 1 ) % box.axis( 2 ) > dot_epsilon )
00458     {
00459         ++non_ortho_count;
00460         print( node, "Box axes are not orthogonal" );
00461     }
00462 
00463     if( !children.empty() )
00464     {
00465         for( int i = 0; i < 2; ++i )
00466         {
00467             OrientedBox other_box;
00468             rval = tool->box( children[i], other_box );
00469             if( MB_SUCCESS != rval )
00470                 return error( children[i], " Error getting oriented box from tree node.  Corrupt tree?" );
00471             //      else if (!box.contained( other_box, epsilon )) {
00472             //        ++child_outside_count;
00473             //        print( children[i], "Parent box does not contain child box." );
00474             //        char string[64];
00475             //        sprintf(string, "     Volume ratio is %f", other_box.volume()/box.volume() );
00476             //        print( children [i], string );
00477             //      }
00478             else
00479             {
00480                 double vol_ratio = other_box.volume() / box.volume();
00481                 if( vol_ratio > 2.0 )
00482                 {
00483                     char string[64];
00484                     sprintf( string, "child/parent volume ratio is %f", vol_ratio );
00485                     print( children[i], string );
00486                     sprintf( string, "   child/parent area ratio is %f", other_box.area() / box.area() );
00487                     print( children[i], string );
00488                 }
00489             }
00490         }
00491     }
00492 
00493     bool bad_element        = false;
00494     bool bad_element_handle = false;
00495     bool bad_element_conn   = false;
00496     bool duplicate_element  = false;
00497     int num_outside         = 0;
00498     bool boundary[6]        = { false, false, false, false, false, false };
00499     for( Range::iterator it = contents.begin(); it != contents.end(); ++it )
00500     {
00501         EntityType type = instance->type_from_handle( *it );
00502         int dim         = CN::Dimension( type );
00503         if( dim != 2 )
00504         {
00505             bad_element = true;
00506             continue;
00507         }
00508 
00509         const EntityHandle* conn;
00510         int conn_len;
00511         rval = instance->get_connectivity( *it, conn, conn_len );
00512         if( MB_SUCCESS != rval )
00513         {
00514             bad_element_handle = true;
00515             continue;
00516         }
00517 
00518         std::vector< CartVect > coords( conn_len );
00519         rval = instance->get_coords( conn, conn_len, coords[0].array() );
00520         if( MB_SUCCESS != rval )
00521         {
00522             bad_element_conn = true;
00523             continue;
00524         }
00525 
00526         bool outside = false;
00527         for( std::vector< CartVect >::iterator j = coords.begin(); j != coords.end(); ++j )
00528         {
00529             if( !box.contained( *j, epsilon ) )
00530                 outside = true;
00531             else
00532                 for( int d = 0; d < 3; ++d )
00533                 {
00534 #if MB_ORIENTED_BOX_UNIT_VECTORS
00535                     double n = box.axis( d ) % ( *j - box.center );
00536                     if( fabs( n - box.length[d] ) <= epsilon ) boundary[2 * d] = true;
00537                     if( fabs( n + box.length[d] ) <= epsilon ) boundary[2 * d + 1] = true;
00538 #else
00539                     double ln   = box.axis( d ).length();
00540                     CartVect v1 = *j - box.center - box.axis[d];
00541                     CartVect v2 = *j - box.center + box.axis[d];
00542                     if( fabs( v1 % box.axis[d] ) <= ln * epsilon ) boundary[2 * d] = true;
00543                     if( fabs( v2 % box.axis[d] ) <= ln * epsilon ) boundary[2 * d + 1] = true;
00544 #endif
00545                 }
00546         }
00547         if( outside ) ++num_outside;
00548 
00549         if( !seen.insert( *it ).second )
00550         {
00551             duplicate_element = true;
00552             ++duplicate_entity_count;
00553         }
00554     }
00555 
00556     CartVect alength( box.axis( 0 ).length(), box.axis( 1 ).length(), box.axis( 2 ).length() );
00557 #if MB_ORIENTED_BOX_UNIT_VECTORS
00558     CartVect length = box.length;
00559 #else
00560     CartVect length = alength;
00561 #endif
00562 
00563     if( length[0] > length[1] || length[0] > length[2] || length[1] > length[2] )
00564     {
00565         ++unsorted_axis_count;
00566         print( node, "Box axes are not ordered from shortest to longest." );
00567     }
00568 
00569 #if MB_ORIENTED_BOX_UNIT_VECTORS
00570     if( fabs( alength[0] - 1.0 ) > epsilon || fabs( alength[1] - 1.0 ) > epsilon || fabs( alength[2] - 1.0 ) > epsilon )
00571     {
00572         ++non_unit_count;
00573         print( node, "Box axes are not unit vectors." );
00574     }
00575 #endif
00576 
00577 #if MB_ORIENTED_BOX_OUTER_RADIUS
00578     if( fabs( length.length() - box.radius ) > tolerance )
00579     {
00580         ++bad_outer_radius_count;
00581         print( node, "Box has incorrect outer radius." );
00582     }
00583 #endif
00584 
00585     if( depth + 1 < settings.max_depth && contents.size() > (unsigned)( 4 * settings.max_leaf_entities ) )
00586     {
00587         char string[64];
00588         sprintf( string, "leaf at depth %d with %u entities", depth, (unsigned)contents.size() );
00589         print( node, string );
00590     }
00591 
00592     bool all_boundaries = true;
00593     for( int f = 0; f < 6; ++f )
00594         all_boundaries = all_boundaries && boundary[f];
00595 
00596     if( bad_element )
00597     {
00598         ++entity_invalid_count;
00599         print( node, "Set contained an entity with an inappropriate dimension." );
00600     }
00601     if( bad_element_handle )
00602     {
00603         ++error_count;
00604         print( node, "Error querying face contained in set." );
00605     }
00606     if( bad_element_conn )
00607     {
00608         ++error_count;
00609         print( node, "Error querying connectivity of element." );
00610     }
00611     if( duplicate_element )
00612     {
00613         print( node, "Elements occur in multiple leaves of tree." );
00614     }
00615     if( num_outside > 0 )
00616     {
00617         ++entity_outside_count;
00618         num_entities_outside += num_outside;
00619         if( printing )
00620             stream << instance->id_from_handle( node ) << ": " << num_outside << " elements outside box." << std::endl;
00621     }
00622     else if( !all_boundaries && !contents.empty() )
00623     {
00624         ++loose_box_count;
00625         print( node, "Box does not fit contained elements tightly." );
00626     }
00627 
00628     return MB_SUCCESS;
00629 }
00630 
00631 class CubitWriter : public OrientedBoxTreeTool::Op
00632 {
00633   public:
00634     CubitWriter( FILE* file_ptr, OrientedBoxTreeTool* tool_ptr ) : file( file_ptr ), tool( tool_ptr ) {}
00635 
00636     ErrorCode visit( EntityHandle node, int depth, bool& descend );
00637     ErrorCode leaf( EntityHandle )
00638     {
00639         return MB_SUCCESS;
00640     }
00641 
00642   private:
00643     FILE* file;
00644     OrientedBoxTreeTool* tool;
00645 };
00646 
00647 ErrorCode CubitWriter::visit( EntityHandle node, int, bool& descend )
00648 {
00649     descend = true;
00650     OrientedBox box;
00651     ErrorCode rval = tool->box( node, box );
00652     if( rval != MB_SUCCESS ) return rval;
00653 
00654     double sign[] = { -1, 1 };
00655     for( int i = 0; i < 2; ++i )
00656         for( int j = 0; j < 2; ++j )
00657             for( int k = 0; k < 2; ++k )
00658             {
00659 #if MB_ORIENTED_BOX_UNIT_VECTORS
00660                 CartVect corner = box.center + box.length[0] * sign[i] * box.axis( 0 ) +
00661                                   box.length[1] * sign[j] * box.axis( 1 ) + box.length[2] * sign[k] * box.axis( 2 );
00662 #else
00663                 CartVect corner =
00664                     box.center + sign[i] * box.axis( 0 ) + sign[j] * box.axis( 1 ) + sign[k] * box.axis( 2 );
00665 #endif
00666                 fprintf( file, "create vertex %f %f %f\n", corner[0], corner[1], corner[2] );
00667             }
00668     fprintf( file, "#{i=Id(\"vertex\")-7}\n" );
00669     fprintf( file, "create surface vertex {i  } {i+1} {i+3} {i+2}\n" );
00670     fprintf( file, "create surface vertex {i+4} {i+5} {i+7} {i+6}\n" );
00671     fprintf( file, "create surface vertex {i+1} {i+0} {i+4} {i+5}\n" );
00672     fprintf( file, "create surface vertex {i+3} {i+2} {i+6} {i+7}\n" );
00673     fprintf( file, "create surface vertex {i+2} {i+0} {i+4} {i+6}\n" );
00674     fprintf( file, "create surface vertex {i+3} {i+1} {i+5} {i+7}\n" );
00675     fprintf( file, "delete vertex {i}-{i+7}\n" );
00676     fprintf( file, "#{s=Id(\"surface\")-5}\n" );
00677     fprintf( file, "create volume surface {s} {s+1} {s+2} {s+3} {s+4} {s+5} noheal\n" );
00678     int id = tool->get_moab_instance()->id_from_handle( node );
00679     fprintf( file, "volume {Id(\"volume\")} name \"cell%d\"\n", id );
00680 
00681     return MB_SUCCESS;
00682 }
00683 
00684 class VtkWriter : public OrientedBoxTreeTool::Op
00685 {
00686   public:
00687     VtkWriter( std::string base_name, Interface* interface ) : baseName( base_name ), instance( interface ) {}
00688 
00689     ErrorCode visit( EntityHandle node, int depth, bool& descend );
00690 
00691     ErrorCode leaf( EntityHandle /*node*/ )
00692     {
00693         return MB_SUCCESS;
00694     }
00695 
00696   private:
00697     std::string baseName;
00698     Interface* instance;
00699 };
00700 
00701 ErrorCode VtkWriter::visit( EntityHandle node, int, bool& descend )
00702 {
00703     descend = true;
00704 
00705     ErrorCode rval;
00706     int count;
00707     rval = instance->get_number_entities_by_handle( node, count );
00708     if( MB_SUCCESS != rval || 0 == count ) return rval;
00709 
00710     int id = instance->id_from_handle( node );
00711     char id_str[32];
00712     sprintf( id_str, "%d", id );
00713 
00714     std::string file_name( baseName );
00715     file_name += ".";
00716     file_name += id_str;
00717     file_name += ".vtk";
00718 
00719     return instance->write_mesh( file_name.c_str(), &node, 1 );
00720 }
00721 
00722 static bool do_ray_fire_test( OrientedBoxTreeTool& tool,
00723                               EntityHandle root_set,
00724                               const char* filename,
00725                               bool have_surface_tree );
00726 
00727 static bool do_closest_point_test( OrientedBoxTreeTool& tool, EntityHandle root_set, bool have_surface_tree );
00728 
00729 static ErrorCode save_tree( Interface* instance, const char* filename, EntityHandle tree_root );
00730 
00731 static bool do_file( const char* filename )
00732 {
00733     ErrorCode rval;
00734     Core instance;
00735     Interface* const iface = &instance;
00736     OrientedBoxTreeTool tool( iface );
00737     bool haveSurfTree = false;
00738 
00739     if( verbosity ) std::cout << filename << std::endl << "------" << std::endl;
00740 
00741     rval = iface->load_mesh( filename );
00742     if( MB_SUCCESS != rval )
00743     {
00744         if( verbosity ) std::cout << "Failed to read file: \"" << filename << '"' << std::endl;
00745         return false;
00746     }
00747 
00748     // IF building from surfaces, get surfaces.
00749     // If AUTO and less than two surfaces, don't build from surfaces.
00750     Range surfaces;
00751     if( surfTree != DISABLE )
00752     {
00753         Tag surftag;
00754         rval = iface->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, surftag );
00755         if( MB_SUCCESS == rval )
00756         {
00757             int dim                 = 2;
00758             const void* tagvalues[] = { &dim };
00759             rval = iface->get_entities_by_type_and_tag( 0, MBENTITYSET, &surftag, tagvalues, 1, surfaces );
00760             if( MB_SUCCESS != rval && MB_ENTITY_NOT_FOUND != rval ) return false;
00761         }
00762         else if( MB_TAG_NOT_FOUND != rval )
00763             return false;
00764 
00765         if( ENABLE == surfTree && surfaces.empty() )
00766         {
00767             std::cerr << "No Surfaces found." << std::endl;
00768             return false;
00769         }
00770 
00771         haveSurfTree = ( ENABLE == surfTree ) || ( surfaces.size() > 1 );
00772     }
00773 
00774     EntityHandle root;
00775     Range entities;
00776     if( !haveSurfTree )
00777     {
00778         rval = iface->get_entities_by_dimension( 0, 2, entities );
00779         if( MB_SUCCESS != rval )
00780         {
00781             std::cerr << "get_entities_by_dimension( 2 ) failed." << std::endl;
00782             return false;
00783         }
00784 
00785         if( entities.empty() )
00786         {
00787             if( verbosity ) std::cout << "File \"" << filename << "\" contains no 2D elements" << std::endl;
00788             return false;
00789         }
00790 
00791         if( verbosity ) std::cout << "Building tree from " << entities.size() << " 2D elements" << std::endl;
00792 
00793         rval = tool.build( entities, root, &settings );
00794         if( MB_SUCCESS != rval )
00795         {
00796             if( verbosity ) std::cout << "Failed to build tree." << std::endl;
00797             return false;
00798         }
00799     }
00800     else
00801     {
00802 
00803         if( verbosity ) std::cout << "Building tree from " << surfaces.size() << " surfaces" << std::endl;
00804 
00805         // Build subtree for each surface, get list of all entities to use later
00806         Range surf_trees, surf_tris;
00807         EntityHandle surf_root;
00808         for( Range::iterator s = surfaces.begin(); s != surfaces.end(); ++s )
00809         {
00810             surf_tris.clear();
00811             rval = iface->get_entities_by_dimension( *s, 2, surf_tris );
00812             if( MB_SUCCESS != rval ) return false;
00813             rval = tool.build( surf_tris, surf_root, &settings );
00814             if( MB_SUCCESS != rval )
00815             {
00816                 if( verbosity ) std::cout << "Failed to build tree for surface." << std::endl;
00817                 return false;
00818             }
00819             surf_trees.insert( surf_root );
00820             entities.merge( surf_tris );
00821             rval = iface->add_entities( surf_root, &*s, 1 );
00822             if( MB_SUCCESS != rval ) return false;
00823         }
00824 
00825         rval = tool.join_trees( surf_trees, root, &settings );
00826         if( MB_SUCCESS != rval )
00827         {
00828             if( verbosity ) std::cout << "Failed to build tree." << std::endl;
00829             return false;
00830         }
00831 
00832         if( verbosity )
00833             std::cout << "Built tree from " << surfaces.size() << " surfaces"
00834                       << " (" << entities.size() - surfaces.size() << " elements)" << std::endl;
00835 
00836         entities.merge( surfaces );
00837     }
00838 
00839     if( write_cubit )
00840     {
00841         std::string name = filename;
00842         name += ".boxes.jou";
00843         FILE* ptr = fopen( name.c_str(), "w+" );
00844         if( !ptr )
00845         {
00846             perror( name.c_str() );
00847         }
00848         else
00849         {
00850             if( verbosity ) std::cout << "Writing: " << name << std::endl;
00851             fprintf( ptr, "graphics off\n" );
00852             CubitWriter op( ptr, &tool );
00853             tool.preorder_traverse( root, op );
00854             fprintf( ptr, "graphics on\n" );
00855             fclose( ptr );
00856         }
00857     }
00858 
00859     if( write_vtk )
00860     {
00861         VtkWriter op( filename, iface );
00862         if( verbosity )
00863             std::cout << "Writing leaf contents as : " << filename << ".xxx.vtk where 'xxx' is the set id" << std::endl;
00864         tool.preorder_traverse( root, op );
00865     }
00866 
00867     bool print_errors = ( verbosity > 2 ), print_contents = ( verbosity > 4 );
00868     if( verbosity > 3 ) tool.print( root, std::cout, print_contents );
00869     if( verbosity > 1 )
00870     {
00871         rval = tool.stats( root, std::cout );
00872         if( MB_SUCCESS != rval ) std::cout << "****** Failed to get tree statistics ******" << std::endl;
00873     }
00874 
00875     TreeValidator op( iface, &tool, print_errors, std::cout, tolerance, haveSurfTree, settings );
00876     rval        = tool.preorder_traverse( root, op );
00877     bool result = op.is_valid();
00878     if( MB_SUCCESS != rval )
00879     {
00880         result = false;
00881         if( verbosity ) std::cout << "Errors traversing tree.  Corrupt tree?" << std::endl;
00882     }
00883 
00884     bool missing = ( op.entity_count != entities.size() );
00885     if( missing ) result = false;
00886 
00887     if( verbosity )
00888     {
00889         if( result )
00890             std::cout << std::endl << "No errors detected." << std::endl;
00891         else
00892             std::cout << std::endl << "*********************** ERROR SUMMARY **********************" << std::endl;
00893         if( op.child_outside_count )
00894             std::cout << "* " << op.child_outside_count << " child boxes not contained in parent." << std::endl;
00895         if( op.entity_outside_count )
00896             std::cout << "* " << op.entity_outside_count << " nodes containing entities outside of box." << std::endl;
00897         if( op.num_entities_outside )
00898             std::cout << "* " << op.num_entities_outside << " entities outside boxes." << std::endl;
00899         if( op.empty_leaf_count ) std::cout << "* " << op.empty_leaf_count << " empty leaf nodes." << std::endl;
00900         if( op.non_empty_non_leaf_count )
00901             std::cout << "* " << op.non_empty_non_leaf_count << " non-leaf nodes containing entities." << std::endl;
00902         if( op.duplicate_entity_count )
00903             std::cout << "* " << op.duplicate_entity_count << " duplicate entities in leaves." << std::endl;
00904         if( op.missing_surface_count )
00905             std::cout << "* " << op.missing_surface_count << " leaves outside surface subtrees." << std::endl;
00906         if( op.multiple_surface_count )
00907             std::cout << "* " << op.multiple_surface_count << " surfaces within other surface subtrees." << std::endl;
00908         if( op.non_ortho_count )
00909             std::cout << "* " << op.non_ortho_count << " boxes with non-orthononal axes." << std::endl;
00910         if( op.non_unit_count ) std::cout << "* " << op.non_unit_count << " boxes with non-unit axes." << std::endl;
00911         if( op.bad_outer_radius_count )
00912             std::cout << "* " << op.bad_outer_radius_count << " boxes incorrect outer radius." << std::endl;
00913         if( op.unsorted_axis_count )
00914             std::cout << "* " << op.unsorted_axis_count << " boxes axes in unsorted order." << std::endl;
00915         if( op.loose_box_count )
00916             std::cout << "* " << op.loose_box_count << " boxes that do not fit the contained entities tightly."
00917                       << std::endl;
00918         if( op.error_count + op.entity_invalid_count )
00919             std::cout << "* " << op.error_count + op.entity_invalid_count << " other errors while traversing tree."
00920                       << std::endl;
00921         if( missing )
00922             std::cout << "* tree built from " << entities.size() << " entities contains " << op.entity_count
00923                       << " entities." << std::endl;
00924         if( !result ) std::cout << "************************************************************" << std::endl;
00925     }
00926 
00927     if( result && save_file_name )
00928     {
00929         if( MB_SUCCESS == save_tree( iface, save_file_name, root ) )
00930             std::cerr << "Wrote '" << save_file_name << "'" << std::endl;
00931         else
00932             std::cerr << "FAILED TO WRITE '" << save_file_name << "'" << std::endl;
00933     }
00934 
00935     if( !do_ray_fire_test( tool, root, filename, haveSurfTree ) )
00936     {
00937         if( verbosity ) std::cout << "Ray fire test failed." << std::endl;
00938         result = false;
00939     }
00940 
00941     if( !do_closest_point_test( tool, root, haveSurfTree ) )
00942     {
00943         if( verbosity ) std::cout << "Closest point test failed" << std::endl;
00944         result = false;
00945     }
00946 
00947     rval = tool.delete_tree( root );
00948     if( MB_SUCCESS != rval )
00949     {
00950         if( verbosity ) std::cout << "delete_tree failed." << std::endl;
00951         result = false;
00952     }
00953 
00954     return result;
00955 }
00956 
00957 static void count_non_tol( std::vector< double > intersections, int& non_tol_count, double& non_tol_dist )
00958 {
00959 
00960     for( size_t i = 0; i < intersections.size(); ++i )
00961     {
00962         if( intersections[i] > tolerance )
00963         {
00964             ++non_tol_count;
00965             if( intersections[i] < non_tol_dist ) non_tol_dist = intersections[i];
00966         }
00967     }
00968 }
00969 
00970 static bool check_ray_intersect_tris( OrientedBoxTreeTool& tool,
00971                                       EntityHandle root_set,
00972                                       RayTest& test,
00973                                       int& non_tol_count,
00974                                       double& non_tol_dist,
00975                                       OrientedBoxTreeTool::TrvStats& stats )
00976 {
00977     ErrorCode rval;
00978     bool result = true;
00979 
00980     non_tol_dist  = std::numeric_limits< double >::max();
00981     non_tol_count = 0;
00982 
00983     std::vector< double > intersections;
00984     std::vector< EntityHandle > facets;
00985     rval = tool.ray_intersect_triangles( intersections, facets, root_set, tolerance, test.point.array(),
00986                                          test.direction.array(), 0, &stats );
00987     if( MB_SUCCESS != rval )
00988     {
00989         if( verbosity ) std::cout << "  Call to OrientedBoxTreeTool::ray_intersect_triangles failed." << std::endl;
00990         result = false;
00991     }
00992     else
00993     {
00994 
00995         if( intersections.size() != test.expected_hits )
00996         {
00997             if( verbosity > 2 )
00998                 std::cout << "  Expected " << test.expected_hits << " and got " << intersections.size()
00999                           << " hits for ray fire of " << test.description << std::endl;
01000             if( verbosity > 3 )
01001             {
01002                 for( unsigned j = 0; j < intersections.size(); ++j )
01003                     std::cout << "  " << intersections[j];
01004                 std::cout << std::endl;
01005             }
01006             result = false;
01007         }
01008 
01009         count_non_tol( intersections, non_tol_count, non_tol_dist );
01010     }
01011     return result;
01012 }
01013 
01014 static bool check_ray_intersect_sets( OrientedBoxTreeTool& tool,
01015                                       EntityHandle root_set,
01016                                       RayTest& test,
01017                                       int& non_tol_count,
01018                                       double& non_tol_dist,
01019                                       OrientedBoxTreeTool::TrvStats& stats )
01020 {
01021 
01022     ErrorCode rval;
01023     bool result = true;
01024 
01025     non_tol_dist  = std::numeric_limits< double >::max();
01026     non_tol_count = 0;
01027 
01028     const int NUM_NON_TOL_INT = 1;
01029 
01030     std::vector< double > intersections;
01031     std::vector< EntityHandle > surfs;
01032     std::vector< EntityHandle > facets;
01033 
01034     OrientedBoxTreeTool::IntersectSearchWindow search_win;
01035     OrientedBoxTreeTool::IntRegCtxt int_reg_ctxt;
01036     rval = tool.ray_intersect_sets( intersections, surfs, facets, root_set, tolerance, test.point.array(),
01037                                     test.direction.array(), search_win, int_reg_ctxt, &stats );
01038 
01039     if( MB_SUCCESS != rval )
01040     {
01041         if( verbosity ) std::cout << "  Call to OrientedBoxTreeTool::ray_intersect_sets failed." << std::endl;
01042         result = false;
01043     }
01044     else
01045     {
01046 
01047         if( surfs.size() != intersections.size() )
01048         {
01049             if( verbosity ) std::cout << "  ray_intersect_sets did not return sets for all intersections." << std::endl;
01050             result = false;
01051         }
01052 
01053         count_non_tol( intersections, non_tol_count, non_tol_dist );
01054 
01055         if( non_tol_count > NUM_NON_TOL_INT )
01056         {
01057             if( verbosity )
01058                 std::cout << "  Requested " << NUM_NON_TOL_INT << "intersections "
01059                           << "  beyond tolerance.  Got " << non_tol_count << std::endl;
01060             result = false;
01061         }
01062     }
01063 
01064     return result;
01065 }
01066 
01067 static bool do_ray_fire_test( OrientedBoxTreeTool& tool,
01068                               EntityHandle root_set,
01069                               const char* filename,
01070                               bool haveSurfTree )
01071 {
01072     if( verbosity > 1 ) std::cout << "beginning ray fire tests" << std::endl;
01073 
01074     OrientedBox box;
01075     ErrorCode rval = tool.box( root_set, box );
01076     if( MB_SUCCESS != rval )
01077     {
01078         if( verbosity ) std::cerr << "Error getting box for tree root set" << std::endl;
01079         return false;
01080     }
01081 
01082     /* Do standard ray fire tests */
01083     std::cout << box << std::endl;
01084 
01085     CartVect origin( 0., 0., 0. );
01086     CartVect unitDiag( 1., 1., 1. );
01087     std::vector< RayTest > tests;
01088     RayTest default_tests[] = { { "large axis through box", 2, box.center - 1.2 * box.scaled_axis( 2 ), box.axis( 2 ) },
01089                                 { "small axis through box", 2, box.center - 1.2 * box.scaled_axis( 0 ), box.axis( 0 ) },
01090                                 { "parallel miss", 0, box.center + 2.0 * box.scaled_axis( 1 ), box.axis( 2 ) },
01091                                 { "skew miss", 0, box.center + box.dimensions(), box.dimensions() * box.axis( 2 ) } };
01092     const size_t num_def_test = sizeof( default_tests ) / sizeof( default_tests[0] );
01093     tests.insert( tests.begin(), &default_tests[0], &default_tests[0] + num_def_test );
01094     tests.insert( tests.end(), default_files_tests[filename].begin(), default_files_tests[filename].end() );
01095 
01096     OrientedBoxTreeTool::TrvStats stats;
01097 
01098     bool result           = true;
01099     const size_t num_test = tests.size();
01100     for( size_t i = 0; i < num_test; ++i )
01101     {
01102         tests[i].direction.normalize();
01103         if( verbosity > 2 )
01104         {
01105             std::cout << ( 0 == i ? "** Common tests\n" : ( num_def_test == i ? "** File-specific tests\n" : "" ) );
01106             std::cout << "  " << tests[i].description << " " << tests[i].point << " " << tests[i].direction
01107                       << std::endl;
01108         }
01109 
01110         int rit_non_tol_count   = 0;
01111         double rit_non_tol_dist = std::numeric_limits< double >::max();
01112         if( !check_ray_intersect_tris( tool, root_set, tests[i], rit_non_tol_count, rit_non_tol_dist, stats ) )
01113         {
01114             result = false;
01115             continue;
01116         }
01117 
01118         if( !haveSurfTree ) continue;
01119 
01120         int ris_non_tol_count   = 0;
01121         double ris_non_tol_dist = std::numeric_limits< double >::max();
01122         if( !check_ray_intersect_sets( tool, root_set, tests[i], ris_non_tol_count, ris_non_tol_dist, stats ) )
01123         {
01124             result = false;
01125             continue;
01126         }
01127 
01128         if( !rit_non_tol_count && ris_non_tol_count )
01129         {
01130             if( verbosity )
01131                 std::cout << "  ray_intersect_sets returned intersection not found by "
01132                              "ray_intersect_triangles"
01133                           << std::endl;
01134             result = false;
01135             continue;
01136         }
01137         else if( rit_non_tol_count && !ris_non_tol_count )
01138         {
01139             if( verbosity )
01140                 std::cout << "  ray_intersect_sets missed intersection found by ray_intersect_triangles" << std::endl;
01141             result = false;
01142             continue;
01143         }
01144         else if( rit_non_tol_count && ris_non_tol_count && fabs( rit_non_tol_dist - ris_non_tol_dist ) > tolerance )
01145         {
01146             if( verbosity )
01147                 std::cout << "  ray_intersect_sets and ray_intersect_triangles did not find same "
01148                              "closest intersection"
01149                           << std::endl;
01150             result = false;
01151         }
01152     }
01153 
01154     /* Do ray fire for any user-specified rays */
01155 
01156     for( size_t i = 0; i < rays.size(); i += 2 )
01157     {
01158         std::cout << rays[i] << "+" << rays[i + 1] << " : ";
01159 
01160         if( !haveSurfTree )
01161         {
01162             Range leaves;
01163             std::vector< double > intersections;
01164             std::vector< EntityHandle > intersection_facets;
01165             rval = tool.ray_intersect_boxes( leaves, root_set, tolerance, rays[i].array(), rays[i + 1].array(), 0,
01166                                              &stats );
01167             if( MB_SUCCESS != rval )
01168             {
01169                 std::cout << "FAILED" << std::endl;
01170                 result = false;
01171                 continue;
01172             }
01173 
01174             if( !leaves.empty() && write_ray_vtk )
01175             {
01176                 std::string num, name( filename );
01177                 std::stringstream s;
01178                 s << ( i / 2 );
01179                 s >> num;
01180                 name += "-ray";
01181                 name += num;
01182                 name += ".vtk";
01183 
01184                 std::vector< EntityHandle > sets( leaves.size() );
01185                 std::copy( leaves.begin(), leaves.end(), sets.begin() );
01186                 tool.get_moab_instance()->write_mesh( name.c_str(), &sets[0], sets.size() );
01187                 if( verbosity ) std::cout << "(Wrote " << name << ") ";
01188             }
01189 
01190             rval = tool.ray_intersect_triangles( intersections, intersection_facets, leaves, tolerance, rays[i].array(),
01191                                                  rays[i + 1].array(), 0 );
01192             if( MB_SUCCESS != rval )
01193             {
01194                 std::cout << "FAILED" << std::endl;
01195                 result = false;
01196                 continue;
01197             }
01198 
01199             if( intersections.empty() )
01200             {
01201                 std::cout << "(none)" << std::endl;
01202                 continue;
01203             }
01204 
01205             std::cout << intersections[0];
01206             for( unsigned j = 1; j < intersections.size(); ++j )
01207                 std::cout << ", " << intersections[j];
01208             std::cout << std::endl;
01209 
01210             if( !leaves.empty() && write_cubit && verbosity > 2 )
01211             {
01212                 std::cout << "  intersected boxes:";
01213                 for( Range::iterator i2 = leaves.begin(); i2 != leaves.end(); ++i2 )
01214                     std::cout << " " << tool.get_moab_instance()->id_from_handle( *i2 );
01215                 std::cout << std::endl;
01216             }
01217         }
01218         else
01219         {
01220             std::vector< double > intersections;
01221             std::vector< EntityHandle > surfaces;
01222             std::vector< EntityHandle > facets;
01223 
01224             OrientedBoxTreeTool::IntersectSearchWindow search_win;
01225             OrientedBoxTreeTool::IntRegCtxt int_reg_ctxt;
01226             rval = tool.ray_intersect_sets( intersections, surfaces, facets, root_set, tolerance, rays[i].array(),
01227                                             rays[i + 1].array(), search_win, int_reg_ctxt, &stats );
01228 
01229             if( MB_SUCCESS != rval )
01230             {
01231                 if( verbosity ) std::cout << "  Call to OrientedBoxTreeTool::ray_intersect_sets failed." << std::endl;
01232                 result = false;
01233                 continue;
01234             }
01235 
01236             if( !surfaces.empty() && write_ray_vtk )
01237             {
01238                 std::string num, name( filename );
01239                 std::stringstream s;
01240                 s << ( i / 2 );
01241                 s >> num;
01242                 name += "-ray";
01243                 name += num;
01244                 name += ".vtk";
01245 
01246                 tool.get_moab_instance()->write_mesh( name.c_str(), &surfaces[0], surfaces.size() );
01247                 if( verbosity ) std::cout << "(Wrote " << name << ") ";
01248             }
01249 
01250             if( intersections.size() != surfaces.size() )
01251             {
01252                 std::cout << "Mismatched output lists." << std::endl;
01253                 result = false;
01254                 continue;
01255             }
01256 
01257             if( intersections.empty() )
01258             {
01259                 std::cout << "(none)" << std::endl;
01260                 continue;
01261             }
01262 
01263             Tag idtag = tool.get_moab_instance()->globalId_tag();
01264             std::vector< int > ids( surfaces.size() );
01265             rval = tool.get_moab_instance()->tag_get_data( idtag, &surfaces[0], surfaces.size(), &ids[0] );
01266             if( MB_SUCCESS != rval )
01267             {
01268                 std::cout << "NO GLOBAL_ID TAG ON SURFACE." << std::endl;
01269                 continue;
01270             }
01271 
01272             // group by surfaces
01273             std::map< int, double > intmap;
01274             for( unsigned j = 0; j < intersections.size(); ++j )
01275                 intmap[ids[j]] = intersections[j];
01276 
01277             std::map< int, double >::iterator it = intmap.begin();
01278             int prevsurf                         = it->first;
01279             std::cout << "Surf" << it->first << " " << it->second;
01280             for( ++it; it != intmap.end(); ++it )
01281             {
01282                 std::cout << ", ";
01283                 if( it->first != prevsurf )
01284                 {
01285                     prevsurf = it->first;
01286                     std::cout << "Surf" << it->first << " ";
01287                 }
01288                 std::cout << it->second;
01289             }
01290             std::cout << std::endl;
01291         }
01292     }
01293 
01294     if( verbosity > 1 )
01295     {
01296         std::cout << "Traversal statistics for ray fire tests: " << std::endl;
01297         stats.print( std::cout );
01298     }
01299 
01300     return result;
01301 }
01302 
01303 static ErrorCode save_tree( Interface* instance, const char* filename, EntityHandle tree_root )
01304 {
01305     ErrorCode rval;
01306     Tag tag;
01307 
01308     rval = instance->tag_get_handle( "OBB_ROOT", 1, MB_TYPE_HANDLE, tag, MB_TAG_SPARSE | MB_TAG_CREAT );
01309     if( MB_SUCCESS != rval ) return rval;
01310 
01311     const EntityHandle root = 0;
01312     rval                    = instance->tag_set_data( tag, &root, 1, &tree_root );
01313     if( MB_SUCCESS != rval ) return rval;
01314 
01315     return instance->write_mesh( filename );
01316 }
01317 
01318 static ErrorCode tri_coords( Interface* moab, EntityHandle tri, CartVect coords[3] )
01319 {
01320     ErrorCode rval;
01321     const EntityHandle* conn;
01322     int len;
01323 
01324     rval = moab->get_connectivity( tri, conn, len, true );
01325     if( MB_SUCCESS != rval ) return rval;
01326     if( len != 3 ) return MB_FAILURE;
01327     rval = moab->get_coords( conn, 3, coords[0].array() );
01328     return rval;
01329 }
01330 
01331 static ErrorCode closest_point_in_triangles( Interface* moab,
01332                                              const CartVect& to_pos,
01333                                              CartVect& result_pos,
01334                                              EntityHandle& result_tri )
01335 {
01336     ErrorCode rval;
01337 
01338     Range triangles;
01339     rval = moab->get_entities_by_type( 0, MBTRI, triangles );
01340     if( MB_SUCCESS != rval ) return rval;
01341 
01342     if( triangles.empty() ) return MB_FAILURE;
01343 
01344     Range::iterator i = triangles.begin();
01345     CartVect coords[3];
01346     rval = tri_coords( moab, *i, coords );
01347     if( MB_SUCCESS != rval ) return rval;
01348     result_tri = *i;
01349     GeomUtil::closest_location_on_tri( to_pos, coords, result_pos );
01350     CartVect diff            = to_pos - result_pos;
01351     double shortest_dist_sqr = diff % diff;
01352 
01353     for( ++i; i != triangles.end(); ++i )
01354     {
01355         rval = tri_coords( moab, *i, coords );
01356         if( MB_SUCCESS != rval ) return rval;
01357         CartVect pos;
01358         GeomUtil::closest_location_on_tri( to_pos, coords, pos );
01359         diff        = to_pos - pos;
01360         double dsqr = diff % diff;
01361         if( dsqr < shortest_dist_sqr )
01362         {
01363             shortest_dist_sqr = dsqr;
01364             result_pos        = pos;
01365             result_tri        = *i;
01366         }
01367     }
01368 
01369     return MB_SUCCESS;
01370 }
01371 
01372 static bool tri_in_set( Interface* moab, EntityHandle set, EntityHandle tri )
01373 {
01374     Range tris;
01375     ErrorCode rval = moab->get_entities_by_type( set, MBTRI, tris );
01376     if( MB_SUCCESS != rval ) return false;
01377     Range::iterator i = tris.find( tri );
01378     return i != tris.end();
01379 }
01380 
01381 static bool do_closest_point_test( OrientedBoxTreeTool& tool, EntityHandle root_set, bool have_surface_tree )
01382 {
01383     if( verbosity > 1 ) std::cout << "beginning closest point tests" << std::endl;
01384 
01385     ErrorCode rval;
01386     Interface* moab = tool.get_moab_instance();
01387     EntityHandle set;
01388     EntityHandle* set_ptr = have_surface_tree ? &set : 0;
01389     bool result           = true;
01390 
01391     // get root box
01392     OrientedBox box;
01393     rval = tool.box( root_set, box );
01394     if( MB_SUCCESS != rval )
01395     {
01396         if( verbosity ) std::cerr << "Invalid tree in do_closest_point_test\n";
01397         return false;
01398     }
01399 
01400     OrientedBoxTreeTool::TrvStats stats;
01401 
01402     // chose some points to test
01403     CartVect points[] = { box.center + box.scaled_axis( 0 ), box.center + 2 * box.scaled_axis( 1 ),
01404                           box.center + 0.5 * box.scaled_axis( 2 ),
01405                           box.center + -2 * box.scaled_axis( 0 ) + -2 * box.scaled_axis( 1 ) +
01406                               -2 * box.scaled_axis( 2 ),
01407                           CartVect( 100, 100, 100 ) };
01408     const int num_pts = sizeof( points ) / sizeof( points[0] );
01409 
01410     // test each point
01411     for( int i = 0; i < num_pts; ++i )
01412     {
01413         if( verbosity >= 3 ) std::cout << "Evaluating closest point to " << points[i] << std::endl;
01414 
01415         CartVect n_result, t_result;
01416         EntityHandle n_tri = 0, t_tri;
01417 
01418         // find closest point the slow way
01419         rval = closest_point_in_triangles( moab, points[i], n_result, n_tri );
01420         if( MB_SUCCESS != rval )
01421         {
01422             std::cerr << "Internal MOAB error in do_closest_point_test" << std::endl;
01423             result = false;
01424             continue;
01425         }
01426 
01427         // find closest point using tree
01428         rval = tool.closest_to_location( points[i].array(), root_set, t_result.array(), t_tri, set_ptr, &stats );
01429         if( MB_SUCCESS != rval )
01430         {
01431             if( verbosity )
01432                 std::cout << "OrientedBoxTreeTool:: closest_to_location( " << points[i] << " ) FAILED!" << std::endl;
01433             result = false;
01434             continue;
01435         }
01436 
01437         CartVect diff = t_result - n_result;
01438         if( diff.length() > tolerance )
01439         {
01440             if( verbosity )
01441                 std::cout << "Closest point to " << points[i] << " INCORRECT! (" << t_result << " != " << n_result
01442                           << ")" << std::endl;
01443             result = false;
01444             continue;
01445         }
01446 
01447         if( t_tri != n_tri )
01448         {
01449             // if result point is triangle, then OK because
01450             // already tested that it is the same location
01451             // as the expected value.  We just have a case where
01452             // the point was on an edge or vertex.
01453             CartVect coords[3];
01454             CartVect diff1( 1, 1, 1 );
01455             rval = tri_coords( moab, t_tri, coords );
01456             if( MB_SUCCESS == rval )
01457             {
01458                 GeomUtil::closest_location_on_tri( points[i], coords, n_result );
01459                 diff1 = n_result - t_result;
01460             }
01461             if( ( diff1 % diff1 ) > tolerance )
01462             {
01463                 if( verbosity )
01464                     std::cout << "Triangle closest to " << points[i] << " INCORRECT! (" << t_tri << " != " << n_tri
01465                               << ")" << std::endl;
01466                 result = false;
01467                 continue;
01468             }
01469         }
01470 
01471         if( set_ptr && !tri_in_set( moab, *set_ptr, t_tri ) )
01472         {
01473             if( verbosity ) std::cout << "Surface closest to " << points[i] << " INCORRECT!" << std::endl;
01474             result = false;
01475             continue;
01476         }
01477     }
01478 
01479     if( verbosity > 1 )
01480     {
01481         std::cout << "Traversal statistics for closest point tests: " << std::endl;
01482         stats.print( std::cout );
01483     }
01484 
01485     return result;
01486 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines