MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /**\file SpatialLocator.hpp 00002 * \class moab::SpatialLocator 00003 * \brief Tool to facilitate spatial location of a point in a mesh 00004 * 00005 * SpatialLocator facilitates searching for points in or performing ray traces on collections of 00006 * mesh entities in 2D or 3D. This searching is facilitated by a tree-based decomposition of the 00007 * mesh. Various types of trees are implemented in MOAB and can be used by this tool, but by 00008 * default it uses AdaptiveKDTree (see child classes of Tree for which others are available). 00009 * Parallel and serial searching are both supported. 00010 * 00011 * SpatialLocator can either cache the search results for points or pass back this information in 00012 * arguments. Cached information is kept in locTable, indexed in the same order as search points 00013 * passed in. This information consists of the entity handle containing the point and the 00014 * parametric coordinates inside that element. Information about the points searched, e.g. the 00015 * entities from which those points are derived, can be stored in the calling application if 00016 * desired. 00017 * 00018 * In parallel, there is a separation between the proc deciding which points to search for (the 00019 * "target" proc), and the proc locating the point in its local mesh (the "source" proc). On the 00020 * source proc, location information is cached in locTable, as in the serial case. By default, this 00021 * location information (handle and parametric coords) is not returned to the target proc, since it 00022 * would be of no use there. Instead, the rank of the source proc locating the point, and the index 00023 * of that location info in the source proc's locTable, is returned; this information is stored on 00024 * the target proc in this class's parLocTable variable. Again, information about the points 00025 * searched should be stored in the calling application, if desired. 00026 * 00027 * This class uses the ElemEvaluator class for specification and evaluation of basis functions (used 00028 * for computing parametric coords within an entity). See documentation and examples for that class 00029 * for usage information. 00030 * 00031 */ 00032 00033 #ifndef MOAB_SPATIALLOCATOR_HPP 00034 #define MOAB_SPATIALLOCATOR_HPP 00035 00036 #include "moab/Types.hpp" 00037 #include "moab/Tree.hpp" 00038 #include "moab/Range.hpp" 00039 #include "moab/TupleList.hpp" 00040 #include "moab/BoundBox.hpp" 00041 #include "moab/SpatialLocatorTimes.hpp" 00042 #include "moab/CpuTimer.hpp" 00043 00044 #include <string> 00045 #include <vector> 00046 #include <map> 00047 #include <cmath> 00048 00049 namespace moab 00050 { 00051 00052 class Interface; 00053 class ElemEvaluator; 00054 class ParallelComm; 00055 00056 class SpatialLocator 00057 { 00058 public: 00059 /* constructor */ 00060 SpatialLocator( Interface* impl, Range& elems, Tree* tree = NULL, ElemEvaluator* eval = NULL ); 00061 00062 /* destructor */ 00063 virtual ~SpatialLocator(); 00064 00065 /* add elements to be searched */ 00066 ErrorCode add_elems( Range& elems ); 00067 00068 /* get bounding box of this locator */ 00069 BoundBox& local_box() 00070 { 00071 return localBox; 00072 } 00073 00074 /* get bounding box of this locator */ 00075 const BoundBox& local_box() const 00076 { 00077 return localBox; 00078 } 00079 00080 /* locate a set of vertices, Range variant */ 00081 ErrorCode locate_points( Range& vertices, EntityHandle* ents, double* params, int* is_inside = NULL, 00082 const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10, 00083 const double inside_tol = 1.0e-6 ); 00084 00085 /* locate a set of points */ 00086 ErrorCode locate_points( const double* pos, int num_points, EntityHandle* ents, double* params, 00087 int* is_inside = NULL, const double rel_iter_tol = 1.0e-10, 00088 const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 ); 00089 00090 /* locate a set of vertices or entity centroids, storing results on TupleList in this class 00091 * Locate a set of vertices or entity centroids, storing the detailed results in member 00092 * variable (TupleList) locTable (see comments on locTable for structure of that tuple). 00093 */ 00094 ErrorCode locate_points( Range& ents, const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10, 00095 const double inside_tol = 1.0e-6 ); 00096 00097 /* locate a set of points, storing results on TupleList in this class 00098 * Locate a set of points, storing the detailed results in member variable (TupleList) locTable 00099 * (see comments on locTable for structure of that tuple). 00100 */ 00101 ErrorCode locate_points( const double* pos, int num_points, const double rel_iter_tol = 1.0e-10, 00102 const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 ); 00103 00104 /* Count the number of located points in locTable 00105 * Return the number of entries in locTable that have non-zero entity handles, which 00106 * represents the number of points in targetEnts that were inside one element in sourceEnts 00107 * 00108 */ 00109 int local_num_located(); 00110 00111 /* Count the number of located points in parLocTable 00112 * Return the number of entries in parLocTable that have a non-negative index in on a remote 00113 * proc in parLocTable, which gives the number of points located in at least one element in a 00114 * remote proc's sourceEnts. 00115 */ 00116 int remote_num_located(); 00117 00118 #ifdef MOAB_HAVE_MPI 00119 /* locate a set of vertices or entity centroids, storing results on TupleList in this class 00120 * Locate a set of vertices or entity centroids, storing the detailed results in member 00121 * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for 00122 * structure of those tuples). 00123 */ 00124 ErrorCode par_locate_points( ParallelComm* pc, Range& vertices, const double rel_iter_tol = 1.0e-10, 00125 const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 ); 00126 00127 /* locate a set of points, storing results on TupleList in this class 00128 * Locate a set of points, storing the detailed results in member 00129 * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for 00130 * structure of those tuples). 00131 */ 00132 ErrorCode par_locate_points( ParallelComm* pc, const double* pos, int num_points, 00133 const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10, 00134 const double inside_tol = 1.0e-6 ); 00135 #endif 00136 00137 /** \brief Return the MOAB interface associated with this locator 00138 */ 00139 Interface* moab() 00140 { 00141 return mbImpl; 00142 } 00143 00144 /* locate a point */ 00145 ErrorCode locate_point( const double* pos, EntityHandle& ent, double* params, int* is_inside = NULL, 00146 const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10, 00147 const double inside_tol = 1.0e-6 ); 00148 00149 /* return the tree */ 00150 Tree* get_tree() 00151 { 00152 return myTree; 00153 } 00154 00155 /* get the locTable 00156 */ 00157 TupleList& loc_table() 00158 { 00159 return locTable; 00160 } 00161 00162 /* get the locTable 00163 */ 00164 const TupleList& loc_table() const 00165 { 00166 return locTable; 00167 } 00168 00169 /* get the parLocTable 00170 */ 00171 TupleList& par_loc_table() 00172 { 00173 return parLocTable; 00174 } 00175 00176 /* get the parLocTable 00177 */ 00178 const TupleList& par_loc_table() const 00179 { 00180 return parLocTable; 00181 } 00182 00183 /* get elemEval */ 00184 ElemEvaluator* elem_eval() 00185 { 00186 return elemEval; 00187 } 00188 00189 /* get elemEval */ 00190 const ElemEvaluator* elem_eval() const 00191 { 00192 return elemEval; 00193 } 00194 00195 /* set elemEval */ 00196 void elem_eval( ElemEvaluator* eval ) 00197 { 00198 elemEval = eval; 00199 if( myTree ) myTree->set_eval( eval ); 00200 } 00201 00202 /** \brief Get spatial locator times object */ 00203 SpatialLocatorTimes& sl_times() 00204 { 00205 return myTimes; 00206 } 00207 00208 /** \brief Get spatial locator times object */ 00209 const SpatialLocatorTimes& sl_times() const 00210 { 00211 return myTimes; 00212 } 00213 00214 private: 00215 #ifdef MOAB_HAVE_MPI 00216 /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box 00217 */ 00218 ErrorCode initialize_intermediate_partition( ParallelComm* pc ); 00219 00220 /* for a given point in space, compute its ijk location in the intermediate decomposition; 00221 * tolerance is used only to test proximity to global box extent, not for local box test 00222 */ 00223 inline ErrorCode get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const; 00224 00225 #if 0 00226 /* given an ijk location in the intermediate partition, return the proc rank for that location 00227 */ 00228 inline int proc_from_ijk(const int *ijk) const; 00229 #endif 00230 00231 /* given a point in space, return the proc responsible for that point from the intermediate 00232 * decomp; no tolerances applied here, so first proc in lexicographic ijk ordering is returned 00233 */ 00234 inline int proc_from_point( const double* pos, const double abs_iter_tol ) const; 00235 00236 /* register my source mesh with intermediate decomposition procs 00237 */ 00238 ErrorCode register_src_with_intermediate_procs( ParallelComm* pc, double abs_iter_tol, TupleList& TLreg_o ); 00239 00240 #endif 00241 00242 /** Create a tree 00243 * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree, 00244 * otherwise creates a BVHTree. 00245 */ 00246 void create_tree(); 00247 00248 /* MOAB instance */ 00249 Interface* mbImpl; 00250 00251 /* elements being located */ 00252 Range myElems; 00253 00254 /* dimension of entities in locator */ 00255 int myDim; 00256 00257 /* tree used for location */ 00258 Tree* myTree; 00259 00260 /* element evaluator */ 00261 ElemEvaluator* elemEval; 00262 00263 /* whether I created the tree or not (determines whether to delete it or not on destruction) */ 00264 bool iCreatedTree; 00265 00266 /* \brief local locations table 00267 * This table stores detailed local location data results from locate_points, that is, location 00268 * data for points located on the local mesh. Data is stored in a TupleList, where each tuple 00269 * consists of (p_i, hs_ul, r[3]_d), where p_i = (int) proc from which request for this point 00270 * was made (0 if serial) hs_ul = (unsigned long) source entity containing the point r[3]_d = 00271 * (double) parametric coordinates of the point in hs 00272 */ 00273 TupleList locTable; 00274 00275 /* \brief parallel locations table 00276 * This table stores information about points located on a local or remote processor. For 00277 * points located on this processor's local mesh, detailed location data is stored in locTable. 00278 * For points located on remote processors, more communication is required to retrieve specific 00279 * location data (usually that information isn't useful on this processor). 00280 * 00281 * The tuple structure of this TupleList is (p_i, ri_i), where: 00282 * p_i = (int) processor rank containing this point 00283 * ri_i = (int) index into locTable on remote proc containing this point's location 00284 * information The indexing of parLocTable corresponds to that of the points/entities passed in. 00285 */ 00286 TupleList parLocTable; 00287 00288 /* \brief Local bounding box, duplicated from tree 00289 */ 00290 BoundBox localBox; 00291 00292 /* \brief Global bounding box, used in parallel spatial location 00293 */ 00294 BoundBox globalBox; 00295 00296 /* \brief Regional delta xyz, used in parallel spatial location 00297 */ 00298 CartVect regDeltaXYZ; 00299 00300 /* \brief Number of regions in each of 3 directions 00301 */ 00302 int regNums[3]; 00303 00304 /* \brief Map from source processor to bounding box of that proc's source mesh 00305 * 00306 */ 00307 std::map< int, BoundBox > srcProcBoxes; 00308 00309 /* \brief Timing object for spatial location 00310 */ 00311 SpatialLocatorTimes myTimes; 00312 00313 /* \brief Timer object to manage overloaded search functions 00314 */ 00315 CpuTimer myTimer; 00316 00317 /* \brief Flag to manage initialization of timer for overloaded search functions 00318 */ 00319 bool timerInitialized; 00320 }; 00321 00322 inline SpatialLocator::~SpatialLocator() 00323 { 00324 if( iCreatedTree && myTree ) delete myTree; 00325 } 00326 00327 inline ErrorCode SpatialLocator::locate_point( const double* pos, EntityHandle& ent, double* params, int* is_inside, 00328 const double rel_iter_tol, const double abs_iter_tol, 00329 const double inside_tol ) 00330 { 00331 return locate_points( pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol ); 00332 } 00333 00334 #ifdef MOAB_HAVE_MPI 00335 inline ErrorCode SpatialLocator::get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const 00336 { 00337 for( int i = 0; i < 3; i++ ) 00338 { 00339 if( point[i] < globalBox.bMin[i] - abs_iter_tol || point[i] > globalBox.bMax[i] + abs_iter_tol ) 00340 ijk[i] = -1; 00341 else 00342 { 00343 ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i]; 00344 if( ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i] + abs_iter_tol ) ijk[i] = regNums[i] - 1; 00345 } 00346 } 00347 00348 return ( ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE ); 00349 ; 00350 } 00351 00352 #if 0 00353 inline int SpatialLocator::proc_from_ijk(const int *ijk) const 00354 { 00355 return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0]; 00356 } 00357 #endif 00358 00359 inline int SpatialLocator::proc_from_point( const double* pos, const double abs_iter_tol ) const 00360 { 00361 int ijk[3]; 00362 ErrorCode rval = get_point_ijk( CartVect( pos ), abs_iter_tol, ijk ); 00363 if( MB_SUCCESS != rval ) return -1; 00364 00365 return ijk[2] * regNums[0] * regNums[1] + ijk[1] * regNums[0] + ijk[0]; 00366 } 00367 #endif 00368 00369 } // namespace moab 00370 00371 #endif