MOAB: Mesh Oriented datABase  (version 5.2.1)
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 <stdlib.h>
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, OrientedBoxTreeTool* tool_ptr, bool print_errors, std::ostream& str,
00360                    double tol, bool surfs, OrientedBoxTreeTool::Settings s )
00361         : instance( instance_ptr ), tool( tool_ptr ), printing( print_errors ), epsilon( tol ), surfaces( surfs ),
00362           stream( str ), settings( s ), entity_count( 0 ), loose_box_count( 0 ), child_outside_count( 0 ),
00363           entity_outside_count( 0 ), num_entities_outside( 0 ), non_ortho_count( 0 ), error_count( 0 ),
00364           empty_leaf_count( 0 ), non_empty_non_leaf_count( 0 ), entity_invalid_count( 0 ), unsorted_axis_count( 0 ),
00365           non_unit_count( 0 ), duplicate_entity_count( 0 ), bad_outer_radius_count( 0 ), missing_surface_count( 0 ),
00366           multiple_surface_count( 0 ), surface_depth( -1 ), surface_handle( 0 )
00367     {
00368     }
00369 
00370     bool is_valid() const
00371     {
00372         return 0 == loose_box_count + child_outside_count + entity_outside_count + num_entities_outside +
00373                         non_ortho_count + error_count + empty_leaf_count + non_empty_non_leaf_count +
00374                         entity_invalid_count + unsorted_axis_count + non_unit_count + duplicate_entity_count +
00375                         bad_outer_radius_count + missing_surface_count + multiple_surface_count;
00376     }
00377 
00378     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
00379 
00380     virtual ErrorCode leaf( EntityHandle )
00381     {
00382         return MB_SUCCESS;
00383     }
00384 };
00385 
00386 ErrorCode TreeValidator::visit( EntityHandle node, int depth, bool& descend )
00387 {
00388     ErrorCode rval;
00389     descend = true;
00390 
00391     Range contents;
00392     rval = instance->get_entities_by_handle( node, contents );
00393     if( MB_SUCCESS != rval ) return error( node, "Error getting contents of tree node.  Corrupt tree?" );
00394     entity_count += contents.size();
00395 
00396     if( surfaces )
00397     {
00398         // if no longer in subtree for previous surface, clear
00399         if( depth <= surface_depth ) surface_depth = -1;
00400 
00401         EntityHandle surface      = 0;
00402         Range::iterator surf_iter = contents.lower_bound( MBENTITYSET );
00403         if( surf_iter != contents.end() )
00404         {
00405             surface = *surf_iter;
00406             contents.erase( surf_iter );
00407         }
00408 
00409         if( surface )
00410         {
00411             if( surface_depth >= 0 )
00412             {
00413                 ++multiple_surface_count;
00414                 print( node, "Multiple surfaces in encountered in same subtree." );
00415             }
00416             else
00417             {
00418                 surface_depth  = depth;
00419                 surface_handle = surface;
00420             }
00421         }
00422     }
00423 
00424     std::vector< EntityHandle > children;
00425     rval = tool->get_moab_instance()->get_child_meshsets( node, children );
00426     if( MB_SUCCESS != rval || ( !children.empty() && children.size() != 2 ) )
00427         return error( node, "Error getting children.  Corrupt tree?" );
00428 
00429     OrientedBox box;
00430     rval = tool->box( node, box );
00431     if( MB_SUCCESS != rval ) return error( node, "Error getting oriented box from tree node.  Corrupt tree?" );
00432 
00433     if( children.empty() && contents.empty() )
00434     {
00435         ++empty_leaf_count;
00436         print( node, "Empty leaf node.\n" );
00437     }
00438     else if( !children.empty() && !contents.empty() )
00439     {
00440         ++non_empty_non_leaf_count;
00441         print( node, "Non-leaf node is not empty." );
00442     }
00443 
00444     if( surfaces && children.empty() && surface_depth < 0 )
00445     {
00446         ++missing_surface_count;
00447         print( node, "Reached leaf node w/out encountering any surface set." );
00448     }
00449 
00450     double dot_epsilon = epsilon * ( box.axis( 0 ) + box.axis( 1 ) + box.axis( 2 ) ).length();
00451     if( box.axis( 0 ) % box.axis( 1 ) > dot_epsilon || box.axis( 0 ) % box.axis( 2 ) > dot_epsilon ||
00452         box.axis( 1 ) % box.axis( 2 ) > dot_epsilon )
00453     {
00454         ++non_ortho_count;
00455         print( node, "Box axes are not orthogonal" );
00456     }
00457 
00458     if( !children.empty() )
00459     {
00460         for( int i = 0; i < 2; ++i )
00461         {
00462             OrientedBox other_box;
00463             rval = tool->box( children[i], other_box );
00464             if( MB_SUCCESS != rval )
00465                 return error( children[i], " Error getting oriented box from tree node.  Corrupt tree?" );
00466             //      else if (!box.contained( other_box, epsilon )) {
00467             //        ++child_outside_count;
00468             //        print( children[i], "Parent box does not contain child box." );
00469             //        char string[64];
00470             //        sprintf(string, "     Volume ratio is %f", other_box.volume()/box.volume() );
00471             //        print( children [i], string );
00472             //      }
00473             else
00474             {
00475                 double vol_ratio = other_box.volume() / box.volume();
00476                 if( vol_ratio > 2.0 )
00477                 {
00478                     char string[64];
00479                     sprintf( string, "child/parent volume ratio is %f", vol_ratio );
00480                     print( children[i], string );
00481                     sprintf( string, "   child/parent area ratio is %f", other_box.area() / box.area() );
00482                     print( children[i], string );
00483                 }
00484             }
00485         }
00486     }
00487 
00488     bool bad_element        = false;
00489     bool bad_element_handle = false;
00490     bool bad_element_conn   = false;
00491     bool duplicate_element  = false;
00492     int num_outside         = 0;
00493     bool boundary[6]        = { false, false, false, false, false, false };
00494     for( Range::iterator it = contents.begin(); it != contents.end(); ++it )
00495     {
00496         EntityType type = instance->type_from_handle( *it );
00497         int dim         = CN::Dimension( type );
00498         if( dim != 2 )
00499         {
00500             bad_element = true;
00501             continue;
00502         }
00503 
00504         const EntityHandle* conn;
00505         int conn_len;
00506         rval = instance->get_connectivity( *it, conn, conn_len );
00507         if( MB_SUCCESS != rval )
00508         {
00509             bad_element_handle = true;
00510             continue;
00511         }
00512 
00513         std::vector< CartVect > coords( conn_len );
00514         rval = instance->get_coords( conn, conn_len, coords[0].array() );
00515         if( MB_SUCCESS != rval )
00516         {
00517             bad_element_conn = true;
00518             continue;
00519         }
00520 
00521         bool outside = false;
00522         for( std::vector< CartVect >::iterator j = coords.begin(); j != coords.end(); ++j )
00523         {
00524             if( !box.contained( *j, epsilon ) )
00525                 outside = true;
00526             else
00527                 for( int d = 0; d < 3; ++d )
00528                 {
00529 #if MB_ORIENTED_BOX_UNIT_VECTORS
00530                     double n = box.axis( d ) % ( *j - box.center );
00531                     if( fabs( n - box.length[d] ) <= epsilon ) boundary[2 * d] = true;
00532                     if( fabs( n + box.length[d] ) <= epsilon ) boundary[2 * d + 1] = true;
00533 #else
00534                     double ln   = box.axis( d ).length();
00535                     CartVect v1 = *j - box.center - box.axis[d];
00536                     CartVect v2 = *j - box.center + box.axis[d];
00537                     if( fabs( v1 % box.axis[d] ) <= ln * epsilon ) boundary[2 * d] = true;
00538                     if( fabs( v2 % box.axis[d] ) <= ln * epsilon ) boundary[2 * d + 1] = true;
00539 #endif
00540                 }
00541         }
00542         if( outside ) ++num_outside;
00543 
00544         if( !seen.insert( *it ).second )
00545         {
00546             duplicate_element = true;
00547             ++duplicate_entity_count;
00548         }
00549     }
00550 
00551     CartVect alength( box.axis( 0 ).length(), box.axis( 1 ).length(), box.axis( 2 ).length() );
00552 #if MB_ORIENTED_BOX_UNIT_VECTORS
00553     CartVect length = box.length;
00554 #else
00555     CartVect length = alength;
00556 #endif
00557 
00558     if( length[0] > length[1] || length[0] > length[2] || length[1] > length[2] )
00559     {
00560         ++unsorted_axis_count;
00561         print( node, "Box axes are not ordered from shortest to longest." );
00562     }
00563 
00564 #if MB_ORIENTED_BOX_UNIT_VECTORS
00565     if( fabs( alength[0] - 1.0 ) > epsilon || fabs( alength[1] - 1.0 ) > epsilon || fabs( alength[2] - 1.0 ) > epsilon )
00566     {
00567         ++non_unit_count;
00568         print( node, "Box axes are not unit vectors." );
00569     }
00570 #endif
00571 
00572 #if MB_ORIENTED_BOX_OUTER_RADIUS
00573     if( fabs( length.length() - box.radius ) > tolerance )
00574     {
00575         ++bad_outer_radius_count;
00576         print( node, "Box has incorrect outer radius." );
00577     }
00578 #endif
00579 
00580     if( depth + 1 < settings.max_depth && contents.size() > (unsigned)( 4 * settings.max_leaf_entities ) )
00581     {
00582         char string[64];
00583         sprintf( string, "leaf at depth %d with %u entities", depth, (unsigned)contents.size() );
00584         print( node, string );
00585     }
00586 
00587     bool all_boundaries = true;
00588     for( int f = 0; f < 6; ++f )
00589         all_boundaries = all_boundaries && boundary[f];
00590 
00591     if( bad_element )
00592     {
00593         ++entity_invalid_count;
00594         print( node, "Set contained an entity with an inappropriate dimension." );
00595     }
00596     if( bad_element_handle )
00597     {
00598         ++error_count;
00599         print( node, "Error querying face contained in set." );
00600     }
00601     if( bad_element_conn )
00602     {
00603         ++error_count;
00604         print( node, "Error querying connectivity of element." );
00605     }
00606     if( duplicate_element ) { print( node, "Elements occur in multiple leaves of tree." ); }
00607     if( num_outside > 0 )
00608     {
00609         ++entity_outside_count;
00610         num_entities_outside += num_outside;
00611         if( printing )
00612             stream << instance->id_from_handle( node ) << ": " << num_outside << " elements outside box." << std::endl;
00613     }
00614     else if( !all_boundaries && !contents.empty() )
00615     {
00616         ++loose_box_count;
00617         print( node, "Box does not fit contained elements tightly." );
00618     }
00619 
00620     return MB_SUCCESS;
00621 }
00622 
00623 class CubitWriter : public OrientedBoxTreeTool::Op
00624 {
00625   public:
00626     CubitWriter( FILE* file_ptr, OrientedBoxTreeTool* tool_ptr ) : file( file_ptr ), tool( tool_ptr ) {}
00627 
00628     ErrorCode visit( EntityHandle node, int depth, bool& descend );
00629     ErrorCode leaf( EntityHandle )
00630     {
00631         return MB_SUCCESS;
00632     }
00633 
00634   private:
00635     FILE* file;
00636     OrientedBoxTreeTool* tool;
00637 };
00638 
00639 ErrorCode CubitWriter::visit( EntityHandle node, int, bool& descend )
00640 {
00641     descend = true;
00642     OrientedBox box;
00643     ErrorCode rval = tool->box( node, box );
00644     if( rval != MB_SUCCESS ) return rval;
00645 
00646     double sign[] = { -1, 1 };
00647     for( int i = 0; i < 2; ++i )
00648         for( int j = 0; j < 2; ++j )
00649             for( int k = 0; k < 2; ++k )
00650             {
00651 #if MB_ORIENTED_BOX_UNIT_VECTORS
00652                 CartVect corner = box.center + box.length[0] * sign[i] * box.axis( 0 ) +
00653                                   box.length[1] * sign[j] * box.axis( 1 ) + box.length[2] * sign[k] * box.axis( 2 );
00654 #else
00655                 CartVect corner =
00656                     box.center + sign[i] * box.axis( 0 ) + sign[j] * box.axis( 1 ) + sign[k] * box.axis( 2 );
00657 #endif
00658                 fprintf( file, "create vertex %f %f %f\n", corner[0], corner[1], corner[2] );
00659             }
00660     fprintf( file, "#{i=Id(\"vertex\")-7}\n" );
00661     fprintf( file, "create surface vertex {i  } {i+1} {i+3} {i+2}\n" );
00662     fprintf( file, "create surface vertex {i+4} {i+5} {i+7} {i+6}\n" );
00663     fprintf( file, "create surface vertex {i+1} {i+0} {i+4} {i+5}\n" );
00664     fprintf( file, "create surface vertex {i+3} {i+2} {i+6} {i+7}\n" );
00665     fprintf( file, "create surface vertex {i+2} {i+0} {i+4} {i+6}\n" );
00666     fprintf( file, "create surface vertex {i+3} {i+1} {i+5} {i+7}\n" );
00667     fprintf( file, "delete vertex {i}-{i+7}\n" );
00668     fprintf( file, "#{s=Id(\"surface\")-5}\n" );
00669     fprintf( file, "create volume surface {s} {s+1} {s+2} {s+3} {s+4} {s+5} noheal\n" );
00670     int id = tool->get_moab_instance()->id_from_handle( node );
00671     fprintf( file, "volume {Id(\"volume\")} name \"cell%d\"\n", id );
00672 
00673     return MB_SUCCESS;
00674 }
00675 
00676 class VtkWriter : public OrientedBoxTreeTool::Op
00677 {
00678   public:
00679     VtkWriter( std::string base_name, Interface* interface ) : baseName( base_name ), instance( interface ) {}
00680 
00681     ErrorCode visit( EntityHandle node, int depth, bool& descend );
00682 
00683     ErrorCode leaf( EntityHandle /*node*/ )
00684     {
00685         return MB_SUCCESS;
00686     }
00687 
00688   private:
00689     std::string baseName;
00690     Interface* instance;
00691 };
00692 
00693 ErrorCode VtkWriter::visit( EntityHandle node, int, bool& descend )
00694 {
00695     descend = true;
00696 
00697     ErrorCode rval;
00698     int count;
00699     rval = instance->get_number_entities_by_handle( node, count );
00700     if( MB_SUCCESS != rval || 0 == count ) return rval;
00701 
00702     int id = instance->id_from_handle( node );
00703     char id_str[32];
00704     sprintf( id_str, "%d", id );
00705 
00706     std::string file_name( baseName );
00707     file_name += ".";
00708     file_name += id_str;
00709     file_name += ".vtk";
00710 
00711     return instance->write_mesh( file_name.c_str(), &node, 1 );
00712 }
00713 
00714 static bool do_ray_fire_test( OrientedBoxTreeTool& tool, EntityHandle root_set, const char* filename,
00715                               bool have_surface_tree );
00716 
00717 static bool do_closest_point_test( OrientedBoxTreeTool& tool, EntityHandle root_set, bool have_surface_tree );
00718 
00719 static ErrorCode save_tree( Interface* instance, const char* filename, EntityHandle tree_root );
00720 
00721 static bool do_file( const char* filename )
00722 {
00723     ErrorCode rval;
00724     Core instance;
00725     Interface* const iface = &instance;
00726     OrientedBoxTreeTool tool( iface );
00727     bool haveSurfTree = false;
00728 
00729     if( verbosity ) std::cout << filename << std::endl << "------" << std::endl;
00730 
00731     rval = iface->load_mesh( filename );
00732     if( MB_SUCCESS != rval )
00733     {
00734         if( verbosity ) std::cout << "Failed to read file: \"" << filename << '"' << std::endl;
00735         return false;
00736     }
00737 
00738     // IF building from surfaces, get surfaces.
00739     // If AUTO and less than two surfaces, don't build from surfaces.
00740     Range surfaces;
00741     if( surfTree != DISABLE )
00742     {
00743         Tag surftag;
00744         rval = iface->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, surftag );
00745         if( MB_SUCCESS == rval )
00746         {
00747             int dim                 = 2;
00748             const void* tagvalues[] = { &dim };
00749             rval = iface->get_entities_by_type_and_tag( 0, MBENTITYSET, &surftag, tagvalues, 1, surfaces );
00750             if( MB_SUCCESS != rval && MB_ENTITY_NOT_FOUND != rval ) return false;
00751         }
00752         else if( MB_TAG_NOT_FOUND != rval )
00753             return false;
00754 
00755         if( ENABLE == surfTree && surfaces.empty() )
00756         {
00757             std::cerr << "No Surfaces found." << std::endl;
00758             return false;
00759         }
00760 
00761         haveSurfTree = ( ENABLE == surfTree ) || ( surfaces.size() > 1 );
00762     }
00763 
00764     EntityHandle root;
00765     Range entities;
00766     if( !haveSurfTree )
00767     {
00768         rval = iface->get_entities_by_dimension( 0, 2, entities );
00769         if( MB_SUCCESS != rval )
00770         {
00771             std::cerr << "get_entities_by_dimension( 2 ) failed." << std::endl;
00772             return false;
00773         }
00774 
00775         if( entities.empty() )
00776         {
00777             if( verbosity ) std::cout << "File \"" << filename << "\" contains no 2D elements" << std::endl;
00778             return false;
00779         }
00780 
00781         if( verbosity ) std::cout << "Building tree from " << entities.size() << " 2D elements" << std::endl;
00782 
00783         rval = tool.build( entities, root, &settings );
00784         if( MB_SUCCESS != rval )
00785         {
00786             if( verbosity ) std::cout << "Failed to build tree." << std::endl;
00787             return false;
00788         }
00789     }
00790     else
00791     {
00792 
00793         if( verbosity ) std::cout << "Building tree from " << surfaces.size() << " surfaces" << std::endl;
00794 
00795         // Build subtree for each surface, get list of all entities to use later
00796         Range surf_trees, surf_tris;
00797         EntityHandle surf_root;
00798         for( Range::iterator s = surfaces.begin(); s != surfaces.end(); ++s )
00799         {
00800             surf_tris.clear();
00801             rval = iface->get_entities_by_dimension( *s, 2, surf_tris );
00802             if( MB_SUCCESS != rval ) return false;
00803             rval = tool.build( surf_tris, surf_root, &settings );
00804             if( MB_SUCCESS != rval )
00805             {
00806                 if( verbosity ) std::cout << "Failed to build tree for surface." << std::endl;
00807                 return false;
00808             }
00809             surf_trees.insert( surf_root );
00810             entities.merge( surf_tris );
00811             rval = iface->add_entities( surf_root, &*s, 1 );
00812             if( MB_SUCCESS != rval ) return false;
00813         }
00814 
00815         rval = tool.join_trees( surf_trees, root, &settings );
00816         if( MB_SUCCESS != rval )
00817         {
00818             if( verbosity ) std::cout << "Failed to build tree." << std::endl;
00819             return false;
00820         }
00821 
00822         if( verbosity )
00823             std::cout << "Built tree from " << surfaces.size() << " surfaces"
00824                       << " (" << entities.size() - surfaces.size() << " elements)" << std::endl;
00825 
00826         entities.merge( surfaces );
00827     }
00828 
00829     if( write_cubit )
00830     {
00831         std::string name = filename;
00832         name += ".boxes.jou";
00833         FILE* ptr = fopen( name.c_str(), "w+" );
00834         if( !ptr ) { perror( name.c_str() ); }
00835         else
00836         {
00837             if( verbosity ) std::cout << "Writing: " << name << std::endl;
00838             fprintf( ptr, "graphics off\n" );
00839             CubitWriter op( ptr, &tool );
00840             tool.preorder_traverse( root, op );
00841             fprintf( ptr, "graphics on\n" );
00842             fclose( ptr );
00843         }
00844     }
00845 
00846     if( write_vtk )
00847     {
00848         VtkWriter op( filename, iface );
00849         if( verbosity )
00850             std::cout << "Writing leaf contents as : " << filename << ".xxx.vtk where 'xxx' is the set id" << std::endl;
00851         tool.preorder_traverse( root, op );
00852     }
00853 
00854     bool print_errors = false, print_contents = false;
00855     switch( verbosity )
00856     {
00857         default:
00858             print_contents = true;
00859         case 4:
00860             tool.print( root, std::cout, print_contents );
00861         case 3:
00862             print_errors = true;
00863         case 2:
00864             rval = tool.stats( root, std::cout );
00865             if( MB_SUCCESS != rval ) std::cout << "****** Failed to get tree statistics ******" << std::endl;
00866         case 1:
00867         case 0:;
00868     }
00869 
00870     TreeValidator op( iface, &tool, print_errors, std::cout, tolerance, haveSurfTree, settings );
00871     rval        = tool.preorder_traverse( root, op );
00872     bool result = op.is_valid();
00873     if( MB_SUCCESS != rval )
00874     {
00875         result = false;
00876         if( verbosity ) std::cout << "Errors traversing tree.  Corrupt tree?" << std::endl;
00877     }
00878 
00879     bool missing = ( op.entity_count != entities.size() );
00880     if( missing ) result = false;
00881 
00882     if( verbosity )
00883     {
00884         if( result )
00885             std::cout << std::endl << "No errors detected." << std::endl;
00886         else
00887             std::cout << std::endl << "*********************** ERROR SUMMARY **********************" << std::endl;
00888         if( op.child_outside_count )
00889             std::cout << "* " << op.child_outside_count << " child boxes not contained in parent." << std::endl;
00890         if( op.entity_outside_count )
00891             std::cout << "* " << op.entity_outside_count << " nodes containing entities outside of box." << std::endl;
00892         if( op.num_entities_outside )
00893             std::cout << "* " << op.num_entities_outside << " entities outside boxes." << std::endl;
00894         if( op.empty_leaf_count ) std::cout << "* " << op.empty_leaf_count << " empty leaf nodes." << std::endl;
00895         if( op.non_empty_non_leaf_count )
00896             std::cout << "* " << op.non_empty_non_leaf_count << " non-leaf nodes containing entities." << std::endl;
00897         if( op.duplicate_entity_count )
00898             std::cout << "* " << op.duplicate_entity_count << " duplicate entities in leaves." << std::endl;
00899         if( op.missing_surface_count )
00900             std::cout << "* " << op.missing_surface_count << " leaves outside surface subtrees." << std::endl;
00901         if( op.multiple_surface_count )
00902             std::cout << "* " << op.multiple_surface_count << " surfaces within other surface subtrees." << std::endl;
00903         if( op.non_ortho_count )
00904             std::cout << "* " << op.non_ortho_count << " boxes with non-orthononal axes." << std::endl;
00905         if( op.non_unit_count ) std::cout << "* " << op.non_unit_count << " boxes with non-unit axes." << std::endl;
00906         if( op.bad_outer_radius_count )
00907             std::cout << "* " << op.bad_outer_radius_count << " boxes incorrect outer radius." << std::endl;
00908         if( op.unsorted_axis_count )
00909             std::cout << "* " << op.unsorted_axis_count << " boxes axes in unsorted order." << std::endl;
00910         if( op.loose_box_count )
00911             std::cout << "* " << op.loose_box_count << " boxes that do not fit the contained entities tightly."
00912                       << std::endl;
00913         if( op.error_count + op.entity_invalid_count )
00914             std::cout << "* " << op.error_count + op.entity_invalid_count << " other errors while traversing tree."
00915                       << std::endl;
00916         if( missing )
00917             std::cout << "* tree built from " << entities.size() << " entities contains " << op.entity_count
00918                       << " entities." << std::endl;
00919         if( !result ) std::cout << "************************************************************" << std::endl;
00920     }
00921 
00922     if( result && save_file_name )
00923     {
00924         if( MB_SUCCESS == save_tree( iface, save_file_name, root ) )
00925             std::cerr << "Wrote '" << save_file_name << "'" << std::endl;
00926         else
00927             std::cerr << "FAILED TO WRITE '" << save_file_name << "'" << std::endl;
00928     }
00929 
00930     if( !do_ray_fire_test( tool, root, filename, haveSurfTree ) )
00931     {
00932         if( verbosity ) std::cout << "Ray fire test failed." << std::endl;
00933         result = false;
00934     }
00935 
00936     if( !do_closest_point_test( tool, root, haveSurfTree ) )
00937     {
00938         if( verbosity ) std::cout << "Closest point test failed" << std::endl;
00939         result = false;
00940     }
00941 
00942     rval = tool.delete_tree( root );
00943     if( MB_SUCCESS != rval )
00944     {
00945         if( verbosity ) std::cout << "delete_tree failed." << std::endl;
00946         result = false;
00947     }
00948 
00949     return result;
00950 }
00951 
00952 static void count_non_tol( std::vector< double > intersections, int& non_tol_count, double& non_tol_dist )
00953 {
00954 
00955     for( size_t i = 0; i < intersections.size(); ++i )
00956     {
00957         if( intersections[i] > tolerance )
00958         {
00959             ++non_tol_count;
00960             if( intersections[i] < non_tol_dist ) non_tol_dist = intersections[i];
00961         }
00962     }
00963 }
00964 
00965 static bool check_ray_intersect_tris( OrientedBoxTreeTool& tool, EntityHandle root_set, RayTest& test,
00966                                       int& non_tol_count, double& non_tol_dist, OrientedBoxTreeTool::TrvStats& stats )
00967 {
00968     ErrorCode rval;
00969     bool result = true;
00970 
00971     non_tol_dist  = std::numeric_limits< double >::max();
00972     non_tol_count = 0;
00973 
00974     std::vector< double > intersections;
00975     std::vector< EntityHandle > facets;
00976     rval = tool.ray_intersect_triangles( intersections, facets, root_set, tolerance, test.point.array(),
00977                                          test.direction.array(), 0, &stats );
00978     if( MB_SUCCESS != rval )
00979     {
00980         if( verbosity ) std::cout << "  Call to OrientedBoxTreeTool::ray_intersect_triangles failed." << std::endl;
00981         result = false;
00982     }
00983     else
00984     {
00985 
00986         if( intersections.size() != test.expected_hits )
00987         {
00988             if( verbosity > 2 )
00989                 std::cout << "  Expected " << test.expected_hits << " and got " << intersections.size()
00990                           << " hits for ray fire of " << test.description << std::endl;
00991             if( verbosity > 3 )
00992             {
00993                 for( unsigned j = 0; j < intersections.size(); ++j )
00994                     std::cout << "  " << intersections[j];
00995                 std::cout << std::endl;
00996             }
00997             result = false;
00998         }
00999 
01000         count_non_tol( intersections, non_tol_count, non_tol_dist );
01001     }
01002     return result;
01003 }
01004 
01005 static bool check_ray_intersect_sets( OrientedBoxTreeTool& tool, EntityHandle root_set, RayTest& test,
01006                                       int& non_tol_count, double& non_tol_dist, OrientedBoxTreeTool::TrvStats& stats )
01007 {
01008 
01009     ErrorCode rval;
01010     bool result = true;
01011 
01012     non_tol_dist  = std::numeric_limits< double >::max();
01013     non_tol_count = 0;
01014 
01015     const int NUM_NON_TOL_INT = 1;
01016 
01017     std::vector< double > intersections;
01018     std::vector< EntityHandle > surfs;
01019     std::vector< EntityHandle > facets;
01020 
01021     OrientedBoxTreeTool::IntersectSearchWindow search_win;
01022     OrientedBoxTreeTool::IntRegCtxt int_reg_ctxt;
01023     rval = tool.ray_intersect_sets( intersections, surfs, facets, root_set, tolerance, test.point.array(),
01024                                     test.direction.array(), search_win, int_reg_ctxt, &stats );
01025 
01026     if( MB_SUCCESS != rval )
01027     {
01028         if( verbosity ) std::cout << "  Call to OrientedBoxTreeTool::ray_intersect_sets failed." << std::endl;
01029         result = false;
01030     }
01031     else
01032     {
01033 
01034         if( surfs.size() != intersections.size() )
01035         {
01036             if( verbosity ) std::cout << "  ray_intersect_sets did not return sets for all intersections." << std::endl;
01037             result = false;
01038         }
01039 
01040         count_non_tol( intersections, non_tol_count, non_tol_dist );
01041 
01042         if( non_tol_count > NUM_NON_TOL_INT )
01043         {
01044             if( verbosity )
01045                 std::cout << "  Requested " << NUM_NON_TOL_INT << "intersections "
01046                           << "  beyond tolerance.  Got " << non_tol_count << std::endl;
01047             result = false;
01048         }
01049     }
01050 
01051     return result;
01052 }
01053 
01054 static bool do_ray_fire_test( OrientedBoxTreeTool& tool, EntityHandle root_set, const char* filename,
01055                               bool haveSurfTree )
01056 {
01057     if( verbosity > 1 ) std::cout << "beginning ray fire tests" << std::endl;
01058 
01059     OrientedBox box;
01060     ErrorCode rval = tool.box( root_set, box );
01061     if( MB_SUCCESS != rval )
01062     {
01063         if( verbosity ) std::cerr << "Error getting box for tree root set" << std::endl;
01064         return false;
01065     }
01066 
01067     /* Do standard ray fire tests */
01068     std::cout << box << std::endl;
01069 
01070     CartVect origin( 0., 0., 0. );
01071     CartVect unitDiag( 1., 1., 1. );
01072     std::vector< RayTest > tests;
01073     RayTest default_tests[] = { { "large axis through box", 2, box.center - 1.2 * box.scaled_axis( 2 ), box.axis( 2 ) },
01074                                 { "small axis through box", 2, box.center - 1.2 * box.scaled_axis( 0 ), box.axis( 0 ) },
01075                                 { "parallel miss", 0, box.center + 2.0 * box.scaled_axis( 1 ), box.axis( 2 ) },
01076                                 { "skew miss", 0, box.center + box.dimensions(), box.dimensions() * box.axis( 2 ) } };
01077     const size_t num_def_test = sizeof( default_tests ) / sizeof( default_tests[0] );
01078     tests.insert( tests.begin(), &default_tests[0], &default_tests[0] + num_def_test );
01079     tests.insert( tests.end(), default_files_tests[filename].begin(), default_files_tests[filename].end() );
01080 
01081     OrientedBoxTreeTool::TrvStats stats;
01082 
01083     bool result           = true;
01084     const size_t num_test = tests.size();
01085     for( size_t i = 0; i < num_test; ++i )
01086     {
01087         tests[i].direction.normalize();
01088         if( verbosity > 2 )
01089         {
01090             std::cout << ( 0 == i ? "** Common tests\n" : ( num_def_test == i ? "** File-specific tests\n" : "" ) );
01091             std::cout << "  " << tests[i].description << " " << tests[i].point << " " << tests[i].direction
01092                       << std::endl;
01093         }
01094 
01095         int rit_non_tol_count   = 0;
01096         double rit_non_tol_dist = std::numeric_limits< double >::max();
01097         if( !check_ray_intersect_tris( tool, root_set, tests[i], rit_non_tol_count, rit_non_tol_dist, stats ) )
01098         {
01099             result = false;
01100             continue;
01101         }
01102 
01103         if( !haveSurfTree ) continue;
01104 
01105         int ris_non_tol_count   = 0;
01106         double ris_non_tol_dist = std::numeric_limits< double >::max();
01107         if( !check_ray_intersect_sets( tool, root_set, tests[i], ris_non_tol_count, ris_non_tol_dist, stats ) )
01108         {
01109             result = false;
01110             continue;
01111         }
01112 
01113         if( !rit_non_tol_count && ris_non_tol_count )
01114         {
01115             if( verbosity )
01116                 std::cout << "  ray_intersect_sets returned intersection not found by "
01117                              "ray_intersect_triangles"
01118                           << std::endl;
01119             result = false;
01120             continue;
01121         }
01122         else if( rit_non_tol_count && !ris_non_tol_count )
01123         {
01124             if( verbosity )
01125                 std::cout << "  ray_intersect_sets missed intersection found by ray_intersect_triangles" << std::endl;
01126             result = false;
01127             continue;
01128         }
01129         else if( rit_non_tol_count && ris_non_tol_count && fabs( rit_non_tol_dist - ris_non_tol_dist ) > tolerance )
01130         {
01131             if( verbosity )
01132                 std::cout << "  ray_intersect_sets and ray_intersect_triangles did not find same "
01133                              "closest intersection"
01134                           << std::endl;
01135             result = false;
01136         }
01137     }
01138 
01139     /* Do ray fire for any user-specified rays */
01140 
01141     for( size_t i = 0; i < rays.size(); i += 2 )
01142     {
01143         std::cout << rays[i] << "+" << rays[i + 1] << " : ";
01144 
01145         if( !haveSurfTree )
01146         {
01147             Range leaves;
01148             std::vector< double > intersections;
01149             std::vector< EntityHandle > intersection_facets;
01150             rval = tool.ray_intersect_boxes( leaves, root_set, tolerance, rays[i].array(), rays[i + 1].array(), 0,
01151                                              &stats );
01152             if( MB_SUCCESS != rval )
01153             {
01154                 std::cout << "FAILED" << std::endl;
01155                 result = false;
01156                 continue;
01157             }
01158 
01159             if( !leaves.empty() && write_ray_vtk )
01160             {
01161                 std::string num, name( filename );
01162                 std::stringstream s;
01163                 s << ( i / 2 );
01164                 s >> num;
01165                 name += "-ray";
01166                 name += num;
01167                 name += ".vtk";
01168 
01169                 std::vector< EntityHandle > sets( leaves.size() );
01170                 std::copy( leaves.begin(), leaves.end(), sets.begin() );
01171                 tool.get_moab_instance()->write_mesh( name.c_str(), &sets[0], sets.size() );
01172                 if( verbosity ) std::cout << "(Wrote " << name << ") ";
01173             }
01174 
01175             rval = tool.ray_intersect_triangles( intersections, intersection_facets, leaves, tolerance, rays[i].array(),
01176                                                  rays[i + 1].array(), 0 );
01177             if( MB_SUCCESS != rval )
01178             {
01179                 std::cout << "FAILED" << std::endl;
01180                 result = false;
01181                 continue;
01182             }
01183 
01184             if( intersections.empty() )
01185             {
01186                 std::cout << "(none)" << std::endl;
01187                 continue;
01188             }
01189 
01190             std::cout << intersections[0];
01191             for( unsigned j = 1; j < intersections.size(); ++j )
01192                 std::cout << ", " << intersections[j];
01193             std::cout << std::endl;
01194 
01195             if( !leaves.empty() && write_cubit && verbosity > 2 )
01196             {
01197                 std::cout << "  intersected boxes:";
01198                 for( Range::iterator i2 = leaves.begin(); i2 != leaves.end(); ++i2 )
01199                     std::cout << " " << tool.get_moab_instance()->id_from_handle( *i2 );
01200                 std::cout << std::endl;
01201             }
01202         }
01203         else
01204         {
01205             std::vector< double > intersections;
01206             std::vector< EntityHandle > surfaces;
01207             std::vector< EntityHandle > facets;
01208 
01209             OrientedBoxTreeTool::IntersectSearchWindow search_win;
01210             OrientedBoxTreeTool::IntRegCtxt int_reg_ctxt;
01211             rval = tool.ray_intersect_sets( intersections, surfaces, facets, root_set, tolerance, rays[i].array(),
01212                                             rays[i + 1].array(), search_win, int_reg_ctxt, &stats );
01213 
01214             if( MB_SUCCESS != rval )
01215             {
01216                 if( verbosity ) std::cout << "  Call to OrientedBoxTreeTool::ray_intersect_sets failed." << std::endl;
01217                 result = false;
01218                 continue;
01219             }
01220 
01221             if( !surfaces.empty() && write_ray_vtk )
01222             {
01223                 std::string num, name( filename );
01224                 std::stringstream s;
01225                 s << ( i / 2 );
01226                 s >> num;
01227                 name += "-ray";
01228                 name += num;
01229                 name += ".vtk";
01230 
01231                 tool.get_moab_instance()->write_mesh( name.c_str(), &surfaces[0], surfaces.size() );
01232                 if( verbosity ) std::cout << "(Wrote " << name << ") ";
01233             }
01234 
01235             if( intersections.size() != surfaces.size() )
01236             {
01237                 std::cout << "Mismatched output lists." << std::endl;
01238                 result = false;
01239                 continue;
01240             }
01241 
01242             if( intersections.empty() )
01243             {
01244                 std::cout << "(none)" << std::endl;
01245                 continue;
01246             }
01247 
01248             Tag idtag = tool.get_moab_instance()->globalId_tag();
01249             std::vector< int > ids( surfaces.size() );
01250             rval = tool.get_moab_instance()->tag_get_data( idtag, &surfaces[0], surfaces.size(), &ids[0] );
01251             if( MB_SUCCESS != rval )
01252             {
01253                 std::cout << "NO GLOBAL_ID TAG ON SURFACE." << std::endl;
01254                 continue;
01255             }
01256 
01257             // group by surfaces
01258             std::map< int, double > intmap;
01259             for( unsigned j = 0; j < intersections.size(); ++j )
01260                 intmap[ids[j]] = intersections[j];
01261 
01262             std::map< int, double >::iterator it = intmap.begin();
01263             int prevsurf                         = it->first;
01264             std::cout << "Surf" << it->first << " " << it->second;
01265             for( ++it; it != intmap.end(); ++it )
01266             {
01267                 std::cout << ", ";
01268                 if( it->first != prevsurf )
01269                 {
01270                     prevsurf = it->first;
01271                     std::cout << "Surf" << it->first << " ";
01272                 }
01273                 std::cout << it->second;
01274             }
01275             std::cout << std::endl;
01276         }
01277     }
01278 
01279     if( verbosity > 1 )
01280     {
01281         std::cout << "Traversal statistics for ray fire tests: " << std::endl;
01282         stats.print( std::cout );
01283     }
01284 
01285     return result;
01286 }
01287 
01288 static ErrorCode save_tree( Interface* instance, const char* filename, EntityHandle tree_root )
01289 {
01290     ErrorCode rval;
01291     Tag tag;
01292 
01293     rval = instance->tag_get_handle( "OBB_ROOT", 1, MB_TYPE_HANDLE, tag, MB_TAG_SPARSE | MB_TAG_CREAT );
01294     if( MB_SUCCESS != rval ) return rval;
01295 
01296     const EntityHandle root = 0;
01297     rval                    = instance->tag_set_data( tag, &root, 1, &tree_root );
01298     if( MB_SUCCESS != rval ) return rval;
01299 
01300     return instance->write_mesh( filename );
01301 }
01302 
01303 static ErrorCode tri_coords( Interface* moab, EntityHandle tri, CartVect coords[3] )
01304 {
01305     ErrorCode rval;
01306     const EntityHandle* conn;
01307     int len;
01308 
01309     rval = moab->get_connectivity( tri, conn, len, true );
01310     if( MB_SUCCESS != rval ) return rval;
01311     if( len != 3 ) return MB_FAILURE;
01312     rval = moab->get_coords( conn, 3, coords[0].array() );
01313     return rval;
01314 }
01315 
01316 static ErrorCode closest_point_in_triangles( Interface* moab, const CartVect& to_pos, CartVect& result_pos,
01317                                              EntityHandle& result_tri )
01318 {
01319     ErrorCode rval;
01320 
01321     Range triangles;
01322     rval = moab->get_entities_by_type( 0, MBTRI, triangles );
01323     if( MB_SUCCESS != rval ) return rval;
01324 
01325     if( triangles.empty() ) return MB_FAILURE;
01326 
01327     Range::iterator i = triangles.begin();
01328     CartVect coords[3];
01329     rval = tri_coords( moab, *i, coords );
01330     if( MB_SUCCESS != rval ) return rval;
01331     result_tri = *i;
01332     GeomUtil::closest_location_on_tri( to_pos, coords, result_pos );
01333     CartVect diff            = to_pos - result_pos;
01334     double shortest_dist_sqr = diff % diff;
01335 
01336     for( ++i; i != triangles.end(); ++i )
01337     {
01338         rval = tri_coords( moab, *i, coords );
01339         if( MB_SUCCESS != rval ) return rval;
01340         CartVect pos;
01341         GeomUtil::closest_location_on_tri( to_pos, coords, pos );
01342         diff        = to_pos - pos;
01343         double dsqr = diff % diff;
01344         if( dsqr < shortest_dist_sqr )
01345         {
01346             shortest_dist_sqr = dsqr;
01347             result_pos        = pos;
01348             result_tri        = *i;
01349         }
01350     }
01351 
01352     return MB_SUCCESS;
01353 }
01354 
01355 static bool tri_in_set( Interface* moab, EntityHandle set, EntityHandle tri )
01356 {
01357     Range tris;
01358     ErrorCode rval = moab->get_entities_by_type( set, MBTRI, tris );
01359     if( MB_SUCCESS != rval ) return false;
01360     Range::iterator i = tris.find( tri );
01361     return i != tris.end();
01362 }
01363 
01364 static bool do_closest_point_test( OrientedBoxTreeTool& tool, EntityHandle root_set, bool have_surface_tree )
01365 {
01366     if( verbosity > 1 ) std::cout << "beginning closest point tests" << std::endl;
01367 
01368     ErrorCode rval;
01369     Interface* moab = tool.get_moab_instance();
01370     EntityHandle set;
01371     EntityHandle* set_ptr = have_surface_tree ? &set : 0;
01372     bool result           = true;
01373 
01374     // get root box
01375     OrientedBox box;
01376     rval = tool.box( root_set, box );
01377     if( MB_SUCCESS != rval )
01378     {
01379         if( verbosity ) std::cerr << "Invalid tree in do_closest_point_test\n";
01380         return false;
01381     }
01382 
01383     OrientedBoxTreeTool::TrvStats stats;
01384 
01385     // chose some points to test
01386     CartVect points[] = { box.center + box.scaled_axis( 0 ), box.center + 2 * box.scaled_axis( 1 ),
01387                           box.center + 0.5 * box.scaled_axis( 2 ),
01388                           box.center + -2 * box.scaled_axis( 0 ) + -2 * box.scaled_axis( 1 ) +
01389                               -2 * box.scaled_axis( 2 ),
01390                           CartVect( 100, 100, 100 ) };
01391     const int num_pts = sizeof( points ) / sizeof( points[0] );
01392 
01393     // test each point
01394     for( int i = 0; i < num_pts; ++i )
01395     {
01396         if( verbosity >= 3 ) std::cout << "Evaluating closest point to " << points[i] << std::endl;
01397 
01398         CartVect n_result, t_result;
01399         EntityHandle n_tri = 0, t_tri;
01400 
01401         // find closest point the slow way
01402         rval = closest_point_in_triangles( moab, points[i], n_result, n_tri );
01403         if( MB_SUCCESS != rval )
01404         {
01405             std::cerr << "Internal MOAB error in do_closest_point_test" << std::endl;
01406             result = false;
01407             continue;
01408         }
01409 
01410         // find closest point using tree
01411         rval = tool.closest_to_location( points[i].array(), root_set, t_result.array(), t_tri, set_ptr, &stats );
01412         if( MB_SUCCESS != rval )
01413         {
01414             if( verbosity )
01415                 std::cout << "OrientedBoxTreeTool:: closest_to_location( " << points[i] << " ) FAILED!" << std::endl;
01416             result = false;
01417             continue;
01418         }
01419 
01420         CartVect diff = t_result - n_result;
01421         if( diff.length() > tolerance )
01422         {
01423             if( verbosity )
01424                 std::cout << "Closest point to " << points[i] << " INCORRECT! (" << t_result << " != " << n_result
01425                           << ")" << std::endl;
01426             result = false;
01427             continue;
01428         }
01429 
01430         if( t_tri != n_tri )
01431         {
01432             // if result point is triangle, then OK because
01433             // already tested that it is the same location
01434             // as the expected value.  We just have a case where
01435             // the point was on an edge or vertex.
01436             CartVect coords[3];
01437             CartVect diff1( 1, 1, 1 );
01438             rval = tri_coords( moab, t_tri, coords );
01439             if( MB_SUCCESS == rval )
01440             {
01441                 GeomUtil::closest_location_on_tri( points[i], coords, n_result );
01442                 diff1 = n_result - t_result;
01443             }
01444             if( ( diff1 % diff1 ) > tolerance )
01445             {
01446                 if( verbosity )
01447                     std::cout << "Triangle closest to " << points[i] << " INCORRECT! (" << t_tri << " != " << n_tri
01448                               << ")" << std::endl;
01449                 result = false;
01450                 continue;
01451             }
01452         }
01453 
01454         if( set_ptr && !tri_in_set( moab, *set_ptr, t_tri ) )
01455         {
01456             if( verbosity ) std::cout << "Surface closest to " << points[i] << " INCORRECT!" << std::endl;
01457             result = false;
01458             continue;
01459         }
01460     }
01461 
01462     if( verbosity > 1 )
01463     {
01464         std::cout << "Traversal statistics for closest point tests: " << std::endl;
01465         stats.print( std::cout );
01466     }
01467 
01468     return result;
01469 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines