MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/Core.hpp" 00002 #include "moab/AdaptiveKDTree.hpp" 00003 #include "moab/CartVect.hpp" 00004 #include "moab/GeomUtil.hpp" 00005 #include "moab/Range.hpp" 00006 #include "TestUtil.hpp" 00007 00008 using namespace moab; 00009 00010 #include <iostream> 00011 #include <algorithm> 00012 #include <sstream> 00013 00014 #ifdef MOAB_HAVE_MPI 00015 #include "moab_mpi.h" 00016 #endif 00017 00018 /* Utility method - compare two range boxes */ 00019 bool box_equal( const AdaptiveKDTreeIter& iter, double x_min, double y_min, double z_min, double x_max, double y_max, 00020 double z_max ) 00021 { 00022 return iter.box_min()[0] == x_min && iter.box_min()[1] == y_min && iter.box_min()[2] == z_min && 00023 iter.box_max()[0] == x_max && iter.box_max()[1] == y_max && iter.box_max()[2] == z_max; 00024 } 00025 00026 void build_triangles( Interface* moab, Range& tris, int num_vert, const double* coords, int num_tri, 00027 const unsigned* conn ) 00028 { 00029 std::vector< EntityHandle > verts( num_vert ); 00030 for( int i = 0; i < num_vert; ++i ) 00031 moab->create_vertex( coords + 3 * i, verts[i] ); 00032 00033 EntityHandle tri, tri_verts[3]; 00034 for( int i = 0; i < num_tri; ++i ) 00035 { 00036 tri_verts[0] = verts[conn[3 * i]]; 00037 tri_verts[1] = verts[conn[3 * i + 1]]; 00038 tri_verts[2] = verts[conn[3 * i + 2]]; 00039 moab->create_element( MBTRI, tri_verts, 3, tri ); 00040 tris.insert( tri ); 00041 } 00042 } 00043 00044 /* Utility method - build a 2x2x2 box composed of two triagles on a side */ 00045 void build_triangle_box_small( Interface* moab, Range& tris ) 00046 { 00047 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 }; 00048 00049 const unsigned conn[] = { 0, 1, 5, 0, 5, 4, 2, 6, 7, 2, 7, 3, 1, 2, 6, 1, 6, 5, 00050 0, 3, 4, 3, 7, 4, 0, 1, 3, 1, 2, 3, 4, 5, 6, 4, 6, 7 }; 00051 00052 build_triangles( moab, tris, 8, coords, 12, conn ); 00053 } 00054 00055 /* build 6x6x6 box centered at origin and composed of 18 triangles per side, 00056 Each side of the box is composed of a 3x3 grid of quads, where each 00057 quad is split diagonally to form two triangles. */ 00058 void build_triangle_box_large( Interface* moab, Range& tris ) 00059 { 00060 const double coords[] = { // corners 00061 -3, -3, -3, 3, -3, -3, 3, 3, -3, -3, 3, -3, -3, -3, 3, 3, -3, 3, 3, 3, 3, -3, 3, 3, 00062 00063 // edges 00064 -1, -3, -3, 1, -3, -3, 3, -1, -3, 3, 1, -3, 1, 3, -3, -1, 3, -3, -3, 1, -3, -3, -1, -3, 00065 00066 -1, -3, 3, 1, -3, 3, 3, -1, 3, 3, 1, 3, 1, 3, 3, -1, 3, 3, -3, 1, 3, -3, -1, 3, 00067 00068 -3, -3, -1, -3, -3, 1, 3, -3, -1, 3, -3, 1, 3, 3, -1, 3, 3, 1, -3, 3, -1, -3, 3, 1, 00069 00070 // faces 00071 -1, -3, -1, 1, -3, -1, 1, -3, 1, -1, -3, 1, 00072 00073 3, -1, -1, 3, 1, -1, 3, 1, 1, 3, -1, 1, 00074 00075 1, 3, -1, -1, 3, -1, -1, 3, 1, 1, 3, 1, 00076 00077 -3, 1, -1, -3, -1, -1, -3, -1, 1, -3, 1, 1, 00078 00079 -1, -1, -3, 1, -1, -3, 1, 1, -3, -1, 1, -3, 00080 00081 -1, -1, 3, 1, -1, 3, 1, 1, 3, -1, 1, 3 00082 }; 00083 const unsigned conn[] = { 00084 // face 0 00085 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, 00086 33, 27, 34, 33, 26, 27, 35, 4, 25, 35, 16, 4, 35, 17, 16, 35, 34, 17, 27, 17, 34, 27, 5, 17, 00087 // face 1 00088 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, 00089 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, 00090 // face 2 00091 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, 00092 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, 00093 // face 3 00094 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, 00095 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, 00096 // face 4 00097 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, 00098 11, 50, 49, 11, 49, 10, 51, 3, 14, 51, 13, 3, 51, 12, 13, 51, 50, 12, 11, 12, 50, 11, 2, 12, 00099 // face 5 00100 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, 00101 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 00102 }; 00103 00104 build_triangles( moab, tris, 56, coords, 108, conn ); 00105 } 00106 00107 /* Utility method - build 2x2x2 octahedron (3D diamond)*/ 00108 void build_triangle_octahedron( Interface* moab, Range& tris ) 00109 { 00110 const double coords[] = { 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1 }; 00111 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 }; 00112 00113 build_triangles( moab, tris, 6, coords, 8, conn ); 00114 } 00115 00116 void test_valid_tree( AdaptiveKDTree* tool, EntityHandle root, FileOptions& options, const Range& expected_tris ) 00117 { 00118 Range all_tris; 00119 ErrorCode rval; 00120 AdaptiveKDTreeIter iter; 00121 CHECK( MB_SUCCESS == tool->get_tree_iterator( root, iter ) ); 00122 do 00123 { 00124 double dep; 00125 rval = options.get_real_option( "MAX_DEPTH", dep ); 00126 CHECK( MB_ENTITY_NOT_FOUND == rval || iter.depth() <= tool->get_max_depth() ); 00127 00128 Range tris; 00129 CHECK( tool->moab()->get_entities_by_type( iter.handle(), MBTRI, tris ) == MB_SUCCESS ); 00130 // CHECK( !tris.empty() ); 00131 all_tris.merge( tris ); 00132 00133 const CartVect min( iter.box_min() ), max( iter.box_max() ); 00134 const CartVect cen( 0.5 * ( min + max ) ), hdim( 0.5 * ( max - min ) ); 00135 00136 for( Range::iterator i = tris.begin(); i != tris.end(); ++i ) 00137 { 00138 EntityHandle tri = *i; 00139 const EntityHandle* conn; 00140 int conn_len; 00141 CHECK( tool->moab()->get_connectivity( tri, conn, conn_len ) == MB_SUCCESS ); 00142 CHECK( conn_len == 3 ); 00143 CartVect coords[3]; 00144 CHECK( tool->moab()->get_coords( conn, 3, coords[0].array() ) == MB_SUCCESS ); 00145 CHECK( GeomUtil::box_tri_overlap( coords, cen, hdim ) ); 00146 } 00147 } while( MB_SUCCESS == ( rval = iter.step() ) ); 00148 00149 CHECK( MB_ENTITY_NOT_FOUND == rval ); 00150 CHECK( expected_tris == all_tris ); 00151 } 00152 00153 /* utility method - check that all tris share a vertex */ 00154 void check_common_vertex( Interface* moab, const EntityHandle* tris, unsigned num_tri, CartVect point ) 00155 { 00156 for( unsigned i = 0; i < num_tri; ++i ) 00157 CHECK( MBTRI == moab->type_from_handle( tris[i] ) ); 00158 00159 ErrorCode rval; 00160 CartVect tri_coords[3]; 00161 const EntityHandle* conn; 00162 int conn_len; 00163 rval = moab->get_connectivity( tris[0], conn, conn_len ); 00164 CHECK( MB_SUCCESS == rval ); 00165 rval = moab->get_coords( conn, 3, tri_coords[0].array() ); 00166 CHECK( MB_SUCCESS == rval ); 00167 tri_coords[0] -= point; 00168 tri_coords[1] -= point; 00169 tri_coords[2] -= point; 00170 EntityHandle vertex = 0; 00171 if( tri_coords[0] % tri_coords[0] < 1e-6 ) 00172 vertex = conn[0]; 00173 else if( tri_coords[1] % tri_coords[1] < 1e-6 ) 00174 vertex = conn[1]; 00175 else if( tri_coords[2] % tri_coords[2] < 1e-6 ) 00176 vertex = conn[2]; 00177 else 00178 CHECK( false ); 00179 for( unsigned j = 1; j < num_tri; ++j ) 00180 { 00181 rval = moab->get_connectivity( tris[j], conn, conn_len ); 00182 CHECK( conn[0] == vertex || conn[1] == vertex || conn[2] == vertex ); 00183 } 00184 } 00185 00186 /* utility method - check that all tris share a vertex */ 00187 void check_point_in_triangles( Interface* moab, const EntityHandle* tris, unsigned num_tris, CartVect point ) 00188 { 00189 ErrorCode rval; 00190 CartVect tri_coords[3], tript; 00191 const EntityHandle* conn; 00192 int conn_len; 00193 00194 for( unsigned i = 0; i < num_tris; ++i ) 00195 { 00196 CHECK( MBTRI == moab->type_from_handle( tris[i] ) ); 00197 00198 rval = moab->get_connectivity( tris[i], conn, conn_len ); 00199 CHECK( MB_SUCCESS == rval ); 00200 rval = moab->get_coords( conn, 3, tri_coords[0].array() ); 00201 CHECK( MB_SUCCESS == rval ); 00202 00203 GeomUtil::closest_location_on_tri( point, tri_coords, tript ); 00204 tript -= point; 00205 CHECK( fabs( tript[0] ) < 1e-6 ); 00206 CHECK( fabs( tript[1] ) < 1e-6 ); 00207 CHECK( fabs( tript[2] ) < 1e-6 ); 00208 } 00209 } 00210 00211 /* Create the following 2D tree (no Z splits) 00212 00213 (6) X = -3 (5) X = 0 00214 / | 00215 / | [8] 00216 / |---------------- (7) Y = 3 00217 [5] / [6] | 00218 / | [7] 00219 / | 00220 / | 00221 -------------------------------------- (1) Y = 0 00222 | | 00223 [2] | | 00224 (3)----------------| | 00225 Y = -2 | [3] | [4] 00226 | | 00227 [1] | | 00228 | | 00229 (2) X = -1 (4) X = 4 00230 */ 00231 void create_simple_2d_tree( AdaptiveKDTree& tool, EntityHandle& root, EntityHandle leaves[9] = 0 ) 00232 { 00233 ErrorCode rval; 00234 AdaptiveKDTree::Plane plane; 00235 00236 // create a single-node tree 00237 double min[3] = { -5, -4, -1 }; 00238 double max[3] = { 5, 4, 1 }; 00239 rval = tool.create_root( min, max, root ); 00240 CHECK( MB_SUCCESS == rval ); 00241 00242 // get iterator for tree 00243 AdaptiveKDTreeIter iter; 00244 rval = tool.get_tree_iterator( root, iter ); 00245 CHECK( MB_SUCCESS == rval ); 00246 CHECK( box_equal( iter, -5, -4, -1, 5, 4, 1 ) ); 00247 00248 // split plane (1) 00249 plane.norm = AdaptiveKDTree::Y; 00250 plane.coord = 0.0; 00251 rval = tool.split_leaf( iter, plane ); 00252 CHECK( MB_SUCCESS == rval ); 00253 CHECK( box_equal( iter, -5, -4, -1, 5, 0, 1 ) ); 00254 00255 // split plane (2) 00256 plane.norm = AdaptiveKDTree::X; 00257 plane.coord = -1.0; 00258 rval = tool.split_leaf( iter, plane ); 00259 CHECK( MB_SUCCESS == rval ); 00260 CHECK( box_equal( iter, -5, -4, -1, -1, 0, 1 ) ); 00261 00262 // split plane (3), leaf [1] 00263 plane.norm = AdaptiveKDTree::Y; 00264 plane.coord = -2.0; 00265 rval = tool.split_leaf( iter, plane ); 00266 CHECK( MB_SUCCESS == rval ); 00267 CHECK( box_equal( iter, -5, -4, -1, -1, -2, 1 ) ); 00268 if( leaves ) leaves[1] = iter.handle(); 00269 00270 // leaf [2] 00271 rval = iter.step(); 00272 CHECK( MB_SUCCESS == rval ); 00273 CHECK( box_equal( iter, -5, -2, -1, -1, 0, 1 ) ); 00274 if( leaves ) leaves[2] = iter.handle(); 00275 00276 // non-leaf right of split plane (2) 00277 rval = iter.step(); 00278 CHECK( MB_SUCCESS == rval ); 00279 CHECK( box_equal( iter, -1, -4, -1, 5, 0, 1 ) ); 00280 00281 // split plane (4) and leaf [3] 00282 plane.norm = AdaptiveKDTree::X; 00283 plane.coord = 4.0; 00284 rval = tool.split_leaf( iter, plane ); 00285 CHECK( MB_SUCCESS == rval ); 00286 CHECK( box_equal( iter, -1, -4, -1, 4, 0, 1 ) ); 00287 if( leaves ) leaves[3] = iter.handle(); 00288 00289 // leaf [4] 00290 rval = iter.step(); 00291 CHECK( MB_SUCCESS == rval ); 00292 CHECK( box_equal( iter, 4, -4, -1, 5, 0, 1 ) ); 00293 if( leaves ) leaves[4] = iter.handle(); 00294 00295 // non-leaf above split plane (1) 00296 rval = iter.step(); 00297 CHECK( MB_SUCCESS == rval ); 00298 CHECK( box_equal( iter, -5, 0, -1, 5, 4, 1 ) ); 00299 00300 // split plane (5) 00301 plane.norm = AdaptiveKDTree::X; 00302 plane.coord = 0.0; 00303 rval = tool.split_leaf( iter, plane ); 00304 CHECK( MB_SUCCESS == rval ); 00305 CHECK( box_equal( iter, -5, 0, -1, 0, 4, 1 ) ); 00306 00307 // split plane (6) and leaf [5] 00308 plane.norm = AdaptiveKDTree::X; 00309 plane.coord = -3.0; 00310 rval = tool.split_leaf( iter, plane ); 00311 CHECK( MB_SUCCESS == rval ); 00312 CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) ); 00313 if( leaves ) leaves[5] = iter.handle(); 00314 00315 // leaf [6]; 00316 rval = iter.step(); 00317 CHECK( MB_SUCCESS == rval ); 00318 CHECK( box_equal( iter, -3, 0, -1, 0, 4, 1 ) ); 00319 if( leaves ) leaves[6] = iter.handle(); 00320 00321 // non-leaf right of split plane (5) 00322 rval = iter.step(); 00323 CHECK( MB_SUCCESS == rval ); 00324 CHECK( box_equal( iter, 0, 0, -1, 5, 4, 1 ) ); 00325 00326 // split plane (7) and leaf [7] 00327 plane.norm = AdaptiveKDTree::Y; 00328 plane.coord = 3.0; 00329 rval = tool.split_leaf( iter, plane ); 00330 CHECK( MB_SUCCESS == rval ); 00331 CHECK( box_equal( iter, 0, 0, -1, 5, 3, 1 ) ); 00332 if( leaves ) leaves[7] = iter.handle(); 00333 00334 // leaf [8]; 00335 rval = iter.step(); 00336 CHECK( MB_SUCCESS == rval ); 00337 CHECK( box_equal( iter, 0, 3, -1, 5, 4, 1 ) ); 00338 if( leaves ) leaves[8] = iter.handle(); 00339 00340 // the end 00341 rval = iter.step(); 00342 CHECK( MB_ENTITY_NOT_FOUND == rval ); 00343 } 00344 00345 void leaf_iterator_test() 00346 { 00347 ErrorCode rval; 00348 Core moab; 00349 Interface* mb = &moab; 00350 AdaptiveKDTree tool( mb ); 00351 00352 // create a single-node tree 00353 EntityHandle root; 00354 double min[3] = { -3, -2, -1 }; 00355 double max[3] = { 1, 2, 3 }; 00356 rval = tool.create_root( min, max, root ); 00357 CHECK( MB_SUCCESS == rval ); 00358 00359 // get iterator for tree 00360 AdaptiveKDTreeIter iter; 00361 rval = tool.get_tree_iterator( root, iter ); 00362 CHECK( MB_SUCCESS == rval ); 00363 CHECK( box_equal( iter, -3, -2, -1, 1, 2, 3 ) ); 00364 00365 // check that steping the iterator fails correctly. 00366 rval = iter.step(); 00367 CHECK( MB_ENTITY_NOT_FOUND == rval ); 00368 rval = iter.step(); 00369 CHECK( MB_SUCCESS != rval ); 00370 00371 rval = tool.reset_tree(); 00372 CHECK( MB_SUCCESS == rval ); 00373 root = 0; 00374 00375 // construct a simple 2D tree for remaining tests 00376 EntityHandle leaves[9]; 00377 create_simple_2d_tree( tool, root, leaves ); 00378 00379 /**** Now traverse tree again, and check neighbors of each leaf *****/ 00380 std::vector< AdaptiveKDTreeIter > list; 00381 00382 // leaf [1] 00383 00384 rval = tool.get_tree_iterator( root, iter ); 00385 CHECK( MB_SUCCESS == rval ); 00386 CHECK( box_equal( iter, -5, -4, -1, -1, -2, 1 ) ); 00387 CHECK( iter.handle() == leaves[1] ); 00388 00389 list.clear(); 00390 rval = iter.get_neighbors( AdaptiveKDTree::X, true, list ); 00391 CHECK( MB_SUCCESS == rval ); 00392 CHECK( list.empty() ); 00393 00394 list.clear(); 00395 rval = iter.get_neighbors( AdaptiveKDTree::X, false, list ); 00396 CHECK( MB_SUCCESS == rval ); 00397 CHECK( list.size() == 1 ); 00398 // should be leaf [3] 00399 CHECK( box_equal( list.front(), -1, -4, -1, 4, 0, 1 ) ); 00400 CHECK( list.front().handle() == leaves[3] ); 00401 00402 list.clear(); 00403 rval = iter.get_neighbors( AdaptiveKDTree::Y, true, list ); 00404 CHECK( MB_SUCCESS == rval ); 00405 CHECK( list.empty() ); 00406 00407 list.clear(); 00408 rval = iter.get_neighbors( AdaptiveKDTree::Y, false, list ); 00409 CHECK( MB_SUCCESS == rval ); 00410 CHECK( list.size() == 1 ); 00411 // should be leaf [2] 00412 CHECK( box_equal( list.front(), -5, -2, -1, -1, 0, 1 ) ); 00413 CHECK( list.front().handle() == leaves[2] ); 00414 00415 list.clear(); 00416 rval = iter.get_neighbors( AdaptiveKDTree::Z, true, list ); 00417 CHECK( MB_SUCCESS == rval ); 00418 CHECK( list.empty() ); 00419 00420 list.clear(); 00421 rval = iter.get_neighbors( AdaptiveKDTree::Z, false, list ); 00422 CHECK( MB_SUCCESS == rval ); 00423 CHECK( list.empty() ); 00424 00425 // leaf [2] 00426 00427 rval = iter.step(); 00428 CHECK( MB_SUCCESS == rval ); 00429 CHECK( box_equal( iter, -5, -2, -1, -1, 0, 1 ) ); 00430 CHECK( iter.handle() == leaves[2] ); 00431 00432 list.clear(); 00433 rval = iter.get_neighbors( AdaptiveKDTree::X, true, list ); 00434 CHECK( MB_SUCCESS == rval ); 00435 CHECK( list.empty() ); 00436 00437 list.clear(); 00438 rval = iter.get_neighbors( AdaptiveKDTree::X, false, list ); 00439 CHECK( MB_SUCCESS == rval ); 00440 CHECK( list.size() == 1 ); 00441 // should be leaf [3] 00442 CHECK( box_equal( list.front(), -1, -4, -1, 4, 0, 1 ) ); 00443 CHECK( list.front().handle() == leaves[3] ); 00444 00445 list.clear(); 00446 rval = iter.get_neighbors( AdaptiveKDTree::Y, true, list ); 00447 CHECK( MB_SUCCESS == rval ); 00448 // should be leaf [1] 00449 CHECK( list.size() == 1 ); 00450 CHECK( box_equal( list.front(), -5, -4, -1, -1, -2, 1 ) ); 00451 CHECK( list.front().handle() == leaves[1] ); 00452 00453 list.clear(); 00454 rval = iter.get_neighbors( AdaptiveKDTree::Y, false, list ); 00455 CHECK( MB_SUCCESS == rval ); 00456 CHECK( list.size() == 2 ); 00457 // should be leaf [5] and leaf[6] 00458 if( list[0].handle() == leaves[6] ) std::swap( list[0], list[1] ); 00459 CHECK( list[0].handle() == leaves[5] ); 00460 CHECK( list[1].handle() == leaves[6] ); 00461 CHECK( box_equal( list[0], -5, 0, -1, -3, 4, 1 ) ); 00462 CHECK( box_equal( list[1], -3, 0, -1, 0, 4, 1 ) ); 00463 00464 list.clear(); 00465 rval = iter.get_neighbors( AdaptiveKDTree::Z, true, list ); 00466 CHECK( MB_SUCCESS == rval ); 00467 CHECK( list.empty() ); 00468 00469 list.clear(); 00470 rval = iter.get_neighbors( AdaptiveKDTree::Z, false, list ); 00471 CHECK( MB_SUCCESS == rval ); 00472 CHECK( list.empty() ); 00473 00474 // leaf [3] 00475 rval = iter.step(); 00476 CHECK( MB_SUCCESS == rval ); 00477 CHECK( box_equal( iter, -1, -4, -1, 4, 0, 1 ) ); 00478 CHECK( iter.handle() == leaves[3] ); 00479 00480 // leaf [4] 00481 rval = iter.step(); 00482 CHECK( MB_SUCCESS == rval ); 00483 CHECK( box_equal( iter, 4, -4, -1, 5, 0, 1 ) ); 00484 CHECK( iter.handle() == leaves[4] ); 00485 00486 // leaf [5] 00487 rval = iter.step(); 00488 CHECK( MB_SUCCESS == rval ); 00489 CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) ); 00490 CHECK( iter.handle() == leaves[5] ); 00491 00492 // leaf [6] 00493 rval = iter.step(); 00494 CHECK( MB_SUCCESS == rval ); 00495 CHECK( box_equal( iter, -3, 0, -1, 0, 4, 1 ) ); 00496 CHECK( iter.handle() == leaves[6] ); 00497 00498 // leaf [7] 00499 rval = iter.step(); 00500 CHECK( MB_SUCCESS == rval ); 00501 CHECK( box_equal( iter, 0, 0, -1, 5, 3, 1 ) ); 00502 CHECK( iter.handle() == leaves[7] ); 00503 00504 // leaf [8] 00505 rval = iter.step(); 00506 CHECK( MB_SUCCESS == rval ); 00507 CHECK( box_equal( iter, 0, 3, -1, 5, 4, 1 ) ); 00508 CHECK( iter.handle() == leaves[8] ); 00509 } 00510 00511 void test_tree_merge_nodes() 00512 { 00513 // build simple tree for tests 00514 ErrorCode rval; 00515 Core moab; 00516 Interface* mb = &moab; 00517 AdaptiveKDTree tool( mb ); 00518 EntityHandle root; 00519 EntityHandle leaves[9]; 00520 create_simple_2d_tree( tool, root, leaves ); 00521 00522 // get iterator for tree 00523 AdaptiveKDTreeIter iter; 00524 rval = tool.get_tree_iterator( root, iter ); 00525 CHECK( MB_SUCCESS == rval ); 00526 00527 // merge leaves 1 and 2 00528 rval = tool.merge_leaf( iter ); 00529 CHECK( MB_SUCCESS == rval ); 00530 CHECK( box_equal( iter, -5, -4, -1, -1, 0, 1 ) ); 00531 00532 // merge leaf 1,2 with 3,4 (implicity merges 3 and 4) 00533 rval = tool.merge_leaf( iter ); 00534 CHECK( MB_SUCCESS == rval ); 00535 CHECK( box_equal( iter, -5, -4, -1, 5, 0, 1 ) ); 00536 00537 // make sure iterator remains valid 00538 rval = iter.step(); 00539 CHECK( MB_SUCCESS == rval ); 00540 // leaf 5 00541 CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) ); 00542 CHECK( iter.handle() == leaves[5] ); 00543 } 00544 00545 void test_build_tree_bisect_triangles() 00546 { 00547 Core moab; 00548 AdaptiveKDTree tool( &moab ); 00549 Range box_tris; 00550 00551 std::ostringstream options; 00552 options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;"; 00553 FileOptions opts( options.str().c_str() ); 00554 EntityHandle root; 00555 ErrorCode rval; 00556 00557 moab.delete_mesh(); 00558 box_tris.clear(); 00559 build_triangle_box_small( &moab, box_tris ); 00560 rval = tool.build_tree( box_tris, &root, &opts ); 00561 CHECK( MB_SUCCESS == rval ); 00562 test_valid_tree( &tool, root, opts, box_tris ); 00563 tool.reset_tree(); 00564 00565 moab.delete_mesh(); 00566 box_tris.clear(); 00567 build_triangle_octahedron( &moab, box_tris ); 00568 rval = tool.build_tree( box_tris, &root, &opts ); 00569 CHECK( MB_SUCCESS == rval ); 00570 test_valid_tree( &tool, root, opts, box_tris ); 00571 tool.reset_tree(); 00572 00573 moab.delete_mesh(); 00574 box_tris.clear(); 00575 build_triangle_box_large( &moab, box_tris ); 00576 rval = tool.build_tree( box_tris, &root, &opts ); 00577 CHECK( MB_SUCCESS == rval ); 00578 test_valid_tree( &tool, root, opts, box_tris ); 00579 tool.reset_tree(); 00580 00581 options << "MAX_DEPTH=2;"; 00582 opts = FileOptions( options.str().c_str() ); 00583 build_triangle_box_large( &moab, box_tris ); 00584 rval = tool.build_tree( box_tris, &root, &opts ); 00585 CHECK( MB_SUCCESS == rval ); 00586 test_valid_tree( &tool, root, opts, box_tris ); 00587 tool.reset_tree(); 00588 } 00589 00590 void test_closest_triangle() 00591 { 00592 Core moab; 00593 AdaptiveKDTree tool( &moab ); 00594 Range box_tris; 00595 00596 std::ostringstream options; 00597 options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;"; 00598 FileOptions opts( options.str().c_str() ); 00599 EntityHandle root; 00600 ErrorCode rval; 00601 EntityHandle tri; 00602 00603 moab.delete_mesh(); 00604 box_tris.clear(); 00605 build_triangle_box_small( &moab, box_tris ); 00606 rval = tool.build_tree( box_tris, &root, &opts ); 00607 CHECK( MB_SUCCESS == rval ); 00608 test_valid_tree( &tool, root, opts, box_tris ); 00609 00610 // test closest to each corner of the box. 00611 for( unsigned i = 0; i < 8; ++i ) 00612 { 00613 int x = i & 1 ? 1 : -1; 00614 int y = i & 2 ? 1 : -1; 00615 int z = i & 4 ? 1 : -1; 00616 double coords[] = { 2.0 * x, 2.0 * y, 2.0 * z }; 00617 CartVect point; 00618 rval = tool.closest_triangle( root, coords, point.array(), tri ); 00619 CHECK( MB_SUCCESS == rval ); 00620 // check result position 00621 CHECK( fabs( x - point[0] ) < 1e-6 ); 00622 CHECK( fabs( y - point[1] ) < 1e-6 ); 00623 CHECK( fabs( z - point[2] ) < 1e-6 ); 00624 // check point is at triangle vertex 00625 check_common_vertex( tool.moab(), &tri, 1, point ); 00626 } 00627 00628 // test location on each face 00629 const CartVect facepts[] = { CartVect( -1.0, 0.5, 0.5 ), CartVect( 1.0, 0.5, 0.5 ), CartVect( 0.5, -1.0, 0.5 ), 00630 CartVect( 0.5, 1.0, 0.5 ), CartVect( 0.5, 0.5, -1.0 ), CartVect( 0.5, 0.5, 1.0 ) }; 00631 for( unsigned i = 0; i < 6; ++i ) 00632 { 00633 CartVect point; 00634 rval = tool.closest_triangle( root, facepts[i].array(), point.array(), tri ); 00635 CHECK( MB_SUCCESS == rval ); 00636 // check result position 00637 const CartVect diff = facepts[i] - point; 00638 CHECK( fabs( diff[0] ) < 1e-6 ); 00639 CHECK( fabs( diff[1] ) < 1e-6 ); 00640 CHECK( fabs( diff[2] ) < 1e-6 ); 00641 // check that point is contained within triangle 00642 check_point_in_triangles( tool.moab(), &tri, 1, point ); 00643 } 00644 00645 // test a location really far from the tree 00646 { 00647 const double far[] = { 0.75, 0.75, 200 }; 00648 CartVect point; 00649 rval = tool.closest_triangle( root, far, point.array(), tri ); 00650 CHECK( MB_SUCCESS == rval ); 00651 CHECK( fabs( point[0] - 0.75 ) < 1e-6 ); 00652 CHECK( fabs( point[1] - 0.75 ) < 1e-6 ); 00653 CHECK( fabs( point[2] - 1.00 ) < 1e-6 ); 00654 // check that point is contained within triangle 00655 check_point_in_triangles( tool.moab(), &tri, 1, point ); 00656 } 00657 00658 // now do it all again with a lot more triangles 00659 tool.reset_tree(); 00660 moab.delete_mesh(); 00661 box_tris.clear(); 00662 00663 build_triangle_box_large( &moab, box_tris ); 00664 rval = tool.build_tree( box_tris, &root, &opts ); 00665 CHECK( MB_SUCCESS == rval ); 00666 test_valid_tree( &tool, root, opts, box_tris ); 00667 00668 // test closest to each corner of the box. 00669 for( unsigned i = 0; i < 8; ++i ) 00670 { 00671 int x = i & 1 ? 1 : -1; 00672 int y = i & 2 ? 1 : -1; 00673 int z = i & 4 ? 1 : -1; 00674 double coords[] = { 4.0 * x, 4.0 * y, 4.0 * z }; 00675 CartVect point; 00676 rval = tool.closest_triangle( root, coords, point.array(), tri ); 00677 CHECK( MB_SUCCESS == rval ); 00678 // check result position 00679 CHECK( fabs( 3.0 * x - point[0] ) < 1e-6 ); 00680 CHECK( fabs( 3.0 * y - point[1] ) < 1e-6 ); 00681 CHECK( fabs( 3.0 * z - point[2] ) < 1e-6 ); 00682 // check point is at triangle vertex 00683 check_common_vertex( tool.moab(), &tri, 1, point ); 00684 } 00685 00686 // test location on each face 00687 const CartVect facepts2[] = { CartVect( -3.0, 0.5, 0.5 ), CartVect( 3.0, 0.5, 0.5 ), CartVect( 0.5, -3.0, 0.5 ), 00688 CartVect( 0.5, 3.0, 0.5 ), CartVect( 0.5, 0.5, -3.0 ), CartVect( 0.5, 0.5, 3.0 ) }; 00689 for( unsigned i = 0; i < 6; ++i ) 00690 { 00691 CartVect point; 00692 rval = tool.closest_triangle( root, facepts2[i].array(), point.array(), tri ); 00693 CHECK( MB_SUCCESS == rval ); 00694 // check result position 00695 const CartVect diff = facepts2[i] - point; 00696 CHECK( fabs( diff[0] ) < 1e-6 ); 00697 CHECK( fabs( diff[1] ) < 1e-6 ); 00698 CHECK( fabs( diff[2] ) < 1e-6 ); 00699 // check that point is contained within triangle 00700 check_point_in_triangles( tool.moab(), &tri, 1, point ); 00701 } 00702 00703 // test a location really far from the tree 00704 { 00705 const double far[] = { 2.75, 2.75, 200 }; 00706 CartVect point; 00707 rval = tool.closest_triangle( root, far, point.array(), tri ); 00708 CHECK( MB_SUCCESS == rval ); 00709 CHECK( fabs( point[0] - 2.75 ) < 1e-6 ); 00710 CHECK( fabs( point[1] - 2.75 ) < 1e-6 ); 00711 CHECK( fabs( point[2] - 3.00 ) < 1e-6 ); 00712 // check that point is contained within triangle 00713 check_point_in_triangles( tool.moab(), &tri, 1, point ); 00714 } 00715 } 00716 00717 void test_sphere_intersect_triangles() 00718 { 00719 Core moab; 00720 AdaptiveKDTree tool( &moab ); 00721 Range box_tris; 00722 00723 std::ostringstream options; 00724 options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;"; 00725 FileOptions opts( options.str().c_str() ); 00726 EntityHandle root; 00727 ErrorCode rval; 00728 std::vector< EntityHandle > triangles; 00729 00730 moab.delete_mesh(); 00731 box_tris.clear(); 00732 build_triangle_box_small( &moab, box_tris ); 00733 rval = tool.build_tree( box_tris, &root, &opts ); 00734 CHECK( MB_SUCCESS == rval ); 00735 test_valid_tree( &tool, root, opts, box_tris ); 00736 00737 // test closest to each corner of the box. 00738 for( unsigned i = 0; i < 8; ++i ) 00739 { 00740 double x = i & 1 ? 1 : -1; 00741 double y = i & 2 ? 1 : -1; 00742 double z = i & 4 ? 1 : -1; 00743 double center[] = { x, y, z }; 00744 triangles.clear(); 00745 rval = tool.sphere_intersect_triangles( root, center, 0.26, triangles ); 00746 CHECK( MB_SUCCESS == rval ); 00747 // expect 3 to 6 triangles, depending on the corner 00748 CHECK( triangles.size() >= 3 ); 00749 CHECK( triangles.size() <= 6 ); 00750 // check point is at the same vertex for each triangle 00751 check_common_vertex( tool.moab(), &triangles[0], triangles.size(), CartVect( center ) ); 00752 } 00753 00754 // now do it all again with a lot more triangles 00755 00756 tool.reset_tree(); 00757 moab.delete_mesh(); 00758 box_tris.clear(); 00759 build_triangle_box_large( &moab, box_tris ); 00760 00761 rval = tool.build_tree( box_tris, &root, &opts ); 00762 CHECK( MB_SUCCESS == rval ); 00763 test_valid_tree( &tool, root, opts, box_tris ); 00764 00765 // test closest to each corner of the box. 00766 for( unsigned i = 0; i < 8; ++i ) 00767 { 00768 int x = i & 1 ? 1 : -1; 00769 int y = i & 2 ? 1 : -1; 00770 int z = i & 4 ? 1 : -1; 00771 double center[] = { 3.0 * x, 3.0 * y, 3.0 * z }; 00772 triangles.clear(); 00773 rval = tool.sphere_intersect_triangles( root, center, 0.26, triangles ); 00774 CHECK( MB_SUCCESS == rval ); 00775 // expect 3 to 6 triangles, depending on the corner 00776 CHECK( triangles.size() >= 3 ); 00777 CHECK( triangles.size() <= 6 ); 00778 // check point is at the same vertex for each triangle 00779 check_common_vertex( tool.moab(), &triangles[0], triangles.size(), CartVect( center ) ); 00780 } 00781 } 00782 00783 void test_ray_intersect_triangles() 00784 { 00785 Core moab; 00786 AdaptiveKDTree tool( &moab ); 00787 Range box_tris; 00788 00789 std::ostringstream options; 00790 options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;"; 00791 FileOptions opts( options.str().c_str() ); 00792 EntityHandle root; 00793 ErrorCode rval; 00794 std::vector< EntityHandle > tris; 00795 std::vector< double > dists; 00796 00797 moab.delete_mesh(); 00798 box_tris.clear(); 00799 build_triangle_box_small( &moab, box_tris ); 00800 rval = tool.build_tree( box_tris, &root, &opts ); 00801 CHECK( MB_SUCCESS == rval ); 00802 test_valid_tree( &tool, root, opts, box_tris ); 00803 00804 // test ray through box parallel to X axis 00805 CartVect dir( 1, 0, 0 ); 00806 CartVect pt( -2, 0.75, 0.75 ); 00807 tris.clear(); 00808 dists.clear(); 00809 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists ); 00810 CHECK( MB_SUCCESS == rval ); 00811 CHECK( tris.size() == 3 ); 00812 CHECK( dists.size() == tris.size() ); 00813 for( unsigned i = 0; i < dists.size(); ++i ) 00814 { 00815 CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 ); 00816 CartVect tript = pt + dists[i] * dir; 00817 check_point_in_triangles( &moab, &tris[i], 1, tript ); 00818 } 00819 00820 // test ray through box parallel to X axis, closest tri only 00821 tris.clear(); 00822 dists.clear(); 00823 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 1 ); 00824 CHECK( MB_SUCCESS == rval ); 00825 CHECK( tris.size() == 1 ); 00826 CHECK( dists.size() == tris.size() ); 00827 CHECK( fabs( dists[0] - 1 ) < 1e-6 ); 00828 check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir ); 00829 00830 // test ray ending within box, parallel to X axis 00831 tris.clear(); 00832 dists.clear(); 00833 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 0, 2.0 ); 00834 CHECK( MB_SUCCESS == rval ); 00835 CHECK( tris.size() == 1 ); 00836 CHECK( dists.size() == tris.size() ); 00837 CHECK( fabs( dists[0] - 1 ) < 1e-6 ); 00838 check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir ); 00839 00840 // test ray starting within box parallel to X axis 00841 tris.clear(); 00842 dists.clear(); 00843 pt = CartVect( 0, .75, .75 ); 00844 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists ); 00845 CHECK( MB_SUCCESS == rval ); 00846 CHECK( tris.size() == 2 ); 00847 CHECK( dists.size() == tris.size() ); 00848 for( unsigned i = 0; i < dists.size(); ++i ) 00849 { 00850 CHECK( fabs( dists[i] - 1 ) < 1e-6 ); 00851 CartVect tript = pt + dists[i] * dir; 00852 check_point_in_triangles( &moab, &tris[i], 1, tript ); 00853 } 00854 00855 // test skew ray through box 00856 dir = CartVect( 0.5 * sqrt( 2.0 ), 0.5 * sqrt( 2.0 ), 0.0 ); 00857 pt = CartVect( 0, -1.5, 0 ); 00858 tris.clear(); 00859 dists.clear(); 00860 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists ); 00861 CHECK( MB_SUCCESS == rval ); 00862 CHECK( tris.size() == 2 ); 00863 CHECK( dists.size() == tris.size() ); 00864 if( dists[0] < dists[1] ) 00865 { 00866 CHECK( fabs( dists[0] - 0.5 * sqrt( 2.0 ) ) < 1e-6 ); 00867 CHECK( fabs( dists[1] - sqrt( 2.0 ) ) < 1e-6 ); 00868 } 00869 check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir ); 00870 check_point_in_triangles( &moab, &tris[1], 1, pt + dists[1] * dir ); 00871 00872 // test ray through box parallel to -Y axis 00873 dir = CartVect( 0, -1, 0 ); 00874 pt = CartVect( 0, 2, 0 ); 00875 tris.clear(); 00876 dists.clear(); 00877 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists ); 00878 CHECK( MB_SUCCESS == rval ); 00879 CHECK( tris.size() == 4 ); 00880 CHECK( dists.size() == tris.size() ); 00881 for( unsigned i = 0; i < dists.size(); ++i ) 00882 { 00883 CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 ); 00884 CartVect tript = pt + dists[i] * dir; 00885 check_point_in_triangles( &moab, &tris[i], 1, tript ); 00886 } 00887 00888 // test ray through box parallel to Z axis 00889 dir = CartVect( 0, 0, 1 ); 00890 pt = CartVect( 0, 0, -2 ); 00891 tris.clear(); 00892 dists.clear(); 00893 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists ); 00894 CHECK( MB_SUCCESS == rval ); 00895 CHECK( tris.size() == 4 ); 00896 CHECK( dists.size() == tris.size() ); 00897 for( unsigned i = 0; i < dists.size(); ++i ) 00898 { 00899 CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 ); 00900 CartVect tript = pt + dists[i] * dir; 00901 check_point_in_triangles( &moab, &tris[i], 1, tript ); 00902 } 00903 00904 // test ray through box parallel to Z axis, limit 2 first 2 intersections 00905 tris.clear(); 00906 dists.clear(); 00907 rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 2 ); 00908 CHECK( MB_SUCCESS == rval ); 00909 CHECK( tris.size() == 2 ); 00910 CHECK( dists.size() == tris.size() ); 00911 for( unsigned i = 0; i < dists.size(); ++i ) 00912 { 00913 CHECK( fabs( dists[i] - 1 ) < 1e-6 ); 00914 CartVect tript = pt + dists[i] * dir; 00915 check_point_in_triangles( &moab, &tris[i], 1, tript ); 00916 } 00917 } 00918 00919 void test_leaf_volume() 00920 { 00921 Core moab; 00922 AdaptiveKDTree tool( &moab ); 00923 EntityHandle root; 00924 00925 create_simple_2d_tree( tool, root ); 00926 AdaptiveKDTreeIter iter; 00927 ErrorCode rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval ); 00928 00929 CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 ); // 1 00930 CHECK_ERR( iter.step() ); 00931 CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 ); // 2 00932 CHECK_ERR( iter.step() ); 00933 CHECK_REAL_EQUAL( 40.0, iter.volume(), 1e-12 ); // 3 00934 CHECK_ERR( iter.step() ); 00935 CHECK_REAL_EQUAL( 8.0, iter.volume(), 1e-12 ); // 4 00936 CHECK_ERR( iter.step() ); 00937 CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 ); // 5 00938 CHECK_ERR( iter.step() ); 00939 CHECK_REAL_EQUAL( 24.0, iter.volume(), 1e-12 ); // 6 00940 CHECK_ERR( iter.step() ); 00941 CHECK_REAL_EQUAL( 30.0, iter.volume(), 1e-12 ); // 7 00942 CHECK_ERR( iter.step() ); 00943 CHECK_REAL_EQUAL( 10.0, iter.volume(), 1e-12 ); // 8 00944 } 00945 00946 void test_leaf_sibling() 00947 { 00948 ErrorCode rval; 00949 Core moab; 00950 AdaptiveKDTree tool( &moab ); 00951 EntityHandle root; 00952 00953 create_simple_2d_tree( tool, root ); 00954 AdaptiveKDTreeIter iter1, iter2; 00955 rval = tool.get_tree_iterator( root, iter1 );CHECK_ERR( rval ); 00956 rval = tool.get_tree_iterator( root, iter2 );CHECK_ERR( rval ); 00957 00958 // iter1 == 1, iter2 == 1 00959 CHECK( !iter1.is_sibling( iter1 ) ); 00960 CHECK( !iter1.is_sibling( iter1.handle() ) ); 00961 CHECK( !iter1.is_sibling( iter2 ) ); 00962 CHECK( iter1.sibling_is_forward() ); 00963 00964 // iter1 == 1, iter2 == 2 00965 rval = iter2.step();CHECK_ERR( rval ); 00966 CHECK( iter1.is_sibling( iter2 ) ); 00967 CHECK( iter1.is_sibling( iter2.handle() ) ); 00968 CHECK( iter2.is_sibling( iter1 ) ); 00969 CHECK( iter2.is_sibling( iter1.handle() ) ); 00970 CHECK( !iter2.sibling_is_forward() ); 00971 00972 // iter1 == 1, iter2 == 3 00973 rval = iter2.step();CHECK_ERR( rval ); 00974 CHECK( !iter1.is_sibling( iter2 ) ); 00975 CHECK( !iter1.is_sibling( iter2.handle() ) ); 00976 CHECK( !iter2.is_sibling( iter1 ) ); 00977 CHECK( !iter2.is_sibling( iter1.handle() ) ); 00978 CHECK( iter2.sibling_is_forward() ); 00979 00980 // iter1 == 2, iter2 == 3 00981 rval = iter1.step();CHECK_ERR( rval ); 00982 CHECK( !iter1.is_sibling( iter2 ) ); 00983 CHECK( !iter1.is_sibling( iter2.handle() ) ); 00984 CHECK( !iter2.is_sibling( iter1 ) ); 00985 CHECK( !iter2.is_sibling( iter1.handle() ) ); 00986 CHECK( !iter1.sibling_is_forward() ); 00987 00988 // iter1 == 4, iter2 == 3 00989 rval = iter1.step();CHECK_ERR( rval ); 00990 CHECK_EQUAL( iter1.handle(), iter2.handle() ); 00991 rval = iter1.step();CHECK_ERR( rval ); 00992 CHECK( iter1.is_sibling( iter2 ) ); 00993 CHECK( iter1.is_sibling( iter2.handle() ) ); 00994 CHECK( iter2.is_sibling( iter1 ) ); 00995 CHECK( iter2.is_sibling( iter1.handle() ) ); 00996 CHECK( !iter1.sibling_is_forward() ); 00997 } 00998 00999 void test_leaf_intersects_plane() 01000 { 01001 ErrorCode rval; 01002 Core moab; 01003 AdaptiveKDTree tool( &moab ); 01004 01005 EntityHandle root; 01006 const double min[3] = { -5, -4, -1 }; 01007 const double max[3] = { 1, 2, 3 }; 01008 rval = tool.create_root( min, max, root );CHECK_ERR( rval ); 01009 01010 AdaptiveKDTreeIter iter; 01011 rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval ); 01012 01013 AdaptiveKDTree::Plane x1 = { min[0] - 1, AdaptiveKDTree::X }; 01014 CHECK( !iter.intersects( x1 ) ); 01015 AdaptiveKDTree::Plane x2 = { min[0], AdaptiveKDTree::X }; 01016 CHECK( iter.intersects( x2 ) ); 01017 AdaptiveKDTree::Plane x3 = { min[0] + 1, AdaptiveKDTree::X }; 01018 CHECK( iter.intersects( x3 ) ); 01019 AdaptiveKDTree::Plane x4 = { max[0] - 1, AdaptiveKDTree::X }; 01020 CHECK( iter.intersects( x4 ) ); 01021 AdaptiveKDTree::Plane x5 = { max[0], AdaptiveKDTree::X }; 01022 CHECK( iter.intersects( x5 ) ); 01023 AdaptiveKDTree::Plane x6 = { max[0] + 1, AdaptiveKDTree::X }; 01024 CHECK( !iter.intersects( x6 ) ); 01025 01026 AdaptiveKDTree::Plane y1 = { min[1] - 1, AdaptiveKDTree::Y }; 01027 CHECK( !iter.intersects( y1 ) ); 01028 AdaptiveKDTree::Plane y2 = { min[1], AdaptiveKDTree::Y }; 01029 CHECK( iter.intersects( y2 ) ); 01030 AdaptiveKDTree::Plane y3 = { min[1] + 1, AdaptiveKDTree::Y }; 01031 CHECK( iter.intersects( y3 ) ); 01032 AdaptiveKDTree::Plane y4 = { max[1] - 1, AdaptiveKDTree::Y }; 01033 CHECK( iter.intersects( y4 ) ); 01034 AdaptiveKDTree::Plane y5 = { max[1], AdaptiveKDTree::Y }; 01035 CHECK( iter.intersects( y5 ) ); 01036 AdaptiveKDTree::Plane y6 = { max[1] + 1, AdaptiveKDTree::Y }; 01037 CHECK( !iter.intersects( y6 ) ); 01038 01039 AdaptiveKDTree::Plane z1 = { min[2] - 1, AdaptiveKDTree::Z }; 01040 CHECK( !iter.intersects( z1 ) ); 01041 AdaptiveKDTree::Plane z2 = { min[2], AdaptiveKDTree::Z }; 01042 CHECK( iter.intersects( z2 ) ); 01043 AdaptiveKDTree::Plane z3 = { min[2] + 1, AdaptiveKDTree::Z }; 01044 CHECK( iter.intersects( z3 ) ); 01045 AdaptiveKDTree::Plane z4 = { max[2] - 1, AdaptiveKDTree::Z }; 01046 CHECK( iter.intersects( z4 ) ); 01047 AdaptiveKDTree::Plane z5 = { max[2], AdaptiveKDTree::Z }; 01048 CHECK( iter.intersects( z5 ) ); 01049 AdaptiveKDTree::Plane z6 = { max[2] + 1, AdaptiveKDTree::Z }; 01050 CHECK( !iter.intersects( z6 ) ); 01051 } 01052 01053 #define CHECK_RAY_XSECTS( PT, DIR, T_IN, T_OUT ) \ 01054 do \ 01055 { \ 01056 CHECK( iter.intersect_ray( ( PT ), ( DIR ), t_in, t_out ) ); \ 01057 CHECK_REAL_EQUAL( ( T_IN ), t_in, 1e-6 ); \ 01058 CHECK_REAL_EQUAL( ( T_OUT ), t_out, 1e-6 ); \ 01059 } while( false ) 01060 01061 void test_leaf_intersects_ray() 01062 { 01063 ErrorCode rval; 01064 Core moab; 01065 AdaptiveKDTree tool( &moab ); 01066 double t_in, t_out; 01067 01068 EntityHandle root; 01069 const double min[3] = { -5, -4, -1 }; 01070 const double max[3] = { 1, 2, 3 }; 01071 rval = tool.create_root( min, max, root );CHECK_ERR( rval ); 01072 01073 AdaptiveKDTreeIter iter; 01074 rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval ); 01075 01076 // start with point inside box 01077 const double pt1[] = { 0, 0, 0 }; 01078 const double dir1[] = { 0.1, 0.1, 0.1 }; 01079 CHECK_RAY_XSECTS( pt1, dir1, 0, 10 ); 01080 const double dir2[] = { 5, 5, 5 }; 01081 CHECK_RAY_XSECTS( pt1, dir2, 0, 0.2 ); 01082 const double pxdir[] = { 1, 0, 0 }; 01083 CHECK_RAY_XSECTS( pt1, pxdir, 0, 1 ); 01084 const double nxdir[] = { -1, 0, 0 }; 01085 CHECK_RAY_XSECTS( pt1, nxdir, 0, 5 ); 01086 const double pydir[] = { 0, 1, 0 }; 01087 CHECK_RAY_XSECTS( pt1, pydir, 0, 2 ); 01088 const double nydir[] = { 0, -1, 0 }; 01089 CHECK_RAY_XSECTS( pt1, nydir, 0, 4 ); 01090 const double pzdir[] = { 0, 0, 1 }; 01091 CHECK_RAY_XSECTS( pt1, pzdir, 0, 3 ); 01092 const double nzdir[] = { 0, 0, -1 }; 01093 CHECK_RAY_XSECTS( pt1, nzdir, 0, 1 ); 01094 01095 // point below box 01096 const double pt2[] = { 0, 0, -2 }; 01097 CHECK_RAY_XSECTS( pt2, dir1, 10, 10 ); 01098 CHECK_RAY_XSECTS( pt2, dir2, 0.2, 0.2 ); 01099 CHECK( !iter.intersect_ray( pt2, pxdir, t_in, t_out ) ); 01100 CHECK( !iter.intersect_ray( pt2, nxdir, t_in, t_out ) ); 01101 CHECK( !iter.intersect_ray( pt2, pydir, t_in, t_out ) ); 01102 CHECK( !iter.intersect_ray( pt2, nydir, t_in, t_out ) ); 01103 CHECK_RAY_XSECTS( pt2, pzdir, 1, 5 ); 01104 CHECK( !iter.intersect_ray( pt2, nzdir, t_in, t_out ) ); 01105 01106 // point right of box 01107 const double pt3[] = { 3, 0, 0 }; 01108 CHECK( !iter.intersect_ray( pt3, dir1, t_in, t_out ) ); 01109 CHECK( !iter.intersect_ray( pt3, dir2, t_in, t_out ) ); 01110 CHECK( !iter.intersect_ray( pt3, pxdir, t_in, t_out ) ); 01111 CHECK_RAY_XSECTS( pt3, nxdir, 2, 8 ); 01112 CHECK( !iter.intersect_ray( pt3, pydir, t_in, t_out ) ); 01113 CHECK( !iter.intersect_ray( pt3, nydir, t_in, t_out ) ); 01114 CHECK( !iter.intersect_ray( pt3, pzdir, t_in, t_out ) ); 01115 CHECK( !iter.intersect_ray( pt3, nzdir, t_in, t_out ) ); 01116 01117 // a few more complex test cases 01118 const double dira[] = { -3, 0, 3 }; 01119 CHECK_RAY_XSECTS( pt3, dira, 2. / 3., 1.0 ); 01120 const double dirb[] = { -2, 0, 3.1 }; 01121 CHECK( !iter.intersect_ray( pt3, dirb, t_in, t_out ) ); 01122 } 01123 01124 int main() 01125 { 01126 #ifdef MOAB_HAVE_MPI 01127 int fail = MPI_Init( 0, 0 ); 01128 if( fail ) return fail; 01129 #endif 01130 01131 int error_count = 0; 01132 01133 error_count += RUN_TEST( leaf_iterator_test ); 01134 error_count += RUN_TEST( test_build_tree_bisect_triangles ); 01135 error_count += RUN_TEST( test_closest_triangle ); 01136 error_count += RUN_TEST( test_sphere_intersect_triangles ); 01137 error_count += RUN_TEST( test_ray_intersect_triangles ); 01138 error_count += RUN_TEST( test_tree_merge_nodes ); 01139 error_count += RUN_TEST( test_leaf_volume ); 01140 error_count += RUN_TEST( test_leaf_sibling ); 01141 error_count += RUN_TEST( test_leaf_intersects_plane ); 01142 error_count += RUN_TEST( test_leaf_intersects_ray ); 01143 01144 #ifdef MOAB_HAVE_MPI 01145 fail = MPI_Finalize(); 01146 if( fail ) return fail; 01147 #endif 01148 01149 return error_count; 01150 }