SpatialLocator.cpp

Go to the documentation of this file.

#include "moab/SpatialLocator.hpp"
#include "moab/Interface.hpp"
#include "moab/ElemEvaluator.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/BVHTree.hpp"
 
// include ScdInterface for box partitioning
#include "moab/ScdInterface.hpp"
 
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#endif
 
namespace moab
{
static bool debug = false;
 
SpatialLocator::SpatialLocator( Interface* impl, Range& elems, Tree* tree, ElemEvaluator* eval )
 : mbImpl( impl ), myElems( elems ), myDim( -1 ), myTree( tree ), elemEval( eval ), iCreatedTree( false ),
 timerInitialized( false )
{
 create_tree();
 
 if( !elems.empty() )
 {
 myDim = mbImpl->dimension_from_handle( *elems.rbegin() );
 ErrorCode rval = myTree->build_tree( myElems );
 if( MB_SUCCESS != rval ) throw rval;
 
 rval = myTree->get_bounding_box( localBox );
 if( MB_SUCCESS != rval ) throw rval;
 }
 
 regNums[0] = regNums[1] = regNums[2] = 0;
}
 
void SpatialLocator::create_tree()
{
 if( myTree ) return;
 
 if( myElems.empty() || mbImpl->type_from_handle( *myElems.rbegin() ) == MBVERTEX )
 // create a kdtree if only vertices
 myTree = new AdaptiveKDTree( mbImpl );
 else
 // otherwise a BVHtree, since it performs better for elements
 myTree = new BVHTree( mbImpl );
 
 iCreatedTree = true;
}
 
ErrorCode SpatialLocator::add_elems( Range& elems )
{
 if( elems.empty() ||
 mbImpl->dimension_from_handle( *elems.begin() ) != mbImpl->dimension_from_handle( *elems.rbegin() ) )
 return MB_FAILURE;
 
 myDim = mbImpl->dimension_from_handle( *elems.begin() );
 myElems = elems;
 
 ErrorCode rval = myTree->build_tree( myElems );
 return rval;
}
 
#ifdef MOAB_HAVE_MPI
ErrorCode SpatialLocator::initialize_intermediate_partition( ParallelComm* pc )
{
 if( !pc ) return MB_FAILURE;
 
 BoundBox gbox;
 
 // step 2
 // get the global bounding box
 double sendbuffer[6];
 double rcvbuffer[6];
 
 localBox.get( sendbuffer ); // fill sendbuffer with local box, max values in [0:2] min values in [3:5]
 sendbuffer[0] *= -1;
 sendbuffer[1] *= -1; // negate Xmin,Ymin,Zmin to get their minimum using MPI_MAX
 sendbuffer[2] *= -1; // to avoid calling MPI_Allreduce again with MPI_MIN
 
 int mpi_err = MPI_Allreduce( sendbuffer, rcvbuffer, 6, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
 if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
 
 rcvbuffer[0] *= -1;
 rcvbuffer[1] *= -1; // negate Xmin,Ymin,Zmin again to get original values
 rcvbuffer[2] *= -1;
 
 globalBox.update_max( &rcvbuffer[3] ); // saving values in globalBox
 globalBox.update_min( &rcvbuffer[0] );
 
 // compute the alternate decomposition; use ScdInterface::compute_partition_sqijk for this
 ScdParData spd;
 spd.partMethod = ScdParData::SQIJK;
 spd.gPeriodic[0] = spd.gPeriodic[1] = spd.gPeriodic[2] = 0;
 double lg = log10( ( localBox.bMax - localBox.bMin ).length() );
 double mfactor = pow( 10.0, 6 - lg );
 int ldims[6], lper[3];
 double dgijk[6];
 localBox.get( dgijk );
 for( int i = 0; i < 6; i++ )
 spd.gDims[i] = dgijk[i] * mfactor;
 ErrorCode rval = ScdInterface::compute_partition( pc->size(), pc->rank(), spd, ldims, lper, regNums );
 if( MB_SUCCESS != rval ) return rval;
 // all we're really interested in is regNums[i], #procs in each direction
 
 for( int i = 0; i < 3; i++ )
 regDeltaXYZ[i] = ( globalBox.bMax[i] - globalBox.bMin[i] ) / double( regNums[i] ); // size of each region
 
 return MB_SUCCESS;
}
 
// this function sets up the TupleList TLreg_o containing the registration messages
// and sends it
ErrorCode SpatialLocator::register_src_with_intermediate_procs( ParallelComm* pc,
 double abs_iter_tol,
 TupleList& TLreg_o )
{
 int corner_ijk[6];
 
 // step 3: compute ijks of local box corners in intermediate partition
 // get corner ijk values for my box
 ErrorCode rval = get_point_ijk( localBox.bMin - CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk );
 if( MB_SUCCESS != rval ) return rval;
 rval = get_point_ijk( localBox.bMax + CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk + 3 );
 if( MB_SUCCESS != rval ) return rval;
 
 // step 4
 // set up TLreg_o
 TLreg_o.initialize( 1, 0, 0, 6, 0 );
 // TLreg_o (int destProc, real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax)
 
 int dest;
 double boxtosend[6];
 
 localBox.get( boxtosend );
 
 // iterate over all regions overlapping with my bounding box using the computerd corner IDs
 for( int k = corner_ijk[2]; k <= corner_ijk[5]; k++ )
 {
 for( int j = corner_ijk[1]; j <= corner_ijk[4]; j++ )
 {
 for( int i = corner_ijk[0]; i <= corner_ijk[3]; i++ )
 {
 dest = k * regNums[0] * regNums[1] + j * regNums[0] + i;
 TLreg_o.push_back( &dest, NULL, NULL, boxtosend );
 }
 }
 }
 
 // step 5
 // send TLreg_o, receive TLrequests_i
 if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLreg_o, 0 );
 
 // step 6
 // Read registration requests from TLreg_o and add to list of procs to forward to
 // get number of tuples sent to me
 
 // read tuples and fill processor list;
 int NN = TLreg_o.get_n();
 for( int i = 0; i < NN; i++ )
 // TLreg_o is now TLrequests_i
 srcProcBoxes[TLreg_o.vi_rd[i]] = BoundBox( TLreg_o.vr_rd + 6 * i );
 
 return MB_SUCCESS;
}
 
ErrorCode SpatialLocator::par_locate_points( ParallelComm* /*pc*/,
 Range& /*vertices*/,
 const double /*rel_iter_tol*/,
 const double /*abs_iter_tol*/,
 const double /*inside_tol*/ )
{
 return MB_UNSUPPORTED_OPERATION;
}
 
bool is_neg( int i )
{
 return ( i == -1 );
}
 
ErrorCode SpatialLocator::par_locate_points( ParallelComm* pc,
 const double* pos,
 int num_points,
 const double rel_iter_tol,
 const double abs_iter_tol,
 const double inside_tol )
{
 ErrorCode rval;
 // TUpleList used for communication
 TupleList TLreg_o; // TLregister_outbound
 TupleList TLquery_o; // TLquery_outbound
 TupleList TLforward_o; // TLforward_outbound
 TupleList TLsearch_results_o; // TLsearch_results_outbound
 
 // initialize timer
 myTimer.time_elapsed();
 timerInitialized = true;
 
 // steps 1-2 - initialize the alternative decomposition box from global box
 rval = initialize_intermediate_partition( pc );
 if( rval != MB_SUCCESS ) return rval;
 
 // steps 3-6 - set up TLreg_o, gs_transfer, gather registrations
 rval = register_src_with_intermediate_procs( pc, abs_iter_tol, TLreg_o );
 if( rval != MB_SUCCESS ) return rval;
 
 myTimes.slTimes[SpatialLocatorTimes::INTMED_INIT] = myTimer.time_elapsed();
 
 // actual parallel point location using intermediate partition
 
 // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
 // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
 // proc point location, ready to send
 // back to tgt procs; src_index of -1 indicates point not located (arguably not
 // useful...)
 
 unsigned int my_rank = ( pc ? pc->proc_config().proc_rank() : 0 );
 
 // TLquery_o: Tuples sent to forwarder proc
 // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
 
 // TLforw_req_i: Tuples to forward to corresponding procs (forwarding requests)
 // TL (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
 
 TLquery_o.initialize( 3, 0, 0, 3, 0 );
 
 int iargs[3];
 
 for( int pnt = 0; pnt < 3 * num_points; pnt += 3 )
 {
 int forw_id =
 proc_from_point( pos + pnt, abs_iter_tol ); // get ID of proc resonsible of the region the proc is in
 
 iargs[0] = forw_id; // toProc
 iargs[1] = my_rank; // originalSourceProc
 iargs[2] = pnt / 3; // targetIndex
 
 TLquery_o.push_back( iargs, NULL, NULL, const_cast< double* >( pos + pnt ) );
 }
 
 // send point search queries to forwarders
 if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLquery_o, 0 );
 
