MOAB: Mesh Oriented datABase  (version 5.3.0)
GeomQueryTool.hpp
Go to the documentation of this file.
00001 #ifndef MOAB_GEOM_QUERY_TOOL_HPP
00002 #define MOAB_GEOM_QUERY_TOOL_HPP
00003 
00004 #ifdef _MSC_VER            /* windows */
00005 #define _USE_MATH_DEFINES  // For M_PI
00006 #endif
00007 
00008 #include "MBTagConventions.hpp"
00009 #include "moab/CartVect.hpp"
00010 #include "moab/Range.hpp"
00011 #include "moab/Core.hpp"
00012 #include "moab/GeomUtil.hpp"
00013 #include "moab/FileOptions.hpp"
00014 #include "moab/EntityHandle.hpp"
00015 #include "moab/GeomTopoTool.hpp"
00016 #include "moab/OrientedBoxTreeTool.hpp"
00017 
00018 #include <vector>
00019 #include <map>
00020 #include <string>
00021 #include <cassert>
00022 
00023 namespace moab
00024 {
00025 
00026 /** \class GeomQueryTool
00027  *
00028  * \brief Tool for querying different aspects of geometric topology sets in MOAB
00029  *
00030  * Given the conventions established in GeomTopoTool for representing
00031  * geometric topology through a hierarchy of meshsets, this tool provides a
00032  * set of methods to query different geometric aspects of those geometric
00033  * topology sets including:
00034  *
00035  *   * measures of surface area and volume
00036  *   * distance to a bounding surface from a point within a volume
00037  *   * test the inclusion of a point within a volume
00038  *   * find the angle of a surface at a point
00039  *
00040  * A feature of these queries is that there is some support for overlapping
00041  * geometries.
00042  *
00043  */
00044 
00045 class GeomQueryTool
00046 {
00047   public:
00048     /* \class RayHistory
00049      *
00050      * In many use cases, it is useful to track some of the history of a ray as
00051      * it passes through a geometry, particularly a geometry represented by
00052      * facets.  For example, given round-off erorr in ray-triangle tests
00053      * (GeomUtil::ray_tri_intersect) used as part of a test for ray-surface
00054      * intersection, it is possible for subsequent queries to oscillate between
00055      * adjacent surfaces.  This class stores information about history of a ray
00056      * that can be used to test for such circumstances so that they can be
00057      * accommodated.
00058      *
00059      */
00060 
00061     class RayHistory
00062     {
00063 
00064       public:
00065         /**
00066          * Clear this entire history-- logically equivalent to creating a new history,
00067          * but probably more efficient.
00068          */
00069         void reset();
00070 
00071         /**
00072          * Clear the history up to the most recent intersection.  This should be
00073          * called when a ray changes direction at the site of a surface crossing,
00074          * a situation that most commonly occurs at a reflecting boundary.
00075          */
00076         void reset_to_last_intersection();
00077 
00078         /**
00079          * Remove the most recent intersection.  This allows a subsequent call
00080          * along the same ray to return the same intersection.
00081          */
00082         void rollback_last_intersection();
00083 
00084         /**
00085          * Get the last intersection in the RayHistory. This will return a null
00086          * EntityHandle (0) if the history is empty.
00087          */
00088         ErrorCode get_last_intersection( EntityHandle& last_facet_hit ) const;
00089 
00090         /**
00091          * @return the number of surface crossings currently represented by this ray history
00092          */
00093         int size() const
00094         {
00095             return prev_facets.size();
00096         }
00097 
00098         /**
00099          * @return Boolean indicating if this entity is in the RayHistory
00100          */
00101         bool in_history( EntityHandle ent ) const;
00102 
00103         /**
00104          * Add entity to the RayHistory
00105          */
00106         void add_entity( EntityHandle ent );
00107 
00108       private:
00109         std::vector< EntityHandle > prev_facets;
00110 
00111         friend class GeomQueryTool;
00112     };
00113 
00114     // Constructors
00115 
00116     GeomQueryTool( Interface* impl, bool find_geoments = false, EntityHandle modelRootSet = 0,
00117                    bool p_rootSets_vector = true, bool restore_rootSets = true, bool trace_counting = false,
00118                    double overlap_thickness = 0., double numerical_precision = 0.001 );
00119 
00120     GeomQueryTool( GeomTopoTool* geomtopotool, bool trace_counting = false, double overlap_thickness = 0.,
00121                    double numerical_precision = 0.001 );
00122 
00123     // Destructor
00124     ~GeomQueryTool();
00125 
00126     ErrorCode initialize();
00127 
00128     /**\brief find the next surface crossing from a given point in a given direction
00129      *
00130      * This is the primary method to enable ray tracing through a geometry.
00131      * Given a volume and a ray, it determines the distance to the nearest intersection
00132      * with a bounding surface of that volume and returns that distance and the
00133      * EntityHandle indicating on which surface that intersection occurs.
00134      * The caller can compute the location of the intersection by adding the
00135      * distance to the ray.
00136      *
00137      * When a series of calls to this function are made along the same ray (e.g. for
00138      * the purpose of tracking a ray through several volumes), the optional history
00139      * argument should be given.  The history prevents previously intersected facets
00140      * from being intersected again.  A single history should be used as long as a
00141      * ray is proceeding forward without changing direction.  This situation is
00142      * sometimes referred to as "streaming."
00143      *
00144      * If a ray changes direction at an intersection site, the caller should call
00145      * reset_to_last_intersection() on the history object before the next ray fire.
00146      *
00147      * @param volume The volume at which to fire the ray
00148      * @param ray_start An array of x,y,z coordinates from which to start the ray.
00149      * @param ray_dir An array of x,y,z coordinates indicating the direction of the ray.
00150      *                Must be of unit length.
00151      * @param next_surf Output parameter indicating the next surface intersected by the ray.
00152      *                If no intersection is found, will be set to 0.
00153      * @param next_surf_dist Output parameter indicating distance to next_surf.  If next_surf is
00154      *                0, this value is undefined and should not be used.
00155      * @param history Optional RayHistory object.  If provided, the facets in the history are
00156      *                assumed to not intersect with the given ray.  The facet intersected
00157      *                by this query will also be added to the history.
00158      * @param dist_limit Optional distance limit.  If provided and > 0, no intersections at a
00159      *                distance further than this value will be returned.
00160      * @param ray_orientation Optional ray orientation. If provided determines intersections
00161      *                along the normal provided, e.g. if -1 allows intersections back along the
00162      *                the ray direction, Default is 1, i.e. exit intersections
00163      * @param stats Optional TrvStats object used to measure performance of underlying OBB
00164      *              ray-firing query.  See OrientedBoxTreeTool.hpp for details.
00165      *
00166      */
00167     ErrorCode ray_fire( const EntityHandle volume, const double ray_start[3], const double ray_dir[3],
00168                         EntityHandle& next_surf, double& next_surf_dist, RayHistory* history = NULL,
00169                         double dist_limit = 0, int ray_orientation = 1, OrientedBoxTreeTool::TrvStats* stats = NULL );
00170 
00171     /**\brief Test if a point is inside or outside a volume
00172      *
00173      * This method finds the point on the boundary of the volume that is nearest
00174      * the test point (x,y,z).  If that point is "close" to a surface, a boundary test
00175      * is performed based on the normal of the surface at that point and the
00176      * optional ray direction (u,v,w).
00177      * @param volume The volume to test
00178      * @param xyz The location to test for volume containment
00179      * @param result Set to 0 if xyz it outside volume, 1 if inside, and -1 if on boundary.
00180      * @param Optional direction to use for underlying ray fire query.  Used to ensure
00181      *        consistent results when a ray direction is known.  If NULL or {0,0,0} is
00182      *        given, a random direction will be used.
00183      * @param history Optional RayHistory object to pass to underlying ray fire query.
00184      *        The history is not modified by this call.
00185      */
00186     ErrorCode point_in_volume( const EntityHandle volume, const double xyz[3], int& result, const double* uvw = NULL,
00187                                const RayHistory* history = NULL );
00188 
00189     /**\brief Robust test if a point is inside or outside a volume using unit sphere area method
00190      *
00191      * This test may be more robust that the standard point_in_volume, but is much slower.
00192      * It does not detect 'on boundary' situations as point_in_volume does.
00193      * @param volume The volume to test
00194      * @param xyz The location to test for volume containment
00195      * @param result Set to 0 if xyz it outside volume, 1 if inside.
00196      */
00197     ErrorCode point_in_volume_slow( const EntityHandle volume, const double xyz[3], int& result );
00198 
00199     /**\brief Find volume for a given location.
00200      *
00201      * Determines which volume contains the point if possible. If no volume is found,
00202      * a null EntityHandle is returned along with a MB_ENTITY_NOT_FOUND ErrorCode.
00203      * @param xyz The location to test
00204      * @param volume Set to volume handle containing the location if found
00205      * @param dir Optional direction to use for underlying ray fire query.  Used to ensure
00206      *        consistent results when a ray direction is known.  If NULL or {0,0,0} is
00207      *        given, a random direction will be used.
00208      */
00209     ErrorCode find_volume( const double xyz[3], EntityHandle& volume, const double* dir = NULL );
00210 
00211     /**\brief Find volume for a given location using loop. (slow)
00212      *
00213      * Loops over all volumes in the model, checking for point containment
00214      * @param xyz The location to test
00215      * @param volume Set to volume handle containing the location if found
00216      * @param dir Optional direction to use for underlying ray fire query.  Used to ensure
00217      *        consistent results when a ray direction is known.  If NULL or {0,0,0} is
00218      *        given, a random direction will be used.
00219      */
00220     ErrorCode find_volume_slow( const double xyz[3], EntityHandle& volume, const double* dir = NULL );
00221 
00222     /**\brief Test if a point is inside or outsize a volume's axis-aligned bounding box
00223      *
00224      * This is used as a fast way to determine whether a point is inside or outside of a volume
00225      * before performing a point_in_volume test which involves firing a ray.
00226      * @param volume The volume to test
00227      * @param point The location to test for bounding box containment
00228      * @param inside set to 0 if point is outside the box, 1 if inside
00229      */
00230     ErrorCode point_in_box( const EntityHandle volume, const double point[3], int& inside );
00231 
00232     /** \brief Given a ray starting at a surface of a volume, check whether the ray enters or exits
00233      * the volume
00234      *
00235      * This function is most useful for rays that change directions at a surface crossing.
00236      * It can be used to check whether a direction change redirects the ray back into the
00237      * originating volume.
00238      *
00239      * @param volume The volume to test
00240      * @param surface A surface on volume
00241      * @param xyz A point location on surface
00242      * @param uvw A (unit) direction vector
00243      * @param result Set to 1 if ray is entering volume, or 0 if it is leaving
00244      * @param history Optional ray history object from a previous call to ray_fire.  If present and
00245      * non-empty, the history is used to look up the surface facet at which the ray begins.  Absent
00246      * a history, the facet nearest to xyz will be looked up.  The history should always be provided
00247      * if available, as it avoids the computational expense of a nearest-facet query.
00248      */
00249     ErrorCode test_volume_boundary( const EntityHandle volume, const EntityHandle surface, const double xyz[3],
00250                                     const double uvw[3], int& result, const RayHistory* history = NULL );
00251 
00252     /**\brief Find the distance to the point on the boundary of the volume closest to the test point
00253      *
00254      * @param volume Volume to query
00255      * @param point Coordinates of test point
00256      * @param result Set to the minimum distance from point to a surface in volume
00257      */
00258     ErrorCode closest_to_location( EntityHandle volume, const double point[3], double& result,
00259                                    EntityHandle* closest_surface = 0 );
00260 
00261     /** Calculate the volume contained in a 'volume' */
00262     ErrorCode measure_volume( EntityHandle volume, double& result );
00263 
00264     /** Calculate sum of area of triangles */
00265     ErrorCode measure_area( EntityHandle surface, double& result );
00266 
00267     /** Get the normal to a given surface at the point on the surface closest to a given point
00268      *
00269      * This method first identifies which facet contains this point and then
00270      * calculates the unit outward normal of that facet.  The facet of the
00271      * provided volume that is nearest the provided point is used for this
00272      * calculation.  The search for that facet can be circumvented by providing
00273      * a RayHistory, in which case the last facet of the history will be used.
00274      *
00275      * @param surf Surface on which to get normal
00276      * @param xyz Point on surf
00277      * @param angle Set to coordinates of surface normal nearest xyz
00278      * @param history Optional ray history from a previous call to ray_fire().
00279      *        If present and non-empty, return the normal
00280      *        of the most recently intersected facet, ignoring xyz.
00281      */
00282     ErrorCode get_normal( EntityHandle surf, const double xyz[3], double angle[3], const RayHistory* history = NULL );
00283 
00284   private:
00285     /**\brief determine the point membership when the point is effectively on the boundary
00286      *
00287      * Called by point_in_volume when the point is with tolerance of the boundary. Compares the
00288      * ray direction with the surface normal to determine a volume membership.
00289      */
00290     ErrorCode boundary_case( EntityHandle volume, int& result, double u, double v, double w, EntityHandle facet,
00291                              EntityHandle surface );
00292 
00293     /** get the solid angle projected by a facet on a unit sphere around a point
00294      *  - used by point_in_volume_slow
00295      */
00296     ErrorCode poly_solid_angle( EntityHandle face, const CartVect& point, double& area );
00297 
00298     /**\brief State object used in calls to ray_fire()
00299      *
00300      * Storage for the "history" of a ray.  This represents the surface facets
00301      * that the ray is known to have crossed, which cannot be crossed again
00302      * as long as the ray does not change direction.  It is intended to be used
00303      * with a series of consecutive calls to ray_fire(), in which a ray passes
00304      * over potentially many surfaces.
00305      */
00306 
00307   public:
00308     /*
00309      Overlap Thickness:
00310      This tolerance is the maximum distance across an overlap. It should be zero
00311      unless the geometry has overlaps. The overlap thickness is set using the dagmc
00312      card. Overlaps must be small enough not to significantly affect physics.
00313        Performance: increasing tolerance decreases performance
00314        Robustness:  increasing tolerance increases robustness
00315        Knowledge:   user must have intuition of overlap thickness
00316     */
00317 
00318     /** Attempt to set a new overlap thickness tolerance, first checking for sanity */
00319 
00320     void set_overlap_thickness( double new_overlap_thickness );
00321 
00322     /*
00323      Numerical Precision:
00324      This tolerance is used for obb.intersect_ray, finding neighborhood of
00325      adjacent triangle for edge/node intersections, and error in advancing
00326      geometric position of particle (x' ~= x + d*u). When determining the
00327      neighborhood of adjacent triangles for edge/node intersections, the facet
00328      based model is expected to be watertight.
00329        Performance: increasing tolerance decreases performance (but not very much)
00330        Robustness:  increasing tolerance increases robustness
00331        Knowledge:   user should not change this tolerance
00332     */
00333 
00334     /** Attempt to set a new numerical precision , first checking for sanity
00335      *  Use of this function is discouraged.
00336      */
00337     void set_numerical_precision( double new_precision );
00338 
00339     double get_numerical_precision()
00340     {
00341         return numericalPrecision;
00342     }
00343 
00344     double get_overlap_thickness()
00345     {
00346         return overlapThickness;
00347     }
00348 
00349     GeomTopoTool* gttool()
00350     {
00351         return geomTopoTool;
00352     }
00353 
00354     Interface* moab_instance()
00355     {
00356         return MBI;
00357     }
00358 
00359   private:
00360     GeomTopoTool* geomTopoTool;
00361     bool owns_gtt;
00362     Interface* MBI;
00363     OrientedBoxTreeTool* obbTreeTool;
00364     bool counting;
00365     long long int n_pt_in_vol_calls;
00366     long long int n_ray_fire_calls;
00367     double overlapThickness, numericalPrecision;
00368     Tag senseTag;
00369 };
00370 
00371 }  // namespace moab
00372 
00373 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines