 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 <cmath>
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 BoundBox& local_box()
70 {
71 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,
82 EntityHandle* ents,
83 double* params,
84 int* is_inside = NULL,
85 const double rel_iter_tol = 1.0e-10,
86 const double abs_iter_tol = 1.0e-10,
87 const double inside_tol = 1.0e-6 );
88
89 /* locate a set of points */
90 ErrorCode locate_points( const double* pos,
91 int num_points,
92 EntityHandle* ents,
93 double* params,
94 int* is_inside = NULL,
95 const double rel_iter_tol = 1.0e-10,
96 const double abs_iter_tol = 1.0e-10,
97 const double inside_tol = 1.0e-6 );
98
99 /* locate a set of vertices or entity centroids, storing results on TupleList in this class
100 * Locate a set of vertices or entity centroids, storing the detailed results in member
101 * variable (TupleList) locTable (see comments on locTable for structure of that tuple).
102 */
103 ErrorCode locate_points( Range& ents,
104 const double rel_iter_tol = 1.0e-10,
105 const double abs_iter_tol = 1.0e-10,
106 const double inside_tol = 1.0e-6 );
107
108 /* locate a set of points, storing results on TupleList in this class
109 * Locate a set of points, storing the detailed results in member variable (TupleList) locTable
110 * (see comments on locTable for structure of that tuple).
111 */
112 ErrorCode locate_points( const double* pos,
113 int num_points,
114 const double rel_iter_tol = 1.0e-10,
115 const double abs_iter_tol = 1.0e-10,
116 const double inside_tol = 1.0e-6 );
117
118 /* Count the number of located points in locTable
119 * Return the number of entries in locTable that have non-zero entity handles, which
120 * represents the number of points in targetEnts that were inside one element in sourceEnts
121 *
122 */
123 int local_num_located();
124
125 /* Count the number of located points in parLocTable
126 * Return the number of entries in parLocTable that have a non-negative index in on a remote
127 * proc in parLocTable, which gives the number of points located in at least one element in a
128 * remote proc's sourceEnts.
129 */
130 int remote_num_located();
131
132 #ifdef MOAB_HAVE_MPI
133 /* locate a set of vertices or entity centroids, storing results on TupleList in this class
134 * Locate a set of vertices or entity centroids, storing the detailed results in member
135 * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
136 * structure of those tuples).
137 */
138 ErrorCode par_locate_points( ParallelComm* pc,
139 Range& vertices,
140 const double rel_iter_tol = 1.0e-10,
141 const double abs_iter_tol = 1.0e-10,
142 const double inside_tol = 1.0e-6 );
143
144 /* locate a set of points, storing results on TupleList in this class
145 * Locate a set of points, storing the detailed results in member
146 * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
147 * structure of those tuples).
148 */
149 ErrorCode par_locate_points( ParallelComm* pc,
150 const double* pos,
151 int num_points,
152 const double rel_iter_tol = 1.0e-10,
153 const double abs_iter_tol = 1.0e-10,
154 const double inside_tol = 1.0e-6 );
155 #endif
156
157 /** \brief Return the MOAB interface associated with this locator
158 */
159 Interface* moab()
160 {
161 return mbImpl;
162 }
163
164 /* locate a point */
165 ErrorCode locate_point( const double* pos,
166 EntityHandle& ent,
167 double* params,
168 int* is_inside = NULL,
169 const double rel_iter_tol = 1.0e-10,
170 const double abs_iter_tol = 1.0e-10,
171 const double inside_tol = 1.0e-6 );
172
173 /* return the tree */
174 Tree* get_tree()
175 {
176 return myTree;
177 }
178
179 /* get the locTable
180 */
181 TupleList& loc_table()
182 {
183 return locTable;
184 }
185
186 /* get the locTable
187 */
188 const TupleList& loc_table() const
189 {
190 return locTable;
191 }
192
193 /* get the parLocTable
194 */
195 TupleList& par_loc_table()
196 {
197 return parLocTable;
198 }
199
200 /* get the parLocTable
201 */
202 const TupleList& par_loc_table() const
203 {
204 return parLocTable;
205 }
206
207 /* get elemEval */
208 ElemEvaluator* elem_eval()
209 {
210 return elemEval;
211 }
212
213 /* get elemEval */
214 const ElemEvaluator* elem_eval() const
215 {
216 return elemEval;
217 }
218
219 /* set elemEval */
220 void elem_eval( ElemEvaluator* eval )
221 {
222 elemEval = eval;
223 if( myTree ) myTree->set_eval( eval );
224 }
225
226 /** \brief Get spatial locator times object */
227 SpatialLocatorTimes& sl_times()
228 {
229 return myTimes;
230 }
231
232 /** \brief Get spatial locator times object */
233 const SpatialLocatorTimes& sl_times() const
234 {
235 return myTimes;
236 }
237
238 private:
239 #ifdef MOAB_HAVE_MPI
240 /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
241 */
242 ErrorCode initialize_intermediate_partition( ParallelComm* pc );
243
244 /* for a given point in space, compute its ijk location in the intermediate decomposition;
245 * tolerance is used only to test proximity to global box extent, not for local box test
246 */
247 inline ErrorCode get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const;
248
249 #if 0
250 /* given an ijk location in the intermediate partition, return the proc rank for that location
251 */
252 inline int proc_from_ijk(const int *ijk) const;
253 #endif
254
255 /* given a point in space, return the proc responsible for that point from the intermediate
256 * decomp; no tolerances applied here, so first proc in lexicographic ijk ordering is returned
257 */
258 inline int proc_from_point( const double* pos, const double abs_iter_tol ) const;
259
260 /* register my source mesh with intermediate decomposition procs
261 */
262 ErrorCode register_src_with_intermediate_procs( ParallelComm* pc, double abs_iter_tol, TupleList& TLreg_o );
263
264 #endif
265
266 /** Create a tree
267 * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree,
268 * otherwise creates a BVHTree.
269 */
270 void create_tree();
271
272 /* MOAB instance */
273 Interface* mbImpl;
274
275 /* elements being located */
276 Range myElems;
277
278 /* dimension of entities in locator */
279 int myDim;
280
281 /* tree used for location */
282 Tree* myTree;
283
284 /* element evaluator */
285 ElemEvaluator* elemEval;
286
287 /* whether I created the tree or not (determines whether to delete it or not on destruction) */
288 bool iCreatedTree;
289
290 /* \brief local locations table
291 * This table stores detailed local location data results from locate_points, that is, location
292 * data for points located on the local mesh. Data is stored in a TupleList, where each tuple
293 * consists of (p_i, hs_ul, r[3]_d), where p_i = (int) proc from which request for this point
294 * was made (0 if serial) hs_ul = (unsigned long) source entity containing the point r[3]_d =
295 * (double) parametric coordinates of the point in hs
296 */
297 TupleList locTable;
298
299 /* \brief parallel locations table
300 * This table stores information about points located on a local or remote processor. For
301 * points located on this processor's local mesh, detailed location data is stored in locTable.
302 * For points located on remote processors, more communication is required to retrieve specific
303 * location data (usually that information isn't useful on this processor).
304 *
305 * The tuple structure of this TupleList is (p_i, ri_i), where:
306 * p_i = (int) processor rank containing this point
307 * ri_i = (int) index into locTable on remote proc containing this point's location
308 * information The indexing of parLocTable corresponds to that of the points/entities passed in.
309 */
310 TupleList parLocTable;
311
312 /* \brief Local bounding box, duplicated from tree
313 */
314 BoundBox localBox;
315
316 /* \brief Global bounding box, used in parallel spatial location
317 */
318 BoundBox globalBox;
319
320 /* \brief Regional delta xyz, used in parallel spatial location
321 */
322 CartVect regDeltaXYZ;
323
324 /* \brief Number of regions in each of 3 directions
325 */
326 int regNums[3];
327
328 /* \brief Map from source processor to bounding box of that proc's source mesh
329 *
330 */
331 std::map< int, BoundBox > srcProcBoxes;
332
333 /* \brief Timing object for spatial location
334 */
335 SpatialLocatorTimes myTimes;
336
337 /* \brief Timer object to manage overloaded search functions
338 */
339 CpuTimer myTimer;
340
341 /* \brief Flag to manage initialization of timer for overloaded search functions
342 */
343 bool timerInitialized;
344 };
345
346 inline SpatialLocator::~SpatialLocator()
347 {
348 if( iCreatedTree && myTree ) delete myTree;
349 }
350
351 inline ErrorCode SpatialLocator::locate_point( const double* pos,
352 EntityHandle& ent,
353 double* params,
354 int* is_inside,
355 const double rel_iter_tol,
356 const double abs_iter_tol,
357 const double inside_tol )
358 {
359 return locate_points( pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
360 }
361
362 #ifdef MOAB_HAVE_MPI
363 inline ErrorCode SpatialLocator::get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const
364 {
365 for( int i = 0; i < 3; i++ )
366 {
367 if( point[i] < globalBox.bMin[i] - abs_iter_tol || point[i] > globalBox.bMax[i] + abs_iter_tol )
368 ijk[i] = -1;
369 else
370 {
371 ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i];
372 if( ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i] + abs_iter_tol ) ijk[i] = regNums[i] - 1;
373 }
374 }
375
376 return ( ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE );
377 ;
378 }
379
380 #if 0
381 inline int SpatialLocator::proc_from_ijk(const int *ijk) const
382 {
383 return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
384 }
385 #endif
386
387 inline int SpatialLocator::proc_from_point( const double* pos, const double abs_iter_tol ) const
388 {
389 int ijk[3];
390 ErrorCode rval = get_point_ijk( CartVect( pos ), abs_iter_tol, ijk );
391 if( MB_SUCCESS != rval ) return -1;
392
393 return ijk[2] * regNums[0] * regNums[1] + ijk[1] * regNums[0] + ijk[0];
394 }
395 #endif
396
397 } // namespace moab
398
399 #endif
1.9.1.