 myTimes.slTimes[SpatialLocatorTimes::INTMED_SEND] = myTimer.time_elapsed();
 
 // now read forwarding requests and forward to corresponding procs
 // TLquery_o is now TLforw_req_i
 
 // TLforward_o: query messages forwarded to corresponding procs
 // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
 
 TLforward_o.initialize( 3, 0, 0, 3, 0 );
 
 int NN = TLquery_o.get_n();
 
 for( int i = 0; i < NN; i++ )
 {
 iargs[1] = TLquery_o.vi_rd[3 * i + 1]; // get OriginalSourceProc
 iargs[2] = TLquery_o.vi_rd[3 * i + 2]; // targetIndex
 CartVect tmp_pnt( TLquery_o.vr_rd + 3 * i );
 
 // compare coordinates to list of bounding boxes
 for( std::map< int, BoundBox >::iterator mit = srcProcBoxes.begin(); mit != srcProcBoxes.end(); ++mit )
 {
 if( ( *mit ).second.contains_point( tmp_pnt.array(), abs_iter_tol ) )
 {
 iargs[0] = ( *mit ).first;
 TLforward_o.push_back( iargs, NULL, NULL, tmp_pnt.array() );
 }
 }
 }
 
 myTimes.slTimes[SpatialLocatorTimes::INTMED_SEARCH] = myTimer.time_elapsed();
 
