moab
|
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], ¶ms[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, ¶ms[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