MOAB: Mesh Oriented datABase  (version 5.4.1)
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,
00117                    bool find_geoments         = false,
00118                    EntityHandle modelRootSet  = 0,
00119                    bool p_rootSets_vector     = true,
00120                    bool restore_rootSets      = true,
00121                    bool trace_counting        = false,
00122                    double overlap_thickness   = 0.,
00123                    double numerical_precision = 0.001 );
00124
00125     GeomQueryTool( GeomTopoTool* geomtopotool,
00126                    bool trace_counting        = false,
00127                    double overlap_thickness   = 0.,
00128                    double numerical_precision = 0.001 );
00129
00130     // Destructor
00131     ~GeomQueryTool();
00132
00133     ErrorCode initialize();
00134
00135     /**\brief find the next surface crossing from a given point in a given direction
00136      *
00137      * This is the primary method to enable ray tracing through a geometry.
00138      * Given a volume and a ray, it determines the distance to the nearest intersection
00139      * with a bounding surface of that volume and returns that distance and the
00140      * EntityHandle indicating on which surface that intersection occurs.
00141      * The caller can compute the location of the intersection by adding the
00142      * distance to the ray.
00143      *
00144      * When a series of calls to this function are made along the same ray (e.g. for
00145      * the purpose of tracking a ray through several volumes), the optional history
00146      * argument should be given.  The history prevents previously intersected facets
00147      * from being intersected again.  A single history should be used as long as a
00148      * ray is proceeding forward without changing direction.  This situation is
00149      * sometimes referred to as "streaming."
00150      *
00151      * If a ray changes direction at an intersection site, the caller should call
00152      * reset_to_last_intersection() on the history object before the next ray fire.
00153      *
00154      * @param volume The volume at which to fire the ray
00155      * @param ray_start An array of x,y,z coordinates from which to start the ray.
00156      * @param ray_dir An array of x,y,z coordinates indicating the direction of the ray.
00157      *                Must be of unit length.
00158      * @param next_surf Output parameter indicating the next surface intersected by the ray.
00159      *                If no intersection is found, will be set to 0.
00160      * @param next_surf_dist Output parameter indicating distance to next_surf.  If next_surf is
00161      *                0, this value is undefined and should not be used.
00162      * @param history Optional RayHistory object.  If provided, the facets in the history are
00163      *                assumed to not intersect with the given ray.  The facet intersected
00164      *                by this query will also be added to the history.
00165      * @param dist_limit Optional distance limit.  If provided and > 0, no intersections at a
00166      *                distance further than this value will be returned.
00167      * @param ray_orientation Optional ray orientation. If provided determines intersections
00168      *                along the normal provided, e.g. if -1 allows intersections back along the
00169      *                the ray direction, Default is 1, i.e. exit intersections
00170      * @param stats Optional TrvStats object used to measure performance of underlying OBB
00171      *              ray-firing query.  See OrientedBoxTreeTool.hpp for details.
00172      *
00173      */
00174     ErrorCode ray_fire( const EntityHandle volume,
00175                         const double ray_start[3],
00176                         const double ray_dir[3],
00177                         EntityHandle& next_surf,
00178                         double& next_surf_dist,
00179                         RayHistory* history                  = NULL,
00180                         double dist_limit                    = 0,
00181                         int ray_orientation                  = 1,
00182                         OrientedBoxTreeTool::TrvStats* stats = NULL );
00183
00184     /**\brief Test if a point is inside or outside a volume
00185      *
00186      * This method finds the point on the boundary of the volume that is nearest
00187      * the test point (x,y,z).  If that point is "close" to a surface, a boundary test
00188      * is performed based on the normal of the surface at that point and the
00189      * optional ray direction (u,v,w).
00190      * @param volume The volume to test
00191      * @param xyz The location to test for volume containment
00192      * @param result Set to 0 if xyz it outside volume, 1 if inside, and -1 if on boundary.
00193      * @param Optional direction to use for underlying ray fire query.  Used to ensure
00194      *        consistent results when a ray direction is known.  If NULL or {0,0,0} is
00195      *        given, a random direction will be used.
00196      * @param history Optional RayHistory object to pass to underlying ray fire query.
00197      *        The history is not modified by this call.
00198      */
00199     ErrorCode point_in_volume( const EntityHandle volume,
00200                                const double xyz[3],
00201                                int& result,
00202                                const double* uvw         = NULL,
00203                                const RayHistory* history = NULL );
00204
00205     /**\brief Robust test if a point is inside or outside a volume using unit sphere area method
00206      *
00207      * This test may be more robust that the standard point_in_volume, but is much slower.
00208      * It does not detect 'on boundary' situations as point_in_volume does.
00209      * @param volume The volume to test
00210      * @param xyz The location to test for volume containment
00211      * @param result Set to 0 if xyz it outside volume, 1 if inside.
00212      */
00213     ErrorCode point_in_volume_slow( const EntityHandle volume, const double xyz[3], int& result );
00214
00215     /**\brief Find volume for a given location.
00216      *
00217      * Determines which volume contains the point if possible. If no volume is found,
00218      * a null EntityHandle is returned along with a MB_ENTITY_NOT_FOUND ErrorCode.
00219      * @param xyz The location to test
00220      * @param volume Set to volume handle containing the location if found
00221      * @param dir Optional direction to use for underlying ray fire query.  Used to ensure
00222      *        consistent results when a ray direction is known.  If NULL or {0,0,0} is
00223      *        given, a random direction will be used.
00224      */
00225     ErrorCode find_volume( const double xyz[3], EntityHandle& volume, const double* dir = NULL );
00226
00227     /**\brief Find volume for a given location using loop. (slow)
00228      *
00229      * Loops over all volumes in the model, checking for point containment
00230      * @param xyz The location to test
00231      * @param volume Set to volume handle containing the location if found
00232      * @param dir Optional direction to use for underlying ray fire query.  Used to ensure
00233      *        consistent results when a ray direction is known.  If NULL or {0,0,0} is
00234      *        given, a random direction will be used.
00235      */
00236     ErrorCode find_volume_slow( const double xyz[3], EntityHandle& volume, const double* dir = NULL );
00237
00238     /**\brief Test if a point is inside or outsize a volume's axis-aligned bounding box
00239      *
00240      * This is used as a fast way to determine whether a point is inside or outside of a volume
00241      * before performing a point_in_volume test which involves firing a ray.
00242      * @param volume The volume to test
00243      * @param point The location to test for bounding box containment
00244      * @param inside set to 0 if point is outside the box, 1 if inside
00245      */
00246     ErrorCode point_in_box( const EntityHandle volume, const double point[3], int& inside );
00247
00248     /** \brief Given a ray starting at a surface of a volume, check whether the ray enters or exits
00249      * the volume
00250      *
00251      * This function is most useful for rays that change directions at a surface crossing.
00252      * It can be used to check whether a direction change redirects the ray back into the
00253      * originating volume.
00254      *
00255      * @param volume The volume to test
00256      * @param surface A surface on volume
00257      * @param xyz A point location on surface
00258      * @param uvw A (unit) direction vector
00259      * @param result Set to 1 if ray is entering volume, or 0 if it is leaving
00260      * @param history Optional ray history object from a previous call to ray_fire.  If present and
00261      * non-empty, the history is used to look up the surface facet at which the ray begins.  Absent
00262      * a history, the facet nearest to xyz will be looked up.  The history should always be provided
00263      * if available, as it avoids the computational expense of a nearest-facet query.
00264      */
00265     ErrorCode test_volume_boundary( const EntityHandle volume,
00266                                     const EntityHandle surface,
00267                                     const double xyz[3],
00268                                     const double uvw[3],
00269                                     int& result,
00270                                     const RayHistory* history = NULL );
00271
00272     /**\brief Find the distance to the point on the boundary of the volume closest to the test point
00273      *
00274      * @param volume Volume to query
00275      * @param point Coordinates of test point
00276      * @param result Set to the minimum distance from point to a surface in volume
00277      */
00278     ErrorCode closest_to_location( EntityHandle volume,
00279                                    const double point[3],
00280                                    double& result,
00281                                    EntityHandle* closest_surface = 0 );
00282
00283     /** Calculate the volume contained in a 'volume' */
00284     ErrorCode measure_volume( EntityHandle volume, double& result );
00285
00286     /** Calculate sum of area of triangles */
00287     ErrorCode measure_area( EntityHandle surface, double& result );
00288
00289     /** Get the normal to a given surface at the point on the surface closest to a given point
00290      *
00291      * This method first identifies which facet contains this point and then
00292      * calculates the unit outward normal of that facet.  The facet of the
00293      * provided volume that is nearest the provided point is used for this
00294      * calculation.  The search for that facet can be circumvented by providing
00295      * a RayHistory, in which case the last facet of the history will be used.
00296      *
00297      * @param surf Surface on which to get normal
00298      * @param xyz Point on surf
00299      * @param angle Set to coordinates of surface normal nearest xyz
00300      * @param history Optional ray history from a previous call to ray_fire().
00301      *        If present and non-empty, return the normal
00302      *        of the most recently intersected facet, ignoring xyz.
00303      */
00304     ErrorCode get_normal( EntityHandle surf, const double xyz[3], double angle[3], const RayHistory* history = NULL );
00305
00306   private:
00307     /**\brief determine the point membership when the point is effectively on the boundary
00308      *
00309      * Called by point_in_volume when the point is with tolerance of the boundary. Compares the
00310      * ray direction with the surface normal to determine a volume membership.
00311      */
00312     ErrorCode boundary_case( EntityHandle volume,
00313                              int& result,
00314                              double u,
00315                              double v,
00316                              double w,
00317                              EntityHandle facet,
00318                              EntityHandle surface );
00319
00320     /** get the solid angle projected by a facet on a unit sphere around a point
00321      *  - used by point_in_volume_slow
00322      */
00323     ErrorCode poly_solid_angle( EntityHandle face, const CartVect& point, double& area );
00324
00325     /**\brief State object used in calls to ray_fire()
00326      *
00327      * Storage for the "history" of a ray.  This represents the surface facets
00328      * that the ray is known to have crossed, which cannot be crossed again
00329      * as long as the ray does not change direction.  It is intended to be used
00330      * with a series of consecutive calls to ray_fire(), in which a ray passes
00331      * over potentially many surfaces.
00332      */
00333
00334   public:
00335     /*
00336      Overlap Thickness:
00337      This tolerance is the maximum distance across an overlap. It should be zero
00338      unless the geometry has overlaps. The overlap thickness is set using the dagmc
00339      card. Overlaps must be small enough not to significantly affect physics.
00340        Performance: increasing tolerance decreases performance
00341        Robustness:  increasing tolerance increases robustness
00342        Knowledge:   user must have intuition of overlap thickness
00343     */
00344
00345     /** Attempt to set a new overlap thickness tolerance, first checking for sanity */
00346
00347     void set_overlap_thickness( double new_overlap_thickness );
00348
00349     /*
00350      Numerical Precision:
00351      This tolerance is used for obb.intersect_ray, finding neighborhood of
00353      geometric position of particle (x' ~= x + d*u). When determining the
00354      neighborhood of adjacent triangles for edge/node intersections, the facet
00355      based model is expected to be watertight.
00356        Performance: increasing tolerance decreases performance (but not very much)
00357        Robustness:  increasing tolerance increases robustness
00358        Knowledge:   user should not change this tolerance
00359     */
00360
00361     /** Attempt to set a new numerical precision , first checking for sanity
00362      *  Use of this function is discouraged.
00363      */
00364     void set_numerical_precision( double new_precision );
00365
00366     double get_numerical_precision()
00367     {
00368         return numericalPrecision;
00369     }
00370
00371     double get_overlap_thickness()
00372     {
00373         return overlapThickness;
00374     }
00375
00376     GeomTopoTool* gttool()
00377     {
00378         return geomTopoTool;
00379     }
00380
00381     Interface* moab_instance()
00382     {
00383         return MBI;
00384     }
00385
00386   private:
00387     GeomTopoTool* geomTopoTool;
00388     bool owns_gtt;
00389     Interface* MBI;
00390     OrientedBoxTreeTool* obbTreeTool;
00391     bool counting;
00392     long long int n_pt_in_vol_calls;
00393     long long int n_ray_fire_calls;
00394     double overlapThickness, numericalPrecision;
00395     Tag senseTag;
00396 };
00397
00398 }  // namespace moab
00399
00400 #endif