MOAB: Mesh Oriented datABase  (version 5.3.1)
sploc_searching_perf.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/AdaptiveKDTree.hpp"
00008 #include "moab/BVHTree.hpp"
00009 #include "moab/ProgOptions.hpp"
00010 #include "moab/CpuTimer.hpp"
00011 #include "moab/ElemEvaluator.hpp"
00012 
00013 #ifdef MOAB_HAVE_MPI
00014 #include "moab_mpi.h"
00015 #endif
00016 
00017 #include <cstdlib>
00018 #include <sstream>
00019 
00020 using namespace moab;
00021 
00022 ErrorCode test_locator( SpatialLocator& sl, int npoints, double rtol, double& cpu_time, double& percent_outside );
00023 ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim );
00024 
00025 int main( int argc, char** argv )
00026 {
00027 #ifdef MOAB_HAVE_MPI
00028     int fail = MPI_Init( &argc, &argv );
00029     if( fail ) return fail;
00030 #else
00031     // silence the warning of parameters not used, in serial; there should be a smarter way :(
00032     argv[0] = argv[argc - argc];
00033 #endif
00034 
00035     int npoints = 100, dim = 3;
00036     int dints = 1, dleafs = 1, ddeps = 1;
00037     bool eval   = false;
00038     double rtol = 1.0e-10;
00039 
00040     ProgOptions po( "tree_searching_perf options" );
00041     po.addOpt< void >( ",e", "Use ElemEvaluator in tree search", &eval );
00042     po.addOpt< int >( "ints,i", "Number of doublings of intervals on each side of scd mesh", &dints );
00043     po.addOpt< int >( "leaf,l", "Number of doublings of maximum number of elements per leaf", &dleafs );
00044     po.addOpt< int >( "max_depth,m", "Number of 5-intervals on maximum depth of tree", &ddeps );
00045     po.addOpt< int >( "npoints,n", "Number of query points", &npoints );
00046     po.addOpt< int >( "dim,d", "Dimension of the mesh", &dim );
00047     po.addOpt< double >( "tol,t", "Relative tolerance of point search", &rtol );
00048     //  po.addOpt<void>( "print,p", "Print tree details", &print_tree);
00049     po.parseCommandLine( argc, argv );
00050 
00051     std::vector< int > ints, deps, leafs;
00052     ints.push_back( 10 );
00053     for( int i = 1; i < dints; i++ )
00054         ints.push_back( 2 * ints[i - 1] );
00055     deps.push_back( 30 );
00056     for( int i = 1; i < ddeps; i++ )
00057         deps.push_back( deps[i - 1] - 5 );
00058     leafs.push_back( 6 );
00059     for( int i = 1; i < dleafs; i++ )
00060         leafs.push_back( 2 * leafs[i - 1] );
00061 
00062     ErrorCode rval = MB_SUCCESS;
00063     std::cout << "Tree_type"
00064               << " "
00065               << "Elems_per_leaf"
00066               << " "
00067               << "Tree_depth"
00068               << " "
00069               << "Ints_per_side"
00070               << " "
00071               << "N_elements"
00072               << " "
00073               << "search_time"
00074               << " "
00075               << "perc_outside"
00076               << " "
00077               << "initTime"
00078               << " "
00079               << "nodesVisited"
00080               << " "
00081               << "leavesVisited"
00082               << " "
00083               << "numTraversals"
00084               << " "
00085               << "leafObjectTests" << std::endl;
00086 
00087     // outermost iteration: # elements
00088     for( std::vector< int >::iterator int_it = ints.begin(); int_it != ints.end(); ++int_it )
00089     {
00090         Core mb;
00091         Range elems;
00092         rval = create_hex_mesh( mb, elems, *int_it, dim );
00093         if( MB_SUCCESS != rval ) return rval;
00094 
00095         // iteration: tree depth
00096         for( std::vector< int >::iterator dep_it = deps.begin(); dep_it != deps.end(); ++dep_it )
00097         {
00098 
00099             // iteration: tree max elems/leaf
00100             for( std::vector< int >::iterator leafs_it = leafs.begin(); leafs_it != leafs.end(); ++leafs_it )
00101             {
00102 
00103                 // iteration: tree type
00104                 for( int tree_tp = 0; tree_tp < 2; tree_tp++ )
00105                 {
00106                     // create tree
00107                     Tree* tree;
00108                     if( 0 == tree_tp )
00109                         tree = new BVHTree( &mb );
00110                     else
00111                         tree = new AdaptiveKDTree( &mb );
00112 
00113                     ElemEvaluator* eeval = NULL;
00114                     if( eval )
00115                     {
00116                         // create an element evaluator
00117                         eeval = new ElemEvaluator( &mb );
00118                         rval  = eeval->set_eval_set( *elems.begin() );
00119                         if( MB_SUCCESS != rval ) return rval;
00120                     }
00121 
00122                     std::ostringstream opts;
00123                     opts << "MAX_DEPTH=" << *dep_it << ";MAX_PER_LEAF=" << *leafs_it;
00124                     FileOptions fo( opts.str().c_str() );
00125                     rval = tree->parse_options( fo );
00126                     if( MB_SUCCESS != rval ) return rval;
00127                     SpatialLocator sl( &mb, elems, tree, eeval );
00128 
00129                     // call evaluation
00130                     double cpu_time, perc_outside;
00131                     rval = test_locator( sl, npoints, rtol, cpu_time, perc_outside );
00132                     if( MB_SUCCESS != rval ) return rval;
00133 
00134                     std::cout << ( tree_tp == 0 ? "BVH" : "KD" ) << " " << *leafs_it << " " << *dep_it << " " << *int_it
00135                               << " " << ( *int_it ) * ( *int_it ) * ( *int_it ) << " " << cpu_time << " "
00136                               << perc_outside << " ";
00137 
00138                     tree->tree_stats().output_all_stats();
00139 
00140                     if( eeval ) delete eeval;
00141 
00142                 }  // tree_tp
00143 
00144             }  // max elems/leaf
00145 
00146         }  // max depth
00147 
00148     }  // # elements
00149 
00150 #ifdef MOAB_HAVE_MPI
00151     fail = MPI_Finalize();
00152     if( fail ) return fail;
00153 #endif
00154 
00155     return 0;
00156 }
00157 
00158 ErrorCode test_locator( SpatialLocator& sl, int npoints, double rtol, double& cpu_time, double& percent_outside )
00159 {
00160     BoundBox box     = sl.local_box();
00161     CartVect box_del = box.bMax - box.bMin;
00162 
00163     std::vector< CartVect > test_pts( npoints ), test_res( npoints );
00164     std::vector< EntityHandle > ents( npoints );
00165     int* is_in = new int[npoints];
00166 
00167     double denom = 1.0 / (double)RAND_MAX;
00168     for( int i = 0; i < npoints; i++ )
00169     {
00170         // generate a small number of random point to test
00171         double rx = (double)rand() * denom, ry = (double)rand() * denom, rz = (double)rand() * denom;
00172         test_pts[i] = box.bMin + CartVect( rx * box_del[0], ry * box_del[1], rz * box_del[2] );
00173     }
00174 
00175     CpuTimer ct;
00176 
00177     // call spatial locator to locate points
00178     ErrorCode rval =
00179         sl.locate_points( test_pts[0].array(), npoints, &ents[0], test_res[0].array(), &is_in[0], rtol, 0.0 );
00180     if( MB_SUCCESS != rval )
00181     {
00182         delete[] is_in;
00183         return rval;
00184     }
00185 
00186     cpu_time = ct.time_elapsed();
00187 
00188     int num_out     = std::count( is_in, is_in + npoints, false );
00189     percent_outside = ( (double)num_out ) / npoints;
00190     delete[] is_in;
00191 
00192     return rval;
00193 }
00194 
00195 ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim )
00196 {
00197     ScdInterface* scdi;
00198     ErrorCode rval = mb.query_interface( scdi );
00199     if( MB_SUCCESS != rval ) return rval;
00200 
00201     HomCoord high( n - 1, -1, -1 );
00202     if( dim > 1 ) high[1] = n - 1;
00203     if( dim > 2 ) high[2] = n - 1;
00204     ScdBox* new_box;
00205     rval = scdi->construct_box( HomCoord( 0, 0, 0 ), high, NULL, 0, new_box );
00206     if( MB_SUCCESS != rval ) return rval;
00207     rval = mb.release_interface( scdi );
00208     if( MB_SUCCESS != rval ) return rval;
00209 
00210     rval = mb.get_entities_by_dimension( 0, dim, elems );
00211     if( MB_SUCCESS != rval ) return rval;
00212 
00213     return rval;
00214 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines