MOAB: Mesh Oriented datABase  (version 5.3.0)
par_spatial_locator_test.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/SpatialLocator.hpp"
00003 #include "moab/Tree.hpp"
00004 #include "moab/HomXform.hpp"
00005 #include "moab/ScdInterface.hpp"
00006 #include "moab/CartVect.hpp"
00007 #include "moab/BVHTree.hpp"
00008 #include "moab/ProgOptions.hpp"
00009 #include "moab/CpuTimer.hpp"
00010 #include "moab/ParallelComm.hpp"
00011 #include "moab_mpi.h"
00012 
00013 #include "TestUtil.hpp"
00014 
00015 #include <cstdlib>
00016 #include <sstream>
00017 #include <string>
00018 
00019 using namespace moab;
00020 
00021 void test_kd_tree();
00022 void test_bvh_tree();
00023 void test_locator( SpatialLocator* sl );
00024 
00025 ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim );
00026 ErrorCode load_file( Interface& mb, std::string& fn, Range& elems );
00027 
00028 int max_depth   = 30;
00029 int npoints     = 1000;
00030 int leaf        = 6;
00031 int tree        = -1;
00032 bool print_tree = false;
00033 int ints        = 10;
00034 std::string fname;
00035 
00036 int main( int argc, char** argv )
00037 {
00038     int fail = MPI_Init( &argc, &argv );
00039     if( fail ) return fail;
00040 
00041     ProgOptions po( "spatial_locator_test options" );
00042     po.addOpt< std::string >( "file,f", "Input filename", &fname );
00043     po.addOpt< int >( "ints,i", "Number of intervals on each side of scd mesh", &ints );
00044     po.addOpt< int >( "leaf,l", "Maximum number of elements per leaf", &leaf );
00045     po.addOpt< int >( "max_depth,m", "Maximum depth of tree", &max_depth );
00046     po.addOpt< int >( "npoints,n", "Number of query points", &npoints );
00047     po.addOpt< void >( "print,p", "Print tree details", &print_tree );
00048     po.addOpt< int >( "tree,t", "Tree type (-1=all (default), 0=AdaptiveKD, 1=BVH", &tree );
00049     po.parseCommandLine( argc, argv );
00050 
00051     if( -1 == tree || 0 == tree ) RUN_TEST( test_kd_tree );
00052     if( -1 == tree || 1 == tree ) RUN_TEST( test_bvh_tree );
00053 
00054     fail = MPI_Finalize();
00055     if( fail ) return fail;
00056 
00057     return 0;
00058 }
00059 
00060 void test_kd_tree()
00061 {
00062     ErrorCode rval;
00063     Core mb;
00064 
00065     // create a simple mesh to test
00066     Range elems;
00067     if( fname.empty() )
00068     {
00069         rval = create_hex_mesh( mb, elems, ints, 3 );CHECK_ERR( rval );
00070     }
00071     else
00072     {
00073         rval = load_file( mb, fname, elems );CHECK_ERR( rval );
00074     }
00075 
00076     // initialize spatial locator with the elements and the default tree type
00077     SpatialLocator* sl = new SpatialLocator( &mb, elems );
00078 
00079     test_locator( sl );
00080 
00081     // destroy spatial locator, and tree along with it
00082     delete sl;
00083 }
00084 
00085 void test_bvh_tree()
00086 {
00087     ErrorCode rval;
00088     Core mb;
00089 
00090     // create a simple mesh to test
00091     Range elems;
00092     if( fname.empty() )
00093     {
00094         rval = create_hex_mesh( mb, elems, ints, 3 );CHECK_ERR( rval );
00095     }
00096     else
00097     {
00098         rval = load_file( mb, fname, elems );CHECK_ERR( rval );
00099     }
00100 
00101     // initialize spatial locator with the elements and a BVH tree
00102     BVHTree bvh( &mb );
00103     std::ostringstream opts;
00104     opts << "MAX_DEPTH=" << max_depth << ";MAX_PER_LEAF=" << leaf;
00105     FileOptions fo( opts.str().c_str() );
00106     rval               = bvh.parse_options( fo );
00107     SpatialLocator* sl = new SpatialLocator( &mb, elems, &bvh );
00108     test_locator( sl );
00109 
00110     // destroy spatial locator, and tree along with it
00111     delete sl;
00112 }
00113 
00114 bool is_neg( int is_neg )
00115 {
00116     return ( is_neg < 0 );
00117 }
00118 
00119 void test_locator( SpatialLocator* sl )
00120 {
00121     CartVect box_del;
00122     BoundBox box = sl->local_box();
00123     box_del      = box.bMax - box.bMin;
00124 
00125     double denom = 1.0 / (double)RAND_MAX;
00126     std::vector< CartVect > test_pts( npoints );
00127 
00128     for( int i = 0; i < npoints; i++ )
00129     {
00130         // generate a small number of random point to test
00131         double rx = (double)rand() * denom, ry = (double)rand() * denom, rz = (double)rand() * denom;
00132         test_pts[i] = box.bMin + CartVect( rx * box_del[0], ry * box_del[1], rz * box_del[2] );
00133     }
00134 
00135     // call spatial locator to locate points
00136     ParallelComm* pc = ParallelComm::get_pcomm( sl->moab(), 0 );
00137     CHECK( pc != NULL );
00138 
00139     ErrorCode rval = sl->par_locate_points( pc, test_pts[0].array(), npoints );CHECK_ERR( rval );
00140     if( pc->rank() == 0 )
00141     {
00142         int num_out = std::count_if( sl->par_loc_table().vi_rd, sl->par_loc_table().vi_rd + 2 * npoints, is_neg );
00143         num_out /= 2;
00144 
00145         std::cout << "Number of points inside an element = " << npoints - num_out << "/" << npoints << " ("
00146                   << 100.0 * ( (double)npoints - num_out ) / npoints << "%)" << std::endl;
00147         std::cout << "Traversal stats:" << std::endl;
00148         sl->get_tree()->tree_stats().output_trav_stats();
00149 
00150         if( print_tree )
00151         {
00152             std::cout << "Tree information: " << std::endl;
00153             rval = sl->get_tree()->print();CHECK_ERR( rval );
00154         }
00155     }
00156 }
00157 
00158 ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim )
00159 {
00160     ScdInterface* scdi;
00161     ErrorCode rval = mb.query_interface( scdi );CHECK_ERR( rval );
00162     ScdParData spd;
00163     spd.gDims[0] = spd.gDims[1] = spd.gDims[2] = 0;
00164     spd.gDims[3]                               = n;
00165     spd.partMethod                             = ScdParData::SQIJK;
00166     if( dim > 1 ) spd.gDims[4] = n;
00167     if( dim > 2 ) spd.gDims[5] = n;
00168     ScdBox* new_box;
00169     rval = scdi->construct_box( HomCoord( 0, 0, 0 ), HomCoord( 0, 0, 0 ), NULL, 0, new_box, NULL, &spd, false, 0 );CHECK_ERR( rval );
00170 
00171     rval = mb.get_entities_by_dimension( 0, dim, elems );CHECK_ERR( rval );
00172 
00173     return rval;
00174 }
00175 
00176 ErrorCode load_file( Interface& mb, std::string& fn, Range& elems )
00177 {
00178     std::string options;
00179     ErrorCode rval;
00180 
00181     options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;";
00182     rval    = mb.load_file( fn.c_str(), NULL, options.c_str() );
00183     if( MB_SUCCESS != rval ) return rval;
00184 
00185     rval = mb.get_entities_by_dimension( 0, 3, elems );
00186     if( MB_SUCCESS != rval ) return rval;
00187 
00188     if( elems.empty() )
00189     {
00190         rval = mb.get_entities_by_dimension( 0, 2, elems );
00191         if( MB_SUCCESS != rval ) return rval;
00192     }
00193 
00194     return MB_SUCCESS;
00195 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines