MOAB: Mesh Oriented datABase  (version 5.5.0)
obb_test.cpp
Go to the documentation of this file.
1 #include "moab/Core.hpp"
2 #include "moab/Range.hpp"
3 #include "moab/OrientedBoxTreeTool.hpp"
4 #include "moab/OrientedBox.hpp"
5 #include "MBTagConventions.hpp"
6 #include "moab/GeomUtil.hpp"
7 #include "moab/CN.hpp"
8 
9 #ifdef MOAB_HAVE_MPI
10 #include "moab_mpi.h"
11 #endif
12 
13 #include "../TestUtil.hpp"
14 #include <iostream>
15 #include <sstream>
16 #include <cstdlib>
17 #include <limits>
18 #include <cstdio>
19 #include <set>
20 #include <map>
21 
22 using namespace moab;
23 
24 struct RayTest
25 {
26  const char* description;
27  unsigned expected_hits;
28  CartVect point, direction;
29 };
30 
31 static const char* NAME = "obb_test";
32 
33 std::map< std::string, std::vector< RayTest > > default_files_tests;
34 typedef std::map< std::string, std::vector< RayTest > >::iterator ray_test_itr;
35 
36 void initialize_default_files()
37 {
38  size_t num_tests;
39  std::vector< RayTest > tests;
40  std::string file = STRINGIFY( MESHDIR ) "/3k-tri-sphere.vtk";
41  RayTest set1[] = { { "triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 99.8792, -5, 0.121729 ) },
42  { "triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( 4.99167, 0, 99.7502 ) },
43  { "triangle node ", 6, CartVect( 0, 0, 0 ), CartVect( 0, 0, 100 ) } };
44 
45  num_tests = sizeof( set1 ) / sizeof( set1[0] );
46  tests.insert( tests.begin(), &set1[0], &set1[num_tests] );
47  default_files_tests[file] = tests;
48  tests.clear();
49 
50 #ifdef MOAB_HAVE_HDF5
51  file = STRINGIFY( MESHDIR ) "/3k-tri-cube.h5m";
52 #else
53  file = STRINGIFY( MESHDIR ) "/3k-tri-cube.vtk";
54 #endif
55 
56  RayTest set2[] = { { "interior triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 0, 0, 5 ) },
57  { "interior triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( -0.25, -0.05, 5 ) },
58  { "interior triangle node ", 5, CartVect( 0, 0, 0 ), CartVect( -0.3, -0.3, 5 ) },
59  { "edge triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 5, 2.9, 2.9 ) },
60  { "edge triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( 5, 5, 2.9 ) },
61  { "edge triangle node ", 6, CartVect( 0, 0, 0 ), CartVect( 5, 5, 3 ) },
62  { "corner triangle interior ", 1, CartVect( 0, 0, 0 ), CartVect( 5, 4.9, 4.9 ) },
63  { "corner triangle edge ", 2, CartVect( 0, 0, 0 ), CartVect( 5, 5, 4.9 ) },
64  { "corner triangle node ", 3, CartVect( 0, 0, 0 ), CartVect( 5, 5, 5 ) } };
65 
66  num_tests = sizeof( set2 ) / sizeof( set2[0] );
67  tests.insert( tests.begin(), &set2[0], &set2[num_tests] );
68  default_files_tests[file] = tests;
69  tests.clear();
70 }
71 
72 // Another possible file
73 // STRINGIFY(MESHDIR) "/4k-tri-plane.vtk",
74 
75 static void usage( const char* error, const char* opt )
76 {
77  const char* default_message = "Invalid option";
78  if( opt && !error ) error = default_message;
79 
80  std::ostream& str = error ? std::cerr : std::cout;
81  if( error )
82  {
83  str << error;
84  if( opt ) str << ": " << opt;
85  str << std::endl;
86  }
87 
88  str << "Usage: " << NAME << " [output opts.] [settings] [file ...]" << std::endl;
89  str << " " << NAME << " -h" << std::endl;
90  if( !error )
91  str << " If no file is specified the defautl test files will be used" << std::endl
92  << " -h print help text. " << std::endl
93  << " -v verbose output (may be specified multiple times) " << std::endl
94  << " -q quiet (minimal output) " << std::endl
95  << " -f <x>:<y>:<z>:<i>:<j>:<k> Do ray fire" << std::endl
96  << " -c write box geometry to Cubit journal file." << std::endl
97  << " -k write leaf contents to vtk files." << std::endl
98  << " -K write contents of leaf boxes intersected by rays to vtk file." << std::endl
99  << " -o <name> Save file containing tree and triangles. Mesh tag \"OBB_ROOT\"." << std::endl
100  << " -t <real> specify tolerance" << std::endl
101  << " -n <int> specify max entities per leaf node " << std::endl
102  << " -l <int> specify max tree levels" << std::endl
103  << " -r <real> specify worst cell split ratio" << std::endl
104  << " -R <real> specify best cell split ratio" << std::endl
105  << " -s force construction of surface tree" << std::endl
106  << " -S do not build surface tree." << std::endl
107  << " (Default: surface tree if file contains multiple surfaces" << std::endl
108  << " -u use unordered (Range) meshsets for tree nodes" << std::endl
109  << " -U use ordered (vector) meshsets for tree nodes" << std::endl
110  << " Verbosity (-q sets to 0, each -v increments, default is 1):" << std::endl
111  << " 0 - no output" << std::endl
112  << " 1 - status messages and error summary" << std::endl
113  << " 2 - print tree statistics " << std::endl
114  << " 3 - output errors for each node" << std::endl
115  << " 4 - print tree" << std::endl
116  << " 5 - print tree w/ contents of each node" << std::endl
117  << " See documentation for OrientedBoxTreeTool::Settings for " << std::endl
118  << " a description of tree generation settings." << std::endl;
119  exit( !!error );
120 }
121 
122 static const char* get_option( int& i, int argc, char* argv[] )
123 {
124  ++i;
125  if( i == argc ) usage( "Expected argument following option", argv[i - 1] );
126  return argv[i];
127 }
128 
129 static int get_int_option( int& i, int argc, char* argv[] )
130 {
131  const char* str = get_option( i, argc, argv );
132  char* end_ptr;
133  long val = strtol( str, &end_ptr, 0 );
134  if( !*str || *end_ptr ) usage( "Expected integer following option", argv[i - 1] );
135  return val;
136 }
137 
138 static double get_double_option( int& i, int argc, char* argv[] )
139 {
140  const char* str = get_option( i, argc, argv );
141  char* end_ptr;
142  double val = strtod( str, &end_ptr );
143  if( !*str || *end_ptr ) usage( "Expected real number following option", argv[i - 1] );
144  return val;
145 }
146 
147 static bool do_file( const char* filename );
148 
149 enum TriOption
150 {
151  DISABLE,
152  ENABLE,
153  AUTO
154 };
155 
156 static int verbosity = 1;
157 static OrientedBoxTreeTool::Settings settings;
158 static double tolerance = 1e-6;
159 static bool write_cubit = false;
160 static bool write_vtk = false;
161 static bool write_ray_vtk = false;
162 static std::vector< CartVect > rays;
163 static const char* save_file_name = 0;
164 static TriOption surfTree = AUTO;
165 
166 static void parse_ray( int& i, int argc, char* argv[] );
167 
168 int main( int argc, char* argv[] )
169 {
170 #ifdef MOAB_HAVE_MPI
171  int fail = MPI_Init( &argc, &argv );
172  if( fail ) return fail;
173 #endif
174 
175  initialize_default_files();
176 
177  std::vector< const char* > file_names;
178  bool flags = true;
179  for( int i = 1; i < argc; ++i )
180  {
181  if( flags && argv[i][0] == '-' )
182  {
183  if( !argv[i][1] || argv[i][2] ) usage( 0, argv[i] );
184  switch( argv[i][1] )
185  {
186  default:
187  usage( 0, argv[i] );
188  break;
189  case '-':
190  flags = false;
191  break;
192  case 'v':
193  ++verbosity;
194  break;
195  case 'q':
196  verbosity = 0;
197  break;
198  case 'h':
199  usage( 0, 0 );
200  break;
201  case 'c':
202  write_cubit = true;
203  break;
204  case 'k':
205  write_vtk = true;
206  break;
207  case 'K':
208  write_ray_vtk = true;
209  break;
210  case 's':
211  surfTree = ENABLE;
212  break;
213  case 'S':
214  surfTree = DISABLE;
215  break;
216  case 'u':
217  settings.set_options = MESHSET_SET;
218  break;
219  case 'U':
220  settings.set_options = MESHSET_ORDERED;
221  break;
222  case 'o':
223  save_file_name = get_option( i, argc, argv );
224  break;
225  case 'n':
226  settings.max_leaf_entities = get_int_option( i, argc, argv );
227  break;
228  case 'l':
229  settings.max_depth = get_int_option( i, argc, argv );
230  break;
231  case 'r':
232  settings.worst_split_ratio = get_double_option( i, argc, argv );
233  break;
234  case 'R':
235  settings.best_split_ratio = get_double_option( i, argc, argv );
236  break;
237  case 't':
238  tolerance = get_double_option( i, argc, argv );
239  break;
240  case 'f':
241  parse_ray( i, argc, argv );
242  break;
243  }
244  }
245  else
246  {
247  file_names.push_back( argv[i] );
248  }
249  }
250 
251  if( verbosity )
252  {
253  Core core;
254  std::string version;
255  core.impl_version( &version );
256  std::cout << version << std::endl;
257  if( verbosity > 1 )
258  std::cout << "max_leaf_entities: " << settings.max_leaf_entities << std::endl
259  << "max_depth: " << settings.max_depth << std::endl
260  << "worst_split_ratio: " << settings.worst_split_ratio << std::endl
261  << "best_split_ratio: " << settings.best_split_ratio << std::endl
262  << "tolerance: " << tolerance << std::endl
263  << "set type: "
264  << ( ( settings.set_options & MESHSET_ORDERED ) ? "ordered" : "set" ) << std::endl
265  << std::endl;
266  }
267 
268  if( !settings.valid() || tolerance < 0.0 )
269  {
270  std::cerr << "Invalid settings specified." << std::endl;
271  return 2;
272  }
273 
274  if( file_names.empty() )
275  {
276  std::cerr << "No file(s) specified." << std::endl;
277  for( ray_test_itr file = default_files_tests.begin(); file != default_files_tests.end(); file++ )
278  {
279  std::cerr << "Using default file \"" << file->first << '"' << std::endl;
280  file_names.push_back( file->first.c_str() );
281  }
282  }
283 
284  if( save_file_name && file_names.size() != 1 )
285  {
286  std::cout << "Only one input file allowed if \"-o\" option is specified." << std::endl;
287  std::cout << "Only testing with single file " << file_names[0] << std::endl;
288  file_names.erase( ++file_names.begin(), file_names.end() );
289  }
290 
291  int exit_val = 0;
292  for( unsigned j = 0; j < file_names.size(); ++j )
293  if( !do_file( file_names[j] ) ) ++exit_val;
294 
295 #ifdef MOAB_HAVE_MPI
296  fail = MPI_Finalize();
297  if( fail ) return fail;
298 #endif
299 
300  return exit_val ? exit_val + 2 : 0;
301 }
302 
303 static void parse_ray( int& i, int argc, char* argv[] )
304 {
305  CartVect point, direction;
306  if( 6 != sscanf( get_option( i, argc, argv ), "%lf:%lf:%lf:%lf:%lf:%lf", &point[0], &point[1], &point[2],
307  &direction[0], &direction[1], &direction[2] ) )
308  usage( "Expected ray specified as <x>:<y>:<z>:<i>:<j>:<k>", 0 );
309  direction.normalize();
310  rays.push_back( point );
311  rays.push_back( direction );
312 }
313 
314 class TreeValidator : public OrientedBoxTreeTool::Op
315 {
316  private:
317  Interface* const instance;
318  OrientedBoxTreeTool* const tool;
319  const bool printing;
320  const double epsilon;
321  bool surfaces;
322  std::ostream& stream;
323  OrientedBoxTreeTool::Settings settings;
324 
325  void print( EntityHandle handle, const char* string )
326  {
327  if( printing ) stream << instance->id_from_handle( handle ) << ": " << string << std::endl;
328  }
329 
330  ErrorCode error( EntityHandle handle, const char* string )
331  {
332  ++error_count;
333  print( handle, string );
334  return MB_SUCCESS;
335  }
336 
337  public:
338  unsigned entity_count;
339 
340  unsigned loose_box_count;
341  unsigned child_outside_count;
342  unsigned entity_outside_count;
343  unsigned num_entities_outside;
344  unsigned non_ortho_count;
345  unsigned error_count;
346  unsigned empty_leaf_count;
347  unsigned non_empty_non_leaf_count;
348  unsigned entity_invalid_count;
349  unsigned unsorted_axis_count;
350  unsigned non_unit_count;
351  unsigned duplicate_entity_count;
352  unsigned bad_outer_radius_count;
353  unsigned missing_surface_count;
354  unsigned multiple_surface_count;
355  std::set< EntityHandle > seen;
356  int surface_depth;
357  EntityHandle surface_handle;
358 
359  TreeValidator( Interface* instance_ptr,
360  OrientedBoxTreeTool* tool_ptr,
361  bool print_errors,
362  std::ostream& str,
363  double tol,
364  bool surfs,
365  OrientedBoxTreeTool::Settings s )
366  : instance( instance_ptr ), tool( tool_ptr ), printing( print_errors ), epsilon( tol ), surfaces( surfs ),
367  stream( str ), settings( s ), entity_count( 0 ), loose_box_count( 0 ), child_outside_count( 0 ),
368  entity_outside_count( 0 ), num_entities_outside( 0 ), non_ortho_count( 0 ), error_count( 0 ),
369  empty_leaf_count( 0 ), non_empty_non_leaf_count( 0 ), entity_invalid_count( 0 ), unsorted_axis_count( 0 ),
370  non_unit_count( 0 ), duplicate_entity_count( 0 ), bad_outer_radius_count( 0 ), missing_surface_count( 0 ),
371  multiple_surface_count( 0 ), surface_depth( -1 ), surface_handle( 0 )
372  {
373  }
374 
375  bool is_valid() const
376  {
377  return 0 == loose_box_count + child_outside_count + entity_outside_count + num_entities_outside +
378  non_ortho_count + error_count + empty_leaf_count + non_empty_non_leaf_count +
379  entity_invalid_count + unsorted_axis_count + non_unit_count + duplicate_entity_count +
380  bad_outer_radius_count + missing_surface_count + multiple_surface_count;
381  }
382 
383  virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
384 
385  virtual ErrorCode leaf( EntityHandle )
386  {
387  return MB_SUCCESS;
388  }
389 };
390 
391 ErrorCode TreeValidator::visit( EntityHandle node, int depth, bool& descend )
392 {
393  ErrorCode rval;
394  descend = true;
395 
396  Range contents;
397  rval = instance->get_entities_by_handle( node, contents );
398  if( MB_SUCCESS != rval ) return error( node, "Error getting contents of tree node. Corrupt tree?" );
399  entity_count += contents.size();
400 
401  if( surfaces )
402  {
403  // if no longer in subtree for previous surface, clear
404  if( depth <= surface_depth ) surface_depth = -1;
405 
406  EntityHandle surface = 0;
407  Range::iterator surf_iter = contents.lower_bound( MBENTITYSET );
408  if( surf_iter != contents.end() )
409  {
410  surface = *surf_iter;
411  contents.erase( surf_iter );
412  }
413 
414  if( surface )
415  {
416  if( surface_depth >= 0 )
417  {
418  ++multiple_surface_count;
419  print( node, "Multiple surfaces in encountered in same subtree." );
420  }
421  else
422  {
423  surface_depth = depth;
424  surface_handle = surface;
425  }
426  }
427  }
428 
429  std::vector< EntityHandle > children;
430  rval = tool->get_moab_instance()->get_child_meshsets( node, children );
431  if( MB_SUCCESS != rval || ( !children.empty() && children.size() != 2 ) )
432  return error( node, "Error getting children. Corrupt tree?" );
433 
434  OrientedBox box;
435  rval = tool->box( node, box );
436  if( MB_SUCCESS != rval ) return error( node, "Error getting oriented box from tree node. Corrupt tree?" );
437 
438  if( children.empty() && contents.empty() )
439  {
440  ++empty_leaf_count;
441  print( node, "Empty leaf node.\n" );
442  }
443  else if( !children.empty() && !contents.empty() )
444  {
445  ++non_empty_non_leaf_count;
446  print( node, "Non-leaf node is not empty." );
447  }
448 
449  if( surfaces && children.empty() && surface_depth < 0 )
450  {
451  ++missing_surface_count;
452  print( node, "Reached leaf node w/out encountering any surface set." );
453  }
454 
455  double dot_epsilon = epsilon * ( box.axis( 0 ) + box.axis( 1 ) + box.axis( 2 ) ).length();
456  if( box.axis( 0 ) % box.axis( 1 ) > dot_epsilon || box.axis( 0 ) % box.axis( 2 ) > dot_epsilon ||
457  box.axis( 1 ) % box.axis( 2 ) > dot_epsilon )
458  {
459  ++non_ortho_count;
460  print( node, "Box axes are not orthogonal" );
461  }
462 
463  if( !children.empty() )
464  {
465  for( int i = 0; i < 2; ++i )
466  {
467  OrientedBox other_box;
468  rval = tool->box( children[i], other_box );
469  if( MB_SUCCESS != rval )
470  return error( children[i], " Error getting oriented box from tree node. Corrupt tree?" );
471  // else if (!box.contained( other_box, epsilon )) {
472  // ++child_outside_count;
473  // print( children[i], "Parent box does not contain child box." );
474  // char string[64];
475  // sprintf(string, " Volume ratio is %f", other_box.volume()/box.volume() );
476  // print( children [i], string );
477  // }
478  else
479  {
480  double vol_ratio = other_box.volume() / box.volume();
481  if( vol_ratio > 2.0 )
482  {
483  char string[64];
484  sprintf( string, "child/parent volume ratio is %f", vol_ratio );
485  print( children[i], string );
486  sprintf( string, " child/parent area ratio is %f", other_box.area() / box.area() );
487  print( children[i], string );
488  }
489  }
490  }
491  }
492 
493  bool bad_element = false;
494  bool bad_element_handle = false;
495  bool bad_element_conn = false;
496  bool duplicate_element = false;
497  int num_outside = 0;
498  bool boundary[6] = { false, false, false, false, false, false };
499  for( Range::iterator it = contents.begin(); it != contents.end(); ++it )
500  {
501  EntityType type = instance->type_from_handle( *it );
502  int dim = CN::Dimension( type );
503  if( dim != 2 )
504  {
505  bad_element = true;
506  continue;
507  }
508 
509  const EntityHandle* conn;
510  int conn_len;
511  rval = instance->get_connectivity( *it, conn, conn_len );
512  if( MB_SUCCESS != rval )
513  {
514  bad_element_handle = true;
515  continue;
516  }
517 
518  std::vector< CartVect > coords( conn_len );
519  rval = instance->get_coords( conn, conn_len, coords[0].array() );
520  if( MB_SUCCESS != rval )
521  {
522  bad_element_conn = true;
523  continue;
524  }
525 
526  bool outside = false;
527  for( std::vector< CartVect >::iterator j = coords.begin(); j != coords.end(); ++j )
528  {
529  if( !box.contained( *j, epsilon ) )
530  outside = true;
531  else
532  for( int d = 0; d < 3; ++d )
533  {
534 #if MB_ORIENTED_BOX_UNIT_VECTORS
535  double n = box.axis( d ) % ( *j - box.center );
536  if( fabs( n - box.length[d] ) <= epsilon ) boundary[2 * d] = true;
537  if( fabs( n + box.length[d] ) <= epsilon ) boundary[2 * d + 1] = true;
538 #else
539  double ln = box.axis( d ).length();
540  CartVect v1 = *j - box.center - box.axis[d];
541  CartVect v2 = *j - box.center + box.axis[d];
542  if( fabs( v1 % box.axis[d] ) <= ln * epsilon ) boundary[2 * d] = true;
543  if( fabs( v2 % box.axis[d] ) <= ln * epsilon ) boundary[2 * d + 1] = true;
544 #endif
545  }
546  }
547  if( outside ) ++num_outside;
548 
549  if( !seen.insert( *it ).second )
550  {
551  duplicate_element = true;
552  ++duplicate_entity_count;
553  }
554  }
555 
556  CartVect alength( box.axis( 0 ).length(), box.axis( 1 ).length(), box.axis( 2 ).length() );
557 #if MB_ORIENTED_BOX_UNIT_VECTORS
558  CartVect length = box.length;
559 #else
560  CartVect length = alength;
561 #endif
562 
563  if( length[0] > length[1] || length[0] > length[2] || length[1] > length[2] )
564  {
565  ++unsorted_axis_count;
566  print( node, "Box axes are not ordered from shortest to longest." );
567  }
568 
569 #if MB_ORIENTED_BOX_UNIT_VECTORS
570  if( fabs( alength[0] - 1.0 ) > epsilon || fabs( alength[1] - 1.0 ) > epsilon || fabs( alength[2] - 1.0 ) > epsilon )
571  {
572  ++non_unit_count;
573  print( node, "Box axes are not unit vectors." );
574  }
575 #endif
576 
577 #if MB_ORIENTED_BOX_OUTER_RADIUS
578  if( fabs( length.length() - box.radius ) > tolerance )
579  {
580  ++bad_outer_radius_count;
581  print( node, "Box has incorrect outer radius." );
582  }
583 #endif
584 
585  if( depth + 1 < settings.max_depth && contents.size() > (unsigned)( 4 * settings.max_leaf_entities ) )
586  {
587  char string[64];
588  sprintf( string, "leaf at depth %d with %u entities", depth, (unsigned)contents.size() );
589  print( node, string );
590  }
591 
592  bool all_boundaries = true;
593  for( int f = 0; f < 6; ++f )
594  all_boundaries = all_boundaries && boundary[f];
595 
596  if( bad_element )
597  {
598  ++entity_invalid_count;
599  print( node, "Set contained an entity with an inappropriate dimension." );
600  }
601  if( bad_element_handle )
602  {
603  ++error_count;
604  print( node, "Error querying face contained in set." );
605  }
606  if( bad_element_conn )
607  {
608  ++error_count;
609  print( node, "Error querying connectivity of element." );
610  }
611  if( duplicate_element )
612  {
613  print( node, "Elements occur in multiple leaves of tree." );
614  }
615  if( num_outside > 0 )
616  {
617  ++entity_outside_count;
618  num_entities_outside += num_outside;
619  if( printing )
620  stream << instance->id_from_handle( node ) << ": " << num_outside << " elements outside box." << std::endl;
621  }
622  else if( !all_boundaries && !contents.empty() )
623  {
624  ++loose_box_count;
625  print( node, "Box does not fit contained elements tightly." );
626  }
627 
628  return MB_SUCCESS;
629 }
630 
631 class CubitWriter : public OrientedBoxTreeTool::Op
632 {
633  public:
634  CubitWriter( FILE* file_ptr, OrientedBoxTreeTool* tool_ptr ) : file( file_ptr ), tool( tool_ptr ) {}
635 
636  ErrorCode visit( EntityHandle node, int depth, bool& descend );
637  ErrorCode leaf( EntityHandle )
638  {
639  return MB_SUCCESS;
640  }
641 
642  private:
643  FILE* file;
644  OrientedBoxTreeTool* tool;
645 };
646 
647 ErrorCode CubitWriter::visit( EntityHandle node, int, bool& descend )
648 {
649  descend = true;
650  OrientedBox box;
651  ErrorCode rval = tool->box( node, box );
652  if( rval != MB_SUCCESS ) return rval;
653 
654  double sign[] = { -1, 1 };
655  for( int i = 0; i < 2; ++i )
656  for( int j = 0; j < 2; ++j )
657  for( int k = 0; k < 2; ++k )
658  {
659 #if MB_ORIENTED_BOX_UNIT_VECTORS
660  CartVect corner = box.center + box.length[0] * sign[i] * box.axis( 0 ) +
661  box.length[1] * sign[j] * box.axis( 1 ) + box.length[2] * sign[k] * box.axis( 2 );
662 #else
663  CartVect corner =
664  box.center + sign[i] * box.axis( 0 ) + sign[j] * box.axis( 1 ) + sign[k] * box.axis( 2 );
665 #endif
666  fprintf( file, "create vertex %f %f %f\n", corner[0], corner[1], corner[2] );
667  }
668  fprintf( file, "#{i=Id(\"vertex\")-7}\n" );
669  fprintf( file, "create surface vertex {i } {i+1} {i+3} {i+2}\n" );
670  fprintf( file, "create surface vertex {i+4} {i+5} {i+7} {i+6}\n" );
671  fprintf( file, "create surface vertex {i+1} {i+0} {i+4} {i+5}\n" );
672  fprintf( file, "create surface vertex {i+3} {i+2} {i+6} {i+7}\n" );
673  fprintf( file, "create surface vertex {i+2} {i+0} {i+4} {i+6}\n" );
674  fprintf( file, "create surface vertex {i+3} {i+1} {i+5} {i+7}\n" );
675  fprintf( file, "delete vertex {i}-{i+7}\n" );
676  fprintf( file, "#{s=Id(\"surface\")-5}\n" );
677  fprintf( file, "create volume surface {s} {s+1} {s+2} {s+3} {s+4} {s+5} noheal\n" );
678  int id = tool->get_moab_instance()->id_from_handle( node );
679  fprintf( file, "volume {Id(\"volume\")} name \"cell%d\"\n", id );
680 
681  return MB_SUCCESS;
682 }
683 
684 class VtkWriter : public OrientedBoxTreeTool::Op
685 {
686  public:
687  VtkWriter( std::string base_name, Interface* interface ) : baseName( base_name ), instance( interface ) {}
688 
689  ErrorCode visit( EntityHandle node, int depth, bool& descend );
690 
691  ErrorCode leaf( EntityHandle /*node*/ )
692  {
693  return MB_SUCCESS;
694  }
695 
696  private:
697  std::string baseName;
698  Interface* instance;
699 };
700 
701 ErrorCode VtkWriter::visit( EntityHandle node, int, bool& descend )
702 {
703  descend = true;
704 
705  ErrorCode rval;
706  int count;
707  rval = instance->get_number_entities_by_handle( node, count );
708  if( MB_SUCCESS != rval || 0 == count ) return rval;
709 
710  int id = instance->id_from_handle( node );
711  char id_str[32];
712  sprintf( id_str, "%d", id );
713 
714  std::string file_name( baseName );
715  file_name += ".";
716  file_name += id_str;
717  file_name += ".vtk";
718 
719  return instance->write_mesh( file_name.c_str(), &node, 1 );
720 }
721 
722 static bool do_ray_fire_test( OrientedBoxTreeTool& tool,
723  EntityHandle root_set,
724  const char* filename,
725  bool have_surface_tree );
726 
727 static bool do_closest_point_test( OrientedBoxTreeTool& tool, EntityHandle root_set, bool have_surface_tree );
728 
729 static ErrorCode save_tree( Interface* instance, const char* filename, EntityHandle tree_root );
730 
731 static bool do_file( const char* filename )
732 {
733  ErrorCode rval;
734  Core instance;
735  Interface* const iface = &instance;
736  OrientedBoxTreeTool tool( iface );
737  bool haveSurfTree = false;
738 
739  if( verbosity ) std::cout << filename << std::endl << "------" << std::endl;
740 
741  rval = iface->load_mesh( filename );
742  if( MB_SUCCESS != rval )
743  {
744  if( verbosity ) std::cout << "Failed to read file: \"" << filename << '"' << std::endl;
745  return false;
746  }
747 
748  // IF building from surfaces, get surfaces.
749  // If AUTO and less than two surfaces, don't build from surfaces.
750  Range surfaces;
751  if( surfTree != DISABLE )
752  {
753  Tag surftag;
754  rval = iface->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, surftag );
755  if( MB_SUCCESS == rval )
756  {
757  int dim = 2;
758  const void* tagvalues[] = { &dim };
759  rval = iface->get_entities_by_type_and_tag( 0, MBENTITYSET, &surftag, tagvalues, 1, surfaces );
760  if( MB_SUCCESS != rval && MB_ENTITY_NOT_FOUND != rval ) return false;
761  }
762  else if( MB_TAG_NOT_FOUND != rval )
763  return false;
764 
765  if( ENABLE == surfTree && surfaces.empty() )
766  {
767  std::cerr << "No Surfaces found." << std::endl;
768  return false;
769  }
770 
771  haveSurfTree = ( ENABLE == surfTree ) || ( surfaces.size() > 1 );
772  }
773 
774  EntityHandle root;
775  Range entities;
776  if( !haveSurfTree )
777  {
778  rval = iface->get_entities_by_dimension( 0, 2, entities );
779  if( MB_SUCCESS != rval )
780  {
781  std::cerr << "get_entities_by_dimension( 2 ) failed." << std::endl;
782  return false;
783  }
784 
785  if( entities.empty() )
786  {
787  if( verbosity ) std::cout << "File \"" << filename << "\" contains no 2D elements" << std::endl;
788  return false;
789  }
790 
791  if( verbosity ) std::cout << "Building tree from " << entities.size() << " 2D elements" << std::endl;
792 
793  rval = tool.build( entities, root, &settings );
794  if( MB_SUCCESS != rval )
795  {
796  if( verbosity ) std::cout << "Failed to build tree." << std::endl;
797  return false;
798  }
799  }
800  else
801  {
802 
803  if( verbosity ) std::cout << "Building tree from " << surfaces.size() << " surfaces" << std::endl;
804 
805  // Build subtree for each surface, get list of all entities to use later
806  Range surf_trees, surf_tris;
807  EntityHandle surf_root;
808  for( Range::iterator s = surfaces.begin(); s != surfaces.end(); ++s )
809  {
810  surf_tris.clear();
811  rval = iface->get_entities_by_dimension( *s, 2, surf_tris );
812  if( MB_SUCCESS != rval ) return false;
813  rval = tool.build( surf_tris, surf_root, &settings );
814  if( MB_SUCCESS != rval )
815  {
816  if( verbosity ) std::cout << "Failed to build tree for surface." << std::endl;
817  return false;
818  }
819  surf_trees.insert( surf_root );
820  entities.merge( surf_tris );
821  rval = iface->add_entities( surf_root, &*s, 1 );
822  if( MB_SUCCESS != rval ) return false;
823  }
824 
825  rval = tool.join_trees( surf_trees, root, &settings );
826  if( MB_SUCCESS != rval )
827  {
828  if( verbosity ) std::cout << "Failed to build tree." << std::endl;
829  return false;
830  }
831 
832  if( verbosity )
833  std::cout << "Built tree from " << surfaces.size() << " surfaces"
834  << " (" << entities.size() - surfaces.size() << " elements)" << std::endl;
835 
836  entities.merge( surfaces );
837  }
838 
839  if( write_cubit )
840  {
841  std::string name = filename;
842  name += ".boxes.jou";
843  FILE* ptr = fopen( name.c_str(), "w+" );
844  if( !ptr )
845  {
846  perror( name.c_str() );
847  }
848  else
849  {
850  if( verbosity ) std::cout << "Writing: " << name << std::endl;
851  fprintf( ptr, "graphics off\n" );
852  CubitWriter op( ptr, &tool );
853  tool.preorder_traverse( root, op );
854  fprintf( ptr, "graphics on\n" );
855  fclose( ptr );
856  }
857  }
858 
859  if( write_vtk )
860  {
861  VtkWriter op( filename, iface );
862  if( verbosity )
863  std::cout << "Writing leaf contents as : " << filename << ".xxx.vtk where 'xxx' is the set id" << std::endl;
864  tool.preorder_traverse( root, op );
865  }
866 
867  bool print_errors = ( verbosity > 2 ), print_contents = ( verbosity > 4 );
868  if( verbosity > 3 ) tool.print( root, std::cout, print_contents );
869  if( verbosity > 1 )
870  {
871  rval = tool.stats( root, std::cout );
872  if( MB_SUCCESS != rval ) std::cout << "****** Failed to get tree statistics ******" << std::endl;
873  }
874 
875  TreeValidator op( iface, &tool, print_errors, std::cout, tolerance, haveSurfTree, settings );
876  rval = tool.preorder_traverse( root, op );
877  bool result = op.is_valid();
878  if( MB_SUCCESS != rval )
879  {
880  result = false;
881  if( verbosity ) std::cout << "Errors traversing tree. Corrupt tree?" << std::endl;
882  }
883 
884  bool missing = ( op.entity_count != entities.size() );
885  if( missing ) result = false;
886 
887  if( verbosity )
888  {
889  if( result )
890  std::cout << std::endl << "No errors detected." << std::endl;
891  else
892  std::cout << std::endl << "*********************** ERROR SUMMARY **********************" << std::endl;
893  if( op.child_outside_count )
894  std::cout << "* " << op.child_outside_count << " child boxes not contained in parent." << std::endl;
895  if( op.entity_outside_count )
896  std::cout << "* " << op.entity_outside_count << " nodes containing entities outside of box." << std::endl;
897  if( op.num_entities_outside )
898  std::cout << "* " << op.num_entities_outside << " entities outside boxes." << std::endl;
899  if( op.empty_leaf_count ) std::cout << "* " << op.empty_leaf_count << " empty leaf nodes." << std::endl;
900  if( op.non_empty_non_leaf_count )
901  std::cout << "* " << op.non_empty_non_leaf_count << " non-leaf nodes containing entities." << std::endl;
902  if( op.duplicate_entity_count )
903  std::cout << "* " << op.duplicate_entity_count << " duplicate entities in leaves." << std::endl;
904  if( op.missing_surface_count )
905  std::cout << "* " << op.missing_surface_count << " leaves outside surface subtrees." << std::endl;
906  if( op.multiple_surface_count )
907  std::cout << "* " << op.multiple_surface_count << " surfaces within other surface subtrees." << std::endl;
908  if( op.non_ortho_count )
909  std::cout << "* " << op.non_ortho_count << " boxes with non-orthononal axes." << std::endl;
910  if( op.non_unit_count ) std::cout << "* " << op.non_unit_count << " boxes with non-unit axes." << std::endl;
911  if( op.bad_outer_radius_count )
912  std::cout << "* " << op.bad_outer_radius_count << " boxes incorrect outer radius." << std::endl;
913  if( op.unsorted_axis_count )
914  std::cout << "* " << op.unsorted_axis_count << " boxes axes in unsorted order." << std::endl;
915  if( op.loose_box_count )
916  std::cout << "* " << op.loose_box_count << " boxes that do not fit the contained entities tightly."
917  << std::endl;
918  if( op.error_count + op.entity_invalid_count )
919  std::cout << "* " << op.error_count + op.entity_invalid_count << " other errors while traversing tree."
920  << std::endl;
921  if( missing )
922  std::cout << "* tree built from " << entities.size() << " entities contains " << op.entity_count
923  << " entities." << std::endl;
924  if( !result ) std::cout << "************************************************************" << std::endl;
925  }
926 
927  if( result && save_file_name )
928  {
929  if( MB_SUCCESS == save_tree( iface, save_file_name, root ) )
930  std::cerr << "Wrote '" << save_file_name << "'" << std::endl;
931  else
932  std::cerr << "FAILED TO WRITE '" << save_file_name << "'" << std::endl;
933  }
934 
935  if( !do_ray_fire_test( tool, root, filename, haveSurfTree ) )
936  {
937  if( verbosity ) std::cout << "Ray fire test failed." << std::endl;
938  result = false;
939  }
940 
941  if( !do_closest_point_test( tool, root, haveSurfTree ) )
942  {
943  if( verbosity ) std::cout << "Closest point test failed" << std::endl;
944  result = false;
945  }
946 
947  rval = tool.delete_tree( root );
948  if( MB_SUCCESS != rval )
949  {
950  if( verbosity ) std::cout << "delete_tree failed." << std::endl;
951  result = false;
952  }
953 
954  return result;
955 }
956 
957 static void count_non_tol( std::vector< double > intersections, int& non_tol_count, double& non_tol_dist )
958 {
959 
960  for( size_t i = 0; i < intersections.size(); ++i )
961  {
962  if( intersections[i] > tolerance )
963  {
964  ++non_tol_count;
965  if( intersections[i] < non_tol_dist ) non_tol_dist = intersections[i];
966  }
967  }
968 }
969 
970 static bool check_ray_intersect_tris( OrientedBoxTreeTool& tool,
971  EntityHandle root_set,
972  RayTest& test,
973  int& non_tol_count,
974  double& non_tol_dist,
975  OrientedBoxTreeTool::TrvStats& stats )
976 {
977  ErrorCode rval;
978  bool result = true;
979 
980  non_tol_dist = std::numeric_limits< double >::max();
981  non_tol_count = 0;
982 
983  std::vector< double > intersections;
984  std::vector< EntityHandle > facets;
985  rval = tool.ray_intersect_triangles( intersections, facets, root_set, tolerance, test.point.array(),
986  test.direction.array(), 0, &stats );
987  if( MB_SUCCESS != rval )
988  {
989  if( verbosity ) std::cout << " Call to OrientedBoxTreeTool::ray_intersect_triangles failed." << std::endl;
990  result = false;
991  }
992  else
993  {
994 
995  if( intersections.size() != test.expected_hits )
996  {
997  if( verbosity > 2 )
998  std::cout << " Expected " << test.expected_hits << " and got " << intersections.size()
999  << " hits for ray fire of " << test.description << std::endl;
1000  if( verbosity > 3 )
1001  {
1002  for( unsigned j = 0; j < intersections.size(); ++j )
1003  std::cout << " " << intersections[j];
1004  std::cout << std::endl;
1005  }
1006  result = false;
1007  }
1008 
1009  count_non_tol( intersections, non_tol_count, non_tol_dist );
1010  }
1011  return result;
1012 }
1013 
1014 static bool check_ray_intersect_sets( OrientedBoxTreeTool& tool,
1015  EntityHandle root_set,
1016  RayTest& test,
1017  int& non_tol_count,
1018  double& non_tol_dist,
1019  OrientedBoxTreeTool::TrvStats& stats )
1020 {
1021 
1022  ErrorCode rval;
1023  bool result = true;
1024 
1025  non_tol_dist = std::numeric_limits< double >::max();
1026  non_tol_count = 0;
1027 
1028  const int NUM_NON_TOL_INT = 1;
1029 
1030  std::vector< double > intersections;
1031  std::vector< EntityHandle > surfs;
1032  std::vector< EntityHandle > facets;
1033 
1034  OrientedBoxTreeTool::IntersectSearchWindow search_win;
1035  OrientedBoxTreeTool::IntRegCtxt int_reg_ctxt;
1036  rval = tool.ray_intersect_sets( intersections, surfs, facets, root_set, tolerance, test.point.array(),
1037  test.direction.array(), search_win, int_reg_ctxt, &stats );
1038 
1039  if( MB_SUCCESS != rval )
1040  {
1041  if( verbosity ) std::cout << " Call to OrientedBoxTreeTool::ray_intersect_sets failed." << std::endl;
1042  result = false;
1043  }
1044  else
1045  {
1046 
1047  if( surfs.size() != intersections.size() )
1048  {
1049  if( verbosity ) std::cout << " ray_intersect_sets did not return sets for all intersections." << std::endl;
1050  result = false;
1051  }
1052 
1053  count_non_tol( intersections, non_tol_count, non_tol_dist );
1054 
1055  if( non_tol_count > NUM_NON_TOL_INT )
1056  {
1057  if( verbosity )
1058  std::cout << " Requested " << NUM_NON_TOL_INT << "intersections "
1059  << " beyond tolerance. Got " << non_tol_count << std::endl;
1060  result = false;
1061  }
1062  }
1063 
1064  return result;
1065 }
1066 
1067 static bool do_ray_fire_test( OrientedBoxTreeTool& tool,
1068  EntityHandle root_set,
1069  const char* filename,
1070  bool haveSurfTree )
1071 {
1072  if( verbosity > 1 ) std::cout << "beginning ray fire tests" << std::endl;
1073 
1074  OrientedBox box;
1075  ErrorCode rval = tool.box( root_set, box );
1076  if( MB_SUCCESS != rval )
1077  {
1078  if( verbosity ) std::cerr << "Error getting box for tree root set" << std::endl;
1079  return false;
1080  }
1081 
1082  /* Do standard ray fire tests */
1083  std::cout << box << std::endl;
1084 
1085  CartVect origin( 0., 0., 0. );
1086  CartVect unitDiag( 1., 1., 1. );
1087  std::vector< RayTest > tests;
1088  RayTest default_tests[] = { { "large axis through box", 2, box.center - 1.2 * box.scaled_axis( 2 ), box.axis( 2 ) },
1089  { "small axis through box", 2, box.center - 1.2 * box.scaled_axis( 0 ), box.axis( 0 ) },
1090  { "parallel miss", 0, box.center + 2.0 * box.scaled_axis( 1 ), box.axis( 2 ) },
1091  { "skew miss", 0, box.center + box.dimensions(), box.dimensions() * box.axis( 2 ) } };
1092  const size_t num_def_test = sizeof( default_tests ) / sizeof( default_tests[0] );
1093  tests.insert( tests.begin(), &default_tests[0], &default_tests[0] + num_def_test );
1094  tests.insert( tests.end(), default_files_tests[filename].begin(), default_files_tests[filename].end() );
1095 
1096  OrientedBoxTreeTool::TrvStats stats;
1097 
1098  bool result = true;
1099  const size_t num_test = tests.size();
1100  for( size_t i = 0; i < num_test; ++i )
1101  {
1102  tests[i].direction.normalize();
1103  if( verbosity > 2 )
1104  {
1105  std::cout << ( 0 == i ? "** Common tests\n" : ( num_def_test == i ? "** File-specific tests\n" : "" ) );
1106  std::cout << " " << tests[i].description << " " << tests[i].point << " " << tests[i].direction
1107  << std::endl;
1108  }
1109 
1110  int rit_non_tol_count = 0;
1111  double rit_non_tol_dist = std::numeric_limits< double >::max();
1112  if( !check_ray_intersect_tris( tool, root_set, tests[i], rit_non_tol_count, rit_non_tol_dist, stats ) )
1113  {
1114  result = false;
1115  continue;
1116  }
1117 
1118  if( !haveSurfTree ) continue;
1119 
1120  int ris_non_tol_count = 0;
1121  double ris_non_tol_dist = std::numeric_limits< double >::max();
1122  if( !check_ray_intersect_sets( tool, root_set, tests[i], ris_non_tol_count, ris_non_tol_dist, stats ) )
1123  {
1124  result = false;
1125  continue;
1126  }
1127 
1128  if( !rit_non_tol_count && ris_non_tol_count )
1129  {
1130  if( verbosity )
1131  std::cout << " ray_intersect_sets returned intersection not found by "
1132  "ray_intersect_triangles"
1133  << std::endl;
1134  result = false;
1135  continue;
1136  }
1137  else if( rit_non_tol_count && !ris_non_tol_count )
1138  {
1139  if( verbosity )
1140  std::cout << " ray_intersect_sets missed intersection found by ray_intersect_triangles" << std::endl;
1141  result = false;
1142  continue;
1143  }
1144  else if( rit_non_tol_count && ris_non_tol_count && fabs( rit_non_tol_dist - ris_non_tol_dist ) > tolerance )
1145  {
1146  if( verbosity )
1147  std::cout << " ray_intersect_sets and ray_intersect_triangles did not find same "
1148  "closest intersection"
1149  << std::endl;
1150  result = false;
1151  }
1152  }
1153 
1154  /* Do ray fire for any user-specified rays */
1155 
1156  for( size_t i = 0; i < rays.size(); i += 2 )
1157  {
1158  std::cout << rays[i] << "+" << rays[i + 1] << " : ";
1159 
1160  if( !haveSurfTree )
1161  {
1162  Range leaves;
1163  std::vector< double > intersections;
1164  std::vector< EntityHandle > intersection_facets;
1165  rval = tool.ray_intersect_boxes( leaves, root_set, tolerance, rays[i].array(), rays[i + 1].array(), 0,
1166  &stats );
1167  if( MB_SUCCESS != rval )
1168  {
1169  std::cout << "FAILED" << std::endl;
1170  result = false;
1171  continue;
1172  }
1173 
1174  if( !leaves.empty() && write_ray_vtk )
1175  {
1176  std::string num, name( filename );
1177  std::stringstream s;
1178  s << ( i / 2 );
1179  s >> num;
1180  name += "-ray";
1181  name += num;
1182  name += ".vtk";
1183 
1184  std::vector< EntityHandle > sets( leaves.size() );
1185  std::copy( leaves.begin(), leaves.end(), sets.begin() );
1186  tool.get_moab_instance()->write_mesh( name.c_str(), &sets[0], sets.size() );
1187  if( verbosity ) std::cout << "(Wrote " << name << ") ";
1188  }
1189 
1190  rval = tool.ray_intersect_triangles( intersections, intersection_facets, leaves, tolerance, rays[i].array(),
1191  rays[i + 1].array(), 0 );
1192  if( MB_SUCCESS != rval )
1193  {
1194  std::cout << "FAILED" << std::endl;
1195  result = false;
1196  continue;
1197  }
1198 
1199  if( intersections.empty() )
1200  {
1201  std::cout << "(none)" << std::endl;
1202  continue;
1203  }
1204 
1205  std::cout << intersections[0];
1206  for( unsigned j = 1; j < intersections.size(); ++j )
1207  std::cout << ", " << intersections[j];
1208  std::cout << std::endl;
1209 
1210  if( !leaves.empty() && write_cubit && verbosity > 2 )
1211  {
1212  std::cout << " intersected boxes:";
1213  for( Range::iterator i2 = leaves.begin(); i2 != leaves.end(); ++i2 )
1214  std::cout << " " << tool.get_moab_instance()->id_from_handle( *i2 );
1215  std::cout << std::endl;
1216  }
1217  }
1218  else
1219  {
1220  std::vector< double > intersections;
1221  std::vector< EntityHandle > surfaces;
1222  std::vector< EntityHandle > facets;
1223 
1224  OrientedBoxTreeTool::IntersectSearchWindow search_win;
1225  OrientedBoxTreeTool::IntRegCtxt int_reg_ctxt;
1226  rval = tool.ray_intersect_sets( intersections, surfaces, facets, root_set, tolerance, rays[i].array(),
1227  rays[i + 1].array(), search_win, int_reg_ctxt, &stats );
1228 
1229  if( MB_SUCCESS != rval )
1230  {
1231  if( verbosity ) std::cout << " Call to OrientedBoxTreeTool::ray_intersect_sets failed." << std::endl;
1232  result = false;
1233  continue;
1234  }
1235 
1236  if( !surfaces.empty() && write_ray_vtk )
1237  {
1238  std::string num, name( filename );
1239  std::stringstream s;
1240  s << ( i / 2 );
1241  s >> num;
1242  name += "-ray";
1243  name += num;
1244  name += ".vtk";
1245 
1246  tool.get_moab_instance()->write_mesh( name.c_str(), &surfaces[0], surfaces.size() );
1247  if( verbosity ) std::cout << "(Wrote " << name << ") ";
1248  }
1249 
1250  if( intersections.size() != surfaces.size() )
1251  {
1252  std::cout << "Mismatched output lists." << std::endl;
1253  result = false;
1254  continue;
1255  }
1256 
1257  if( intersections.empty() )
1258  {
1259  std::cout << "(none)" << std::endl;
1260  continue;
1261  }
1262 
1263  Tag idtag = tool.get_moab_instance()->globalId_tag();
1264  std::vector< int > ids( surfaces.size() );
1265  rval = tool.get_moab_instance()->tag_get_data( idtag, &surfaces[0], surfaces.size(), &ids[0] );
1266  if( MB_SUCCESS != rval )
1267  {
1268  std::cout << "NO GLOBAL_ID TAG ON SURFACE." << std::endl;
1269  continue;
1270  }
1271 
1272  // group by surfaces
1273  std::map< int, double > intmap;
1274  for( unsigned j = 0; j < intersections.size(); ++j )
1275  intmap[ids[j]] = intersections[j];
1276 
1277  std::map< int, double >::iterator it = intmap.begin();
1278  int prevsurf = it->first;
1279  std::cout << "Surf" << it->first << " " << it->second;
1280  for( ++it; it != intmap.end(); ++it )
1281  {
1282  std::cout << ", ";
1283  if( it->first != prevsurf )
1284  {
1285  prevsurf = it->first;
1286  std::cout << "Surf" << it->first << " ";
1287  }
1288  std::cout << it->second;
1289  }
1290  std::cout << std::endl;
1291  }
1292  }
1293 
1294  if( verbosity > 1 )
1295  {
1296  std::cout << "Traversal statistics for ray fire tests: " << std::endl;
1297  stats.print( std::cout );
1298  }
1299 
1300  return result;
1301 }
1302 
1303 static ErrorCode save_tree( Interface* instance, const char* filename, EntityHandle tree_root )
1304 {
1305  ErrorCode rval;
1306  Tag tag;
1307 
1308  rval = instance->tag_get_handle( "OBB_ROOT", 1, MB_TYPE_HANDLE, tag, MB_TAG_SPARSE | MB_TAG_CREAT );
1309  if( MB_SUCCESS != rval ) return rval;
1310 
1311  const EntityHandle root = 0;
1312  rval = instance->tag_set_data( tag, &root, 1, &tree_root );
1313  if( MB_SUCCESS != rval ) return rval;
1314 
1315  return instance->write_mesh( filename );
1316 }
1317 
1318 static ErrorCode tri_coords( Interface* moab, EntityHandle tri, CartVect coords[3] )
1319 {
1320  ErrorCode rval;
1321  const EntityHandle* conn;
1322  int len;
1323 
1324  rval = moab->get_connectivity( tri, conn, len, true );
1325  if( MB_SUCCESS != rval ) return rval;
1326  if( len != 3 ) return MB_FAILURE;
1327  rval = moab->get_coords( conn, 3, coords[0].array() );
1328  return rval;
1329 }
1330 
1331 static ErrorCode closest_point_in_triangles( Interface* moab,
1332  const CartVect& to_pos,
1333  CartVect& result_pos,
1334  EntityHandle& result_tri )
1335 {
1336  ErrorCode rval;
1337 
1338  Range triangles;
1339  rval = moab->get_entities_by_type( 0, MBTRI, triangles );
1340  if( MB_SUCCESS != rval ) return rval;
1341 
1342  if( triangles.empty() ) return MB_FAILURE;
1343 
1344  Range::iterator i = triangles.begin();
1345  CartVect coords[3];
1346  rval = tri_coords( moab, *i, coords );
1347  if( MB_SUCCESS != rval ) return rval;
1348  result_tri = *i;
1349  GeomUtil::closest_location_on_tri( to_pos, coords, result_pos );
1350  CartVect diff = to_pos - result_pos;
1351  double shortest_dist_sqr = diff % diff;
1352 
1353  for( ++i; i != triangles.end(); ++i )
1354  {
1355  rval = tri_coords( moab, *i, coords );
1356  if( MB_SUCCESS != rval ) return rval;
1357  CartVect pos;
1358  GeomUtil::closest_location_on_tri( to_pos, coords, pos );
1359  diff = to_pos - pos;
1360  double dsqr = diff % diff;
1361  if( dsqr < shortest_dist_sqr )
1362  {
1363  shortest_dist_sqr = dsqr;
1364  result_pos = pos;
1365  result_tri = *i;
1366  }
1367  }
1368 
1369  return MB_SUCCESS;
1370 }
1371 
1372 static bool tri_in_set( Interface* moab, EntityHandle set, EntityHandle tri )
1373 {
1374  Range tris;
1375  ErrorCode rval = moab->get_entities_by_type( set, MBTRI, tris );
1376  if( MB_SUCCESS != rval ) return false;
1377  Range::iterator i = tris.find( tri );
1378  return i != tris.end();
1379 }
1380 
1381 static bool do_closest_point_test( OrientedBoxTreeTool& tool, EntityHandle root_set, bool have_surface_tree )
1382 {
1383  if( verbosity > 1 ) std::cout << "beginning closest point tests" << std::endl;
1384 
1385  ErrorCode rval;
1386  Interface* moab = tool.get_moab_instance();
1387  EntityHandle set;
1388  EntityHandle* set_ptr = have_surface_tree ? &set : 0;
1389  bool result = true;
1390 
1391  // get root box
1392  OrientedBox box;
1393  rval = tool.box( root_set, box );
1394  if( MB_SUCCESS != rval )
1395  {
1396  if( verbosity ) std::cerr << "Invalid tree in do_closest_point_test\n";
1397  return false;
1398  }
1399 
1400  OrientedBoxTreeTool::TrvStats stats;
1401 
1402  // chose some points to test
1403  CartVect points[] = { box.center + box.scaled_axis( 0 ), box.center + 2 * box.scaled_axis( 1 ),
1404  box.center + 0.5 * box.scaled_axis( 2 ),
1405  box.center + -2 * box.scaled_axis( 0 ) + -2 * box.scaled_axis( 1 ) +
1406  -2 * box.scaled_axis( 2 ),
1407  CartVect( 100, 100, 100 ) };
1408  const int num_pts = sizeof( points ) / sizeof( points[0] );
1409 
1410  // test each point
1411  for( int i = 0; i < num_pts; ++i )
1412  {
1413  if( verbosity >= 3 ) std::cout << "Evaluating closest point to " << points[i] << std::endl;
1414 
1415  CartVect n_result, t_result;
1416  EntityHandle n_tri = 0, t_tri;
1417 
1418  // find closest point the slow way
1419  rval = closest_point_in_triangles( moab, points[i], n_result, n_tri );
1420  if( MB_SUCCESS != rval )
1421  {
1422  std::cerr << "Internal MOAB error in do_closest_point_test" << std::endl;
1423  result = false;
1424  continue;
1425  }
1426 
1427  // find closest point using tree
1428  rval = tool.closest_to_location( points[i].array(), root_set, t_result.array(), t_tri, set_ptr, &stats );
1429  if( MB_SUCCESS != rval )
1430  {
1431  if( verbosity )
1432  std::cout << "OrientedBoxTreeTool:: closest_to_location( " << points[i] << " ) FAILED!" << std::endl;
1433  result = false;
1434  continue;
1435  }
1436 
1437  CartVect diff = t_result - n_result;
1438  if( diff.length() > tolerance )
1439  {
1440  if( verbosity )
1441  std::cout << "Closest point to " << points[i] << " INCORRECT! (" << t_result << " != " << n_result
1442  << ")" << std::endl;
1443  result = false;
1444  continue;
1445  }
1446 
1447  if( t_tri != n_tri )
1448  {
1449  // if result point is triangle, then OK because
1450  // already tested that it is the same location
1451  // as the expected value. We just have a case where
1452  // the point was on an edge or vertex.
1453  CartVect coords[3];
1454  CartVect diff1( 1, 1, 1 );
1455  rval = tri_coords( moab, t_tri, coords );
1456  if( MB_SUCCESS == rval )
1457  {
1458  GeomUtil::closest_location_on_tri( points[i], coords, n_result );
1459  diff1 = n_result - t_result;
1460  }
1461  if( ( diff1 % diff1 ) > tolerance )
1462  {
1463  if( verbosity )
1464  std::cout << "Triangle closest to " << points[i] << " INCORRECT! (" << t_tri << " != " << n_tri
1465  << ")" << std::endl;
1466  result = false;
1467  continue;
1468  }
1469  }
1470 
1471  if( set_ptr && !tri_in_set( moab, *set_ptr, t_tri ) )
1472  {
1473  if( verbosity ) std::cout << "Surface closest to " << points[i] << " INCORRECT!" << std::endl;
1474  result = false;
1475  continue;
1476  }
1477  }
1478 
1479  if( verbosity > 1 )
1480  {
1481  std::cout << "Traversal statistics for closest point tests: " << std::endl;
1482  stats.print( std::cout );
1483  }
1484 
1485  return result;
1486 }