Branch data Line data Source code
1 : : #ifndef MOAB_GEOM_QUERY_TOOL_HPP
2 : : #define MOAB_GEOM_QUERY_TOOL_HPP
3 : :
4 : : #include "MBTagConventions.hpp"
5 : : #include "moab/CartVect.hpp"
6 : : #include "moab/Range.hpp"
7 : : #include "moab/Core.hpp"
8 : : #include "moab/GeomUtil.hpp"
9 : : #include "moab/FileOptions.hpp"
10 : : #include "moab/EntityHandle.hpp"
11 : : #include "moab/GeomTopoTool.hpp"
12 : : #include "moab/OrientedBoxTreeTool.hpp"
13 : :
14 : : #include <vector>
15 : : #include <map>
16 : : #include <string>
17 : : #include <assert.h>
18 : :
19 : : namespace moab
20 : : {
21 : :
22 : : /** \class GeomQueryTool
23 : : *
24 : : * \brief Tool for querying different aspects of geometric topology sets in MOAB
25 : : *
26 : : * Given the conventions established in GeomTopoTool for representing
27 : : * geometric topology through a hierarchy of meshsets, this tool provides a
28 : : * set of methods to query different geometric aspects of those geometric
29 : : * topology sets including:
30 : : *
31 : : * * measures of surface area and volume
32 : : * * distance to a bounding surface from a point within a volume
33 : : * * test the inclusion of a point within a volume
34 : : * * find the angle of a surface at a point
35 : : *
36 : : * A feature of these queries is that there is some support for overlapping
37 : : * geometries.
38 : : *
39 : : */
40 : :
41 : : class GeomQueryTool
42 : : {
43 : : public:
44 : : /* \class RayHistory
45 : : *
46 : : * In many use cases, it is useful to track some of the history of a ray as
47 : : * it passes through a geometry, particularly a geometry represented by
48 : : * facets. For example, given round-off erorr in ray-triangle tests
49 : : * (GeomUtil::ray_tri_intersect) used as part of a test for ray-surface
50 : : * intersection, it is possible for subsequent queries to oscillate between
51 : : * adjacent surfaces. This class stores information about history of a ray
52 : : * that can be used to test for such circumstances so that they can be
53 : : * accommodated.
54 : : *
55 : : */
56 : :
57 : 152 : class RayHistory
58 : : {
59 : :
60 : : public:
61 : : /**
62 : : * Clear this entire history-- logically equivalent to creating a new history,
63 : : * but probably more efficient.
64 : : */
65 : : void reset();
66 : :
67 : : /**
68 : : * Clear the history up to the most recent intersection. This should be
69 : : * called when a ray changes direction at the site of a surface crossing,
70 : : * a situation that most commonly occurs at a reflecting boundary.
71 : : */
72 : : void reset_to_last_intersection();
73 : :
74 : : /**
75 : : * Remove the most recent intersection. This allows a subsequent call
76 : : * along the same ray to return the same intersection.
77 : : */
78 : : void rollback_last_intersection();
79 : :
80 : : /**
81 : : * Get the last intersection in the RayHistory. This will return a null
82 : : * EntityHandle (0) if the history is empty.
83 : : */
84 : : ErrorCode get_last_intersection( EntityHandle& last_facet_hit ) const;
85 : :
86 : : /**
87 : : * @return the number of surface crossings currently represented by this ray history
88 : : */
89 : 1 : int size() const
90 : : {
91 : 1 : return prev_facets.size();
92 : : }
93 : :
94 : : /**
95 : : * @return Boolean indicating if this entity is in the RayHistory
96 : : */
97 : : bool in_history( EntityHandle ent ) const;
98 : :
99 : : /**
100 : : * Add entity to the RayHistory
101 : : */
102 : : void add_entity( EntityHandle ent );
103 : :
104 : : private:
105 : : std::vector< EntityHandle > prev_facets;
106 : :
107 : : friend class GeomQueryTool;
108 : : };
109 : :
110 : : // Constructors
111 : :
112 : : GeomQueryTool( Interface* impl, bool find_geoments = false, EntityHandle modelRootSet = 0,
113 : : bool p_rootSets_vector = true, bool restore_rootSets = true, bool trace_counting = false,
114 : : double overlap_thickness = 0., double numerical_precision = 0.001 );
115 : :
116 : : GeomQueryTool( GeomTopoTool* geomtopotool, bool trace_counting = false, double overlap_thickness = 0.,
117 : : double numerical_precision = 0.001 );
118 : :
119 : : // Destructor
120 : : ~GeomQueryTool();
121 : :
122 : : ErrorCode initialize();
123 : :
124 : : /**\brief find the next surface crossing from a given point in a given direction
125 : : *
126 : : * This is the primary method to enable ray tracing through a geometry.
127 : : * Given a volume and a ray, it determines the distance to the nearest intersection
128 : : * with a bounding surface of that volume and returns that distance and the
129 : : * EntityHandle indicating on which surface that intersection occurs.
130 : : * The caller can compute the location of the intersection by adding the
131 : : * distance to the ray.
132 : : *
133 : : * When a series of calls to this function are made along the same ray (e.g. for
134 : : * the purpose of tracking a ray through several volumes), the optional history
135 : : * argument should be given. The history prevents previously intersected facets
136 : : * from being intersected again. A single history should be used as long as a
137 : : * ray is proceeding forward without changing direction. This situation is
138 : : * sometimes referred to as "streaming."
139 : : *
140 : : * If a ray changes direction at an intersection site, the caller should call
141 : : * reset_to_last_intersection() on the history object before the next ray fire.
142 : : *
143 : : * @param volume The volume at which to fire the ray
144 : : * @param ray_start An array of x,y,z coordinates from which to start the ray.
145 : : * @param ray_dir An array of x,y,z coordinates indicating the direction of the ray.
146 : : * Must be of unit length.
147 : : * @param next_surf Output parameter indicating the next surface intersected by the ray.
148 : : * If no intersection is found, will be set to 0.
149 : : * @param next_surf_dist Output parameter indicating distance to next_surf. If next_surf is
150 : : * 0, this value is undefined and should not be used.
151 : : * @param history Optional RayHistory object. If provided, the facets in the history are
152 : : * assumed to not intersect with the given ray. The facet intersected
153 : : * by this query will also be added to the history.
154 : : * @param dist_limit Optional distance limit. If provided and > 0, no intersections at a
155 : : * distance further than this value will be returned.
156 : : * @param ray_orientation Optional ray orientation. If provided determines intersections
157 : : * along the normal provided, e.g. if -1 allows intersections back along the
158 : : * the ray direction, Default is 1, i.e. exit intersections
159 : : * @param stats Optional TrvStats object used to measure performance of underlying OBB
160 : : * ray-firing query. See OrientedBoxTreeTool.hpp for details.
161 : : *
162 : : */
163 : : ErrorCode ray_fire( const EntityHandle volume, const double ray_start[3], const double ray_dir[3],
164 : : EntityHandle& next_surf, double& next_surf_dist, RayHistory* history = NULL,
165 : : double dist_limit = 0, int ray_orientation = 1, OrientedBoxTreeTool::TrvStats* stats = NULL );
166 : :
167 : : /**\brief Test if a point is inside or outside a volume
168 : : *
169 : : * This method finds the point on the boundary of the volume that is nearest
170 : : * the test point (x,y,z). If that point is "close" to a surface, a boundary test
171 : : * is performed based on the normal of the surface at that point and the
172 : : * optional ray direction (u,v,w).
173 : : * @param volume The volume to test
174 : : * @param xyz The location to test for volume containment
175 : : * @param result Set to 0 if xyz it outside volume, 1 if inside, and -1 if on boundary.
176 : : * @param Optional direction to use for underlying ray fire query. Used to ensure
177 : : * consistent results when a ray direction is known. If NULL or {0,0,0} is
178 : : * given, a random direction will be used.
179 : : * @param history Optional RayHistory object to pass to underlying ray fire query.
180 : : * The history is not modified by this call.
181 : : */
182 : : ErrorCode point_in_volume( const EntityHandle volume, const double xyz[3], int& result, const double* uvw = NULL,
183 : : const RayHistory* history = NULL );
184 : :
185 : : /**\brief Robust test if a point is inside or outside a volume using unit sphere area method
186 : : *
187 : : * This test may be more robust that the standard point_in_volume, but is much slower.
188 : : * It does not detect 'on boundary' situations as point_in_volume does.
189 : : * @param volume The volume to test
190 : : * @param xyz The location to test for volume containment
191 : : * @param result Set to 0 if xyz it outside volume, 1 if inside.
192 : : */
193 : : ErrorCode point_in_volume_slow( const EntityHandle volume, const double xyz[3], int& result );
194 : :
195 : : /**\brief Find volume for a given location.
196 : : *
197 : : * Determines which volume contains the point if possible. If no volume is found,
198 : : * a null EntityHandle is returned along with a MB_ENTITY_NOT_FOUND ErrorCode.
199 : : * @param xyz The location to test
200 : : * @param volume Set to volume handle containing the location if found
201 : : * @param dir Optional direction to use for underlying ray fire query. Used to ensure
202 : : * consistent results when a ray direction is known. If NULL or {0,0,0} is
203 : : * given, a random direction will be used.
204 : : */
205 : : ErrorCode find_volume( const double xyz[3], EntityHandle& volume, const double* dir = NULL );
206 : :
207 : : /**\brief Find volume for a given location using loop. (slow)
208 : : *
209 : : * Loops over all volumes in the model, checking for point containment
210 : : * @param xyz The location to test
211 : : * @param volume Set to volume handle containing the location if found
212 : : * @param dir Optional direction to use for underlying ray fire query. Used to ensure
213 : : * consistent results when a ray direction is known. If NULL or {0,0,0} is
214 : : * given, a random direction will be used.
215 : : */
216 : : ErrorCode find_volume_slow( const double xyz[3], EntityHandle& volume, const double* dir = NULL );
217 : :
218 : : /**\brief Test if a point is inside or outsize a volume's axis-aligned bounding box
219 : : *
220 : : * This is used as a fast way to determine whether a point is inside or outside of a volume
221 : : * before performing a point_in_volume test which involves firing a ray.
222 : : * @param volume The volume to test
223 : : * @param point The location to test for bounding box containment
224 : : * @param inside set to 0 if point is outside the box, 1 if inside
225 : : */
226 : : ErrorCode point_in_box( const EntityHandle volume, const double point[3], int& inside );
227 : :
228 : : /** \brief Given a ray starting at a surface of a volume, check whether the ray enters or exits
229 : : * the volume
230 : : *
231 : : * This function is most useful for rays that change directions at a surface crossing.
232 : : * It can be used to check whether a direction change redirects the ray back into the
233 : : * originating volume.
234 : : *
235 : : * @param volume The volume to test
236 : : * @param surface A surface on volume
237 : : * @param xyz A point location on surface
238 : : * @param uvw A (unit) direction vector
239 : : * @param result Set to 1 if ray is entering volume, or 0 if it is leaving
240 : : * @param history Optional ray history object from a previous call to ray_fire. If present and
241 : : * non-empty, the history is used to look up the surface facet at which the ray begins. Absent
242 : : * a history, the facet nearest to xyz will be looked up. The history should always be provided
243 : : * if available, as it avoids the computational expense of a nearest-facet query.
244 : : */
245 : : ErrorCode test_volume_boundary( const EntityHandle volume, const EntityHandle surface, const double xyz[3],
246 : : const double uvw[3], int& result, const RayHistory* history = NULL );
247 : :
248 : : /**\brief Find the distance to the point on the boundary of the volume closest to the test point
249 : : *
250 : : * @param volume Volume to query
251 : : * @param point Coordinates of test point
252 : : * @param result Set to the minimum distance from point to a surface in volume
253 : : */
254 : : ErrorCode closest_to_location( EntityHandle volume, const double point[3], double& result,
255 : : EntityHandle* closest_surface = 0 );
256 : :
257 : : /** Calculate the volume contained in a 'volume' */
258 : : ErrorCode measure_volume( EntityHandle volume, double& result );
259 : :
260 : : /** Calculate sum of area of triangles */
261 : : ErrorCode measure_area( EntityHandle surface, double& result );
262 : :
263 : : /** Get the normal to a given surface at the point on the surface closest to a given point
264 : : *
265 : : * This method first identifies which facet contains this point and then
266 : : * calculates the unit outward normal of that facet. The facet of the
267 : : * provided volume that is nearest the provided point is used for this
268 : : * calculation. The search for that facet can be circumvented by providing
269 : : * a RayHistory, in which case the last facet of the history will be used.
270 : : *
271 : : * @param surf Surface on which to get normal
272 : : * @param xyz Point on surf
273 : : * @param angle Set to coordinates of surface normal nearest xyz
274 : : * @param history Optional ray history from a previous call to ray_fire().
275 : : * If present and non-empty, return the normal
276 : : * of the most recently intersected facet, ignoring xyz.
277 : : */
278 : : ErrorCode get_normal( EntityHandle surf, const double xyz[3], double angle[3], const RayHistory* history = NULL );
279 : :
280 : : private:
281 : : /**\brief determine the point membership when the point is effectively on the boundary
282 : : *
283 : : * Called by point_in_volume when the point is with tolerance of the boundary. Compares the
284 : : * ray direction with the surface normal to determine a volume membership.
285 : : */
286 : : ErrorCode boundary_case( EntityHandle volume, int& result, double u, double v, double w, EntityHandle facet,
287 : : EntityHandle surface );
288 : :
289 : : /** get the solid angle projected by a facet on a unit sphere around a point
290 : : * - used by point_in_volume_slow
291 : : */
292 : : ErrorCode poly_solid_angle( EntityHandle face, const CartVect& point, double& area );
293 : :
294 : : /**\brief State object used in calls to ray_fire()
295 : : *
296 : : * Storage for the "history" of a ray. This represents the surface facets
297 : : * that the ray is known to have crossed, which cannot be crossed again
298 : : * as long as the ray does not change direction. It is intended to be used
299 : : * with a series of consecutive calls to ray_fire(), in which a ray passes
300 : : * over potentially many surfaces.
301 : : */
302 : :
303 : : public:
304 : : /*
305 : : Overlap Thickness:
306 : : This tolerance is the maximum distance across an overlap. It should be zero
307 : : unless the geometry has overlaps. The overlap thickness is set using the dagmc
308 : : card. Overlaps must be small enough not to significantly affect physics.
309 : : Performance: increasing tolerance decreases performance
310 : : Robustness: increasing tolerance increases robustness
311 : : Knowledge: user must have intuition of overlap thickness
312 : : */
313 : :
314 : : /** Attempt to set a new overlap thickness tolerance, first checking for sanity */
315 : :
316 : : void set_overlap_thickness( double new_overlap_thickness );
317 : :
318 : : /*
319 : : Numerical Precision:
320 : : This tolerance is used for obb.intersect_ray, finding neighborhood of
321 : : adjacent triangle for edge/node intersections, and error in advancing
322 : : geometric position of particle (x' ~= x + d*u). When determining the
323 : : neighborhood of adjacent triangles for edge/node intersections, the facet
324 : : based model is expected to be watertight.
325 : : Performance: increasing tolerance decreases performance (but not very much)
326 : : Robustness: increasing tolerance increases robustness
327 : : Knowledge: user should not change this tolerance
328 : : */
329 : :
330 : : /** Attempt to set a new numerical precision , first checking for sanity
331 : : * Use of this function is discouraged.
332 : : */
333 : : void set_numerical_precision( double new_precision );
334 : :
335 : : double get_numerical_precision()
336 : : {
337 : : return numericalPrecision;
338 : : }
339 : :
340 : : double get_overlap_thickness()
341 : : {
342 : : return overlapThickness;
343 : : }
344 : :
345 : 108 : GeomTopoTool* gttool()
346 : : {
347 : 108 : return geomTopoTool;
348 : : }
349 : :
350 : 28 : Interface* moab_instance()
351 : : {
352 : 28 : return MBI;
353 : : }
354 : :
355 : : private:
356 : : GeomTopoTool* geomTopoTool;
357 : : bool owns_gtt;
358 : : Interface* MBI;
359 : : OrientedBoxTreeTool* obbTreeTool;
360 : : bool counting;
361 : : long long int n_pt_in_vol_calls;
362 : : long long int n_ray_fire_calls;
363 : : double overlapThickness, numericalPrecision;
364 : : Tag senseTag;
365 : : };
366 : :
367 : : } // namespace moab
368 : :
369 : : #endif
|
