MOAB: Mesh Oriented datABase  (version 5.5.0)
adaptive_kd_tree_tests.cpp
Go to the documentation of this file.
1 #include "moab/Core.hpp"
2 #include "moab/AdaptiveKDTree.hpp"
3 #include "moab/CartVect.hpp"
4 #include "moab/GeomUtil.hpp"
5 #include "moab/Range.hpp"
6 #include "TestUtil.hpp"
7 
8 using namespace moab;
9 
10 #include <iostream>
11 #include <algorithm>
12 #include <sstream>
13 
14 #ifdef MOAB_HAVE_MPI
15 #include "moab_mpi.h"
16 #endif
17 
18 /* Utility method - compare two range boxes */
19 bool box_equal( const AdaptiveKDTreeIter& iter,
20  double x_min,
21  double y_min,
22  double z_min,
23  double x_max,
24  double y_max,
25  double z_max )
26 {
27  return iter.box_min()[0] == x_min && iter.box_min()[1] == y_min && iter.box_min()[2] == z_min &&
28  iter.box_max()[0] == x_max && iter.box_max()[1] == y_max && iter.box_max()[2] == z_max;
29 }
30 
31 void build_triangles( Interface* moab,
32  Range& tris,
33  int num_vert,
34  const double* coords,
35  int num_tri,
36  const unsigned* conn )
37 {
38  std::vector< EntityHandle > verts( num_vert );
39  for( int i = 0; i < num_vert; ++i )
40  moab->create_vertex( coords + 3 * i, verts[i] );
41 
42  EntityHandle tri, tri_verts[3];
43  for( int i = 0; i < num_tri; ++i )
44  {
45  tri_verts[0] = verts[conn[3 * i]];
46  tri_verts[1] = verts[conn[3 * i + 1]];
47  tri_verts[2] = verts[conn[3 * i + 2]];
48  moab->create_element( MBTRI, tri_verts, 3, tri );
49  tris.insert( tri );
50  }
51 }
52 
53 /* Utility method - build a 2x2x2 box composed of two triagles on a side */
54 void build_triangle_box_small( Interface* moab, Range& tris )
55 {
56  const double coords[] = { -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1 };
57 
58  const unsigned conn[] = { 0, 1, 5, 0, 5, 4, 2, 6, 7, 2, 7, 3, 1, 2, 6, 1, 6, 5,
59  0, 3, 4, 3, 7, 4, 0, 1, 3, 1, 2, 3, 4, 5, 6, 4, 6, 7 };
60 
61  build_triangles( moab, tris, 8, coords, 12, conn );
62 }
63 
64 /* build 6x6x6 box centered at origin and composed of 18 triangles per side,
65  Each side of the box is composed of a 3x3 grid of quads, where each
66  quad is split diagonally to form two triangles. */
67 void build_triangle_box_large( Interface* moab, Range& tris )
68 {
69  const double coords[] = { // corners
70  -3, -3, -3, 3, -3, -3, 3, 3, -3, -3, 3, -3, -3, -3, 3, 3, -3, 3, 3, 3, 3, -3, 3, 3,
71 
72  // edges
73  -1, -3, -3, 1, -3, -3, 3, -1, -3, 3, 1, -3, 1, 3, -3, -1, 3, -3, -3, 1, -3, -3, -1, -3,
74 
75  -1, -3, 3, 1, -3, 3, 3, -1, 3, 3, 1, 3, 1, 3, 3, -1, 3, 3, -3, 1, 3, -3, -1, 3,
76 
77  -3, -3, -1, -3, -3, 1, 3, -3, -1, 3, -3, 1, 3, 3, -1, 3, 3, 1, -3, 3, -1, -3, 3, 1,
78 
79  // faces
80  -1, -3, -1, 1, -3, -1, 1, -3, 1, -1, -3, 1,
81 
82  3, -1, -1, 3, 1, -1, 3, 1, 1, 3, -1, 1,
83 
84  1, 3, -1, -1, 3, -1, -1, 3, 1, 1, 3, 1,
85 
86  -3, 1, -1, -3, -1, -1, -3, -1, 1, -3, 1, 1,
87 
88  -1, -1, -3, 1, -1, -3, 1, 1, -3, -1, 1, -3,
89 
90  -1, -1, 3, 1, -1, 3, 1, 1, 3, -1, 1, 3 };
91  const unsigned conn[] = {
92  // face 0
93  0, 8, 24, 8, 32, 24, 8, 33, 32, 8, 9, 33, 9, 1, 33, 1, 26, 33, 24, 35, 25, 24, 32, 35, 32, 33, 35, 33, 34, 35,
94  33, 27, 34, 33, 26, 27, 35, 4, 25, 35, 16, 4, 35, 17, 16, 35, 34, 17, 27, 17, 34, 27, 5, 17,
95  // face 1
96  36, 26, 1, 36, 1, 10, 36, 10, 11, 36, 11, 37, 11, 28, 37, 11, 2, 28, 36, 27, 26, 36, 39, 27, 36, 38, 39, 36, 37,
97  38, 37, 28, 38, 28, 29, 38, 18, 5, 27, 18, 27, 39, 18, 39, 38, 18, 38, 19, 6, 19, 38, 6, 38, 29,
98  // face 2
99  12, 28, 2, 12, 40, 28, 12, 41, 40, 12, 13, 41, 3, 41, 13, 3, 30, 41, 43, 29, 28, 43, 28, 40, 43, 40, 41, 43, 41,
100  42, 41, 31, 42, 41, 30, 31, 43, 6, 29, 43, 20, 6, 43, 21, 20, 43, 42, 21, 21, 42, 31, 21, 31, 7,
101  // face 3
102  44, 30, 3, 44, 3, 14, 44, 14, 15, 44, 15, 45, 15, 24, 45, 15, 0, 24, 44, 31, 30, 44, 47, 31, 44, 46, 47, 44, 45,
103  46, 46, 45, 24, 46, 24, 25, 31, 22, 7, 31, 47, 22, 46, 22, 47, 46, 23, 22, 46, 4, 23, 46, 25, 4,
104  // face 4
105  8, 15, 0, 8, 48, 15, 8, 49, 48, 8, 9, 49, 1, 49, 9, 1, 10, 49, 51, 14, 15, 51, 15, 48, 51, 48, 49, 51, 49, 50,
106  11, 50, 49, 11, 49, 10, 51, 3, 14, 51, 13, 3, 51, 12, 13, 51, 50, 12, 11, 12, 50, 11, 2, 12,
107  // face 5
108  4, 52, 16, 4, 23, 52, 22, 52, 23, 22, 55, 52, 22, 21, 55, 22, 7, 21, 17, 16, 52, 17, 52, 53, 54, 53, 52, 54, 52,
109  55, 54, 55, 21, 54, 21, 20, 18, 5, 17, 18, 17, 53, 18, 53, 54, 18, 54, 19, 6, 19, 54, 6, 54, 20 };
110 
111  build_triangles( moab, tris, 56, coords, 108, conn );
112 }
113 
114 /* Utility method - build 2x2x2 octahedron (3D diamond)*/
115 void build_triangle_octahedron( Interface* moab, Range& tris )
116 {
117  const double coords[] = { 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1 };
118  const unsigned conn[] = { 0, 1, 4, 1, 2, 4, 2, 3, 4, 3, 0, 4, 1, 0, 5, 2, 1, 5, 3, 2, 5, 0, 3, 5 };
119 
120  build_triangles( moab, tris, 6, coords, 8, conn );
121 }
122 
123 void test_valid_tree( AdaptiveKDTree* tool, EntityHandle root, FileOptions& options, const Range& expected_tris )
124 {
125  Range all_tris;
126  ErrorCode rval;
127  AdaptiveKDTreeIter iter;
128  CHECK( MB_SUCCESS == tool->get_tree_iterator( root, iter ) );
129  do
130  {
131  double dep;
132  rval = options.get_real_option( "MAX_DEPTH", dep );
133  CHECK( MB_ENTITY_NOT_FOUND == rval || iter.depth() <= tool->get_max_depth() );
134 
135  Range tris;
136  CHECK( tool->moab()->get_entities_by_type( iter.handle(), MBTRI, tris ) == MB_SUCCESS );
137  // CHECK( !tris.empty() );
138  all_tris.merge( tris );
139 
140  const CartVect min( iter.box_min() ), max( iter.box_max() );
141  const CartVect cen( 0.5 * ( min + max ) ), hdim( 0.5 * ( max - min ) );
142 
143  for( Range::iterator i = tris.begin(); i != tris.end(); ++i )
144  {
145  EntityHandle tri = *i;
146  const EntityHandle* conn;
147  int conn_len;
148  CHECK( tool->moab()->get_connectivity( tri, conn, conn_len ) == MB_SUCCESS );
149  CHECK( conn_len == 3 );
150  CartVect coords[3];
151  CHECK( tool->moab()->get_coords( conn, 3, coords[0].array() ) == MB_SUCCESS );
152  CHECK( GeomUtil::box_tri_overlap( coords, cen, hdim ) );
153  }
154  } while( MB_SUCCESS == ( rval = iter.step() ) );
155 
156  CHECK( MB_ENTITY_NOT_FOUND == rval );
157  CHECK( expected_tris == all_tris );
158 }
159 
160 /* utility method - check that all tris share a vertex */
161 void check_common_vertex( Interface* moab, const EntityHandle* tris, unsigned num_tri, CartVect point )
162 {
163  for( unsigned i = 0; i < num_tri; ++i )
164  CHECK( MBTRI == moab->type_from_handle( tris[i] ) );
165 
166  ErrorCode rval;
167  CartVect tri_coords[3];
168  const EntityHandle* conn;
169  int conn_len;
170  rval = moab->get_connectivity( tris[0], conn, conn_len );
171  CHECK( MB_SUCCESS == rval );
172  rval = moab->get_coords( conn, 3, tri_coords[0].array() );
173  CHECK( MB_SUCCESS == rval );
174  tri_coords[0] -= point;
175  tri_coords[1] -= point;
176  tri_coords[2] -= point;
177  EntityHandle vertex = 0;
178  if( tri_coords[0] % tri_coords[0] < 1e-6 )
179  vertex = conn[0];
180  else if( tri_coords[1] % tri_coords[1] < 1e-6 )
181  vertex = conn[1];
182  else if( tri_coords[2] % tri_coords[2] < 1e-6 )
183  vertex = conn[2];
184  else
185  CHECK( false );
186  for( unsigned j = 1; j < num_tri; ++j )
187  {
188  rval = moab->get_connectivity( tris[j], conn, conn_len );
189  CHECK( conn[0] == vertex || conn[1] == vertex || conn[2] == vertex );
190  }
191 }
192 
193 /* utility method - check that all tris share a vertex */
194 void check_point_in_triangles( Interface* moab, const EntityHandle* tris, unsigned num_tris, CartVect point )
195 {
196  ErrorCode rval;
197  CartVect tri_coords[3], tript;
198  const EntityHandle* conn;
199  int conn_len;
200 
201  for( unsigned i = 0; i < num_tris; ++i )
202  {
203  CHECK( MBTRI == moab->type_from_handle( tris[i] ) );
204 
205  rval = moab->get_connectivity( tris[i], conn, conn_len );
206  CHECK( MB_SUCCESS == rval );
207  rval = moab->get_coords( conn, 3, tri_coords[0].array() );
208  CHECK( MB_SUCCESS == rval );
209 
210  GeomUtil::closest_location_on_tri( point, tri_coords, tript );
211  tript -= point;
212  CHECK( fabs( tript[0] ) < 1e-6 );
213  CHECK( fabs( tript[1] ) < 1e-6 );
214  CHECK( fabs( tript[2] ) < 1e-6 );
215  }
216 }
217 
218 /* Create the following 2D tree (no Z splits)
219 
220  (6) X = -3 (5) X = 0
221  / |
222  / | [8]
223  / |---------------- (7) Y = 3
224  [5] / [6] |
225  / | [7]
226  / |
227  / |
228  -------------------------------------- (1) Y = 0
229  | |
230  [2] | |
231 (3)----------------| |
232 Y = -2 | [3] | [4]
233  | |
234  [1] | |
235  | |
236  (2) X = -1 (4) X = 4
237  */
238 void create_simple_2d_tree( AdaptiveKDTree& tool, EntityHandle& root, EntityHandle leaves[9] = 0 )
239 {
240  ErrorCode rval;
241  AdaptiveKDTree::Plane plane;
242 
243  // create a single-node tree
244  double min[3] = { -5, -4, -1 };
245  double max[3] = { 5, 4, 1 };
246  rval = tool.create_root( min, max, root );
247  CHECK( MB_SUCCESS == rval );
248 
249  // get iterator for tree
250  AdaptiveKDTreeIter iter;
251  rval = tool.get_tree_iterator( root, iter );
252  CHECK( MB_SUCCESS == rval );
253  CHECK( box_equal( iter, -5, -4, -1, 5, 4, 1 ) );
254 
255  // split plane (1)
256  plane.norm = AdaptiveKDTree::Y;
257  plane.coord = 0.0;
258  rval = tool.split_leaf( iter, plane );
259  CHECK( MB_SUCCESS == rval );
260  CHECK( box_equal( iter, -5, -4, -1, 5, 0, 1 ) );
261 
262  // split plane (2)
263  plane.norm = AdaptiveKDTree::X;
264  plane.coord = -1.0;
265  rval = tool.split_leaf( iter, plane );
266  CHECK( MB_SUCCESS == rval );
267  CHECK( box_equal( iter, -5, -4, -1, -1, 0, 1 ) );
268 
269  // split plane (3), leaf [1]
270  plane.norm = AdaptiveKDTree::Y;
271  plane.coord = -2.0;
272  rval = tool.split_leaf( iter, plane );
273  CHECK( MB_SUCCESS == rval );
274  CHECK( box_equal( iter, -5, -4, -1, -1, -2, 1 ) );
275  if( leaves ) leaves[1] = iter.handle();
276 
277  // leaf [2]
278  rval = iter.step();
279  CHECK( MB_SUCCESS == rval );
280  CHECK( box_equal( iter, -5, -2, -1, -1, 0, 1 ) );
281  if( leaves ) leaves[2] = iter.handle();
282 
283  // non-leaf right of split plane (2)
284  rval = iter.step();
285  CHECK( MB_SUCCESS == rval );
286  CHECK( box_equal( iter, -1, -4, -1, 5, 0, 1 ) );
287 
288  // split plane (4) and leaf [3]
289  plane.norm = AdaptiveKDTree::X;
290  plane.coord = 4.0;
291  rval = tool.split_leaf( iter, plane );
292  CHECK( MB_SUCCESS == rval );
293  CHECK( box_equal( iter, -1, -4, -1, 4, 0, 1 ) );
294  if( leaves ) leaves[3] = iter.handle();
295 
296  // leaf [4]
297  rval = iter.step();
298  CHECK( MB_SUCCESS == rval );
299  CHECK( box_equal( iter, 4, -4, -1, 5, 0, 1 ) );
300  if( leaves ) leaves[4] = iter.handle();
301 
302  // non-leaf above split plane (1)
303  rval = iter.step();
304  CHECK( MB_SUCCESS == rval );
305  CHECK( box_equal( iter, -5, 0, -1, 5, 4, 1 ) );
306 
307  // split plane (5)
308  plane.norm = AdaptiveKDTree::X;
309  plane.coord = 0.0;
310  rval = tool.split_leaf( iter, plane );
311  CHECK( MB_SUCCESS == rval );
312  CHECK( box_equal( iter, -5, 0, -1, 0, 4, 1 ) );
313 
314  // split plane (6) and leaf [5]
315  plane.norm = AdaptiveKDTree::X;
316  plane.coord = -3.0;
317  rval = tool.split_leaf( iter, plane );
318  CHECK( MB_SUCCESS == rval );
319  CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) );
320  if( leaves ) leaves[5] = iter.handle();
321 
322  // leaf [6];
323  rval = iter.step();
324  CHECK( MB_SUCCESS == rval );
325  CHECK( box_equal( iter, -3, 0, -1, 0, 4, 1 ) );
326  if( leaves ) leaves[6] = iter.handle();
327 
328  // non-leaf right of split plane (5)
329  rval = iter.step();
330  CHECK( MB_SUCCESS == rval );
331  CHECK( box_equal( iter, 0, 0, -1, 5, 4, 1 ) );
332 
333  // split plane (7) and leaf [7]
334  plane.norm = AdaptiveKDTree::Y;
335  plane.coord = 3.0;
336  rval = tool.split_leaf( iter, plane );
337  CHECK( MB_SUCCESS == rval );
338  CHECK( box_equal( iter, 0, 0, -1, 5, 3, 1 ) );
339  if( leaves ) leaves[7] = iter.handle();
340 
341  // leaf [8];
342  rval = iter.step();
343  CHECK( MB_SUCCESS == rval );
344  CHECK( box_equal( iter, 0, 3, -1, 5, 4, 1 ) );
345  if( leaves ) leaves[8] = iter.handle();
346 
347  // the end
348  rval = iter.step();
349  CHECK( MB_ENTITY_NOT_FOUND == rval );
350 }
351 
352 void leaf_iterator_test()
353 {
354  ErrorCode rval;
355  Core moab;
356  Interface* mb = &moab;
357  AdaptiveKDTree tool( mb );
358 
359  // create a single-node tree
360  EntityHandle root;
361  double min[3] = { -3, -2, -1 };
362  double max[3] = { 1, 2, 3 };
363  rval = tool.create_root( min, max, root );
364  CHECK( MB_SUCCESS == rval );
365 
366  // get iterator for tree
367  AdaptiveKDTreeIter iter;
368  rval = tool.get_tree_iterator( root, iter );
369  CHECK( MB_SUCCESS == rval );
370  CHECK( box_equal( iter, -3, -2, -1, 1, 2, 3 ) );
371 
372  // check that steping the iterator fails correctly.
373  rval = iter.step();
374  CHECK( MB_ENTITY_NOT_FOUND == rval );
375  rval = iter.step();
376  CHECK( MB_SUCCESS != rval );
377 
378  rval = tool.reset_tree();
379  CHECK( MB_SUCCESS == rval );
380  root = 0;
381 
382  // construct a simple 2D tree for remaining tests
383  EntityHandle leaves[9];
384  create_simple_2d_tree( tool, root, leaves );
385 
386  /**** Now traverse tree again, and check neighbors of each leaf *****/
387  std::vector< AdaptiveKDTreeIter > list;
388 
389  // leaf [1]
390 
391  rval = tool.get_tree_iterator( root, iter );
392  CHECK( MB_SUCCESS == rval );
393  CHECK( box_equal( iter, -5, -4, -1, -1, -2, 1 ) );
394  CHECK( iter.handle() == leaves[1] );
395 
396  list.clear();
397  rval = iter.get_neighbors( AdaptiveKDTree::X, true, list );
398  CHECK( MB_SUCCESS == rval );
399  CHECK( list.empty() );
400 
401  list.clear();
402  rval = iter.get_neighbors( AdaptiveKDTree::X, false, list );
403  CHECK( MB_SUCCESS == rval );
404  CHECK( list.size() == 1 );
405  // should be leaf [3]
406  CHECK( box_equal( list.front(), -1, -4, -1, 4, 0, 1 ) );
407  CHECK( list.front().handle() == leaves[3] );
408 
409  list.clear();
410  rval = iter.get_neighbors( AdaptiveKDTree::Y, true, list );
411  CHECK( MB_SUCCESS == rval );
412  CHECK( list.empty() );
413 
414  list.clear();
415  rval = iter.get_neighbors( AdaptiveKDTree::Y, false, list );
416  CHECK( MB_SUCCESS == rval );
417  CHECK( list.size() == 1 );
418  // should be leaf [2]
419  CHECK( box_equal( list.front(), -5, -2, -1, -1, 0, 1 ) );
420  CHECK( list.front().handle() == leaves[2] );
421 
422  list.clear();
423  rval = iter.get_neighbors( AdaptiveKDTree::Z, true, list );
424  CHECK( MB_SUCCESS == rval );
425  CHECK( list.empty() );
426 
427  list.clear();
428  rval = iter.get_neighbors( AdaptiveKDTree::Z, false, list );
429  CHECK( MB_SUCCESS == rval );
430  CHECK( list.empty() );
431 
432  // leaf [2]
433 
434  rval = iter.step();
435  CHECK( MB_SUCCESS == rval );
436  CHECK( box_equal( iter, -5, -2, -1, -1, 0, 1 ) );
437  CHECK( iter.handle() == leaves[2] );
438 
439  list.clear();
440  rval = iter.get_neighbors( AdaptiveKDTree::X, true, list );
441  CHECK( MB_SUCCESS == rval );
442  CHECK( list.empty() );
443 
444  list.clear();
445  rval = iter.get_neighbors( AdaptiveKDTree::X, false, list );
446  CHECK( MB_SUCCESS == rval );
447  CHECK( list.size() == 1 );
448  // should be leaf [3]
449  CHECK( box_equal( list.front(), -1, -4, -1, 4, 0, 1 ) );
450  CHECK( list.front().handle() == leaves[3] );
451 
452  list.clear();
453  rval = iter.get_neighbors( AdaptiveKDTree::Y, true, list );
454  CHECK( MB_SUCCESS == rval );
455  // should be leaf [1]
456  CHECK( list.size() == 1 );
457  CHECK( box_equal( list.front(), -5, -4, -1, -1, -2, 1 ) );
458  CHECK( list.front().handle() == leaves[1] );
459 
460  list.clear();
461  rval = iter.get_neighbors( AdaptiveKDTree::Y, false, list );
462  CHECK( MB_SUCCESS == rval );
463  CHECK( list.size() == 2 );
464  // should be leaf [5] and leaf[6]
465  if( list[0].handle() == leaves[6] ) std::swap( list[0], list[1] );
466  CHECK( list[0].handle() == leaves[5] );
467  CHECK( list[1].handle() == leaves[6] );
468  CHECK( box_equal( list[0], -5, 0, -1, -3, 4, 1 ) );
469  CHECK( box_equal( list[1], -3, 0, -1, 0, 4, 1 ) );
470 
471  list.clear();
472  rval = iter.get_neighbors( AdaptiveKDTree::Z, true, list );
473  CHECK( MB_SUCCESS == rval );
474  CHECK( list.empty() );
475 
476  list.clear();
477  rval = iter.get_neighbors( AdaptiveKDTree::Z, false, list );
478  CHECK( MB_SUCCESS == rval );
479  CHECK( list.empty() );
480 
481  // leaf [3]
482  rval = iter.step();
483  CHECK( MB_SUCCESS == rval );
484  CHECK( box_equal( iter, -1, -4, -1, 4, 0, 1 ) );
485  CHECK( iter.handle() == leaves[3] );
486 
487  // leaf [4]
488  rval = iter.step();
489  CHECK( MB_SUCCESS == rval );
490  CHECK( box_equal( iter, 4, -4, -1, 5, 0, 1 ) );
491  CHECK( iter.handle() == leaves[4] );
492 
493  // leaf [5]
494  rval = iter.step();
495  CHECK( MB_SUCCESS == rval );
496  CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) );
497  CHECK( iter.handle() == leaves[5] );
498 
499  // leaf [6]
500  rval = iter.step();
501  CHECK( MB_SUCCESS == rval );
502  CHECK( box_equal( iter, -3, 0, -1, 0, 4, 1 ) );
503  CHECK( iter.handle() == leaves[6] );
504 
505  // leaf [7]
506  rval = iter.step();
507  CHECK( MB_SUCCESS == rval );
508  CHECK( box_equal( iter, 0, 0, -1, 5, 3, 1 ) );
509  CHECK( iter.handle() == leaves[7] );
510 
511  // leaf [8]
512  rval = iter.step();
513  CHECK( MB_SUCCESS == rval );
514  CHECK( box_equal( iter, 0, 3, -1, 5, 4, 1 ) );
515  CHECK( iter.handle() == leaves[8] );
516 }
517 
518 void test_tree_merge_nodes()
519 {
520  // build simple tree for tests
521  ErrorCode rval;
522  Core moab;
523  Interface* mb = &moab;
524  AdaptiveKDTree tool( mb );
525  EntityHandle root;
526  EntityHandle leaves[9];
527  create_simple_2d_tree( tool, root, leaves );
528 
529  // get iterator for tree
530  AdaptiveKDTreeIter iter;
531  rval = tool.get_tree_iterator( root, iter );
532  CHECK( MB_SUCCESS == rval );
533 
534  // merge leaves 1 and 2
535  rval = tool.merge_leaf( iter );
536  CHECK( MB_SUCCESS == rval );
537  CHECK( box_equal( iter, -5, -4, -1, -1, 0, 1 ) );
538 
539  // merge leaf 1,2 with 3,4 (implicity merges 3 and 4)
540  rval = tool.merge_leaf( iter );
541  CHECK( MB_SUCCESS == rval );
542  CHECK( box_equal( iter, -5, -4, -1, 5, 0, 1 ) );
543 
544  // make sure iterator remains valid
545  rval = iter.step();
546  CHECK( MB_SUCCESS == rval );
547  // leaf 5
548  CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) );
549  CHECK( iter.handle() == leaves[5] );
550 }
551 
552 void test_build_tree_bisect_triangles()
553 {
554  Core moab;
555  AdaptiveKDTree tool( &moab );
556  Range box_tris;
557 
558  std::ostringstream options;
559  options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
560  FileOptions opts( options.str().c_str() );
561  EntityHandle root;
562  ErrorCode rval;
563 
564  moab.delete_mesh();
565  box_tris.clear();
566  build_triangle_box_small( &moab, box_tris );
567  rval = tool.build_tree( box_tris, &root, &opts );
568  CHECK( MB_SUCCESS == rval );
569  test_valid_tree( &tool, root, opts, box_tris );
570  tool.reset_tree();
571 
572  moab.delete_mesh();
573  box_tris.clear();
574  build_triangle_octahedron( &moab, box_tris );
575  rval = tool.build_tree( box_tris, &root, &opts );
576  CHECK( MB_SUCCESS == rval );
577  test_valid_tree( &tool, root, opts, box_tris );
578  tool.reset_tree();
579 
580  moab.delete_mesh();
581  box_tris.clear();
582  build_triangle_box_large( &moab, box_tris );
583  rval = tool.build_tree( box_tris, &root, &opts );
584  CHECK( MB_SUCCESS == rval );
585  test_valid_tree( &tool, root, opts, box_tris );
586  tool.reset_tree();
587 
588  options << "MAX_DEPTH=2;";
589  opts = FileOptions( options.str().c_str() );
590  build_triangle_box_large( &moab, box_tris );
591  rval = tool.build_tree( box_tris, &root, &opts );
592  CHECK( MB_SUCCESS == rval );
593  test_valid_tree( &tool, root, opts, box_tris );
594  tool.reset_tree();
595 }
596 
597 void test_closest_triangle()
598 {
599  Core moab;
600  AdaptiveKDTree tool( &moab );
601  Range box_tris;
602 
603  std::ostringstream options;
604  options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
605  FileOptions opts( options.str().c_str() );
606  EntityHandle root;
607  ErrorCode rval;
608  EntityHandle tri;
609 
610  moab.delete_mesh();
611  box_tris.clear();
612  build_triangle_box_small( &moab, box_tris );
613  rval = tool.build_tree( box_tris, &root, &opts );
614  CHECK( MB_SUCCESS == rval );
615  test_valid_tree( &tool, root, opts, box_tris );
616 
617  // test closest to each corner of the box.
618  for( unsigned i = 0; i < 8; ++i )
619  {
620  int x = i & 1 ? 1 : -1;
621  int y = i & 2 ? 1 : -1;
622  int z = i & 4 ? 1 : -1;
623  double coords[] = { 2.0 * x, 2.0 * y, 2.0 * z };
624  CartVect point;
625  rval = tool.closest_triangle( root, coords, point.array(), tri );
626  CHECK( MB_SUCCESS == rval );
627  // check result position
628  CHECK( fabs( x - point[0] ) < 1e-6 );
629  CHECK( fabs( y - point[1] ) < 1e-6 );
630  CHECK( fabs( z - point[2] ) < 1e-6 );
631  // check point is at triangle vertex
632  check_common_vertex( tool.moab(), &tri, 1, point );
633  }
634 
635  // test location on each face
636  const CartVect facepts[] = { CartVect( -1.0, 0.5, 0.5 ), CartVect( 1.0, 0.5, 0.5 ), CartVect( 0.5, -1.0, 0.5 ),
637  CartVect( 0.5, 1.0, 0.5 ), CartVect( 0.5, 0.5, -1.0 ), CartVect( 0.5, 0.5, 1.0 ) };
638  for( unsigned i = 0; i < 6; ++i )
639  {
640  CartVect point;
641  rval = tool.closest_triangle( root, facepts[i].array(), point.array(), tri );
642  CHECK( MB_SUCCESS == rval );
643  // check result position
644  const CartVect diff = facepts[i] - point;
645  CHECK( fabs( diff[0] ) < 1e-6 );
646  CHECK( fabs( diff[1] ) < 1e-6 );
647  CHECK( fabs( diff[2] ) < 1e-6 );
648  // check that point is contained within triangle
649  check_point_in_triangles( tool.moab(), &tri, 1, point );
650  }
651 
652  // test a location really far from the tree
653  {
654  const double far[] = { 0.75, 0.75, 200 };
655  CartVect point;
656  rval = tool.closest_triangle( root, far, point.array(), tri );
657  CHECK( MB_SUCCESS == rval );
658  CHECK( fabs( point[0] - 0.75 ) < 1e-6 );
659  CHECK( fabs( point[1] - 0.75 ) < 1e-6 );
660  CHECK( fabs( point[2] - 1.00 ) < 1e-6 );
661  // check that point is contained within triangle
662  check_point_in_triangles( tool.moab(), &tri, 1, point );
663  }
664 
665  // now do it all again with a lot more triangles
666  tool.reset_tree();
667  moab.delete_mesh();
668  box_tris.clear();
669 
670  build_triangle_box_large( &moab, box_tris );
671  rval = tool.build_tree( box_tris, &root, &opts );
672  CHECK( MB_SUCCESS == rval );
673  test_valid_tree( &tool, root, opts, box_tris );
674 
675  // test closest to each corner of the box.
676  for( unsigned i = 0; i < 8; ++i )
677  {
678  int x = i & 1 ? 1 : -1;
679  int y = i & 2 ? 1 : -1;
680  int z = i & 4 ? 1 : -1;
681  double coords[] = { 4.0 * x, 4.0 * y, 4.0 * z };
682  CartVect point;
683  rval = tool.closest_triangle( root, coords, point.array(), tri );
684  CHECK( MB_SUCCESS == rval );
685  // check result position
686  CHECK( fabs( 3.0 * x - point[0] ) < 1e-6 );
687  CHECK( fabs( 3.0 * y - point[1] ) < 1e-6 );
688  CHECK( fabs( 3.0 * z - point[2] ) < 1e-6 );
689  // check point is at triangle vertex
690  check_common_vertex( tool.moab(), &tri, 1, point );
691  }
692 
693  // test location on each face
694  const CartVect facepts2[] = { CartVect( -3.0, 0.5, 0.5 ), CartVect( 3.0, 0.5, 0.5 ), CartVect( 0.5, -3.0, 0.5 ),
695  CartVect( 0.5, 3.0, 0.5 ), CartVect( 0.5, 0.5, -3.0 ), CartVect( 0.5, 0.5, 3.0 ) };
696  for( unsigned i = 0; i < 6; ++i )
697  {
698  CartVect point;
699  rval = tool.closest_triangle( root, facepts2[i].array(), point.array(), tri );
700  CHECK( MB_SUCCESS == rval );
701  // check result position
702  const CartVect diff = facepts2[i] - point;
703  CHECK( fabs( diff[0] ) < 1e-6 );
704  CHECK( fabs( diff[1] ) < 1e-6 );
705  CHECK( fabs( diff[2] ) < 1e-6 );
706  // check that point is contained within triangle
707  check_point_in_triangles( tool.moab(), &tri, 1, point );
708  }
709 
710  // test a location really far from the tree
711  {
712  const double far[] = { 2.75, 2.75, 200 };
713  CartVect point;
714  rval = tool.closest_triangle( root, far, point.array(), tri );
715  CHECK( MB_SUCCESS == rval );
716  CHECK( fabs( point[0] - 2.75 ) < 1e-6 );
717  CHECK( fabs( point[1] - 2.75 ) < 1e-6 );
718  CHECK( fabs( point[2] - 3.00 ) < 1e-6 );
719  // check that point is contained within triangle
720  check_point_in_triangles( tool.moab(), &tri, 1, point );
721  }
722 }
723 
724 void test_sphere_intersect_triangles()
725 {
726  Core moab;
727  AdaptiveKDTree tool( &moab );
728  Range box_tris;
729 
730  std::ostringstream options;
731  options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
732  FileOptions opts( options.str().c_str() );
733  EntityHandle root;
734  ErrorCode rval;
735  std::vector< EntityHandle > triangles;
736 
737  moab.delete_mesh();
738  box_tris.clear();
739  build_triangle_box_small( &moab, box_tris );
740  rval = tool.build_tree( box_tris, &root, &opts );
741  CHECK( MB_SUCCESS == rval );
742  test_valid_tree( &tool, root, opts, box_tris );
743 
744  // test closest to each corner of the box.
745  for( unsigned i = 0; i < 8; ++i )
746  {
747  double x = i & 1 ? 1 : -1;
748  double y = i & 2 ? 1 : -1;
749  double z = i & 4 ? 1 : -1;
750  double center[] = { x, y, z };
751  triangles.clear();
752  rval = tool.sphere_intersect_triangles( root, center, 0.26, triangles );
753  CHECK( MB_SUCCESS == rval );
754  // expect 3 to 6 triangles, depending on the corner
755  CHECK( triangles.size() >= 3 );
756  CHECK( triangles.size() <= 6 );
757  // check point is at the same vertex for each triangle
758  check_common_vertex( tool.moab(), &triangles[0], triangles.size(), CartVect( center ) );
759  }
760 
761  // now do it all again with a lot more triangles
762 
763  tool.reset_tree();
764  moab.delete_mesh();
765  box_tris.clear();
766  build_triangle_box_large( &moab, box_tris );
767 
768  rval = tool.build_tree( box_tris, &root, &opts );
769  CHECK( MB_SUCCESS == rval );
770  test_valid_tree( &tool, root, opts, box_tris );
771 
772  // test closest to each corner of the box.
773  for( unsigned i = 0; i < 8; ++i )
774  {
775  int x = i & 1 ? 1 : -1;
776  int y = i & 2 ? 1 : -1;
777  int z = i & 4 ? 1 : -1;
778  double center[] = { 3.0 * x, 3.0 * y, 3.0 * z };
779  triangles.clear();
780  rval = tool.sphere_intersect_triangles( root, center, 0.26, triangles );
781  CHECK( MB_SUCCESS == rval );
782  // expect 3 to 6 triangles, depending on the corner
783  CHECK( triangles.size() >= 3 );
784  CHECK( triangles.size() <= 6 );
785  // check point is at the same vertex for each triangle
786  check_common_vertex( tool.moab(), &triangles[0], triangles.size(), CartVect( center ) );
787  }
788 }
789 
790 void test_ray_intersect_triangles()
791 {
792  Core moab;
793  AdaptiveKDTree tool( &moab );
794  Range box_tris;
795 
796  std::ostringstream options;
797  options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
798  FileOptions opts( options.str().c_str() );
799  EntityHandle root;
800  ErrorCode rval;
801  std::vector< EntityHandle > tris;
802  std::vector< double > dists;
803 
804  moab.delete_mesh();
805  box_tris.clear();
806  build_triangle_box_small( &moab, box_tris );
807  rval = tool.build_tree( box_tris, &root, &opts );
808  CHECK( MB_SUCCESS == rval );
809  test_valid_tree( &tool, root, opts, box_tris );
810 
811  // test ray through box parallel to X axis
812  CartVect dir( 1, 0, 0 );
813  CartVect pt( -2, 0.75, 0.75 );
814  tris.clear();
815  dists.clear();
816  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
817  CHECK( MB_SUCCESS == rval );
818  CHECK( tris.size() == 3 );
819  CHECK( dists.size() == tris.size() );
820  for( unsigned i = 0; i < dists.size(); ++i )
821  {
822  CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 );
823  CartVect tript = pt + dists[i] * dir;
824  check_point_in_triangles( &moab, &tris[i], 1, tript );
825  }
826 
827  // test ray through box parallel to X axis, closest tri only
828  tris.clear();
829  dists.clear();
830  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 1 );
831  CHECK( MB_SUCCESS == rval );
832  CHECK( tris.size() == 1 );
833  CHECK( dists.size() == tris.size() );
834  CHECK( fabs( dists[0] - 1 ) < 1e-6 );
835  check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir );
836 
837  // test ray ending within box, parallel to X axis
838  tris.clear();
839  dists.clear();
840  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 0, 2.0 );
841  CHECK( MB_SUCCESS == rval );
842  CHECK( tris.size() == 1 );
843  CHECK( dists.size() == tris.size() );
844  CHECK( fabs( dists[0] - 1 ) < 1e-6 );
845  check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir );
846 
847  // test ray starting within box parallel to X axis
848  tris.clear();
849  dists.clear();
850  pt = CartVect( 0, .75, .75 );
851  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
852  CHECK( MB_SUCCESS == rval );
853  CHECK( tris.size() == 2 );
854  CHECK( dists.size() == tris.size() );
855  for( unsigned i = 0; i < dists.size(); ++i )
856  {
857  CHECK( fabs( dists[i] - 1 ) < 1e-6 );
858  CartVect tript = pt + dists[i] * dir;
859  check_point_in_triangles( &moab, &tris[i], 1, tript );
860  }
861 
862  // test skew ray through box
863  dir = CartVect( 0.5 * sqrt( 2.0 ), 0.5 * sqrt( 2.0 ), 0.0 );
864  pt = CartVect( 0, -1.5, 0 );
865  tris.clear();
866  dists.clear();
867  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
868  CHECK( MB_SUCCESS == rval );
869  CHECK( tris.size() == 2 );
870  CHECK( dists.size() == tris.size() );
871  if( dists[0] < dists[1] )
872  {
873  CHECK( fabs( dists[0] - 0.5 * sqrt( 2.0 ) ) < 1e-6 );
874  CHECK( fabs( dists[1] - sqrt( 2.0 ) ) < 1e-6 );
875  }
876  check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir );
877  check_point_in_triangles( &moab, &tris[1], 1, pt + dists[1] * dir );
878 
879  // test ray through box parallel to -Y axis
880  dir = CartVect( 0, -1, 0 );
881  pt = CartVect( 0, 2, 0 );
882  tris.clear();
883  dists.clear();
884  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
885  CHECK( MB_SUCCESS == rval );
886  CHECK( tris.size() == 4 );
887  CHECK( dists.size() == tris.size() );
888  for( unsigned i = 0; i < dists.size(); ++i )
889  {
890  CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 );
891  CartVect tript = pt + dists[i] * dir;
892  check_point_in_triangles( &moab, &tris[i], 1, tript );
893  }
894 
895  // test ray through box parallel to Z axis
896  dir = CartVect( 0, 0, 1 );
897  pt = CartVect( 0, 0, -2 );
898  tris.clear();
899  dists.clear();
900  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
901  CHECK( MB_SUCCESS == rval );
902  CHECK( tris.size() == 4 );
903  CHECK( dists.size() == tris.size() );
904  for( unsigned i = 0; i < dists.size(); ++i )
905  {
906  CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 );
907  CartVect tript = pt + dists[i] * dir;
908  check_point_in_triangles( &moab, &tris[i], 1, tript );
909  }
910 
911  // test ray through box parallel to Z axis, limit 2 first 2 intersections
912  tris.clear();
913  dists.clear();
914  rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 2 );
915  CHECK( MB_SUCCESS == rval );
916  CHECK( tris.size() == 2 );
917  CHECK( dists.size() == tris.size() );
918  for( unsigned i = 0; i < dists.size(); ++i )
919  {
920  CHECK( fabs( dists[i] - 1 ) < 1e-6 );
921  CartVect tript = pt + dists[i] * dir;
922  check_point_in_triangles( &moab, &tris[i], 1, tript );
923  }
924 }
925 
926 void test_leaf_volume()
927 {
928  Core moab;
929  AdaptiveKDTree tool( &moab );
930  EntityHandle root;
931 
932  create_simple_2d_tree( tool, root );
933  AdaptiveKDTreeIter iter;
934  ErrorCode rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
935 
936  CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 ); // 1
937  CHECK_ERR( iter.step() );
938  CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 ); // 2
939  CHECK_ERR( iter.step() );
940  CHECK_REAL_EQUAL( 40.0, iter.volume(), 1e-12 ); // 3
941  CHECK_ERR( iter.step() );
942  CHECK_REAL_EQUAL( 8.0, iter.volume(), 1e-12 ); // 4
943  CHECK_ERR( iter.step() );
944  CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 ); // 5
945  CHECK_ERR( iter.step() );
946  CHECK_REAL_EQUAL( 24.0, iter.volume(), 1e-12 ); // 6
947  CHECK_ERR( iter.step() );
948  CHECK_REAL_EQUAL( 30.0, iter.volume(), 1e-12 ); // 7
949  CHECK_ERR( iter.step() );
950  CHECK_REAL_EQUAL( 10.0, iter.volume(), 1e-12 ); // 8
951 }
952 
953 void test_leaf_sibling()
954 {
955  ErrorCode rval;
956  Core moab;
957  AdaptiveKDTree tool( &moab );
958  EntityHandle root;
959 
960  create_simple_2d_tree( tool, root );
961  AdaptiveKDTreeIter iter1, iter2;
962  rval = tool.get_tree_iterator( root, iter1 );CHECK_ERR( rval );
963  rval = tool.get_tree_iterator( root, iter2 );CHECK_ERR( rval );
964 
965  // iter1 == 1, iter2 == 1
966  CHECK( !iter1.is_sibling( iter1 ) );
967  CHECK( !iter1.is_sibling( iter1.handle() ) );
968  CHECK( !iter1.is_sibling( iter2 ) );
969  CHECK( iter1.sibling_is_forward() );
970 
971  // iter1 == 1, iter2 == 2
972  rval = iter2.step();CHECK_ERR( rval );
973  CHECK( iter1.is_sibling( iter2 ) );
974  CHECK( iter1.is_sibling( iter2.handle() ) );
975  CHECK( iter2.is_sibling( iter1 ) );
976  CHECK( iter2.is_sibling( iter1.handle() ) );
977  CHECK( !iter2.sibling_is_forward() );
978 
979  // iter1 == 1, iter2 == 3
980  rval = iter2.step();CHECK_ERR( rval );
981  CHECK( !iter1.is_sibling( iter2 ) );
982  CHECK( !iter1.is_sibling( iter2.handle() ) );
983  CHECK( !iter2.is_sibling( iter1 ) );
984  CHECK( !iter2.is_sibling( iter1.handle() ) );
985  CHECK( iter2.sibling_is_forward() );
986 
987  // iter1 == 2, iter2 == 3
988  rval = iter1.step();CHECK_ERR( rval );
989  CHECK( !iter1.is_sibling( iter2 ) );
990  CHECK( !iter1.is_sibling( iter2.handle() ) );
991  CHECK( !iter2.is_sibling( iter1 ) );
992  CHECK( !iter2.is_sibling( iter1.handle() ) );
993  CHECK( !iter1.sibling_is_forward() );
994 
995  // iter1 == 4, iter2 == 3
996  rval = iter1.step();CHECK_ERR( rval );
997  CHECK_EQUAL( iter1.handle(), iter2.handle() );
998  rval = iter1.step();CHECK_ERR( rval );
999  CHECK( iter1.is_sibling( iter2 ) );
1000  CHECK( iter1.is_sibling( iter2.handle() ) );
1001  CHECK( iter2.is_sibling( iter1 ) );
1002  CHECK( iter2.is_sibling( iter1.handle() ) );
1003  CHECK( !iter1.sibling_is_forward() );
1004 }
1005 
1006 void test_leaf_intersects_plane()
1007 {
1008  ErrorCode rval;
1009  Core moab;
1010  AdaptiveKDTree tool( &moab );
1011 
1012  EntityHandle root;
1013  const double min[3] = { -5, -4, -1 };
1014  const double max[3] = { 1, 2, 3 };
1015  rval = tool.create_root( min, max, root );CHECK_ERR( rval );
1016 
1017  AdaptiveKDTreeIter iter;
1018  rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
1019 
1020  AdaptiveKDTree::Plane x1 = { min[0] - 1, AdaptiveKDTree::X };
1021  CHECK( !iter.intersects( x1 ) );
1022  AdaptiveKDTree::Plane x2 = { min[0], AdaptiveKDTree::X };
1023  CHECK( iter.intersects( x2 ) );
1024  AdaptiveKDTree::Plane x3 = { min[0] + 1, AdaptiveKDTree::X };
1025  CHECK( iter.intersects( x3 ) );
1026  AdaptiveKDTree::Plane x4 = { max[0] - 1, AdaptiveKDTree::X };
1027  CHECK( iter.intersects( x4 ) );
1028  AdaptiveKDTree::Plane x5 = { max[0], AdaptiveKDTree::X };
1029  CHECK( iter.intersects( x5 ) );
1030  AdaptiveKDTree::Plane x6 = { max[0] + 1, AdaptiveKDTree::X };
1031  CHECK( !iter.intersects( x6 ) );
1032 
1033  AdaptiveKDTree::Plane y1 = { min[1] - 1, AdaptiveKDTree::Y };
1034  CHECK( !iter.intersects( y1 ) );
1035  AdaptiveKDTree::Plane y2 = { min[1], AdaptiveKDTree::Y };
1036  CHECK( iter.intersects( y2 ) );
1037  AdaptiveKDTree::Plane y3 = { min[1] + 1, AdaptiveKDTree::Y };
1038  CHECK( iter.intersects( y3 ) );
1039  AdaptiveKDTree::Plane y4 = { max[1] - 1, AdaptiveKDTree::Y };
1040  CHECK( iter.intersects( y4 ) );
1041  AdaptiveKDTree::Plane y5 = { max[1], AdaptiveKDTree::Y };
1042  CHECK( iter.intersects( y5 ) );
1043  AdaptiveKDTree::Plane y6 = { max[1] + 1, AdaptiveKDTree::Y };
1044  CHECK( !iter.intersects( y6 ) );
1045 
1046  AdaptiveKDTree::Plane z1 = { min[2] - 1, AdaptiveKDTree::Z };
1047  CHECK( !iter.intersects( z1 ) );
1048  AdaptiveKDTree::Plane z2 = { min[2], AdaptiveKDTree::Z };
1049  CHECK( iter.intersects( z2 ) );
1050  AdaptiveKDTree::Plane z3 = { min[2] + 1, AdaptiveKDTree::Z };
1051  CHECK( iter.intersects( z3 ) );
1052  AdaptiveKDTree::Plane z4 = { max[2] - 1, AdaptiveKDTree::Z };
1053  CHECK( iter.intersects( z4 ) );
1054  AdaptiveKDTree::Plane z5 = { max[2], AdaptiveKDTree::Z };
1055  CHECK( iter.intersects( z5 ) );
1056  AdaptiveKDTree::Plane z6 = { max[2] + 1, AdaptiveKDTree::Z };
1057  CHECK( !iter.intersects( z6 ) );
1058 }
1059 
1060 #define CHECK_RAY_XSECTS( PT, DIR, T_IN, T_OUT ) \
1061  do \
1062  { \
1063  CHECK( iter.intersect_ray( ( PT ), ( DIR ), t_in, t_out ) ); \
1064  CHECK_REAL_EQUAL( ( T_IN ), t_in, 1e-6 ); \
1065  CHECK_REAL_EQUAL( ( T_OUT ), t_out, 1e-6 ); \
1066  } while( false )
1067 
1068 void test_leaf_intersects_ray()
1069 {
1070  ErrorCode rval;
1071  Core moab;
1072  AdaptiveKDTree tool( &moab );
1073  double t_in, t_out;
1074 
1075  EntityHandle root;
1076  const double min[3] = { -5, -4, -1 };
1077  const double max[3] = { 1, 2, 3 };
1078  rval = tool.create_root( min, max, root );CHECK_ERR( rval );
1079 
1080  AdaptiveKDTreeIter iter;
1081  rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
1082 
1083  // start with point inside box
1084  const double pt1[] = { 0, 0, 0 };
1085  const double dir1[] = { 0.1, 0.1, 0.1 };
1086  CHECK_RAY_XSECTS( pt1, dir1, 0, 10 );
1087  const double dir2[] = { 5, 5, 5 };
1088  CHECK_RAY_XSECTS( pt1, dir2, 0, 0.2 );
1089  const double pxdir[] = { 1, 0, 0 };
1090  CHECK_RAY_XSECTS( pt1, pxdir, 0, 1 );
1091  const double nxdir[] = { -1, 0, 0 };
1092  CHECK_RAY_XSECTS( pt1, nxdir, 0, 5 );
1093  const double pydir[] = { 0, 1, 0 };
1094  CHECK_RAY_XSECTS( pt1, pydir, 0, 2 );
1095  const double nydir[] = { 0, -1, 0 };
1096  CHECK_RAY_XSECTS( pt1, nydir, 0, 4 );
1097  const double pzdir[] = { 0, 0, 1 };
1098  CHECK_RAY_XSECTS( pt1, pzdir, 0, 3 );
1099  const double nzdir[] = { 0, 0, -1 };
1100  CHECK_RAY_XSECTS( pt1, nzdir, 0, 1 );
1101 
1102  // point below box
1103  const double pt2[] = { 0, 0, -2 };
1104  CHECK_RAY_XSECTS( pt2, dir1, 10, 10 );
1105  CHECK_RAY_XSECTS( pt2, dir2, 0.2, 0.2 );
1106  CHECK( !iter.intersect_ray( pt2, pxdir, t_in, t_out ) );
1107  CHECK( !iter.intersect_ray( pt2, nxdir, t_in, t_out ) );
1108  CHECK( !iter.intersect_ray( pt2, pydir, t_in, t_out ) );
1109  CHECK( !iter.intersect_ray( pt2, nydir, t_in, t_out ) );
1110  CHECK_RAY_XSECTS( pt2, pzdir, 1, 5 );
1111  CHECK( !iter.intersect_ray( pt2, nzdir, t_in, t_out ) );
1112 
1113  // point right of box
1114  const double pt3[] = { 3, 0, 0 };
1115  CHECK( !iter.intersect_ray( pt3, dir1, t_in, t_out ) );
1116  CHECK( !iter.intersect_ray( pt3, dir2, t_in, t_out ) );
1117  CHECK( !iter.intersect_ray( pt3, pxdir, t_in, t_out ) );
1118  CHECK_RAY_XSECTS( pt3, nxdir, 2, 8 );
1119  CHECK( !iter.intersect_ray( pt3, pydir, t_in, t_out ) );
1120  CHECK( !iter.intersect_ray( pt3, nydir, t_in, t_out ) );
1121  CHECK( !iter.intersect_ray( pt3, pzdir, t_in, t_out ) );
1122  CHECK( !iter.intersect_ray( pt3, nzdir, t_in, t_out ) );
1123 
1124  // a few more complex test cases
1125  const double dira[] = { -3, 0, 3 };
1126  CHECK_RAY_XSECTS( pt3, dira, 2. / 3., 1.0 );
1127  const double dirb[] = { -2, 0, 3.1 };
1128  CHECK( !iter.intersect_ray( pt3, dirb, t_in, t_out ) );
1129 }
1130 
1131 int main()
1132 {
1133 #ifdef MOAB_HAVE_MPI
1134  int fail = MPI_Init( 0, 0 );
1135  if( fail ) return fail;
1136 #endif
1137 
1138  int error_count = 0;
1139 
1140  error_count += RUN_TEST( leaf_iterator_test );
1141  error_count += RUN_TEST( test_build_tree_bisect_triangles );
1142  error_count += RUN_TEST( test_closest_triangle );
1143  error_count += RUN_TEST( test_sphere_intersect_triangles );
1144  error_count += RUN_TEST( test_ray_intersect_triangles );
1145  error_count += RUN_TEST( test_tree_merge_nodes );
1146  error_count += RUN_TEST( test_leaf_volume );
1147  error_count += RUN_TEST( test_leaf_sibling );
1148  error_count += RUN_TEST( test_leaf_intersects_plane );
1149  error_count += RUN_TEST( test_leaf_intersects_ray );
1150 
1151 #ifdef MOAB_HAVE_MPI
1152  fail = MPI_Finalize();
1153  if( fail ) return fail;
1154 #endif
1155 
1156  return error_count;
1157 }