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