 if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLforward_o, 0 );
 
 myTimes.slTimes[SpatialLocatorTimes::SRC_SEND] = myTimer.time_elapsed();
 
 // cache time here, because locate_points also calls elapsed functions and we want to account
 // for tuple list initialization here
 double tstart = myTimer.time_since_birth();
 
 // step 12
 // now read Point Search requests
 // TLforward_o is now TLsearch_req_i
 // TLsearch_req_i: (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
 
 NN = TLforward_o.get_n();
 
 // TLsearch_results_o
 // TL: (OriginalSourceProc, targetIndex, sourceIndex, U,V,W);
 TLsearch_results_o.initialize( 3, 0, 0, 0, 0 );
 
 // step 13 is done in test_local_box
 
 std::vector< double > params( 3 * NN );
 std::vector< int > is_inside( NN, 0 );
 std::vector< EntityHandle > ents( NN, 0 );
 
 rval = locate_points( TLforward_o.vr_rd, TLforward_o.get_n(), &ents[0], &params[0], &is_inside[0], rel_iter_tol,
 abs_iter_tol, inside_tol );
 if( MB_SUCCESS != rval ) return rval;
 
 locTable.initialize( 1, 0, 1, 3, 0 );
 locTable.enableWriteAccess();
 for( int i = 0; i < NN; i++ )
 {
 if( is_inside[i] )
 {
 iargs[0] = TLforward_o.vi_rd[3 * i + 1];
 iargs[1] = TLforward_o.vi_rd[3 * i + 2];
 iargs[2] = locTable.get_n();
 TLsearch_results_o.push_back( iargs, NULL, NULL, NULL );
 Ulong ent_ulong = (Ulong)ents[i];
 sint forward = (sint)TLforward_o.vi_rd[3 * i + 1];
 locTable.push_back( &forward, NULL, &ent_ulong, &params[3 * i] );
 }
 }
 locTable.disableWriteAccess();
 
 myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_since_birth() - tstart;
 myTimer.time_elapsed(); // call this to reset last time called
 
 // step 14: send TLsearch_results_o and receive TLloc_i
 if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLsearch_results_o, 0 );
 
 myTimes.slTimes[SpatialLocatorTimes::TARG_RETURN] = myTimer.time_elapsed();
 
 // store proc/index tuples in parLocTable
 parLocTable.initialize( 2, 0, 0, 0, num_points );
 parLocTable.enableWriteAccess();
 std::fill( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, -1 );
 
 for( unsigned int i = 0; i < TLsearch_results_o.get_n(); i++ )
 {
 int idx = TLsearch_results_o.vi_rd[3 * i + 1];
 parLocTable.vi_wr[2 * idx] = TLsearch_results_o.vi_rd[3 * i];
 parLocTable.vi_wr[2 * idx + 1] = TLsearch_results_o.vi_rd[3 * i + 2];
 }
 
 if( debug )
 {
 int num_found =
 num_points - 0.5 * std::count_if( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, is_neg );
 std::cout << "Points found = " << num_found << "/" << num_points << " ("
 << 100.0 * ( (double)num_found / num_points ) << "%)" << std::endl;
 }
 
 myTimes.slTimes[SpatialLocatorTimes::TARG_STORE] = myTimer.time_elapsed();
 
 return MB_SUCCESS;
}
 
#endif
 
ErrorCode SpatialLocator::locate_points( Range& verts,
 const double rel_iter_tol,
 const double abs_iter_tol,
 const double inside_tol )
{
 bool i_initialized = false;
 if( !timerInitialized )
 {
 myTimer.time_elapsed();
 timerInitialized = true;
 i_initialized = true;
 }
 
 assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
 std::vector< double > pos( 3 * verts.size() );
 ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
 if( MB_SUCCESS != rval ) return rval;
 rval = locate_points( &pos[0], verts.size(), rel_iter_tol, abs_iter_tol, inside_tol );
 if( MB_SUCCESS != rval ) return rval;
 
 // only call this if I'm the top-level function, since it resets the last time called
 if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
 
 return MB_SUCCESS;
}
 
