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