MOAB: Mesh Oriented datABase  (version 5.4.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,
00115                                                                 double abs_iter_tol,
00116                                                                 TupleList& TLreg_o )
00117 {
00118     int corner_ijk[6];
00119 
00120     // step 3: compute ijks of local box corners in intermediate partition
00121     // get corner ijk values for my box
00122     ErrorCode rval = get_point_ijk( localBox.bMin - CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk );
00123     if( MB_SUCCESS != rval ) return rval;
00124     rval = get_point_ijk( localBox.bMax + CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk + 3 );
00125     if( MB_SUCCESS != rval ) return rval;
00126 
00127     // step 4
00128     // set up TLreg_o
00129     TLreg_o.initialize( 1, 0, 0, 6, 0 );
00130     // TLreg_o (int destProc, real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax)
00131 
00132     int dest;
00133     double boxtosend[6];
00134 
00135     localBox.get( boxtosend );
00136 
00137     // iterate over all regions overlapping with my bounding box using the computerd corner IDs
00138     for( int k = corner_ijk[2]; k <= corner_ijk[5]; k++ )
00139     {
00140         for( int j = corner_ijk[1]; j <= corner_ijk[4]; j++ )
00141         {
00142             for( int i = corner_ijk[0]; i <= corner_ijk[3]; i++ )
00143             {
00144                 dest = k * regNums[0] * regNums[1] + j * regNums[0] + i;
00145                 TLreg_o.push_back( &dest, NULL, NULL, boxtosend );
00146             }
00147         }
00148     }
00149 
00150     // step 5
00151     // send TLreg_o, receive TLrequests_i
00152     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLreg_o, 0 );
00153 
00154     // step 6
00155     // Read registration requests from TLreg_o and add to list of procs to forward to
00156     // get number of tuples sent to me
00157 
00158     // read tuples and fill processor list;
00159     int NN = TLreg_o.get_n();
00160     for( int i = 0; i < NN; i++ )
00161         // TLreg_o is now TLrequests_i
00162         srcProcBoxes[TLreg_o.vi_rd[i]] = BoundBox( TLreg_o.vr_rd + 6 * i );
00163 
00164     return MB_SUCCESS;
00165 }
00166 
00167 ErrorCode SpatialLocator::par_locate_points( ParallelComm* /*pc*/,
00168                                              Range& /*vertices*/,
00169                                              const double /*rel_iter_tol*/,
00170                                              const double /*abs_iter_tol*/,
00171                                              const double /*inside_tol*/ )
00172 {
00173     return MB_UNSUPPORTED_OPERATION;
00174 }
00175 
00176 bool is_neg( int i )
00177 {
00178     return ( i == -1 );
00179 }
00180 
00181 ErrorCode SpatialLocator::par_locate_points( ParallelComm* pc,
00182                                              const double* pos,
00183                                              int num_points,
00184                                              const double rel_iter_tol,
00185                                              const double abs_iter_tol,
00186                                              const double inside_tol )
00187 {
00188     ErrorCode rval;
00189     // TUpleList used for communication
00190     TupleList TLreg_o;             // TLregister_outbound
00191     TupleList TLquery_o;           // TLquery_outbound
00192     TupleList TLforward_o;         // TLforward_outbound
00193     TupleList TLsearch_results_o;  // TLsearch_results_outbound
00194 
00195     // initialize timer
00196     myTimer.time_elapsed();
00197     timerInitialized = true;
00198 
00199     // steps 1-2 - initialize the alternative decomposition box from global box
00200     rval = initialize_intermediate_partition( pc );
00201     if( rval != MB_SUCCESS ) return rval;
00202 
00203     // steps 3-6 - set up TLreg_o, gs_transfer, gather registrations
00204     rval = register_src_with_intermediate_procs( pc, abs_iter_tol, TLreg_o );
00205     if( rval != MB_SUCCESS ) return rval;
00206 
00207     myTimes.slTimes[SpatialLocatorTimes::INTMED_INIT] = myTimer.time_elapsed();
00208 
00209     // actual parallel point location using intermediate partition
00210 
00211     // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
00212     // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
00213     // proc point location, ready to send
00214     //             back to tgt procs; src_index of -1 indicates point not located (arguably not
00215     //             useful...)
00216 
00217     unsigned int my_rank = ( pc ? pc->proc_config().proc_rank() : 0 );
00218 
00219     // TLquery_o: Tuples sent to forwarder proc
00220     // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
00221 
00222     // TLforw_req_i: Tuples to forward to corresponding procs (forwarding requests)
00223     // TL (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
00224 
00225     TLquery_o.initialize( 3, 0, 0, 3, 0 );
00226 
00227     int iargs[3];
00228 
00229     for( int pnt = 0; pnt < 3 * num_points; pnt += 3 )
00230     {
00231         int forw_id =
00232             proc_from_point( pos + pnt, abs_iter_tol );  // get ID of proc resonsible of the region the proc is in
00233 
00234         iargs[0] = forw_id;  // toProc
00235         iargs[1] = my_rank;  // originalSourceProc
00236         iargs[2] = pnt / 3;  // targetIndex
00237 
00238         TLquery_o.push_back( iargs, NULL, NULL, const_cast< double* >( pos + pnt ) );
00239     }
00240 
00241     // send point search queries to forwarders
00242     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLquery_o, 0 );
00243 
00244     myTimes.slTimes[SpatialLocatorTimes::INTMED_SEND] = myTimer.time_elapsed();
00245 
00246     // now read forwarding requests and forward to corresponding procs
00247     // TLquery_o is now TLforw_req_i
00248 
00249     // TLforward_o: query messages forwarded to corresponding procs
00250     // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
00251 
00252     TLforward_o.initialize( 3, 0, 0, 3, 0 );
00253 
00254     int NN = TLquery_o.get_n();
00255 
00256     for( int i = 0; i < NN; i++ )
00257     {
00258         iargs[1] = TLquery_o.vi_rd[3 * i + 1];  // get OriginalSourceProc
00259         iargs[2] = TLquery_o.vi_rd[3 * i + 2];  // targetIndex
00260         CartVect tmp_pnt( TLquery_o.vr_rd + 3 * i );
00261 
00262         // compare coordinates to list of bounding boxes
00263         for( std::map< int, BoundBox >::iterator mit = srcProcBoxes.begin(); mit != srcProcBoxes.end(); ++mit )
00264         {
00265             if( ( *mit ).second.contains_point( tmp_pnt.array(), abs_iter_tol ) )
00266             {
00267                 iargs[0] = ( *mit ).first;
00268                 TLforward_o.push_back( iargs, NULL, NULL, tmp_pnt.array() );
00269             }
00270         }
00271     }
00272 
00273     myTimes.slTimes[SpatialLocatorTimes::INTMED_SEARCH] = myTimer.time_elapsed();
00274 
00275     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLforward_o, 0 );
00276 
00277     myTimes.slTimes[SpatialLocatorTimes::SRC_SEND] = myTimer.time_elapsed();
00278 
00279     // cache time here, because locate_points also calls elapsed functions and we want to account
00280     // for tuple list initialization here
00281     double tstart = myTimer.time_since_birth();
00282 
00283     // step 12
00284     // now read Point Search requests
00285     // TLforward_o is now TLsearch_req_i
00286     // TLsearch_req_i: (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
00287 
00288     NN = TLforward_o.get_n();
00289 
00290     // TLsearch_results_o
00291     // TL: (OriginalSourceProc, targetIndex, sourceIndex, U,V,W);
00292     TLsearch_results_o.initialize( 3, 0, 0, 0, 0 );
00293 
00294     // step 13 is done in test_local_box
00295 
00296     std::vector< double > params( 3 * NN );
00297     std::vector< int > is_inside( NN, 0 );
00298     std::vector< EntityHandle > ents( NN, 0 );
00299 
00300     rval = locate_points( TLforward_o.vr_rd, TLforward_o.get_n(), &ents[0], &params[0], &is_inside[0], rel_iter_tol,
00301                           abs_iter_tol, inside_tol );
00302     if( MB_SUCCESS != rval ) return rval;
00303 
00304     locTable.initialize( 1, 0, 1, 3, 0 );
00305     locTable.enableWriteAccess();
00306     for( int i = 0; i < NN; i++ )
00307     {
00308         if( is_inside[i] )
00309         {
00310             iargs[0] = TLforward_o.vi_rd[3 * i + 1];
00311             iargs[1] = TLforward_o.vi_rd[3 * i + 2];
00312             iargs[2] = locTable.get_n();
00313             TLsearch_results_o.push_back( iargs, NULL, NULL, NULL );
00314             Ulong ent_ulong = (Ulong)ents[i];
00315             sint forward    = (sint)TLforward_o.vi_rd[3 * i + 1];
00316             locTable.push_back( &forward, NULL, &ent_ulong, &params[3 * i] );
00317         }
00318     }
00319     locTable.disableWriteAccess();
00320 
00321     myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_since_birth() - tstart;
00322     myTimer.time_elapsed();  // call this to reset last time called
00323 
00324     // step 14: send TLsearch_results_o and receive TLloc_i
00325     if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLsearch_results_o, 0 );
00326 
00327     myTimes.slTimes[SpatialLocatorTimes::TARG_RETURN] = myTimer.time_elapsed();
00328 
00329     // store proc/index tuples in parLocTable
00330     parLocTable.initialize( 2, 0, 0, 0, num_points );
00331     parLocTable.enableWriteAccess();
00332     std::fill( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, -1 );
00333 
00334     for( unsigned int i = 0; i < TLsearch_results_o.get_n(); i++ )
00335     {
00336         int idx                        = TLsearch_results_o.vi_rd[3 * i + 1];
00337         parLocTable.vi_wr[2 * idx]     = TLsearch_results_o.vi_rd[3 * i];
00338         parLocTable.vi_wr[2 * idx + 1] = TLsearch_results_o.vi_rd[3 * i + 2];
00339     }
00340 
00341     if( debug )
00342     {
00343         int num_found =
00344             num_points - 0.5 * std::count_if( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, is_neg );
00345         std::cout << "Points found = " << num_found << "/" << num_points << " ("
00346                   << 100.0 * ( (double)num_found / num_points ) << "%)" << std::endl;
00347     }
00348 
00349     myTimes.slTimes[SpatialLocatorTimes::TARG_STORE] = myTimer.time_elapsed();
00350 
00351     return MB_SUCCESS;
00352 }
00353 
00354 #endif
00355 
00356 ErrorCode SpatialLocator::locate_points( Range& verts,
00357                                          const double rel_iter_tol,
00358                                          const double abs_iter_tol,
00359                                          const double inside_tol )
00360 {
00361     bool i_initialized = false;
00362     if( !timerInitialized )
00363     {
00364         myTimer.time_elapsed();
00365         timerInitialized = true;
00366         i_initialized    = true;
00367     }
00368 
00369     assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
00370     std::vector< double > pos( 3 * verts.size() );
00371     ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
00372     if( MB_SUCCESS != rval ) return rval;
00373     rval = locate_points( &pos[0], verts.size(), rel_iter_tol, abs_iter_tol, inside_tol );
00374     if( MB_SUCCESS != rval ) return rval;
00375 
00376     // only call this if I'm the top-level function, since it resets the last time called
00377     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00378 
00379     return MB_SUCCESS;
00380 }
00381 
00382 ErrorCode SpatialLocator::locate_points( const double* pos,
00383                                          int num_points,
00384                                          const double rel_iter_tol,
00385                                          const double abs_iter_tol,
00386                                          const double inside_tol )
00387 {
00388     bool i_initialized = false;
00389     if( !timerInitialized )
00390     {
00391         myTimer.time_elapsed();
00392         timerInitialized = true;
00393         i_initialized    = true;
00394     }
00395     // initialize to tuple structure (p_ui, hs_ul, r[3]_d) (see header comments for locTable)
00396     locTable.initialize( 1, 0, 1, 3, num_points );
00397     locTable.enableWriteAccess();
00398 
00399     // pass storage directly into locate_points, since we know those arrays are contiguous
00400     ErrorCode rval = locate_points( pos, num_points, (EntityHandle*)locTable.vul_wr, locTable.vr_wr, NULL, rel_iter_tol,
00401                                     abs_iter_tol, inside_tol );
00402     std::fill( locTable.vi_wr, locTable.vi_wr + num_points, 0 );
00403     locTable.set_n( num_points );
00404     if( MB_SUCCESS != rval ) return rval;
00405 
00406     // only call this if I'm the top-level function, since it resets the last time called
00407     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00408 
00409     return MB_SUCCESS;
00410 }
00411 
00412 ErrorCode SpatialLocator::locate_points( Range& verts,
00413                                          EntityHandle* ents,
00414                                          double* params,
00415                                          int* is_inside,
00416                                          const double rel_iter_tol,
00417                                          const double abs_iter_tol,
00418                                          const double inside_tol )
00419 {
00420     bool i_initialized = false;
00421     if( !timerInitialized )
00422     {
00423         myTimer.time_elapsed();
00424         timerInitialized = true;
00425         i_initialized    = true;
00426     }
00427 
00428     assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
00429     std::vector< double > pos( 3 * verts.size() );
00430     ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
00431     if( MB_SUCCESS != rval ) return rval;
00432     rval = locate_points( &pos[0], verts.size(), ents, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
00433 
00434     // only call this if I'm the top-level function, since it resets the last time called
00435     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00436 
00437     return rval;
00438 }
00439 
00440 ErrorCode SpatialLocator::locate_points( const double* pos,
00441                                          int num_points,
00442                                          EntityHandle* ents,
00443                                          double* params,
00444                                          int* is_inside,
00445                                          const double /* rel_iter_tol */,
00446                                          const double abs_iter_tol,
00447                                          const double inside_tol )
00448 {
00449     bool i_initialized = false;
00450     if( !timerInitialized )
00451     {
00452         myTimer.time_elapsed();
00453         timerInitialized = true;
00454         i_initialized    = true;
00455     }
00456 
00457     /*
00458     double tmp_abs_iter_tol = abs_iter_tol;
00459     if (rel_iter_tol && !tmp_abs_iter_tol) {
00460         // relative epsilon given, translate to absolute epsilon using box dimensions
00461       tmp_abs_iter_tol = rel_iter_tol * localBox.diagonal_length();
00462     }
00463     */
00464 
00465     if( elemEval && myTree->get_eval() != elemEval ) myTree->set_eval( elemEval );
00466 
00467     ErrorCode rval = MB_SUCCESS;
00468     for( int i = 0; i < num_points; i++ )
00469     {
00470         int i3 = 3 * i;
00471         ErrorCode tmp_rval =
00472             myTree->point_search( pos + i3, ents[i], abs_iter_tol, inside_tol, NULL, NULL, (CartVect*)( params + i3 ) );
00473         if( MB_SUCCESS != tmp_rval )
00474         {
00475             rval = tmp_rval;
00476             continue;
00477         }
00478 
00479         if( debug && !ents[i] )
00480         {
00481             std::cout << "Point " << i << " not found; point: (" << pos[i3] << "," << pos[i3 + 1] << "," << pos[i3 + 2]
00482                       << ")" << std::endl;
00483         }
00484 
00485         if( is_inside ) is_inside[i] = ( ents[i] ? true : false );
00486     }
00487 
00488     // only call this if I'm the top-level function, since it resets the last time called
00489     if( i_initialized ) myTimes.slTimes[SpatialLocatorTimes::SRC_SEARCH] = myTimer.time_elapsed();
00490 
00491     return rval;
00492 }
00493 
00494 /* Count the number of located points in locTable
00495  * Return the number of entries in locTable that have non-zero entity handles, which
00496  * represents the number of points in targetEnts that were inside one element in sourceEnts
00497  *
00498  */
00499 int SpatialLocator::local_num_located()
00500 {
00501     int num_located = locTable.get_n() - std::count( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
00502     if( num_located != (int)locTable.get_n() )
00503     {
00504         Ulong* nl = std::find( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
00505         if( nl )
00506         {
00507             int idx = nl - locTable.vul_rd;
00508             if( idx )
00509             {
00510             }
00511         }
00512     }
00513     return num_located;
00514 }
00515 
00516 /* Count the number of located points in parLocTable
00517  * Return the number of entries in parLocTable that have a non-negative index in on a remote
00518  * proc in parLocTable, which gives the number of points located in at least one element in a
00519  * remote proc's sourceEnts.
00520  */
00521 int SpatialLocator::remote_num_located()
00522 {
00523     int located = 0;
00524     for( unsigned int i = 0; i < parLocTable.get_n(); i++ )
00525         if( parLocTable.vi_rd[2 * i] != -1 ) located++;
00526     return located;
00527 }
00528 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines