MOAB: Mesh Oriented datABase  (version 5.3.0)
GeomQueryTool.cpp
Go to the documentation of this file.
00001 #include "moab/GeomQueryTool.hpp"
00002 
00003 #ifdef WIN32               /* windows */
00004 #define _USE_MATH_DEFINES  // For M_PI
00005 #endif
00006 #include <string>
00007 #include <iostream>
00008 #include <fstream>
00009 #include <sstream>
00010 #include <limits>
00011 #include <algorithm>
00012 #include <set>
00013 
00014 #include <cctype>
00015 #include <cstring>
00016 #include <cstdlib>
00017 #include <cstdio>
00018 
00019 #include "moab/OrientedBoxTreeTool.hpp"
00020 
00021 const bool debug = false;
00022 #ifdef __DEBUG
00023 debug = true;
00024 #endif
00025 
00026 namespace moab
00027 {
00028 
00029 /** \class FindVolume_IntRegCtxt
00030  *
00031  *
00032  *\brief An intersection context used for finding a volume
00033  *
00034  * This context is used to find the nearest intersection location
00035  * and is intended for use with a global surface tree from
00036  * the GeomTopoTool.
00037  *
00038  * The behavior of this context is relatively simple in that it
00039  * returns only one intersection distance, surface, and facet.
00040  * Intersections of any orientation are accepted. The positive
00041  * value of the search window is limited to the current nearest
00042  * intersection distance.
00043  *
00044  */
00045 
00046 class FindVolumeIntRegCtxt : public OrientedBoxTreeTool::IntRegCtxt
00047 {
00048 
00049   public:
00050     // Constructor
00051     FindVolumeIntRegCtxt()
00052     {
00053         // initialize return vectors
00054         // only one hit is returned in this context
00055         intersections.push_back( std::numeric_limits< double >::max() );
00056         sets.push_back( 0 );
00057         facets.push_back( 0 );
00058     }
00059 
00060     ErrorCode register_intersection( EntityHandle set, EntityHandle tri, double dist,
00061                                      OrientedBoxTreeTool::IntersectSearchWindow& search_win,
00062                                      GeomUtil::intersection_type /*it*/ )
00063     {
00064         // update dist, set, and triangle hit if
00065         // we found a new minimum distance
00066         double abs_dist = fabs( dist );
00067         if( abs_dist < fabs( intersections[0] ) )
00068         {
00069             intersections[0] = dist;
00070             sets[0]          = set;
00071             facets[0]        = tri;
00072 
00073             // narrow search window based on the hit distance
00074             pos               = abs_dist;
00075             neg               = -abs_dist;
00076             search_win.first  = &pos;
00077             search_win.second = &neg;
00078         }
00079 
00080         return MB_SUCCESS;
00081     }
00082 
00083     // storage for updated window values during search
00084     double pos;
00085     double neg;
00086 };
00087 
00088 /** \class GQT_IntRegCtxt
00089  *
00090  *\brief An implementation of an Intersection Registration Context for use GQT ray-firing
00091  *
00092  * This context uses a variety of tests and conditions to confirm whether or
00093  * not to accumulate an intersection, to ensure robustness for ray firing.
00094  *
00095  * This context only accumulates intersections that are oriented parallel to
00096  * the 'desiredOrient', if provided, with respect to 'geomVol', using
00097  * information in the in 'senseTag'.
00098  *
00099  * This context only accumulates a single intersection out of a set of
00100  * multiple intersections that fall in the same 'neighborhood', where a
00101  * 'neighborhood' is defined as facets that share edges or vertices.
00102  *
00103  * This context only accumulates piercing intersections.  This is relevant
00104  * for intersections that are found to be on an edge or vertex by the
00105  * Plucker test.  Such intersections are piercing if the ray has the same
00106  * orientation w.r.t. to all fecets that share that edge or vertex.
00107  *
00108  * This context tests intersections against a list of 'prevFacets' to
00109  * prevent a ray from crossing the same facet more than once.  The user is
00110  * responsible for ensuring that this list is reset when appropriate.
00111  *
00112  * This context accumulates all intersections within 'tol' of the
00113  * start of the ray and if the number of intersections within the
00114  * 'tol' of the ray start point is less than 'minTolInt', the next
00115  * closest intersection. If the desired result is only the closest
00116  * intersection, 'minTolInt' should be 0.  This function will return all
00117  * intersections, regardless of distance from the start of the ray, if
00118  * 'minTolInt' is negative.
00119  *
00120  */
00121 
00122 class GQT_IntRegCtxt : public OrientedBoxTreeTool::IntRegCtxt
00123 {
00124 
00125   private:
00126     // Input
00127     OrientedBoxTreeTool* tool;
00128     const CartVect ray_origin;
00129     const CartVect ray_direction;
00130     const double tol; /* used for box.intersect_ray, radius of
00131                          neighborhood for adjacent triangles,
00132                          and old mode of add_intersection */
00133     const int minTolInt;
00134 
00135     // Optional Input - to screen RTIs by orientation and edge/node intersection
00136     const EntityHandle* rootSet; /* used for sphere_intersect */
00137     const EntityHandle* geomVol; /* used for determining surface sense */
00138     const Tag* senseTag;         /* allows screening by triangle orientation.
00139                                     both geomVol and senseTag must be used together. */
00140     const int* desiredOrient;    /* points to desired orientation of ray with
00141                                     respect to surf normal, if this feature is used.
00142                                     Must point to -1 (reverse) or 1 (forward).
00143                                     geomVol and senseTag are needed for this feature */
00144 
00145     // Optional Input - to avoid returning these as RTIs
00146     const std::vector< EntityHandle >* prevFacets; /* intersections on these triangles
00147                                                       will not be returned */
00148 
00149     // Other Variables
00150     std::vector< std::vector< EntityHandle > > neighborhoods;
00151     std::vector< EntityHandle > neighborhood;
00152 
00153     void add_intersection( EntityHandle set, EntityHandle tri, double dist,
00154                            OrientedBoxTreeTool::IntersectSearchWindow& search_win );
00155     void append_intersection( EntityHandle set, EntityHandle facet, double dist );
00156     void set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist );
00157     void add_mode1_intersection( EntityHandle set, EntityHandle facet, double dist,
00158                                  OrientedBoxTreeTool::IntersectSearchWindow& search_win );
00159     bool edge_node_piercing_intersect( const EntityHandle tri, const CartVect& ray_direction,
00160                                        const GeomUtil::intersection_type int_type,
00161                                        const std::vector< EntityHandle >& close_tris,
00162                                        const std::vector< int >& close_senses, const Interface* MBI,
00163                                        std::vector< EntityHandle >* neighborhood_tris = 0 );
00164 
00165     bool in_prevFacets( const EntityHandle tri );
00166     bool in_neighborhoods( const EntityHandle tri );
00167 
00168   public:
00169     GQT_IntRegCtxt( OrientedBoxTreeTool* obbtool, const double ray_point[3], const double ray_dir[3], double tolerance,
00170                     int min_tolerance_intersections, const EntityHandle* root_set, const EntityHandle* geom_volume,
00171                     const Tag* sense_tag, const int* desired_orient, const std::vector< EntityHandle >* prev_facets )
00172         : tool( obbtool ), ray_origin( ray_point ), ray_direction( ray_dir ), tol( tolerance ),
00173           minTolInt( min_tolerance_intersections ), rootSet( root_set ), geomVol( geom_volume ), senseTag( sense_tag ),
00174           desiredOrient( desired_orient ), prevFacets( prev_facets ){
00175 
00176           };
00177 
00178     virtual ErrorCode register_intersection( EntityHandle set, EntityHandle triangle, double distance,
00179                                              OrientedBoxTreeTool::IntersectSearchWindow&,
00180                                              GeomUtil::intersection_type int_type );
00181 
00182     virtual ErrorCode update_orient( EntityHandle set, int* surfTriOrient );
00183     virtual const int* getDesiredOrient()
00184     {
00185         return desiredOrient;
00186     };
00187 };
00188 
00189 ErrorCode GQT_IntRegCtxt::update_orient( EntityHandle set, int* surfTriOrient )
00190 {
00191 
00192     ErrorCode rval;
00193 
00194     // Get desired orientation of surface wrt volume. Use this to return only
00195     // exit or entrance intersections.
00196     if( geomVol && senseTag && desiredOrient && surfTriOrient )
00197     {
00198         if( 1 != *desiredOrient && -1 != *desiredOrient )
00199         {
00200             std::cerr << "error: desired orientation must be 1 (forward) or -1 (reverse)" << std::endl;
00201         }
00202         EntityHandle vols[2];
00203         rval = tool->get_moab_instance()->tag_get_data( *senseTag, &set, 1, vols );
00204         assert( MB_SUCCESS == rval );
00205         if( MB_SUCCESS != rval ) return rval;
00206         if( vols[0] == vols[1] )
00207         {
00208             std::cerr << "error: surface has positive and negative sense wrt same volume" << std::endl;
00209             return MB_FAILURE;
00210         }
00211         // surfTriOrient will be used by plucker_ray_tri_intersect to avoid
00212         // intersections with wrong orientation.
00213         if( *geomVol == vols[0] ) { *surfTriOrient = *desiredOrient * 1; }
00214         else if( *geomVol == vols[1] )
00215         {
00216             *surfTriOrient = *desiredOrient * ( -1 );
00217         }
00218         else
00219         {
00220             assert( false );
00221             return MB_FAILURE;
00222         }
00223     }
00224 
00225     return MB_SUCCESS;
00226 }
00227 
00228 bool GQT_IntRegCtxt::in_prevFacets( const EntityHandle tri )
00229 {
00230     return ( prevFacets && ( ( *prevFacets ).end() != find( ( *prevFacets ).begin(), ( *prevFacets ).end(), tri ) ) );
00231 }
00232 
00233 bool GQT_IntRegCtxt::in_neighborhoods( const EntityHandle tri )
00234 {
00235     bool same_neighborhood = false;
00236     for( unsigned i = 0; i < neighborhoods.size(); ++i )
00237     {
00238         if( neighborhoods[i].end() != find( neighborhoods[i].begin(), neighborhoods[i].end(), tri ) )
00239         {
00240             same_neighborhood = true;
00241             continue;
00242         }
00243     }
00244     return same_neighborhood;
00245 }
00246 
00247 /**\brief Determine if a ray-edge/node intersection is glancing or piercing.
00248  *        This function avoids asking for upward adjacencies to prevent their
00249  *        creation.
00250  *\param tri           The intersected triangle
00251  *\param ray_dir       The direction of the ray
00252  *\param int_type      The type of intersection (EDGE0, EDGE1, NODE2, ...)
00253  *\param close_tris    Vector of triangles in the proximity of the intersection
00254  *\param close_senses  Vector of surface senses for tris in the proximity of
00255  *                     the intersection
00256  *\param neighborhood  Vector of triangles in the topological neighborhood of the intersection
00257  *\return              True if piercing, false otherwise.
00258  */
00259 bool GQT_IntRegCtxt::edge_node_piercing_intersect( const EntityHandle tri, const CartVect& ray_dir,
00260                                                    const GeomUtil::intersection_type int_type,
00261                                                    const std::vector< EntityHandle >& close_tris,
00262                                                    const std::vector< int >& close_senses, const Interface* MBI,
00263                                                    std::vector< EntityHandle >* neighborhood_tris )
00264 {
00265 
00266     // get the node of the triangle
00267     const EntityHandle* conn = NULL;
00268     int len                  = 0;
00269     ErrorCode rval           = MBI->get_connectivity( tri, conn, len );
00270     if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE;
00271 
00272     // get adjacent tris (and keep their corresponding senses)
00273     std::vector< EntityHandle > adj_tris;
00274     std::vector< int > adj_senses;
00275 
00276     // node intersection
00277     if( GeomUtil::NODE0 == int_type || GeomUtil::NODE1 == int_type || GeomUtil::NODE2 == int_type )
00278     {
00279 
00280         // get the intersected node
00281         EntityHandle node;
00282         if( GeomUtil::NODE0 == int_type )
00283             node = conn[0];
00284         else if( GeomUtil::NODE1 == int_type )
00285             node = conn[1];
00286         else
00287             node = conn[2];
00288 
00289         // get tris adjacent to node
00290         for( unsigned i = 0; i < close_tris.size(); ++i )
00291         {
00292             const EntityHandle* con = NULL;
00293             rval                    = MBI->get_connectivity( close_tris[i], con, len );
00294             if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE;
00295 
00296             if( node == con[0] || node == con[1] || node == con[2] )
00297             {
00298                 adj_tris.push_back( close_tris[i] );
00299                 adj_senses.push_back( close_senses[i] );
00300             }
00301         }
00302         if( adj_tris.empty() )
00303         {
00304             std::cerr << "error: no tris are adjacent to the node" << std::endl;
00305             return MB_FAILURE;
00306         }
00307         // edge intersection
00308     }
00309     else if( GeomUtil::EDGE0 == int_type || GeomUtil::EDGE1 == int_type || GeomUtil::EDGE2 == int_type )
00310     {
00311 
00312         // get the endpoints of the edge
00313         EntityHandle endpts[2];
00314         if( GeomUtil::EDGE0 == int_type )
00315         {
00316             endpts[0] = conn[0];
00317             endpts[1] = conn[1];
00318         }
00319         else if( GeomUtil::EDGE1 == int_type )
00320         {
00321             endpts[0] = conn[1];
00322             endpts[1] = conn[2];
00323         }
00324         else
00325         {
00326             endpts[0] = conn[2];
00327             endpts[1] = conn[0];
00328         }
00329 
00330         // get tris adjacent to edge
00331         for( unsigned i = 0; i < close_tris.size(); ++i )
00332         {
00333             const EntityHandle* con = NULL;
00334             rval                    = MBI->get_connectivity( close_tris[i], con, len );
00335             if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE;
00336 
00337             // check both orientations in case close_tris are not on the same surface
00338             if( ( endpts[0] == con[0] && endpts[1] == con[1] ) || ( endpts[0] == con[1] && endpts[1] == con[0] ) ||
00339                 ( endpts[0] == con[1] && endpts[1] == con[2] ) || ( endpts[0] == con[2] && endpts[1] == con[1] ) ||
00340                 ( endpts[0] == con[2] && endpts[1] == con[0] ) || ( endpts[0] == con[0] && endpts[1] == con[2] ) )
00341             {
00342                 adj_tris.push_back( close_tris[i] );
00343                 adj_senses.push_back( close_senses[i] );
00344             }
00345         }
00346         // In a 2-manifold each edge is adjacent to exactly 2 tris
00347         if( 2 != adj_tris.size() )
00348         {
00349             std::cerr << "error: edge of a manifold must be topologically adjacent to exactly 2 tris" << std::endl;
00350             MBI->list_entities( endpts, 2 );
00351             return true;
00352         }
00353     }
00354     else
00355     {
00356         std::cerr << "error: special case not an node/edge intersection" << std::endl;
00357         return MB_FAILURE;
00358     }
00359 
00360     // The close tris were in proximity to the intersection. The adj_tris are
00361     // topologically adjacent to the intersection (the neighborhood).
00362     if( neighborhood_tris ) ( *neighborhood_tris ).assign( adj_tris.begin(), adj_tris.end() );
00363 
00364     // determine glancing/piercing
00365     // If a desired_orientation was used in this call to ray_intersect_sets,
00366     // the plucker_ray_tri_intersect will have already used it. For a piercing
00367     // intersection, the normal of all tris must have the same orientation.
00368     int sign = 0;
00369     for( unsigned i = 0; i < adj_tris.size(); ++i )
00370     {
00371         const EntityHandle* con = NULL;
00372         rval                    = MBI->get_connectivity( adj_tris[i], con, len );
00373         if( MB_SUCCESS != rval || 3 != len ) return MB_FAILURE;
00374         CartVect coords[3];
00375         rval = MBI->get_coords( con, len, coords[0].array() );
00376         if( MB_SUCCESS != rval ) return MB_FAILURE;
00377 
00378         // get normal of triangle
00379         CartVect v0     = coords[1] - coords[0];
00380         CartVect v1     = coords[2] - coords[0];
00381         CartVect norm   = adj_senses[i] * ( v0 * v1 );
00382         double dot_prod = norm % ray_dir;
00383 
00384         // if the sign has not yet been decided, choose it
00385         if( 0 == sign && 0 != dot_prod )
00386         {
00387             if( 0 < dot_prod )
00388                 sign = 1;
00389             else
00390                 sign = -1;
00391         }
00392 
00393         // intersection is glancing if tri and ray do not point in same direction
00394         // for every triangle
00395         if( 0 != sign && 0 > sign * dot_prod ) return false;
00396     }
00397     return true;
00398 }
00399 
00400 ErrorCode GQT_IntRegCtxt::register_intersection( EntityHandle set, EntityHandle t, double int_dist,
00401                                                  OrientedBoxTreeTool::IntersectSearchWindow& search_win,
00402                                                  GeomUtil::intersection_type int_type )
00403 {
00404     ErrorCode rval;
00405 
00406     // Do not accept intersections if they are in the vector of previously intersected
00407     // facets.
00408     if( in_prevFacets( t ) ) return MB_SUCCESS;
00409 
00410     // Do not accept intersections if they are in the neighborhood of previous
00411     // intersections.
00412     if( in_neighborhoods( t ) ) return MB_SUCCESS;
00413 
00414     neighborhood.clear();
00415 
00416     // Handle special case of edge/node intersection. Accept piercing
00417     // intersections and reject glancing intersections.
00418     // The edge_node_intersection function needs to know surface sense wrt volume.
00419     // A less-robust implementation could work without sense information.
00420     // Would it ever be useful to accept a glancing intersection?
00421     if( GeomUtil::INTERIOR != int_type && rootSet && geomVol && senseTag )
00422     {
00423         // get triangles in the proximity of the intersection
00424         std::vector< EntityHandle > close_tris;
00425         std::vector< int > close_senses;
00426         rval = tool->get_close_tris( ray_origin + int_dist * ray_direction, tol, rootSet, geomVol, senseTag, close_tris,
00427                                      close_senses );
00428 
00429         if( MB_SUCCESS != rval ) return rval;
00430 
00431         if( !edge_node_piercing_intersect( t, ray_direction, int_type, close_tris, close_senses,
00432                                            tool->get_moab_instance(), &neighborhood ) )
00433             return MB_SUCCESS;
00434     }
00435     else
00436     {
00437         neighborhood.push_back( t );
00438     }
00439 
00440     // NOTE: add_intersection may modify the 'neg_ray_len' and 'nonneg_ray_len'
00441     //       members, which will affect subsequent calls to ray_tri_intersect
00442     //       in this loop.
00443     add_intersection( set, t, int_dist, search_win );
00444 
00445     return MB_SUCCESS;
00446 }
00447 
00448 void GQT_IntRegCtxt::append_intersection( EntityHandle set, EntityHandle facet, double dist )
00449 {
00450     intersections.push_back( dist );
00451     sets.push_back( set );
00452     facets.push_back( facet );
00453     neighborhoods.push_back( neighborhood );
00454     return;
00455 }
00456 
00457 void GQT_IntRegCtxt::set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist )
00458 {
00459     intersections[len_idx] = dist;
00460     sets[len_idx]          = set;
00461     facets[len_idx]        = facet;
00462     return;
00463 }
00464 
00465 /* Mode 1: Used if neg_ray_len and nonneg_ray_len are specified
00466    variables used:     nonneg_ray_len, neg_ray_len
00467    variables not used: min_tol_int, tol
00468    1) keep the closest nonneg intersection and one negative intersection, if closer
00469 */
00470 void GQT_IntRegCtxt::add_mode1_intersection( EntityHandle set, EntityHandle facet, double dist,
00471                                              OrientedBoxTreeTool::IntersectSearchWindow& search_win )
00472 {
00473     if( 2 != intersections.size() )
00474     {
00475         intersections.resize( 2, 0 );
00476         sets.resize( 2, 0 );
00477         facets.resize( 2, 0 );
00478         // must initialize this for comparison below
00479         intersections[0] = -std::numeric_limits< double >::max();
00480     }
00481 
00482     // negative case
00483     if( 0.0 > dist )
00484     {
00485         set_intersection( 0, set, facet, dist );
00486         search_win.second = &intersections[0];
00487         // nonnegative case
00488     }
00489     else
00490     {
00491         set_intersection( 1, set, facet, dist );
00492         search_win.first = &intersections[1];
00493         // if the intersection is closer than the negative one, remove the negative one
00494         if( dist < -*( search_win.second ) )
00495         {
00496             set_intersection( 0, 0, 0, -intersections[1] );
00497             search_win.second = &intersections[0];
00498         }
00499     }
00500     //    std::cout << "add_intersection: dist = " << dist << " search_win.second=" <<
00501     //    *search_win.second
00502     //          << " search_win.first=" << *search_win.first << std::endl;
00503     return;
00504 }
00505 
00506 void GQT_IntRegCtxt::add_intersection( EntityHandle set, EntityHandle facet, double dist,
00507                                        OrientedBoxTreeTool::IntersectSearchWindow& search_win )
00508 {
00509 
00510     // Mode 1, detected by non-null neg_ray_len pointer
00511     // keep the closest nonneg intersection and one negative intersection, if closer
00512     if( search_win.second && search_win.first ) { return add_mode1_intersection( set, facet, dist, search_win ); }
00513 
00514     // ---------------------------------------------------------------------------
00515     /*   Mode 2: Used if neg_ray_len is not specified
00516          variables used:     min_tol_int, tol, search_win.first
00517          variables not used: neg_ray_len
00518          1) if(min_tol_int<0) return all intersections
00519          2) otherwise return all inside tolerance and unless there are >min_tol_int
00520          inside of tolerance, return the closest outside of tolerance */
00521     // Mode 2
00522     // If minTolInt is less than zero, return all intersections
00523     if( minTolInt < 0 && dist > -tol )
00524     {
00525         append_intersection( set, facet, dist );
00526         neighborhoods.push_back( neighborhood );
00527         return;
00528     }
00529 
00530     // Check if the 'len' pointer is pointing into the intersection
00531     // list.  If this is the case, then the list contains, at that
00532     // location, an intersection greater than the tolerance away from
00533     // the base point of the ray.
00534     int len_idx = -1;
00535     if( search_win.first && search_win.first >= &intersections[0] &&
00536         search_win.first < &intersections[0] + intersections.size() )
00537         len_idx = search_win.first - &intersections[0];
00538 
00539     // If the intersection is within tol of the ray base point, we
00540     // always add it to the list.
00541     if( dist <= tol )
00542     {
00543         // If the list contains an intersection outside the tolerance...
00544         if( len_idx >= 0 )
00545         {
00546             // If we no longer want an intersection outside the tolerance,
00547             // remove it.
00548             if( (int)intersections.size() >= minTolInt )
00549             {
00550                 set_intersection( len_idx, set, facet, dist );
00551                 // From now on, we want only intersections within the tolerance,
00552                 // so update length accordingly
00553                 search_win.first = &tol;
00554             }
00555             // Otherwise appended to the list and update pointer
00556             else
00557             {
00558                 append_intersection( set, facet, dist );
00559                 search_win.first = &intersections[len_idx];
00560             }
00561         }
00562         // Otherwise just append it
00563         else
00564         {
00565             append_intersection( set, facet, dist );
00566             // If we have all the intersections we want, set
00567             // length such that we will only find further intersections
00568             // within the tolerance
00569             if( (int)intersections.size() >= minTolInt ) search_win.first = &tol;
00570         }
00571     }
00572     // Otherwise the intersection is outside the tolerance
00573     // If we already have an intersection outside the tolerance and
00574     // this one is closer, replace the existing one with this one.
00575     else if( len_idx >= 0 )
00576     {
00577         if( dist <= *search_win.first ) { set_intersection( len_idx, set, facet, dist ); }
00578     }
00579     // Otherwise if we want an intersection outside the tolerance
00580     // and don'thave one yet, add it.
00581     else if( (int)intersections.size() < minTolInt )
00582     {
00583         append_intersection( set, facet, dist );
00584         // update length.  this is currently the closest intersection, so
00585         // only want further intersections that are closer than this one.
00586         search_win.first = &intersections.back();
00587     }
00588 }
00589 
00590 GeomQueryTool::GeomQueryTool( Interface* impl, bool find_geomsets, EntityHandle modelRootSet, bool p_rootSets_vector,
00591                               bool restore_rootSets, bool trace_counting, double overlap_thickness,
00592                               double numerical_precision )
00593     : owns_gtt( true )
00594 {
00595     geomTopoTool = new GeomTopoTool( impl, find_geomsets, modelRootSet, p_rootSets_vector, restore_rootSets );
00596 
00597     senseTag = geomTopoTool->get_sense_tag();
00598 
00599     obbTreeTool = geomTopoTool->obb_tree();
00600     MBI         = geomTopoTool->get_moab_instance();
00601 
00602     counting           = trace_counting;
00603     overlapThickness   = overlap_thickness;
00604     numericalPrecision = numerical_precision;
00605 
00606     // reset query counters
00607     n_pt_in_vol_calls = 0;
00608     n_ray_fire_calls  = 0;
00609 }
00610 
00611 GeomQueryTool::GeomQueryTool( GeomTopoTool* geomtopotool, bool trace_counting, double overlap_thickness,
00612                               double numerical_precision )
00613     : owns_gtt( false )
00614 {
00615 
00616     geomTopoTool = geomtopotool;
00617 
00618     senseTag = geomTopoTool->get_sense_tag();
00619 
00620     obbTreeTool = geomTopoTool->obb_tree();
00621     MBI         = geomTopoTool->get_moab_instance();
00622 
00623     counting           = trace_counting;
00624     overlapThickness   = overlap_thickness;
00625     numericalPrecision = numerical_precision;
00626 
00627     // reset query counters
00628     n_pt_in_vol_calls = 0;
00629     n_ray_fire_calls  = 0;
00630 }
00631 
00632 GeomQueryTool::~GeomQueryTool()
00633 {
00634     if( owns_gtt ) { delete geomTopoTool; }
00635 }
00636 
00637 ErrorCode GeomQueryTool::initialize()
00638 {
00639 
00640     ErrorCode rval;
00641 
00642     rval = geomTopoTool->find_geomsets();MB_CHK_SET_ERR( rval, "Failed to find geometry sets" );
00643 
00644     rval = geomTopoTool->setup_implicit_complement();MB_CHK_SET_ERR( rval, "Couldn't setup the implicit complement" );
00645 
00646     rval = geomTopoTool->construct_obb_trees();MB_CHK_SET_ERR( rval, "Failed to construct OBB trees" );
00647 
00648     return MB_SUCCESS;
00649 }
00650 
00651 void GeomQueryTool::RayHistory::reset()
00652 {
00653     prev_facets.clear();
00654 }
00655 
00656 void GeomQueryTool::RayHistory::reset_to_last_intersection()
00657 {
00658 
00659     if( prev_facets.size() > 1 )
00660     {
00661         prev_facets[0] = prev_facets.back();
00662         prev_facets.resize( 1 );
00663     }
00664 }
00665 
00666 void GeomQueryTool::RayHistory::rollback_last_intersection()
00667 {
00668     if( prev_facets.size() ) prev_facets.pop_back();
00669 }
00670 
00671 ErrorCode GeomQueryTool::RayHistory::get_last_intersection( EntityHandle& last_facet_hit ) const
00672 {
00673     if( prev_facets.size() > 0 )
00674     {
00675         last_facet_hit = prev_facets.back();
00676         return MB_SUCCESS;
00677     }
00678     else
00679     {
00680         return MB_ENTITY_NOT_FOUND;
00681     }
00682 }
00683 
00684 bool GeomQueryTool::RayHistory::in_history( EntityHandle ent ) const
00685 {
00686     return std::find( prev_facets.begin(), prev_facets.end(), ent ) != prev_facets.end();
00687 }
00688 
00689 void GeomQueryTool::RayHistory::add_entity( EntityHandle ent )
00690 {
00691     prev_facets.push_back( ent );
00692 }
00693 
00694 ErrorCode GeomQueryTool::ray_fire( const EntityHandle volume, const double point[3], const double dir[3],
00695                                    EntityHandle& next_surf, double& next_surf_dist, RayHistory* history,
00696                                    double user_dist_limit, int ray_orientation, OrientedBoxTreeTool::TrvStats* stats )
00697 {
00698 
00699     // take some stats that are independent of nps
00700     if( counting )
00701     {
00702         ++n_ray_fire_calls;
00703         if( 0 == n_ray_fire_calls % 10000000 )
00704         {
00705             std::cout << "n_ray_fires=" << n_ray_fire_calls << " n_pt_in_vols=" << n_pt_in_vol_calls << std::endl;
00706         }
00707     }
00708 
00709     if( debug )
00710     {
00711         std::cout << "ray_fire:"
00712                   << " xyz=" << point[0] << " " << point[1] << " " << point[2] << " uvw=" << dir[0] << " " << dir[1]
00713                   << " " << dir[2] << " entity_handle=" << volume << std::endl;
00714     }
00715 
00716     const double huge_val = std::numeric_limits< double >::max();
00717     double dist_limit     = huge_val;
00718     if( user_dist_limit > 0 ) dist_limit = user_dist_limit;
00719 
00720     // don't recreate these every call
00721     std::vector< double > dists;
00722     std::vector< EntityHandle > surfs;
00723     std::vector< EntityHandle > facets;
00724 
00725     EntityHandle root;
00726     ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the obb tree root of the volume" );
00727 
00728     // check behind the ray origin for intersections
00729     double neg_ray_len;
00730     if( 0 == overlapThickness ) { neg_ray_len = -numericalPrecision; }
00731     else
00732     {
00733         neg_ray_len = -overlapThickness;
00734     }
00735 
00736     // optionally, limit the nonneg_ray_len with the distance to next collision.
00737     double nonneg_ray_len = dist_limit;
00738 
00739     // the nonneg_ray_len should not be less than -neg_ray_len, or an overlap
00740     // may be missed due to optimization within ray_intersect_sets
00741     if( nonneg_ray_len < -neg_ray_len ) nonneg_ray_len = -neg_ray_len;
00742     if( 0 > nonneg_ray_len || 0 <= neg_ray_len ) { MB_SET_ERR( MB_FAILURE, "Incorrect ray length provided" ); }
00743 
00744     // min_tolerance_intersections is passed but not used in this call
00745     const int min_tolerance_intersections = 0;
00746 
00747     // numericalPrecision is used for box.intersect_ray and find triangles in the
00748     // neighborhood of edge/node intersections.
00749     GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), point, dir, numericalPrecision, min_tolerance_intersections,
00750                                  &root, &volume, &senseTag, &ray_orientation,
00751                                  history ? &( history->prev_facets ) : NULL );
00752 
00753     OrientedBoxTreeTool::IntersectSearchWindow search_win( &nonneg_ray_len, &neg_ray_len );
00754     rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, point, dir,
00755                                                          search_win, int_reg_ctxt, stats );
00756 
00757     MB_CHK_SET_ERR( rval, "Ray query failed" );
00758 
00759     // If no distances are returned, the particle is lost unless the physics limit
00760     // is being used. If the physics limit is being used, there is no way to tell
00761     // if the particle is lost. To avoid ambiguity, DO NOT use the distance limit
00762     // unless you know lost particles do not occur.
00763     if( dists.empty() )
00764     {
00765         next_surf = 0;
00766         if( debug ) { std::cout << "          next_surf=0 dist=(undef)" << std::endl; }
00767         return MB_SUCCESS;
00768     }
00769 
00770     // Assume that a (neg, nonneg) pair of RTIs could be returned,
00771     // however, only one or the other may exist. dists[] may be populated, but
00772     // intersections are ONLY indicated by nonzero surfs[] and facets[].
00773     if( 2 != dists.size() || 2 != facets.size() ) { MB_SET_ERR( MB_FAILURE, "Incorrect number of facets/distances" ); }
00774     if( 0.0 < dists[0] || 0.0 > dists[1] ) { MB_SET_ERR( MB_FAILURE, "Invalid intersection distance signs" ); }
00775 
00776     // If both negative and nonnegative RTIs are returned, the negative RTI must
00777     // closer to the origin.
00778     if( ( 0 != facets[0] && 0 != facets[1] ) && ( -dists[0] > dists[1] ) )
00779     {
00780         MB_SET_ERR( MB_FAILURE, "Invalid intersection distance values" );
00781     }
00782 
00783     // If an RTI is found at negative distance, perform a PMT to see if the
00784     // particle is inside an overlap.
00785     int exit_idx = -1;
00786     if( 0 != facets[0] )
00787     {
00788         // get the next volume
00789         std::vector< EntityHandle > vols;
00790         EntityHandle nx_vol;
00791         rval = MBI->get_parent_meshsets( surfs[0], vols );MB_CHK_SET_ERR( rval, "Failed to get the parent meshsets" );
00792         if( 2 != vols.size() ) { MB_SET_ERR( MB_FAILURE, "Invaid number of parent volumes found" ); }
00793         if( vols.front() == volume ) { nx_vol = vols.back(); }
00794         else
00795         {
00796             nx_vol = vols.front();
00797         }
00798         // Check to see if the point is actually in the next volume.
00799         // The list of previous facets is used to topologically identify the
00800         // "on_boundary" result of the PMT. This avoids a test that uses proximity
00801         // (a tolerance).
00802         int result;
00803         rval = point_in_volume( nx_vol, point, result, dir, history );MB_CHK_SET_ERR( rval, "Point in volume query failed" );
00804         if( 1 == result ) exit_idx = 0;
00805     }
00806 
00807     // if the negative distance is not the exit, try the nonnegative distance
00808     if( -1 == exit_idx && 0 != facets[1] ) exit_idx = 1;
00809 
00810     // if the exit index is still unknown, the particle is lost
00811     if( -1 == exit_idx )
00812     {
00813         next_surf = 0;
00814         if( debug ) { std::cout << "next surf hit = 0, dist = (undef)" << std::endl; }
00815         return MB_SUCCESS;
00816     }
00817 
00818     // return the intersection
00819     next_surf      = surfs[exit_idx];
00820     next_surf_dist = ( 0 > dists[exit_idx] ? 0 : dists[exit_idx] );
00821 
00822     if( history ) { history->prev_facets.push_back( facets[exit_idx] ); }
00823 
00824     if( debug )
00825     {
00826         if( 0 > dists[exit_idx] ) { std::cout << "          OVERLAP track length=" << dists[exit_idx] << std::endl; }
00827         std::cout << "          next_surf = " << next_surf  // todo: use geomtopotool to get id by entity handle
00828                   << ", dist = " << next_surf_dist << " new_pt=";
00829         for( int i = 0; i < 3; ++i )
00830         {
00831             std::cout << point[i] + dir[i] * next_surf_dist << " ";
00832         }
00833         std::cout << std::endl;
00834     }
00835 
00836     return MB_SUCCESS;
00837 }
00838 
00839 ErrorCode GeomQueryTool::point_in_volume( const EntityHandle volume, const double xyz[3], int& result,
00840                                           const double* uvw, const RayHistory* history )
00841 {
00842     // take some stats that are independent of nps
00843     if( counting ) ++n_pt_in_vol_calls;
00844 
00845     // early fail for piv - see if point inside the root level obb
00846     // if its not even in the box dont bother doing anything else
00847     ErrorCode rval = point_in_box( volume, xyz, result );
00848     if( !result )
00849     {
00850         result = 0;
00851         return MB_SUCCESS;
00852     }
00853 
00854     // get OBB Tree for volume
00855     EntityHandle root;
00856     rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to find the volume's obb tree root" );
00857 
00858     // Don't recreate these every call. These cannot be the same as the ray_fire
00859     // vectors because both are used simultaneously.
00860     std::vector< double > dists;
00861     std::vector< EntityHandle > surfs;
00862     std::vector< EntityHandle > facets;
00863     std::vector< int > dirs;
00864 
00865     // if uvw is not given or is full of zeros, use a random direction
00866     double u = 0, v = 0, w = 0;
00867 
00868     if( uvw )
00869     {
00870         u = uvw[0];
00871         v = uvw[1], w = uvw[2];
00872     }
00873 
00874     if( u == 0 && v == 0 && w == 0 )
00875     {
00876         u                      = rand();
00877         v                      = rand();
00878         w                      = rand();
00879         const double magnitude = sqrt( u * u + v * v + w * w );
00880         u /= magnitude;
00881         v /= magnitude;
00882         w /= magnitude;
00883     }
00884 
00885     const double ray_direction[] = { u, v, w };
00886 
00887     // if overlaps, ray must be cast to infinity and all RTIs must be returned
00888     const double large = 1e15;
00889     double ray_length  = large;
00890 
00891     // If overlaps occur, the pt is inside if traveling along the ray from the
00892     // origin, there are ever more exits than entrances. In lieu of implementing
00893     // that, all intersections to infinity are required if overlaps occur (expensive)
00894     int min_tolerance_intersections;
00895     if( 0 != overlapThickness )
00896     {
00897         min_tolerance_intersections = -1;
00898         // only the first intersection is needed if overlaps do not occur (cheap)
00899     }
00900     else
00901     {
00902         min_tolerance_intersections = 1;
00903     }
00904 
00905     // Get intersection(s) of forward and reverse orientation. Do not return
00906     // glancing intersections or previous facets.
00907     GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), xyz, ray_direction, numericalPrecision,
00908                                  min_tolerance_intersections, &root, &volume, &senseTag, NULL,
00909                                  history ? &( history->prev_facets ) : NULL );
00910 
00911     OrientedBoxTreeTool::IntersectSearchWindow search_win( &ray_length, (double*)NULL );
00912     rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, xyz,
00913                                                          ray_direction, search_win, int_reg_ctxt );MB_CHK_SET_ERR( rval, "Ray fire query failed" );
00914 
00915     // determine orientation of all intersections
00916     // 1 for entering, 0 for leaving, -1 for tangent
00917     // Tangent intersections are not returned from ray_tri_intersect.
00918     dirs.resize( dists.size() );
00919     for( unsigned i = 0; i < dists.size(); ++i )
00920     {
00921         rval = boundary_case( volume, dirs[i], u, v, w, facets[i], surfs[i] );MB_CHK_SET_ERR( rval, "Failed to resolve boundary case" );
00922     }
00923 
00924     // count all crossings
00925     if( 0 != overlapThickness )
00926     {
00927         int sum = 0;
00928         for( unsigned i = 0; i < dirs.size(); ++i )
00929         {
00930             if( 1 == dirs[i] )
00931                 sum += 1;  // +1 for entering
00932             else if( 0 == dirs[i] )
00933                 sum -= 1;  // -1 for leaving
00934             else if( -1 == dirs[i] )
00935             {  //  0 for tangent
00936                 std::cout << "direction==tangent" << std::endl;
00937                 sum += 0;
00938             }
00939             else
00940             {
00941                 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
00942             }
00943         }
00944 
00945         // inside/outside depends on the sum
00946         if( 0 < sum )
00947             result = 0;  // pt is outside (for all vols)
00948         else if( 0 > sum )
00949             result = 1;  // pt is inside  (for all vols)
00950         else if( geomTopoTool->is_implicit_complement( volume ) )
00951             result = 1;  // pt is inside  (for impl_compl_vol)
00952         else
00953             result = 0;  // pt is outside (for all other vols)
00954 
00955         // Only use the first crossing
00956     }
00957     else
00958     {
00959         if( dirs.empty() )
00960         {
00961             result = 0;  // pt is outside
00962         }
00963         else
00964         {
00965             int smallest = std::min_element( dists.begin(), dists.end() ) - dists.begin();
00966             if( 1 == dirs[smallest] )
00967                 result = 0;  // pt is outside
00968             else if( 0 == dirs[smallest] )
00969                 result = 1;  // pt is inside
00970             else if( -1 == dirs[smallest] )
00971             {
00972                 // Should not be here because Plucker ray-triangle test does not
00973                 // return coplanar rays as intersections.
00974                 std::cout << "direction==tangent" << std::endl;
00975                 result = -1;
00976             }
00977             else
00978             {
00979                 MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
00980             }
00981         }
00982     }
00983 
00984     if( debug )
00985         std::cout << "pt_in_vol: result=" << result << " xyz=" << xyz[0] << " " << xyz[1] << " " << xyz[2]
00986                   << " uvw=" << u << " " << v << " " << w << " vol_id=" << volume
00987                   << std::endl;  // todo: use geomtopotool to get id by entity handle
00988 
00989     return MB_SUCCESS;
00990 }
00991 
00992 /**
00993  *  \brief For the volume pointed to and the point wished to be tested, returns
00994  *   whether the point is inside or outside the bounding box of the volume.
00995  * inside = 0, not inside, inside = 1, inside
00996  */
00997 ErrorCode GeomQueryTool::point_in_box( EntityHandle volume, const double point[3], int& inside )
00998 {
00999     double minpt[3];
01000     double maxpt[3];
01001     ErrorCode rval = geomTopoTool->get_bounding_coords( volume, minpt, maxpt );MB_CHK_SET_ERR( rval, "Failed to get the bounding coordinates of the volume" );
01002 
01003     // early exits
01004     if( point[0] > maxpt[0] || point[0] < minpt[0] )
01005     {
01006         inside = 0;
01007         return rval;
01008     }
01009     if( point[1] > maxpt[1] || point[1] < minpt[1] )
01010     {
01011         inside = 0;
01012         return rval;
01013     }
01014     if( point[2] > maxpt[2] || point[2] < minpt[2] )
01015     {
01016         inside = 0;
01017         return rval;
01018     }
01019     inside = 1;
01020     return rval;
01021 }
01022 
01023 ErrorCode GeomQueryTool::test_volume_boundary( const EntityHandle volume, const EntityHandle surface,
01024                                                const double xyz[3], const double uvw[3], int& result,
01025                                                const RayHistory* history )
01026 {
01027     ErrorCode rval;
01028     int dir;
01029 
01030     if( history && history->prev_facets.size() )
01031     {
01032         // the current facet is already available
01033         rval = boundary_case( volume, dir, uvw[0], uvw[1], uvw[2], history->prev_facets.back(), surface );MB_CHK_SET_ERR( rval, "Failed to resolve the boundary case" );
01034     }
01035     else
01036     {
01037         // look up nearest facet
01038 
01039         // Get OBB Tree for surface
01040         EntityHandle root;
01041         rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's OBB tree root" );
01042 
01043         // Get closest triangle on surface
01044         const CartVect point( xyz );
01045         CartVect nearest;
01046         EntityHandle facet_out;
01047         rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out );MB_CHK_SET_ERR( rval, "Failed to find the closest point to location" );
01048 
01049         rval = boundary_case( volume, dir, uvw[0], uvw[1], uvw[2], facet_out, surface );MB_CHK_SET_ERR( rval, "Failed to resolve the boundary case" );
01050     }
01051 
01052     result = dir;
01053 
01054     return MB_SUCCESS;
01055 }
01056 
01057 // use spherical area test to determine inside/outside of a polyhedron.
01058 ErrorCode GeomQueryTool::point_in_volume_slow( EntityHandle volume, const double xyz[3], int& result )
01059 {
01060     ErrorCode rval;
01061     Range faces;
01062     std::vector< EntityHandle > surfs;
01063     std::vector< int > senses;
01064     double sum = 0.0;
01065     const CartVect point( xyz );
01066 
01067     rval = MBI->get_child_meshsets( volume, surfs );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
01068 
01069     senses.resize( surfs.size() );
01070     rval = geomTopoTool->get_surface_senses( volume, surfs.size(), &surfs[0], &senses[0] );MB_CHK_SET_ERR( rval, "Failed to get the volume's surface senses" );
01071 
01072     for( unsigned i = 0; i < surfs.size(); ++i )
01073     {
01074         if( !senses[i] )  // skip non-manifold surfaces
01075             continue;
01076 
01077         double surf_area = 0.0, face_area;
01078         faces.clear();
01079         rval = MBI->get_entities_by_dimension( surfs[i], 2, faces );MB_CHK_SET_ERR( rval, "Failed to get the surface entities by dimension" );
01080 
01081         for( Range::iterator j = faces.begin(); j != faces.end(); ++j )
01082         {
01083             rval = poly_solid_angle( *j, point, face_area );MB_CHK_SET_ERR( rval, "Failed to determin the polygon's solid angle" );
01084 
01085             surf_area += face_area;
01086         }
01087 
01088         sum += senses[i] * surf_area;
01089     }
01090 
01091     result = fabs( sum ) > 2.0 * M_PI;
01092     return MB_SUCCESS;
01093 }
01094 
01095 ErrorCode GeomQueryTool::find_volume( const double xyz[3], EntityHandle& volume, const double* dir )
01096 {
01097     ErrorCode rval;
01098     volume = 0;
01099 
01100     EntityHandle global_surf_tree_root = geomTopoTool->get_one_vol_root();
01101 
01102     // fast check - make sure point is in the implicit complement bounding box
01103     int ic_result;
01104     EntityHandle ic_handle;
01105     rval = geomTopoTool->get_implicit_complement( ic_handle );MB_CHK_SET_ERR( rval, "Failed to get the implicit complement handle" );
01106 
01107     rval = point_in_box( ic_handle, xyz, ic_result );MB_CHK_SET_ERR( rval, "Failed to check implicit complement for containment" );
01108     if( ic_result == 0 )
01109     {
01110         volume = 0;
01111         return MB_ENTITY_NOT_FOUND;
01112     }
01113 
01114     // if geomTopoTool doesn't have a global tree, use a loop over vols (slow)
01115     if( !global_surf_tree_root )
01116     {
01117         rval = find_volume_slow( xyz, volume, dir );
01118         return rval;
01119     }
01120 
01121     moab::CartVect uvw( 0.0 );
01122 
01123     if( dir )
01124     {
01125         uvw[0] = dir[0];
01126         uvw[1] = dir[1];
01127         uvw[2] = dir[2];
01128     }
01129 
01130     if( uvw == 0.0 )
01131     {
01132         uvw[0] = rand();
01133         uvw[1] = rand();
01134         uvw[2] = rand();
01135     }
01136 
01137     // always normalize direction
01138     uvw.normalize();
01139 
01140     // fire a ray along dir and get surface
01141     const double huge_val = std::numeric_limits< double >::max();
01142     double pos_ray_len    = huge_val;
01143     double neg_ray_len    = -huge_val;
01144 
01145     // RIS output data
01146     std::vector< double > dists;
01147     std::vector< EntityHandle > surfs;
01148     std::vector< EntityHandle > facets;
01149 
01150     FindVolumeIntRegCtxt find_vol_reg_ctxt;
01151     OrientedBoxTreeTool::IntersectSearchWindow search_win( &pos_ray_len, &neg_ray_len );
01152     rval =
01153         geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, global_surf_tree_root, numericalPrecision,
01154                                                       xyz, uvw.array(), search_win, find_vol_reg_ctxt );MB_CHK_SET_ERR( rval, "Failed in global tree ray fire" );
01155 
01156     // if there was no intersection, no volume is found
01157     if( surfs.size() == 0 || surfs[0] == 0 )
01158     {
01159         volume = 0;
01160         return MB_ENTITY_NOT_FOUND;
01161     }
01162 
01163     // get the positive distance facet, surface hit
01164     EntityHandle facet = facets[0];
01165     EntityHandle surf  = surfs[0];
01166 
01167     // get these now, we're going to use them no matter what
01168     EntityHandle fwd_vol, bwd_vol;
01169     rval = geomTopoTool->get_surface_senses( surf, fwd_vol, bwd_vol );MB_CHK_SET_ERR( rval, "Failed to get sense data" );
01170     EntityHandle parent_vols[2];
01171     parent_vols[0] = fwd_vol;
01172     parent_vols[1] = bwd_vol;
01173 
01174     // get triangle normal
01175     std::vector< EntityHandle > conn;
01176     CartVect coords[3];
01177     rval = MBI->get_connectivity( &facet, 1, conn );MB_CHK_SET_ERR( rval, "Failed to get triangle connectivity" );
01178 
01179     rval = MBI->get_coords( &conn[0], 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get triangle coordinates" );
01180 
01181     CartVect normal = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] );
01182     normal.normalize();
01183 
01184     // reverse direction if a hit in the negative direction is found
01185     if( dists[0] < 0 ) { uvw *= -1; }
01186 
01187     // if this is a "forward" intersection return the first sense entity
01188     // otherwise return the second, "reverse" sense entity
01189     double dot_prod = uvw % normal;
01190     int idx         = dot_prod > 0.0 ? 0 : 1;
01191 
01192     if( dot_prod == 0.0 )
01193     {
01194         std::cerr << "Tangent dot product in find_volume. Shouldn't be here." << std::endl;
01195         volume = 0;
01196         return MB_FAILURE;
01197     }
01198 
01199     volume = parent_vols[idx];
01200 
01201     return MB_SUCCESS;
01202 }
01203 
01204 ErrorCode GeomQueryTool::find_volume_slow( const double xyz[3], EntityHandle& volume, const double* dir )
01205 {
01206     ErrorCode rval;
01207     volume = 0;
01208     // get all volumes
01209     Range all_vols;
01210     rval = geomTopoTool->get_gsets_by_dimension( 3, all_vols );MB_CHK_SET_ERR( rval, "Failed to get all volumes in the model" );
01211 
01212     Range::iterator it;
01213     int result = 0;
01214     for( it = all_vols.begin(); it != all_vols.end(); it++ )
01215     {
01216         rval = point_in_volume( *it, xyz, result, dir );MB_CHK_SET_ERR( rval, "Failed in point in volume loop" );
01217         if( result )
01218         {
01219             volume = *it;
01220             break;
01221         }
01222     }
01223     return volume ? MB_SUCCESS : MB_ENTITY_NOT_FOUND;
01224 }
01225 
01226 // detemine distance to nearest surface
01227 ErrorCode GeomQueryTool::closest_to_location( EntityHandle volume, const double coords[3], double& result,
01228                                               EntityHandle* closest_surface )
01229 {
01230     // Get OBB Tree for volume
01231     EntityHandle root;
01232     ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's obb tree root" );
01233 
01234     // Get closest triangles in volume
01235     const CartVect point( coords );
01236     CartVect nearest;
01237     EntityHandle facet_out;
01238 
01239     rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out,
01240                                                           closest_surface );MB_CHK_SET_ERR( rval, "Failed to get the closest intersection to location" );
01241     // calculate distance between point and nearest facet
01242     result = ( point - nearest ).length();
01243 
01244     return MB_SUCCESS;
01245 }
01246 
01247 // calculate volume of polyhedron
01248 ErrorCode GeomQueryTool::measure_volume( EntityHandle volume, double& result )
01249 {
01250     ErrorCode rval;
01251     std::vector< EntityHandle > surfaces;
01252     result = 0.0;
01253 
01254     // don't try to calculate volume of implicit complement
01255     if( geomTopoTool->is_implicit_complement( volume ) )
01256     {
01257         result = 1.0;
01258         return MB_SUCCESS;
01259     }
01260 
01261     // get surfaces from volume
01262     rval = MBI->get_child_meshsets( volume, surfaces );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
01263 
01264     // get surface senses
01265     std::vector< int > senses( surfaces.size() );
01266     rval = geomTopoTool->get_surface_senses( volume, surfaces.size(), &surfaces[0], &senses[0] );MB_CHK_SET_ERR( rval, "Failed to retrieve surface-volume sense data. Cannot calculate volume" );
01267 
01268     for( unsigned i = 0; i < surfaces.size(); ++i )
01269     {
01270         // skip non-manifold surfaces
01271         if( !senses[i] ) continue;
01272 
01273         // get triangles in surface
01274         Range triangles;
01275         rval = MBI->get_entities_by_dimension( surfaces[i], 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
01276 
01277         if( !triangles.all_of_type( MBTRI ) )
01278         {
01279             std::cout << "WARNING: Surface " << surfaces[i]  // todo: use geomtopotool to get id by entity handle
01280                       << " contains non-triangle elements. Volume calculation may be incorrect." << std::endl;
01281             triangles.clear();
01282             rval = MBI->get_entities_by_type( surfaces[i], MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
01283         }
01284 
01285         // calculate signed volume beneath surface (x 6.0)
01286         double surf_sum = 0.0;
01287         const EntityHandle* conn;
01288         int len;
01289         CartVect coords[3];
01290         for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
01291         {
01292             rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current triangle" );
01293             if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); }
01294             rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the current triangle's vertices" );
01295 
01296             coords[1] -= coords[0];
01297             coords[2] -= coords[0];
01298             surf_sum += ( coords[0] % ( coords[1] * coords[2] ) );
01299         }
01300         result += senses[i] * surf_sum;
01301     }
01302 
01303     result /= 6.0;
01304     return MB_SUCCESS;
01305 }
01306 
01307 // sum area of elements in surface
01308 ErrorCode GeomQueryTool::measure_area( EntityHandle surface, double& result )
01309 {
01310     // get triangles in surface
01311     Range triangles;
01312     ErrorCode rval = MBI->get_entities_by_dimension( surface, 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" );
01313     if( !triangles.all_of_type( MBTRI ) )
01314     {
01315         std::cout << "WARNING: Surface " << surface  // todo: use geomtopotool to get id by entity handle
01316                   << " contains non-triangle elements. Area calculation may be incorrect." << std::endl;
01317         triangles.clear();
01318         rval = MBI->get_entities_by_type( surface, MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to the surface's triangle entities" );
01319     }
01320 
01321     // calculate sum of area of triangles
01322     result = 0.0;
01323     const EntityHandle* conn;
01324     int len;
01325     CartVect coords[3];
01326     for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
01327     {
01328         rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's connectivity" );
01329         if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); }
01330         rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's vertex coordinates" );
01331 
01332         // calculated area using cross product of triangle edges
01333         CartVect v1 = coords[1] - coords[0];
01334         CartVect v2 = coords[2] - coords[0];
01335         CartVect xp = v1 * v2;
01336         result += xp.length();
01337     }
01338     result *= 0.5;
01339     return MB_SUCCESS;
01340 }
01341 
01342 ErrorCode GeomQueryTool::get_normal( EntityHandle surf, const double in_pt[3], double angle[3],
01343                                      const RayHistory* history )
01344 {
01345     EntityHandle root;
01346     ErrorCode rval = geomTopoTool->get_root( surf, root );MB_CHK_SET_ERR( rval, "Failed to get the surface's obb tree root" );
01347 
01348     std::vector< EntityHandle > facets;
01349 
01350     // if no history or history empty, use nearby facets
01351     if( !history || ( history->prev_facets.size() == 0 ) )
01352     {
01353         rval = geomTopoTool->obb_tree()->closest_to_location( in_pt, root, numericalPrecision, facets );MB_CHK_SET_ERR( rval, "Failed to get closest intersection to location" );
01354     }
01355     // otherwise use most recent facet in history
01356     else
01357     {
01358         facets.push_back( history->prev_facets.back() );
01359     }
01360 
01361     CartVect coords[3], normal( 0.0 );
01362     const EntityHandle* conn;
01363     int len;
01364     for( unsigned i = 0; i < facets.size(); ++i )
01365     {
01366         rval = MBI->get_connectivity( facets[i], conn, len );MB_CHK_SET_ERR( rval, "Failed to get facet connectivity" );
01367         if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); }
01368 
01369         rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
01370 
01371         coords[1] -= coords[0];
01372         coords[2] -= coords[0];
01373         normal += coords[1] * coords[2];
01374     }
01375 
01376     normal.normalize();
01377     normal.get( angle );
01378 
01379     return MB_SUCCESS;
01380 }
01381 
01382 /* SECTION II (private) */
01383 
01384 // If point is on boundary, then this function is called to
01385 // discriminate cases in which the ray is entering or leaving.
01386 // result= 1 -> inside volume or entering volume
01387 // result= 0 -> outside volume or leaving volume
01388 // result=-1 -> on boundary with null or tangent uvw
01389 ErrorCode GeomQueryTool::boundary_case( EntityHandle volume, int& result, double u, double v, double w,
01390                                         EntityHandle facet, EntityHandle surface )
01391 {
01392     ErrorCode rval;
01393 
01394     // test to see if uvw is provided
01395     if( u <= 1.0 && v <= 1.0 && w <= 1.0 )
01396     {
01397 
01398         const CartVect ray_vector( u, v, w );
01399         CartVect coords[3], normal( 0.0 );
01400         const EntityHandle* conn;
01401         int len, sense_out;
01402 
01403         rval = MBI->get_connectivity( facet, conn, len );MB_CHK_SET_ERR( rval, "Failed to get the triangle's connectivity" );
01404         if( 3 != len ) { MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" ); }
01405 
01406         rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
01407 
01408         rval = geomTopoTool->get_sense( surface, volume, sense_out );MB_CHK_SET_ERR( rval, "Failed to get the surface's sense with respect to it's volume" );
01409 
01410         coords[1] -= coords[0];
01411         coords[2] -= coords[0];
01412         normal = sense_out * ( coords[1] * coords[2] );
01413 
01414         double sense = ray_vector % normal;
01415 
01416         if( sense < 0.0 )
01417         {
01418             result = 1;  // inside or entering
01419         }
01420         else if( sense > 0.0 )
01421         {
01422             result = 0;  // outside or leaving
01423         }
01424         else if( sense == 0.0 )
01425         {
01426             result = -1;  // tangent, therefore on boundary
01427         }
01428         else
01429         {
01430             result = -1;  // failure
01431             MB_SET_ERR( MB_FAILURE, "Failed to resolve boundary case" );
01432         }
01433 
01434         // if uvw not provided, return on_boundary.
01435     }
01436     else
01437     {
01438         result = -1;  // on boundary
01439         return MB_SUCCESS;
01440     }
01441 
01442     return MB_SUCCESS;
01443 }
01444 
01445 // point_in_volume_slow, including poly_solid_angle helper subroutine
01446 // are adapted from "Point in Polyhedron Testing Using Spherical Polygons", Paulo Cezar
01447 // Pinto Carvalho and Paulo Roma Cavalcanti, _Graphics Gems V_, pg. 42.  Original algorithm
01448 // was described in "An Efficient Point In Polyhedron Algorithm", Jeff Lane, Bob Magedson,
01449 // and Mike Rarick, _Computer Vision, Graphics, and Image Processing 26_, pg. 118-225, 1984.
01450 
01451 // helper function for point_in_volume_slow.  calculate area of a polygon
01452 // projected into a unit-sphere space
01453 ErrorCode GeomQueryTool::poly_solid_angle( EntityHandle face, const CartVect& point, double& area )
01454 {
01455     ErrorCode rval;
01456 
01457     // Get connectivity
01458     const EntityHandle* conn;
01459     int len;
01460     rval = MBI->get_connectivity( face, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the polygon" );
01461 
01462     // Allocate space to store vertices
01463     CartVect coords_static[4];
01464     std::vector< CartVect > coords_dynamic;
01465     CartVect* coords = coords_static;
01466     if( (unsigned)len > ( sizeof( coords_static ) / sizeof( coords_static[0] ) ) )
01467     {
01468         coords_dynamic.resize( len );
01469         coords = &coords_dynamic[0];
01470     }
01471 
01472     // get coordinates
01473     rval = MBI->get_coords( conn, len, coords->array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the polygon vertices" );
01474 
01475     // calculate normal
01476     CartVect norm( 0.0 ), v1, v0 = coords[1] - coords[0];
01477     for( int i = 2; i < len; ++i )
01478     {
01479         v1 = coords[i] - coords[0];
01480         norm += v0 * v1;
01481         v0 = v1;
01482     }
01483 
01484     // calculate area
01485     double s, ang;
01486     area = 0.0;
01487     CartVect r, n1, n2, b, a = coords[len - 1] - coords[0];
01488     for( int i = 0; i < len; ++i )
01489     {
01490         r   = coords[i] - point;
01491         b   = coords[( i + 1 ) % len] - coords[i];
01492         n1  = a * r;                                          // = norm1 (magnitude is important)
01493         n2  = r * b;                                          // = norm2 (magnitude is important)
01494         s   = ( n1 % n2 ) / ( n1.length() * n2.length() );    // = cos(angle between norm1,norm2)
01495         ang = s <= -1.0 ? M_PI : s >= 1.0 ? 0.0 : acos( s );  // = acos(s)
01496         s   = ( b * a ) % norm;                               // =orientation of triangle wrt point
01497         area += s > 0.0 ? M_PI - ang : M_PI + ang;
01498         a = -b;
01499     }
01500 
01501     area -= M_PI * ( len - 2 );
01502     if( ( norm % r ) > 0 ) area = -area;
01503     return MB_SUCCESS;
01504 }
01505 
01506 void GeomQueryTool::set_overlap_thickness( double new_thickness )
01507 {
01508 
01509     if( new_thickness < 0 || new_thickness > 100 )
01510     {
01511         std::cerr << "Invalid overlap_thickness = " << new_thickness << std::endl;
01512     }
01513     else
01514     {
01515         overlapThickness = new_thickness;
01516     }
01517     std::cout << "Set overlap thickness = " << overlapThickness << std::endl;
01518 }
01519 
01520 void GeomQueryTool::set_numerical_precision( double new_precision )
01521 {
01522 
01523     if( new_precision <= 0 || new_precision > 1 )
01524     {
01525         std::cerr << "Invalid numerical_precision = " << numericalPrecision << std::endl;
01526     }
01527     else
01528     {
01529         numericalPrecision = new_precision;
01530     }
01531 
01532     std::cout << "Set numerical precision = " << numericalPrecision << std::endl;
01533 }
01534 
01535 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines