MOAB: Mesh Oriented datABase  (version 5.5.0)
point_location.cpp
Go to the documentation of this file.
1 #include "moab/Core.hpp"
2 #include "moab/AdaptiveKDTree.hpp"
3 #include "moab/Range.hpp"
4 #include "moab/CartVect.hpp"
5 #include "moab/GeomUtil.hpp"
6 #include <iostream>
7 #include <ctime>
8 #include <cstdlib>
9 #include <cassert>
10 #include <sstream>
11 
12 #define CHK( ErrorCode ) \
13  do \
14  { \
15  if( MB_SUCCESS != ( ErrorCode ) ) fail( ( ErrorCode ), __FILE__, __LINE__ ); \
16  } while( false )
17 
18 using namespace moab;
19 
20 static void fail( ErrorCode error_code, const char* file_name, int line_number );
21 
22 enum TreeType
23 {
24  UseKDTree,
25  UseNoTree,
26  UseDefaultTree = UseKDTree
27 };
28 
29 const char* default_str = "(default)";
30 const char* empty_str = "";
31 inline const char* is_default_tree( TreeType type )
32 {
33  return type == UseDefaultTree ? default_str : empty_str;
34 }
35 
36 const long DEFAULT_NUM_TEST = 100000;
37 const long HARD_MAX_UNIQUE_POINTS = 100000;
38 const long HARD_MIN_UNIQUE_POINTS = 1000;
39 const long FRACTION_UNIQUE_POINTS = 100;
40 
41 static void usage( char* argv0, bool help = false )
42 {
43  const char* usage_str = "[-k|-v] [-n <N>] [-d <N>] [-e <N>] <input_mesh>";
44  if( !help )
45  {
46  std::cerr << "Usage: " << argv0 << " " << usage_str << std::endl;
47  std::cerr << " : " << argv0 << " -h" << std::endl;
48  exit( 1 );
49  }
50 
51  std::cout << "Usage: " << argv0 << " " << usage_str << std::endl
52  << " -k : Use kD Tree " << is_default_tree( UseKDTree ) << std::endl
53  << " -v : Use no tree" << is_default_tree( UseNoTree ) << std::endl
54  << " -n : Specify number of test points (default = " << DEFAULT_NUM_TEST << ")" << std::endl
55  << " -d : Specify maximum tree depth" << std::endl
56  << " -e : Specify maximum elements per leaf" << std::endl
57  << " <input_mesh> : Mesh to build and query." << std::endl
58  << std::endl;
59  exit( 0 );
60 }
61 
62 static void generate_random_points( Interface& mesh,
63  size_t num_points,
64  std::vector< CartVect >& points,
65  std::vector< EntityHandle >& point_elems );
66 
67 static void do_kdtree_test( Interface& mesh,
68  int tree_depth,
69  int elem_per_leaf,
70  long num_test,
71  const std::vector< CartVect >& points,
72  std::vector< EntityHandle >& point_elems,
73  clock_t& build_time,
74  clock_t& test_time,
75  size_t& depth );
76 
77 static void do_linear_test( Interface& mesh,
78  int tree_depth,
79  int elem_per_leaf,
80  long num_test,
81  const std::vector< CartVect >& points,
82  std::vector< EntityHandle >& point_elems,
83  clock_t& build_time,
84  clock_t& test_time,
85  size_t& depth );
86 
87 int main( int argc, char* argv[] )
88 {
89 
90  const char* input_file = 0;
91  long tree_depth = -1;
92  long elem_per_leaf = -1;
93  long num_points = DEFAULT_NUM_TEST;
94  TreeType type = UseDefaultTree;
95  char* endptr;
96 
97  // PARSE ARGUMENTS
98 
99  long* valptr = 0;
100  for( int i = 1; i < argc; ++i )
101  {
102  if( valptr )
103  {
104  *valptr = strtol( argv[i], &endptr, 0 );
105  if( *valptr < 1 || *endptr )
106  {
107  std::cerr << "Invalid value for '" << argv[i - 1] << "' flag: " << argv[i] << std::endl;
108  exit( 1 );
109  }
110  valptr = 0;
111  }
112  else if( argv[i][0] == '-' && argv[i][1] && !argv[i][2] )
113  {
114  switch( argv[i][1] )
115  {
116  case 'h':
117  usage( argv[0], true );
118  break;
119  case 'k':
120  type = UseKDTree;
121  break;
122  case 'v':
123  type = UseNoTree;
124  break;
125  case 'd':
126  valptr = &tree_depth;
127  break;
128  case 'e':
129  valptr = &elem_per_leaf;
130  break;
131  case 'n':
132  valptr = &num_points;
133  break;
134  default:
135  std::cerr << "Invalid flag: " << argv[i] << std::endl;
136  usage( argv[0] );
137  break;
138  }
139  }
140  else if( !input_file )
141  {
142  input_file = argv[i];
143  }
144  else
145  {
146  std::cerr << "Unexpected argument: " << argv[i] << std::endl;
147  usage( argv[0] );
148  }
149  }
150  if( valptr )
151  {
152  std::cerr << "Expected value following '" << argv[argc - 1] << "' flag" << std::endl;
153  usage( argv[0] );
154  }
155  if( !input_file )
156  {
157  std::cerr << "No input file specified." << std::endl;
158  usage( argv[0] );
159  }
160 
161  // LOAD MESH
162 
163  Core moab;
164  Interface& mb = moab;
165  ErrorCode rval;
166  std::string init_msg, msg;
167  mb.get_last_error( init_msg );
168  rval = mb.load_file( input_file );
169  if( MB_SUCCESS != rval )
170  {
171  std::cerr << input_file << " : failed to read file." << std::endl;
172  mb.get_last_error( msg );
173  if( msg != init_msg ) std::cerr << "message : " << msg << std::endl;
174  }
175 
176  // GENERATE TEST POINTS
177 
178  int num_unique = num_points / FRACTION_UNIQUE_POINTS;
179  if( num_unique > HARD_MAX_UNIQUE_POINTS )
180  num_unique = HARD_MAX_UNIQUE_POINTS;
181  else if( num_unique < HARD_MIN_UNIQUE_POINTS )
182  num_unique = num_points;
183  std::vector< CartVect > points;
184  std::vector< EntityHandle > elems;
185  generate_random_points( mb, num_unique, points, elems );
186 
187  // GET MEMORY USE BEFORE BUILDING TREE
188 
189  unsigned long long init_total_storage;
190  mb.estimated_memory_use( 0, 0, &init_total_storage );
191 
192  // RUN TIMING TEST
193  clock_t build_time, test_time;
194  size_t actual_depth = 0;
195  std::vector< EntityHandle > results( points.size() );
196  switch( type )
197  {
198  default:
199  case UseKDTree:
200  do_kdtree_test( mb, tree_depth, elem_per_leaf, num_points, points, results, build_time, test_time,
201  actual_depth );
202  break;
203  case UseNoTree:
204  do_linear_test( mb, tree_depth, elem_per_leaf, num_points, points, results, build_time, test_time,
205  actual_depth );
206  break;
207  }
208 
209  unsigned long long fini_total_storage;
210  mb.estimated_memory_use( 0, 0, &fini_total_storage );
211 
212  // VALIDATE RESULTS
213  int fail = 0;
214  if( results != elems )
215  {
216  std::cout << "WARNING: Test produced invalid results" << std::endl;
217  fail = 1;
218  }
219 
220  // SUMMARIZE RESULTS
221  std::cout << "Number of test queries: " << num_points << std::endl;
222  std::cout << "Tree build time: " << ( (double)build_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
223  std::cout << "Total query time: " << ( (double)test_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
224  std::cout << "Time per query: " << ( (double)test_time ) / CLOCKS_PER_SEC / num_points << " seconds" << std::endl;
225  std::cout << "Tree depth: " << actual_depth << std::endl;
226  if( actual_depth )
227  std::cout << "Total query time/depth: " << ( (double)test_time ) / CLOCKS_PER_SEC / actual_depth << " seconds"
228  << std::endl;
229  std::cout << std::endl;
230  std::cout << "Bytes before tree construction: " << init_total_storage << std::endl;
231  std::cout << "Bytes after tree construction: " << fini_total_storage << std::endl;
232  std::cout << "Difference: " << fini_total_storage - init_total_storage << " bytes" << std::endl;
233  return fail;
234 }
235 
236 void fail( ErrorCode error_code, const char* file, int line )
237 {
238  std::cerr << "Internal error (error code " << error_code << ") at " << file << ":" << line << std::endl;
239  abort();
240 }
241 
242 const int HexSign[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
243  { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
244 
245 static CartVect random_point_in_hex( Interface& mb, EntityHandle hex )
246 {
247  const double f = RAND_MAX / 2;
248  CartVect xi( ( (double)rand() - f ) / f, ( (double)rand() - f ) / f, ( (double)rand() - f ) / f );
249  CartVect coords[8];
250  const EntityHandle* conn;
251  int len;
252  ErrorCode rval = mb.get_connectivity( hex, conn, len, true );
253  if( len != 8 && MB_SUCCESS != rval )
254  {
255  std::cerr << "Invalid element" << std::endl;
256  assert( false );
257  abort();
258  }
259  rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
260  CHK( rval );
261 
262  CartVect point( 0, 0, 0 );
263  for( unsigned i = 0; i < 8; ++i )
264  {
265  double coeff = 0.125;
266  for( unsigned j = 0; j < 3; ++j )
267  coeff *= 1 + HexSign[i][j] * xi[j];
268  point += coeff * coords[i];
269  }
270 
271  return point;
272 }
273 
274 void generate_random_points( Interface& mb,
275  size_t num_points,
276  std::vector< CartVect >& points,
277  std::vector< EntityHandle >& point_elems )
278 {
279  Range elems;
280  ErrorCode rval;
281  rval = mb.get_entities_by_dimension( 0, 3, elems );
282  CHK( rval );
283  if( !elems.all_of_type( MBHEX ) )
284  {
285  std::cerr << "Warning: ignoring non-hexahedral elements." << std::endl;
286  std::pair< Range::iterator, Range::iterator > p = elems.equal_range( MBHEX );
287  elems.erase( p.second, elems.end() );
288  elems.erase( elems.begin(), p.first );
289  }
290  if( elems.empty() )
291  {
292  std::cerr << "Input file contains no hexahedral elements." << std::endl;
293  exit( 1 );
294  }
295 
296  points.resize( num_points );
297  point_elems.resize( num_points );
298  const size_t num_elem = elems.size();
299  for( size_t i = 0; i < num_points; ++i )
300  {
301  size_t offset = 0;
302  for( size_t x = num_elem; x > 0; x /= RAND_MAX )
303  offset += rand();
304  offset %= num_elem;
305  point_elems[i] = elems[offset];
306  points[i] = random_point_in_hex( mb, point_elems[i] );
307  }
308 }
309 
310 void do_kdtree_test( Interface& mb,
311  int tree_depth,
312  int elem_per_leaf,
313  long num_test,
314  const std::vector< CartVect >& points,
315  std::vector< EntityHandle >& point_elems,
316  clock_t& build_time,
317  clock_t& test_time,
318  size_t& depth )
319 {
320  ErrorCode rval;
321  clock_t init = clock();
322  AdaptiveKDTree tool( &mb );
323  EntityHandle root;
324  std::ostringstream options;
325  if( tree_depth > 0 ) options << "MAX_DEPTH=" << tree_depth << ";";
326  if( elem_per_leaf > 0 ) options << "MAX_PER_LEAF=" << elem_per_leaf << ";";
327  Range all_hexes;
328  rval = mb.get_entities_by_type( 0, MBHEX, all_hexes );
329  CHK( rval );
330  FileOptions opt( options.str().c_str() );
331  rval = tool.build_tree( all_hexes, &root, &opt );
332  CHK( rval );
333  all_hexes.clear();
334  build_time = clock() - init;
335 
336  EntityHandle leaf;
337  std::vector< EntityHandle > hexes;
338  std::vector< EntityHandle >::iterator j;
339  CartVect coords[8];
340  const EntityHandle* conn;
341  int len;
342  for( long i = 0; i < num_test; ++i )
343  {
344  const size_t idx = (size_t)i % points.size();
345  rval = tool.point_search( points[idx].array(), leaf, 1.0e-10, 1.0e-6, NULL, &root );
346  CHK( rval );
347  hexes.clear();
348  rval = mb.get_entities_by_handle( leaf, hexes );
349  CHK( rval );
350  for( j = hexes.begin(); j != hexes.end(); ++j )
351  {
352  rval = mb.get_connectivity( *j, conn, len, true );
353  CHK( rval );
354  rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
355  CHK( rval );
356  if( GeomUtil::point_in_trilinear_hex( coords, points[idx], 1e-12 ) )
357  {
358  point_elems[idx] = *j;
359  break;
360  }
361  }
362  if( j == hexes.end() ) point_elems[idx] = 0;
363  }
364 
365  test_time = clock() - build_time;
366  double tmp_box[3];
367  unsigned max_d;
368  tool.get_info( root, tmp_box, tmp_box, max_d );
369  depth = max_d;
370 }
371 
372 void do_linear_test( Interface& mb,
373  int,
374  int,
375  long num_test,
376  const std::vector< CartVect >& points,
377  std::vector< EntityHandle >& point_elems,
378  clock_t& build_time,
379  clock_t& test_time,
380  size_t& depth )
381 {
382  clock_t init = clock();
383  Range hexes;
384  ErrorCode rval = mb.get_entities_by_type( 0, MBHEX, hexes );
385  CHK( rval );
386  depth = 0;
387  point_elems.resize( points.size() );
388  build_time = clock() - init;
389 
390  CartVect coords[8];
391  const EntityHandle* conn;
392  int len;
393  for( long i = 0; i < num_test; ++i )
394  {
395  const size_t idx = (size_t)i % points.size();
396  for( Range::iterator h = hexes.begin(); h != hexes.end(); ++h )
397  {
398  rval = mb.get_connectivity( *h, conn, len, true );
399  CHK( rval );
400  rval = mb.get_coords( conn, 8, reinterpret_cast< double* >( coords ) );
401  CHK( rval );
402  if( GeomUtil::point_in_trilinear_hex( coords, points[idx], 1e-12 ) )
403  {
404  point_elems[idx] = *h;
405  break;
406  }
407  }
408  }
409 
410  test_time = clock() - build_time;
411 }