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