MOAB: Mesh Oriented datABase  (version 5.2.1)
SpatialLocator.cpp
Go to the documentation of this file.
00001 #include "moab/SpatialLocator.hpp"
00002 #include "moab/Interface.hpp"
00003 #include "moab/ElemEvaluator.hpp"
00004 #include "moab/AdaptiveKDTree.hpp"
00005 #include "moab/BVHTree.hpp"
00006 
00007 // include ScdInterface for box partitioning
00008 #include "moab/ScdInterface.hpp"
00009 
00010 #ifdef MOAB_HAVE_MPI
00011 #include "moab/ParallelComm.hpp"
00012 #endif
00013 
00014 namespace moab
00015 {
00016 static bool debug = false;
00017 
00018 SpatialLocator::SpatialLocator( Interface* impl, Range& elems, Tree* tree, ElemEvaluator* eval )
00019     : mbImpl( impl ), myElems( elems ), myDim( -1 ), myTree( tree ), elemEval( eval ), iCreatedTree( false ),
00020       timerInitialized( false )
00021 {
00022     create_tree();
00023 
00024     if( !elems.empty() )
00025     {
00026         myDim          = mbImpl->dimension_from_handle( *elems.rbegin() );
00027         ErrorCode rval = myTree->build_tree( myElems );
00028         if( MB_SUCCESS != rval ) throw rval;
00029 
00030         rval = myTree->get_bounding_box( localBox );
00031         if( MB_SUCCESS != rval ) throw rval;
00032     }
00033 
00034     regNums[0] = regNums[1] = regNums[2] = 0;
00035 }
00036 
00037 void SpatialLocator::create_tree()
00038 {
00039     if( myTree ) return;
00040 
00041     if( myElems.empty() || mbImpl->type_from_handle( *myElems.rbegin() ) == MBVERTEX )
00042         // create a kdtree if only vertices
00043         myTree = new AdaptiveKDTree( mbImpl );
00044     else
00045         // otherwise a BVHtree, since it performs better for elements
00046         myTree = new BVHTree( mbImpl );
00047 
00048     iCreatedTree = true;
00049 }
00050 
00051 ErrorCode SpatialLocator::add_elems( Range& elems )
00052 {
00053     if( elems.empty() ||
00054         mbImpl->dimension_from_handle( *elems.begin() ) != mbImpl->dimension_from_handle( *elems.rbegin() ) )
00055         return MB_FAILURE;
00056 
00057     myDim   = mbImpl->dimension_from_handle( *elems.begin() );
00058     myElems = elems;
00059 
00060     ErrorCode rval = myTree->build_tree( myElems );
00061     return rval;
00062 }
00063 
00064 #ifdef MOAB_HAVE_MPI
00065 ErrorCode SpatialLocator::initialize_intermediate_partition( ParallelComm* pc )
00066 {
00067     if( !pc ) return MB_FAILURE;
00068 
00069     BoundBox gbox;
00070 
00071     // step 2
00072     // get the global bounding box
00073     double sendbuffer[6];
00074     double rcvbuffer[6];
00075 
00076     localBox.get( sendbuffer );  // fill sendbuffer with local box, max values in [0:2] min values in [3:5]
00077     sendbuffer[0] *= -1;
00078     sendbuffer[1] *= -1;  // negate Xmin,Ymin,Zmin to get their minimum using MPI_MAX
00079     sendbuffer[2] *= -1;  // to avoid calling MPI_Allreduce again with MPI_MIN
00080 
00081     int mpi_err = MPI_Allreduce( sendbuffer, rcvbuffer, 6, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
00082     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00083 
00084     rcvbuffer[0] *= -1;
00085     rcvbuffer[1] *= -1;  // negate Xmin,Ymin,Zmin again to get original values
00086     rcvbuffer[2] *= -1;
00087 
00088     globalBox.update_max( &rcvbuffer[3] );  // saving values in globalBox
00089     globalBox.update_min( &rcvbuffer[0] );
00090 
00091     // compute the alternate decomposition; use ScdInterface::compute_partition_sqijk for this
00092     ScdParData spd;
00093     spd.partMethod   = ScdParData::SQIJK;
00094     spd.gPeriodic[0] = spd.gPeriodic[1] = spd.gPeriodic[2] = 0;
00095     double lg                                              = log10( ( localBox.bMax - localBox.bMin ).length() );
00096     double mfactor                                         = pow( 10.0, 6 - lg );
00097     int ldims[6], lper[3];
00098     double dgijk[6];
00099     localBox.get( dgijk );
00100     for( int i = 0; i < 6; i++ )
00101         spd.gDims[i] = dgijk[i] * mfactor;
00102     ErrorCode rval = ScdInterface::compute_partition( pc->size(), pc->rank(), spd, ldims, lper, regNums );
00103     if( MB_SUCCESS != rval ) return rval;
00104     // all we're really interested in is regNums[i], #procs in each direction
00105 
00106     for( int i = 0; i < 3; i++ )
00107         regDeltaXYZ[i] = ( globalBox.bMax[i] - globalBox.bMin[i] ) / double( regNums[i] );  // size of each region
00108 
00109     return MB_SUCCESS;
00110 }
00111 
00112 // this function sets up the TupleList TLreg_o containing the registration messages
00113 // and sends it
00114 ErrorCode SpatialLocator::register_src_with_intermediate_procs( ParallelComm* pc, double abs_iter_tol,
00115                                                                 TupleList& TLreg_o )
00116 {
00117     int corner_ijk[6];
00118 
00119     // step 3: compute ijks of local box corners in intermediate partition
00120     // get corner ijk values for my box
00121     ErrorCode rval = get_point_ijk( localBox.bMin - CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk );
00122     if( MB_SUCCESS != rval ) return rval;
00123     rval = get_point_ijk( localBox.bMax + CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk + 3 );
00124     if( MB_SUCCESS != rval ) return rval;
00125 
00126     // step 4
00127     // set up TLreg_o
00128     TLreg_o.initialize( 1, 0, 0, 6, 0 );
00129     // TLreg_o (int destProc, real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax)
00130 
00131     int dest;
00132     double boxtosend[6];
00133 
00134     localBox.get( boxtosend );
00135 
00136     // iterate over all regions overlapping with my bounding box using the computerd corner IDs
00137     for( int k = corner_ijk[2]; k <= corner_ijk[5]; k++ )
00138     {
00139         for( int j = corner_ijk[1]; j <= corner_ijk[4]; j++ )
00140         {
00141             for( int i = corner_ijk[0]; i <= corner_ijk[3]; i++ )
00142             {
00143                 dest = k * regNums[0] * regNums[1] + j * regNums[0] + i;
00144                 TLreg_o.push_back( &dest, NULL, NULL, boxtosend );
00145             }
00146         }
00147     }
00148 
00149     // step 5
00150     // send TLreg_o, receive TLrequests_i
00151     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLreg_o, 0 );
00152 
00153     // step 6
00154     // Read registration requests from TLreg_o and add to list of procs to forward to
00155     // get number of tuples sent to me
00156 
00157     // read tuples and fill processor list;
00158     int NN = TLreg_o.get_n();
00159     for( int i = 0; i < NN; i++ )
00160         // TLreg_o is now TLrequests_i
00161         srcProcBoxes[TLreg_o.vi_rd[i]] = BoundBox( TLreg_o.vr_rd + 6 * i );
00162 
00163     return MB_SUCCESS;
00164 }
00165 
00166 ErrorCode SpatialLocator::par_locate_points( ParallelComm* /*pc*/, Range& /*vertices*/, const double /*rel_iter_tol*/,
00167                                              const double /*abs_iter_tol*/, const double /*inside_tol*/ )
00168 {
00169     return MB_UNSUPPORTED_OPERATION;
00170 }
00171 
00172 bool is_neg( int i )
00173 {
00174     return ( i == -1 );
00175 }
00176 
00177 ErrorCode SpatialLocator::par_locate_points( ParallelComm* pc, const double* pos, int num_points,
00178                                              const double rel_iter_tol, const double abs_iter_tol,
00179                                              const double inside_tol )
00180 {
00181     ErrorCode rval;
00182     // TUpleList used for communication
00183     TupleList TLreg_o;             // TLregister_outbound
00184     TupleList TLquery_o;           // TLquery_outbound
00185     TupleList TLforward_o;         // TLforward_outbound
00186     TupleList TLsearch_results_o;  // TLsearch_results_outbound
00187 
00188     // initialize timer
00189     myTimer.time_elapsed();
00190     timerInitialized = true;
00191 
00192     // steps 1-2 - initialize the alternative decomposition box from global box
00193     rval = initialize_intermediate_partition( pc );
00194     if( rval != MB_SUCCESS ) return rval;
00195 
00196     // steps 3-6 - set up TLreg_o, gs_transfer, gather registrations
00197     rval = register_src_with_intermediate_procs( pc, abs_iter_tol, TLreg_o );
00198     if( rval != MB_SUCCESS ) return rval;
00199 
00200     myTimes.slTimes[SpatialLocatorTimes::INTMED_INIT] = myTimer.time_elapsed();
00201 
00202     // actual parallel point location using intermediate partition
00203 
00204     // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
00205     // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
00206     // proc point location, ready to send
00207     //             back to tgt procs; src_index of -1 indicates point not located (arguably not
00208     //             useful...)
00209 
00210     unsigned int my_rank = ( pc ? pc->proc_config().proc_rank() : 0 );
00211 
00212     // TLquery_o: Tuples sent to forwarder proc
00213     // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
00214 
00215     // TLforw_req_i: Tuples to forward to corresponding procs (forwarding requests)
00216     // TL (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
00217 
00218     TLquery_o.initialize( 3, 0, 0, 3, 0 );
00219 
00220     int iargs[3];
00221 
00222     for( int pnt = 0; pnt < 3 * num_points; pnt += 3 )
00223     {
00224         int forw_id =
00225             proc_from_point( pos + pnt, abs_iter_tol );  // get ID of proc resonsible of the region the proc is in
00226 
00227         iargs[0] = forw_id;  // toProc
00228         iargs[1] = my_rank;  // originalSourceProc
00229         iargs[2] = pnt / 3;  // targetIndex
00230 
00231         TLquery_o.push_back( iargs, NULL, NULL, const_cast< double* >( pos + pnt ) );
00232     }
00233 
00234     // send point search queries to forwarders
00235     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLquery_o, 0 );
00236 
00237     myTimes.slTimes[SpatialLocatorTimes::INTMED_SEND] = myTimer.time_elapsed();
00238 
00239     // now read forwarding requests and forward to corresponding procs
00240     // TLquery_o is now TLforw_req_i
00241 
00242     // TLforward_o: query messages forwarded to corresponding procs
00243     // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
00244 
00245     TLforward_o.initialize( 3, 0, 0, 3, 0 );
00246 
00247     int NN = TLquery_o.get_n();
00248 
00249     for( int i = 0; i < NN; i++ )
00250     {
00251         iargs[1] = TLquery_o.vi_rd[3 * i + 1];  // get OriginalSourceProc
00252         iargs[2] = TLquery_o.vi_rd[3 * i + 2];  // targetIndex
00253         CartVect tmp_pnt( TLquery_o.vr_rd + 3 * i );
00254 
00255         // compare coordinates to list of bounding boxes
00256         for( std::map< int, BoundBox >::iterator mit = srcProcBoxes.begin(); mit != srcProcBoxes.end(); ++mit )
00257         {
00258             if( ( *mit ).second.contains_point( tmp_pnt.array(), abs_iter_tol ) )
00259             {
00260                 iargs[0] = ( *mit ).first;
00261                 TLforward_o.push_back( iargs, NULL, NULL, tmp_pnt.array() );
00262             }
00263         }
00264     }
00265 
00266     myTimes.slTimes[SpatialLocatorTimes::INTMED_SEARCH] = myTimer.time_elapsed();
00267 
00268     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLforward_o, 0 );
00269 
00270     myTimes.slTimes[SpatialLocatorTimes::SRC_SEND] = myTimer.time_elapsed();
00271 
00272     // cache time here, because locate_points also calls elapsed functions and we want to account
00273     // for tuple list initialization here
00274     double tstart = myTimer.time_since_birth();
00275 
00276     // step 12
00277     // now read Point Search requests
00278     // TLforward_o is now TLsearch_req_i
00279     // TLsearch_req_i: (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
00280 
00281     NN = TLforward_o.get_n();
00282 
00283     // TLsearch_results_o
00284     // TL: (OriginalSourceProc, targetIndex, sourceIndex, U,V,W);
00285     TLsearch_results_o.initialize( 3, 0, 0, 0, 0 );
00286 
00287     // step 13 is done in test_local_box
00288 
00289     std::vector< double > params( 3 * NN );
00290     std::vector< int > is_inside( NN, 0 );
00291     std::vector< EntityHandle > ents( NN, 0 );
00292 
00293     rval = locate_points( TLforward_o.vr_rd, TLforward_o.get_n(), &ents[0], &params[0], &is_inside[0], rel_iter_tol,
00294                           abs_iter_tol, inside_tol );
00295     if( MB_SUCCESS != rval ) return rval;
00296 
00297     locTable.initialize( 1, 0, 1, 3, 0 );
00298     locTable.enableWriteAccess();
00299     for( int i = 0; i < NN; i++ )
00300     {
00301         if( is_inside[i] )
00302         {
00303             iargs[0] = TLforward_o.vi_rd[3 * i + 1];
00304             iargs[1] = TLforward_o.vi_rd[3 * i + 2];
00305             iargs[2] = locTable.get_n();
00306             TLsearch_results_o.push_back( iargs, NULL, NULL, NULL );
00307             Ulong ent_ulong = (Ulong)ents[i];
00308             sint forward    = (sint)TLforward_o.vi_rd[3 * i + 1];
00309             locTable.push_back( &forward, NULL, &ent_ulong, &params[3 * i] );
00310         }
00311     }
00312     locTable.disableWriteAccess();
00313 
00314     myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_since_birth() - tstart;
00315     myTimer.time_elapsed();  // call this to reset last time called
00316 
00317     // step 14: send TLsearch_results_o and receive TLloc_i
00318     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLsearch_results_o, 0 );
00319 
00320     myTimes.slTimes[SpatialLocatorTimes::TARG_RETURN] = myTimer.time_elapsed();
00321 
00322     // store proc/index tuples in parLocTable
00323     parLocTable.initialize( 2, 0, 0, 0, num_points );
00324     parLocTable.enableWriteAccess();
00325     std::fill( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, -1 );
00326 
00327     for( unsigned int i = 0; i < TLsearch_results_o.get_n(); i++ )
00328     {
00329         int idx                        = TLsearch_results_o.vi_rd[3 * i + 1];
00330         parLocTable.vi_wr[2 * idx]     = TLsearch_results_o.vi_rd[3 * i];
00331         parLocTable.vi_wr[2 * idx + 1] = TLsearch_results_o.vi_rd[3 * i + 2];
00332     }
00333 
00334     if( debug )
00335     {
00336         int num_found =
00337             num_points - 0.5 * std::count_if( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, is_neg );
00338         std::cout << "Points found = " << num_found << "/" << num_points << " ("
00339                   << 100.0 * ( (double)num_found / num_points ) << "%)" << std::endl;
00340     }
00341 
00342     myTimes.slTimes[SpatialLocatorTimes::TARG_STORE] = myTimer.time_elapsed();
00343 
00344     return MB_SUCCESS;
00345 }
00346 
00347 #endif
00348 
00349 ErrorCode SpatialLocator::locate_points( Range& verts, const double rel_iter_tol, const double abs_iter_tol,
00350                                          const double inside_tol )
00351 {
00352     bool i_initialized = false;
00353     if( !timerInitialized )
00354     {
00355         myTimer.time_elapsed();
00356         timerInitialized = true;
00357         i_initialized    = true;
00358     }
00359 
00360     assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
00361     std::vector< double > pos( 3 * verts.size() );
00362     ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
00363     if( MB_SUCCESS != rval ) return rval;
00364     rval = locate_points( &pos[0], verts.size(), rel_iter_tol, abs_iter_tol, inside_tol );
00365     if( MB_SUCCESS != rval ) return rval;
00366 
00367     // only call this if I'm the top-level function, since it resets the last time called
00368     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00369 
00370     return MB_SUCCESS;
00371 }
00372 
00373 ErrorCode SpatialLocator::locate_points( const double* pos, int num_points, const double rel_iter_tol,
00374                                          const double abs_iter_tol, const double inside_tol )
00375 {
00376     bool i_initialized = false;
00377     if( !timerInitialized )
00378     {
00379         myTimer.time_elapsed();
00380         timerInitialized = true;
00381         i_initialized    = true;
00382     }
00383     // initialize to tuple structure (p_ui, hs_ul, r[3]_d) (see header comments for locTable)
00384     locTable.initialize( 1, 0, 1, 3, num_points );
00385     locTable.enableWriteAccess();
00386 
00387     // pass storage directly into locate_points, since we know those arrays are contiguous
00388     ErrorCode rval = locate_points( pos, num_points, (EntityHandle*)locTable.vul_wr, locTable.vr_wr, NULL, rel_iter_tol,
00389                                     abs_iter_tol, inside_tol );
00390     std::fill( locTable.vi_wr, locTable.vi_wr + num_points, 0 );
00391     locTable.set_n( num_points );
00392     if( MB_SUCCESS != rval ) return rval;
00393 
00394     // only call this if I'm the top-level function, since it resets the last time called
00395     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00396 
00397     return MB_SUCCESS;
00398 }
00399 
00400 ErrorCode SpatialLocator::locate_points( Range& verts, EntityHandle* ents, double* params, int* is_inside,
00401                                          const double rel_iter_tol, const double abs_iter_tol, const double inside_tol )
00402 {
00403     bool i_initialized = false;
00404     if( !timerInitialized )
00405     {
00406         myTimer.time_elapsed();
00407         timerInitialized = true;
00408         i_initialized    = true;
00409     }
00410 
00411     assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
00412     std::vector< double > pos( 3 * verts.size() );
00413     ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
00414     if( MB_SUCCESS != rval ) return rval;
00415     rval = locate_points( &pos[0], verts.size(), ents, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
00416 
00417     // only call this if I'm the top-level function, since it resets the last time called
00418     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00419 
00420     return rval;
00421 }
00422 
00423 ErrorCode SpatialLocator::locate_points( const double* pos, int num_points, EntityHandle* ents, double* params,
00424                                          int* is_inside, const double /* rel_iter_tol */, const double abs_iter_tol,
00425                                          const double inside_tol )
00426 {
00427     bool i_initialized = false;
00428     if( !timerInitialized )
00429     {
00430         myTimer.time_elapsed();
00431         timerInitialized = true;
00432         i_initialized    = true;
00433     }
00434 
00435     /*
00436     double tmp_abs_iter_tol = abs_iter_tol;
00437     if (rel_iter_tol && !tmp_abs_iter_tol) {
00438         // relative epsilon given, translate to absolute epsilon using box dimensions
00439       tmp_abs_iter_tol = rel_iter_tol * localBox.diagonal_length();
00440     }
00441     */
00442 
00443     if( elemEval && myTree->get_eval() != elemEval ) myTree->set_eval( elemEval );
00444 
00445     ErrorCode rval = MB_SUCCESS;
00446     for( int i = 0; i < num_points; i++ )
00447     {
00448         int i3 = 3 * i;
00449         ErrorCode tmp_rval =
00450             myTree->point_search( pos + i3, ents[i], abs_iter_tol, inside_tol, NULL, NULL, (CartVect*)( params + i3 ) );
00451         if( MB_SUCCESS != tmp_rval )
00452         {
00453             rval = tmp_rval;
00454             continue;
00455         }
00456 
00457         if( debug && !ents[i] )
00458         {
00459             std::cout << "Point " << i << " not found; point: (" << pos[i3] << "," << pos[i3 + 1] << "," << pos[i3 + 2]
00460                       << ")" << std::endl;
00461         }
00462 
00463         if( is_inside ) is_inside[i] = ( ents[i] ? true : false );
00464     }
00465 
00466     // only call this if I'm the top-level function, since it resets the last time called
00467     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00468 
00469     return rval;
00470 }
00471 
00472 /* Count the number of located points in locTable
00473  * Return the number of entries in locTable that have non-zero entity handles, which
00474  * represents the number of points in targetEnts that were inside one element in sourceEnts
00475  *
00476  */
00477 int SpatialLocator::local_num_located()
00478 {
00479     int num_located = locTable.get_n() - std::count( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
00480     if( num_located != (int)locTable.get_n() )
00481     {
00482         Ulong* nl = std::find( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
00483         if( nl )
00484         {
00485             int idx = nl - locTable.vul_rd;
00486             if( idx ) {}
00487         }
00488     }
00489     return num_located;
00490 }
00491 
00492 /* Count the number of located points in parLocTable
00493  * Return the number of entries in parLocTable that have a non-negative index in on a remote
00494  * proc in parLocTable, which gives the number of points located in at least one element in a
00495  * remote proc's sourceEnts.
00496  */
00497 int SpatialLocator::remote_num_located()
00498 {
00499     int located = 0;
00500     for( unsigned int i = 0; i < parLocTable.get_n(); i++ )
00501         if( parLocTable.vi_rd[2 * i] != -1 ) located++;
00502     return located;
00503 }
00504 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines