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