Branch data Line data Source code
1 : : /**\file SpatialLocator.hpp
2 : : * \class moab::SpatialLocator
3 : : * \brief Tool to facilitate spatial location of a point in a mesh
4 : : *
5 : : * SpatialLocator facilitates searching for points in or performing ray traces on collections of
6 : : * mesh entities in 2D or 3D. This searching is facilitated by a tree-based decomposition of the
7 : : * mesh. Various types of trees are implemented in MOAB and can be used by this tool, but by
8 : : * default it uses AdaptiveKDTree (see child classes of Tree for which others are available).
9 : : * Parallel and serial searching are both supported.
10 : : *
11 : : * SpatialLocator can either cache the search results for points or pass back this information in
12 : : * arguments. Cached information is kept in locTable, indexed in the same order as search points
13 : : * passed in. This information consists of the entity handle containing the point and the
14 : : * parametric coordinates inside that element. Information about the points searched, e.g. the
15 : : * entities from which those points are derived, can be stored in the calling application if
16 : : * desired.
17 : : *
18 : : * In parallel, there is a separation between the proc deciding which points to search for (the
19 : : * "target" proc), and the proc locating the point in its local mesh (the "source" proc). On the
20 : : * source proc, location information is cached in locTable, as in the serial case. By default, this
21 : : * location information (handle and parametric coords) is not returned to the target proc, since it
22 : : * would be of no use there. Instead, the rank of the source proc locating the point, and the index
23 : : * of that location info in the source proc's locTable, is returned; this information is stored on
24 : : * the target proc in this class's parLocTable variable. Again, information about the points
25 : : * searched should be stored in the calling application, if desired.
26 : : *
27 : : * This class uses the ElemEvaluator class for specification and evaluation of basis functions (used
28 : : * for computing parametric coords within an entity). See documentation and examples for that class
29 : : * for usage information.
30 : : *
31 : : */
32 : :
33 : : #ifndef MOAB_SPATIALLOCATOR_HPP
34 : : #define MOAB_SPATIALLOCATOR_HPP
35 : :
36 : : #include "moab/Types.hpp"
37 : : #include "moab/Tree.hpp"
38 : : #include "moab/Range.hpp"
39 : : #include "moab/TupleList.hpp"
40 : : #include "moab/BoundBox.hpp"
41 : : #include "moab/SpatialLocatorTimes.hpp"
42 : : #include "moab/CpuTimer.hpp"
43 : :
44 : : #include <string>
45 : : #include <vector>
46 : : #include <map>
47 : : #include <math.h>
48 : :
49 : : namespace moab
50 : : {
51 : :
52 : : class Interface;
53 : : class ElemEvaluator;
54 : : class ParallelComm;
55 : :
56 : : class SpatialLocator
57 : : {
58 : : public:
59 : : /* constructor */
60 : : SpatialLocator( Interface* impl, Range& elems, Tree* tree = NULL, ElemEvaluator* eval = NULL );
61 : :
62 : : /* destructor */
63 : : virtual ~SpatialLocator();
64 : :
65 : : /* add elements to be searched */
66 : : ErrorCode add_elems( Range& elems );
67 : :
68 : : /* get bounding box of this locator */
69 : 4 : BoundBox& local_box()
70 : : {
71 : 4 : return localBox;
72 : : }
73 : :
74 : : /* get bounding box of this locator */
75 : : const BoundBox& local_box() const
76 : : {
77 : : return localBox;
78 : : }
79 : :
80 : : /* locate a set of vertices, Range variant */
81 : : ErrorCode locate_points( Range& vertices, EntityHandle* ents, double* params, int* is_inside = NULL,
82 : : const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
83 : : const double inside_tol = 1.0e-6 );
84 : :
85 : : /* locate a set of points */
86 : : ErrorCode locate_points( const double* pos, int num_points, EntityHandle* ents, double* params,
87 : : int* is_inside = NULL, const double rel_iter_tol = 1.0e-10,
88 : : const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 );
89 : :
90 : : /* locate a set of vertices or entity centroids, storing results on TupleList in this class
91 : : * Locate a set of vertices or entity centroids, storing the detailed results in member
92 : : * variable (TupleList) locTable (see comments on locTable for structure of that tuple).
93 : : */
94 : : ErrorCode locate_points( Range& ents, const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
95 : : const double inside_tol = 1.0e-6 );
96 : :
97 : : /* locate a set of points, storing results on TupleList in this class
98 : : * Locate a set of points, storing the detailed results in member variable (TupleList) locTable
99 : : * (see comments on locTable for structure of that tuple).
100 : : */
101 : : ErrorCode locate_points( const double* pos, int num_points, const double rel_iter_tol = 1.0e-10,
102 : : const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 );
103 : :
104 : : /* Count the number of located points in locTable
105 : : * Return the number of entries in locTable that have non-zero entity handles, which
106 : : * represents the number of points in targetEnts that were inside one element in sourceEnts
107 : : *
108 : : */
109 : : int local_num_located();
110 : :
111 : : /* Count the number of located points in parLocTable
112 : : * Return the number of entries in parLocTable that have a non-negative index in on a remote
113 : : * proc in parLocTable, which gives the number of points located in at least one element in a
114 : : * remote proc's sourceEnts.
115 : : */
116 : : int remote_num_located();
117 : :
118 : : #ifdef MOAB_HAVE_MPI
119 : : /* locate a set of vertices or entity centroids, storing results on TupleList in this class
120 : : * Locate a set of vertices or entity centroids, storing the detailed results in member
121 : : * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
122 : : * structure of those tuples).
123 : : */
124 : : ErrorCode par_locate_points( ParallelComm* pc, Range& vertices, const double rel_iter_tol = 1.0e-10,
125 : : const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 );
126 : :
127 : : /* locate a set of points, storing results on TupleList in this class
128 : : * Locate a set of points, storing the detailed results in member
129 : : * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
130 : : * structure of those tuples).
131 : : */
132 : : ErrorCode par_locate_points( ParallelComm* pc, const double* pos, int num_points,
133 : : const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
134 : : const double inside_tol = 1.0e-6 );
135 : : #endif
136 : :
137 : : /** \brief Return the MOAB interface associated with this locator
138 : : */
139 : : Interface* moab()
140 : : {
141 : : return mbImpl;
142 : : }
143 : :
144 : : /* locate a point */
145 : : ErrorCode locate_point( const double* pos, EntityHandle& ent, double* params, int* is_inside = NULL,
146 : : const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
147 : : const double inside_tol = 1.0e-6 );
148 : :
149 : : /* return the tree */
150 : 4 : Tree* get_tree()
151 : : {
152 : 4 : return myTree;
153 : : }
154 : :
155 : : /* get the locTable
156 : : */
157 : : TupleList& loc_table()
158 : : {
159 : : return locTable;
160 : : }
161 : :
162 : : /* get the locTable
163 : : */
164 : : const TupleList& loc_table() const
165 : : {
166 : : return locTable;
167 : : }
168 : :
169 : : /* get the parLocTable
170 : : */
171 : : TupleList& par_loc_table()
172 : : {
173 : : return parLocTable;
174 : : }
175 : :
176 : : /* get the parLocTable
177 : : */
178 : : const TupleList& par_loc_table() const
179 : : {
180 : : return parLocTable;
181 : : }
182 : :
183 : : /* get elemEval */
184 : : ElemEvaluator* elem_eval()
185 : : {
186 : : return elemEval;
187 : : }
188 : :
189 : : /* get elemEval */
190 : : const ElemEvaluator* elem_eval() const
191 : : {
192 : : return elemEval;
193 : : }
194 : :
195 : : /* set elemEval */
196 : : void elem_eval( ElemEvaluator* eval )
197 : : {
198 : : elemEval = eval;
199 : : if( myTree ) myTree->set_eval( eval );
200 : : }
201 : :
202 : : /** \brief Get spatial locator times object */
203 : : SpatialLocatorTimes& sl_times()
204 : : {
205 : : return myTimes;
206 : : }
207 : :
208 : : /** \brief Get spatial locator times object */
209 : : const SpatialLocatorTimes& sl_times() const
210 : : {
211 : : return myTimes;
212 : : }
213 : :
214 : : private:
215 : : #ifdef MOAB_HAVE_MPI
216 : : /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
217 : : */
218 : : ErrorCode initialize_intermediate_partition( ParallelComm* pc );
219 : :
220 : : /* for a given point in space, compute its ijk location in the intermediate decomposition;
221 : : * tolerance is used only to test proximity to global box extent, not for local box test
222 : : */
223 : : inline ErrorCode get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const;
224 : :
225 : : #if 0
226 : : /* given an ijk location in the intermediate partition, return the proc rank for that location
227 : : */
228 : : inline int proc_from_ijk(const int *ijk) const;
229 : : #endif
230 : :
231 : : /* given a point in space, return the proc responsible for that point from the intermediate
232 : : * decomp; no tolerances applied here, so first proc in lexicographic ijk ordering is returned
233 : : */
234 : : inline int proc_from_point( const double* pos, const double abs_iter_tol ) const;
235 : :
236 : : /* register my source mesh with intermediate decomposition procs
237 : : */
238 : : ErrorCode register_src_with_intermediate_procs( ParallelComm* pc, double abs_iter_tol, TupleList& TLreg_o );
239 : :
240 : : #endif
241 : :
242 : : /** Create a tree
243 : : * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree,
244 : : * otherwise creates a BVHTree.
245 : : */
246 : : void create_tree();
247 : :
248 : : /* MOAB instance */
249 : : Interface* mbImpl;
250 : :
251 : : /* elements being located */
252 : : Range myElems;
253 : :
254 : : /* dimension of entities in locator */
255 : : int myDim;
256 : :
257 : : /* tree used for location */
258 : : Tree* myTree;
259 : :
260 : : /* element evaluator */
261 : : ElemEvaluator* elemEval;
262 : :
263 : : /* whether I created the tree or not (determines whether to delete it or not on destruction) */
264 : : bool iCreatedTree;
265 : :
266 : : /* \brief local locations table
267 : : * This table stores detailed local location data results from locate_points, that is, location
268 : : * data for points located on the local mesh. Data is stored in a TupleList, where each tuple
269 : : * consists of (p_i, hs_ul, r[3]_d), where p_i = (int) proc from which request for this point
270 : : * was made (0 if serial) hs_ul = (unsigned long) source entity containing the point r[3]_d =
271 : : * (double) parametric coordinates of the point in hs
272 : : */
273 : : TupleList locTable;
274 : :
275 : : /* \brief parallel locations table
276 : : * This table stores information about points located on a local or remote processor. For
277 : : * points located on this processor's local mesh, detailed location data is stored in locTable.
278 : : * For points located on remote processors, more communication is required to retrieve specific
279 : : * location data (usually that information isn't useful on this processor).
280 : : *
281 : : * The tuple structure of this TupleList is (p_i, ri_i), where:
282 : : * p_i = (int) processor rank containing this point
283 : : * ri_i = (int) index into locTable on remote proc containing this point's location
284 : : * information The indexing of parLocTable corresponds to that of the points/entities passed in.
285 : : */
286 : : TupleList parLocTable;
287 : :
288 : : /* \brief Local bounding box, duplicated from tree
289 : : */
290 : : BoundBox localBox;
291 : :
292 : : /* \brief Global bounding box, used in parallel spatial location
293 : : */
294 : : BoundBox globalBox;
295 : :
296 : : /* \brief Regional delta xyz, used in parallel spatial location
297 : : */
298 : : CartVect regDeltaXYZ;
299 : :
300 : : /* \brief Number of regions in each of 3 directions
301 : : */
302 : : int regNums[3];
303 : :
304 : : /* \brief Map from source processor to bounding box of that proc's source mesh
305 : : *
306 : : */
307 : : std::map< int, BoundBox > srcProcBoxes;
308 : :
309 : : /* \brief Timing object for spatial location
310 : : */
311 : : SpatialLocatorTimes myTimes;
312 : :
313 : : /* \brief Timer object to manage overloaded search functions
314 : : */
315 : : CpuTimer myTimer;
316 : :
317 : : /* \brief Flag to manage initialization of timer for overloaded search functions
318 : : */
319 : : bool timerInitialized;
320 : : };
321 : :
322 : 6 : inline SpatialLocator::~SpatialLocator()
323 : : {
324 [ - + ][ # # ]: 2 : if( iCreatedTree && myTree ) delete myTree;
[ # # ]
325 [ - + ]: 4 : }
326 : :
327 : : inline ErrorCode SpatialLocator::locate_point( const double* pos, EntityHandle& ent, double* params, int* is_inside,
328 : : const double rel_iter_tol, const double abs_iter_tol,
329 : : const double inside_tol )
330 : : {
331 : : return locate_points( pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
332 : : }
333 : :
334 : : #ifdef MOAB_HAVE_MPI
335 : 0 : inline ErrorCode SpatialLocator::get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const
336 : : {
337 [ # # ]: 0 : for( int i = 0; i < 3; i++ )
338 : : {
339 [ # # ][ # # ]: 0 : if( point[i] < globalBox.bMin[i] - abs_iter_tol || point[i] > globalBox.bMax[i] + abs_iter_tol )
[ # # ]
340 : 0 : ijk[i] = -1;
341 : : else
342 : : {
343 : 0 : ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i];
344 [ # # ][ # # ]: 0 : if( ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i] + abs_iter_tol ) ijk[i] = regNums[i] - 1;
[ # # ]
345 : : }
346 : : }
347 : :
348 [ # # ][ # # ]: 0 : return ( ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE );
[ # # ]
349 : : ;
350 : : }
351 : :
352 : : #if 0
353 : : inline int SpatialLocator::proc_from_ijk(const int *ijk) const
354 : : {
355 : : return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
356 : : }
357 : : #endif
358 : :
359 : 0 : inline int SpatialLocator::proc_from_point( const double* pos, const double abs_iter_tol ) const
360 : : {
361 : : int ijk[3];
362 [ # # ][ # # ]: 0 : ErrorCode rval = get_point_ijk( CartVect( pos ), abs_iter_tol, ijk );
363 [ # # ]: 0 : if( MB_SUCCESS != rval ) return -1;
364 : :
365 : 0 : return ijk[2] * regNums[0] * regNums[1] + ijk[1] * regNums[0] + ijk[0];
366 : : }
367 : : #endif
368 : :
369 : : } // namespace moab
370 : :
371 : : #endif
|