ErrorCode SpatialLocator::locate_points( const double* pos,
 int num_points,
 const double rel_iter_tol,
 const double abs_iter_tol,
 const double inside_tol )
{
 bool i_initialized = false;
 if( !timerInitialized )
 {
 myTimer.time_elapsed();
 timerInitialized = true;
 i_initialized = true;
 }
 // initialize to tuple structure (p_ui, hs_ul, r[3]_d) (see header comments for locTable)
 locTable.initialize( 1, 0, 1, 3, num_points );
 locTable.enableWriteAccess();
 
 // pass storage directly into locate_points, since we know those arrays are contiguous
 ErrorCode rval = locate_points( pos, num_points, (EntityHandle*)locTable.vul_wr, locTable.vr_wr, NULL, rel_iter_tol,
 abs_iter_tol, inside_tol );
 std::fill( locTable.vi_wr, locTable.vi_wr + num_points, 0 );
 locTable.set_n( num_points );
 if( MB_SUCCESS != rval ) return rval;
 
 // only call this if I'm the top-level function, since it resets the last time called
 if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
 
 return MB_SUCCESS;
}
 
ErrorCode SpatialLocator::locate_points( Range& verts,
 EntityHandle* ents,
 double* params,
 int* is_inside,
 const double rel_iter_tol,
 const double abs_iter_tol,
 const double inside_tol )
{
 bool i_initialized = false;
 if( !timerInitialized )
 {
 myTimer.time_elapsed();
 timerInitialized = true;
 i_initialized = true;
 }
 
 assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
 std::vector< double > pos( 3 * verts.size() );
 ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
 if( MB_SUCCESS != rval ) return rval;
 rval = locate_points( &pos[0], verts.size(), ents, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
 
 // only call this if I'm the top-level function, since it resets the last time called
 if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
 
 return rval;
}
 
ErrorCode SpatialLocator::locate_points( const double* pos,
 int num_points,
 EntityHandle* ents,
 double* params,
 int* is_inside,
 const double /* rel_iter_tol */,
 const double abs_iter_tol,
 const double inside_tol )
{
 bool i_initialized = false;
 if( !timerInitialized )
 {
 myTimer.time_elapsed();
 timerInitialized = true;
 i_initialized = true;
 }
 
 /*
 double tmp_abs_iter_tol = abs_iter_tol;
 if (rel_iter_tol && !tmp_abs_iter_tol) {
 // relative epsilon given, translate to absolute epsilon using box dimensions
 tmp_abs_iter_tol = rel_iter_tol * localBox.diagonal_length();
 }
 */
 
 if( elemEval && myTree->get_eval() != elemEval ) myTree->set_eval( elemEval );
 
 ErrorCode rval = MB_SUCCESS;
 for( int i = 0; i < num_points; i++ )
 {
 int i3 = 3 * i;
 ErrorCode tmp_rval =
 myTree->point_search( pos + i3, ents[i], abs_iter_tol, inside_tol, NULL, NULL, (CartVect*)( params + i3 ) );
 if( MB_SUCCESS != tmp_rval )
 {
 rval = tmp_rval;
 continue;
 }
 
 if( debug && !ents[i] )
 {
 std::cout << "Point " << i << " not found; point: (" << pos[i3] << "," << pos[i3 + 1] << "," << pos[i3 + 2]
 << ")" << std::endl;
 }
 
 if( is_inside ) is_inside[i] = ( ents[i] ? true : false );
 }
 
 // only call this if I'm the top-level function, since it resets the last time called
 if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
 
 return rval;
}
 
/* Count the number of located points in locTable
 * Return the number of entries in locTable that have non-zero entity handles, which
 * represents the number of points in targetEnts that were inside one element in sourceEnts
 *
 */
int SpatialLocator::local_num_located()
{
 int num_located = locTable.get_n() - std::count( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
 if( num_located != (int)locTable.get_n() )
 {
 Ulong* nl = std::find( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
 if( nl )
 {
 int idx = nl - locTable.vul_rd;
 if( idx )
 {
 }
 }
 }
 return num_located;
}
 
/* Count the number of located points in parLocTable
 * Return the number of entries in parLocTable that have a non-negative index in on a remote
 * proc in parLocTable, which gives the number of points located in at least one element in a
 * remote proc's sourceEnts.
 */
int SpatialLocator::remote_num_located()
{
 int located = 0;
 for( unsigned int i = 0; i < parLocTable.get_n(); i++ )
 if( parLocTable.vi_rd[2 * i] != -1 ) located++;
 return located;
}
} // namespace moab