Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
GeomQueryTool.cpp
Go to the documentation of this file.
1 #include "moab/GeomQueryTool.hpp"
2 
3 #ifdef WIN32 /* windows */
4 #define _USE_MATH_DEFINES // For M_PI
5 #endif
6 #include <string>
7 #include <iostream>
8 #include <fstream>
9 #include <sstream>
10 #include <limits>
11 #include <algorithm>
12 #include <set>
13 
14 #include <cctype>
15 #include <cstring>
16 #include <cstdlib>
17 #include <cstdio>
18 
20 
21 #ifdef __DEBUG
22 #define MB_GQT_DEBUG
23 #endif
24 
25 namespace moab
26 {
27 
28 /** \class FindVolume_IntRegCtxt
29  *
30  *
31  *\brief An intersection context used for finding a volume
32  *
33  * This context is used to find the nearest intersection location
34  * and is intended for use with a global surface tree from
35  * the GeomTopoTool.
36  *
37  * The behavior of this context is relatively simple in that it
38  * returns only one intersection distance, surface, and facet.
39  * Intersections of any orientation are accepted. The positive
40  * value of the search window is limited to the current nearest
41  * intersection distance.
42  *
43  */
44 
46 {
47 
48  public:
49  // Constructor
51  {
52  // initialize return vectors
53  // only one hit is returned in this context
54  intersections.push_back( std::numeric_limits< double >::max() );
55  sets.push_back( 0 );
56  facets.push_back( 0 );
57  }
58 
60  EntityHandle tri,
61  double dist,
64  {
65  // update dist, set, and triangle hit if
66  // we found a new minimum distance
67  double abs_dist = fabs( dist );
68  if( abs_dist < fabs( intersections[0] ) )
69  {
70  intersections[0] = dist;
71  sets[0] = set;
72  facets[0] = tri;
73 
74  // narrow search window based on the hit distance
75  pos = abs_dist;
76  neg = -abs_dist;
77  search_win.first = &pos;
78  search_win.second = &neg;
79  }
80 
81  return MB_SUCCESS;
82  }
83 
84  // storage for updated window values during search
85  double pos;
86  double neg;
87 };
88 
89 /** \class GQT_IntRegCtxt
90  *
91  *\brief An implementation of an Intersection Registration Context for use GQT ray-firing
92  *
93  * This context uses a variety of tests and conditions to confirm whether or
94  * not to accumulate an intersection, to ensure robustness for ray firing.
95  *
96  * This context only accumulates intersections that are oriented parallel to
97  * the 'desiredOrient', if provided, with respect to 'geomVol', using
98  * information in the in 'senseTag'.
99  *
100  * This context only accumulates a single intersection out of a set of
101  * multiple intersections that fall in the same 'neighborhood', where a
102  * 'neighborhood' is defined as facets that share edges or vertices.
103  *
104  * This context only accumulates piercing intersections. This is relevant
105  * for intersections that are found to be on an edge or vertex by the
106  * Plucker test. Such intersections are piercing if the ray has the same
107  * orientation w.r.t. to all fecets that share that edge or vertex.
108  *
109  * This context tests intersections against a list of 'prevFacets' to
110  * prevent a ray from crossing the same facet more than once. The user is
111  * responsible for ensuring that this list is reset when appropriate.
112  *
113  * This context accumulates all intersections within 'tol' of the
114  * start of the ray and if the number of intersections within the
115  * 'tol' of the ray start point is less than 'minTolInt', the next
116  * closest intersection. If the desired result is only the closest
117  * intersection, 'minTolInt' should be 0. This function will return all
118  * intersections, regardless of distance from the start of the ray, if
119  * 'minTolInt' is negative.
120  *
121  */
122 
124 {
125 
126  private:
127  // Input
131  const double tol; /* used for box.intersect_ray, radius of
132  neighborhood for adjacent triangles,
133  and old mode of add_intersection */
134  const int minTolInt;
135 
136  // Optional Input - to screen RTIs by orientation and edge/node intersection
137  const EntityHandle* rootSet; /* used for sphere_intersect */
138  const EntityHandle* geomVol; /* used for determining surface sense */
139  const Tag* senseTag; /* allows screening by triangle orientation.
140  both geomVol and senseTag must be used together. */
141  const int* desiredOrient; /* points to desired orientation of ray with
142  respect to surf normal, if this feature is used.
143  Must point to -1 (reverse) or 1 (forward).
144  geomVol and senseTag are needed for this feature */
145 
146  // Optional Input - to avoid returning these as RTIs
147  const std::vector< EntityHandle >* prevFacets; /* intersections on these triangles
148  will not be returned */
149 
150  // Other Variables
151  std::vector< std::vector< EntityHandle > > neighborhoods;
152  std::vector< EntityHandle > neighborhood;
153 
154  void add_intersection( EntityHandle set,
155  EntityHandle tri,
156  double dist,
158  void append_intersection( EntityHandle set, EntityHandle facet, double dist );
159  void set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist );
161  EntityHandle facet,
162  double dist,
165  const CartVect& ray_direction,
166  const GeomUtil::intersection_type int_type,
167  const std::vector< EntityHandle >& close_tris,
168  const std::vector< int >& close_senses,
169  const Interface* MBI,
170  std::vector< EntityHandle >* neighborhood_tris = 0 );
171 
172  bool in_prevFacets( const EntityHandle tri );
173  bool in_neighborhoods( const EntityHandle tri );
174 
175  public:
177  const double ray_point[3],
178  const double ray_dir[3],
179  double tolerance,
180  int min_tolerance_intersections,
181  const EntityHandle* root_set,
182  const EntityHandle* geom_volume,
183  const Tag* sense_tag,
184  const int* desired_orient,
185  const std::vector< EntityHandle >* prev_facets )
186  : tool( obbtool ), ray_origin( ray_point ), ray_direction( ray_dir ), tol( tolerance ),
187  minTolInt( min_tolerance_intersections ), rootSet( root_set ), geomVol( geom_volume ), senseTag( sense_tag ),
188  desiredOrient( desired_orient ), prevFacets( prev_facets ){
189 
190  };
191 
193  EntityHandle triangle,
194  double distance,
196  GeomUtil::intersection_type int_type );
197 
198  virtual ErrorCode update_orient( EntityHandle set, int* surfTriOrient );
199  virtual const int* getDesiredOrient()
200  {
201  return desiredOrient;
202  };
203 };
204 
206 {
207 
208  ErrorCode rval;
209 
210  // Get desired orientation of surface wrt volume. Use this to return only
211  // exit or entrance intersections.
212  if( geomVol && senseTag && desiredOrient && surfTriOrient )
213  {
214  if( 1 != *desiredOrient && -1 != *desiredOrient )
215  {
216  std::cerr << "error: desired orientation must be 1 (forward) or -1 (reverse)" << std::endl;
217  }
218  EntityHandle vols[2];
219  rval = tool->get_moab_instance()->tag_get_data( *senseTag, &set, 1, vols );
220  assert( MB_SUCCESS == rval );
221  if( MB_SUCCESS != rval ) return rval;
222  if( vols[0] == vols[1] )
223  {
224  std::cerr << "error: surface has positive and negative sense wrt same volume" << std::endl;
225  return MB_FAILURE;
226  }
227  // surfTriOrient will be used by plucker_ray_tri_intersect to avoid
228  // intersections with wrong orientation.
229  if( *geomVol == vols[0] )
230  {
231  *surfTriOrient = *desiredOrient * 1;
232  }
233  else if( *geomVol == vols[1] )
234  {
235  *surfTriOrient = *desiredOrient * ( -1 );
236  }
237  else
238  {
239  assert( false );
240  return MB_FAILURE;
241  }
242  }
243 
244  return MB_SUCCESS;
245 }
246 
248 {
249  return ( prevFacets && ( ( *prevFacets ).end() != find( ( *prevFacets ).begin(), ( *prevFacets ).end(), tri ) ) );
250 }
251 
253 {
254  bool same_neighborhood = false;
255  for( unsigned i = 0; i < neighborhoods.size(); ++i )
256  {
257  if( neighborhoods[i].end() != find( neighborhoods[i].begin(), neighborhoods[i].end(), tri ) )
258  {
259  same_neighborhood = true;
260  continue;
261  }
262  }
263  return same_neighborhood;
264 }
265 
266 /**\brief Determine if a ray-edge/node intersection is glancing or piercing.
267  * This function avoids asking for upward adjacencies to prevent their
268  * creation.
269  *\param tri The intersected triangle
270  *\param ray_dir The direction of the ray
271  *\param int_type The type of intersection (EDGE0, EDGE1, NODE2, ...)
272  *\param close_tris Vector of triangles in the proximity of the intersection
273  *\param close_senses Vector of surface senses for tris in the proximity of
274  * the intersection
275  *\param neighborhood Vector of triangles in the topological neighborhood of the intersection
276  *\return True if piercing, false otherwise.
277  */
279  const CartVect& ray_dir,
280  const GeomUtil::intersection_type int_type,
281  const std::vector< EntityHandle >& close_tris,
282  const std::vector< int >& close_senses,
283  const Interface* MBI,
284  std::vector< EntityHandle >* neighborhood_tris )
285 {
286 
287  // get the node of the triangle
288  const EntityHandle* conn = NULL;
289  int len = 0;
290  ErrorCode rval = MBI->get_connectivity( tri, conn, len );MB_CHK_ERR_RET_VAL( rval, false );
291  // NOTE: original code is next line, but return type was wrong; returning true to get same net effect
292  if( 3 != len ) return false;
293 
294  // get adjacent tris (and keep their corresponding senses)
295  std::vector< EntityHandle > adj_tris;
296  std::vector< int > adj_senses;
297 
298  // node intersection
299  if( GeomUtil::NODE0 == int_type || GeomUtil::NODE1 == int_type || GeomUtil::NODE2 == int_type )
300  {
301 
302  // get the intersected node
303  EntityHandle node;
304  if( GeomUtil::NODE0 == int_type )
305  node = conn[0];
306  else if( GeomUtil::NODE1 == int_type )
307  node = conn[1];
308  else
309  node = conn[2];
310 
311  // get tris adjacent to node
312  for( unsigned i = 0; i < close_tris.size(); ++i )
313  {
314  const EntityHandle* con = NULL;
315  rval = MBI->get_connectivity( close_tris[i], con, len );MB_CHK_ERR_RET_VAL( rval, false );
316  if( 3 != len ) return false;
317 
318  if( node == con[0] || node == con[1] || node == con[2] )
319  {
320  adj_tris.push_back( close_tris[i] );
321  adj_senses.push_back( close_senses[i] );
322  }
323  }
324  if( adj_tris.empty() )
325  {
326  std::cerr << "error: no tris are adjacent to the node" << std::endl;
327  return false;
328  }
329  // edge intersection
330  }
331  else if( GeomUtil::EDGE0 == int_type || GeomUtil::EDGE1 == int_type || GeomUtil::EDGE2 == int_type )
332  {
333 
334  // get the endpoints of the edge
335  EntityHandle endpts[2];
336  if( GeomUtil::EDGE0 == int_type )
337  {
338  endpts[0] = conn[0];
339  endpts[1] = conn[1];
340  }
341  else if( GeomUtil::EDGE1 == int_type )
342  {
343  endpts[0] = conn[1];
344  endpts[1] = conn[2];
345  }
346  else
347  {
348  endpts[0] = conn[2];
349  endpts[1] = conn[0];
350  }
351 
352  // get tris adjacent to edge
353  for( unsigned i = 0; i < close_tris.size(); ++i )
354  {
355  const EntityHandle* con = NULL;
356  rval = MBI->get_connectivity( close_tris[i], con, len );MB_CHK_ERR_RET_VAL( rval, false );
357  if( 3 != len ) return false;
358 
359  // check both orientations in case close_tris are not on the same surface
360  if( ( endpts[0] == con[0] && endpts[1] == con[1] ) || ( endpts[0] == con[1] && endpts[1] == con[0] ) ||
361  ( endpts[0] == con[1] && endpts[1] == con[2] ) || ( endpts[0] == con[2] && endpts[1] == con[1] ) ||
362  ( endpts[0] == con[2] && endpts[1] == con[0] ) || ( endpts[0] == con[0] && endpts[1] == con[2] ) )
363  {
364  adj_tris.push_back( close_tris[i] );
365  adj_senses.push_back( close_senses[i] );
366  }
367  }
368  // In a 2-manifold each edge is adjacent to exactly 2 tris
369  if( 2 != adj_tris.size() )
370  {
371  std::cerr << "error: edge of a manifold must be topologically adjacent to exactly 2 tris" << std::endl;
372  MBI->list_entities( endpts, 2 );
373  return true;
374  }
375  }
376  else
377  {
378  std::cerr << "error: special case not an node/edge intersection" << std::endl;
379  return false;
380  }
381 
382  // The close tris were in proximity to the intersection. The adj_tris are
383  // topologically adjacent to the intersection (the neighborhood).
384  if( neighborhood_tris ) ( *neighborhood_tris ).assign( adj_tris.begin(), adj_tris.end() );
385 
386  // determine glancing/piercing
387  // If a desired_orientation was used in this call to ray_intersect_sets,
388  // the plucker_ray_tri_intersect will have already used it. For a piercing
389  // intersection, the normal of all tris must have the same orientation.
390  int sign = 0;
391  for( unsigned i = 0; i < adj_tris.size(); ++i )
392  {
393  const EntityHandle* con = NULL;
394  rval = MBI->get_connectivity( adj_tris[i], con, len );MB_CHK_ERR_RET_VAL( rval, false );
395  if( 3 != len ) return false;
396 
397  CartVect coords[3];
398  rval = MBI->get_coords( con, len, coords[0].array() );MB_CHK_ERR_RET_VAL( rval, false );
399 
400  // get normal of triangle
401  CartVect v0 = coords[1] - coords[0];
402  CartVect v1 = coords[2] - coords[0];
403  CartVect norm = adj_senses[i] * ( v0 * v1 );
404  double dot_prod = norm % ray_dir;
405 
406  // if the sign has not yet been decided, choose it
407  if( 0 == sign && 0 != dot_prod )
408  {
409  if( 0 < dot_prod )
410  sign = 1;
411  else
412  sign = -1;
413  }
414 
415  // intersection is glancing if tri and ray do not point in same direction
416  // for every triangle
417  if( 0 != sign && 0 > sign * dot_prod ) return false;
418  }
419  return true;
420 }
421 
423  EntityHandle t,
424  double int_dist,
426  GeomUtil::intersection_type int_type )
427 {
428  ErrorCode rval;
429 
430  // Do not accept intersections if they are in the vector of previously intersected
431  // facets.
432  if( in_prevFacets( t ) ) return MB_SUCCESS;
433 
434  // Do not accept intersections if they are in the neighborhood of previous
435  // intersections.
436  if( in_neighborhoods( t ) ) return MB_SUCCESS;
437 
438  neighborhood.clear();
439 
440  // Handle special case of edge/node intersection. Accept piercing
441  // intersections and reject glancing intersections.
442  // The edge_node_intersection function needs to know surface sense wrt volume.
443  // A less-robust implementation could work without sense information.
444  // Would it ever be useful to accept a glancing intersection?
445  if( GeomUtil::INTERIOR != int_type && rootSet && geomVol && senseTag )
446  {
447  // get triangles in the proximity of the intersection
448  std::vector< EntityHandle > close_tris;
449  std::vector< int > close_senses;
450  rval = tool->get_close_tris( ray_origin + int_dist * ray_direction, tol, rootSet, geomVol, senseTag, close_tris,
451  close_senses );
452 
453  if( MB_SUCCESS != rval ) return rval;
454 
455  if( !edge_node_piercing_intersect( t, ray_direction, int_type, close_tris, close_senses,
457  return MB_SUCCESS;
458  }
459  else
460  {
461  neighborhood.push_back( t );
462  }
463 
464  // NOTE: add_intersection may modify the 'neg_ray_len' and 'nonneg_ray_len'
465  // members, which will affect subsequent calls to ray_tri_intersect
466  // in this loop.
467  add_intersection( set, t, int_dist, search_win );
468 
469  return MB_SUCCESS;
470 }
471 
473 {
474  intersections.push_back( dist );
475  sets.push_back( set );
476  facets.push_back( facet );
477  neighborhoods.push_back( neighborhood );
478  return;
479 }
480 
481 void GQT_IntRegCtxt::set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist )
482 {
483  intersections[len_idx] = dist;
484  sets[len_idx] = set;
485  facets[len_idx] = facet;
486  return;
487 }
488 
489 /* Mode 1: Used if neg_ray_len and nonneg_ray_len are specified
490  variables used: nonneg_ray_len, neg_ray_len
491  variables not used: min_tol_int, tol
492  1) keep the closest nonneg intersection and one negative intersection, if closer
493 */
495  EntityHandle facet,
496  double dist,
498 {
499  if( 2 != intersections.size() )
500  {
501  intersections.resize( 2, 0 );
502  sets.resize( 2, 0 );
503  facets.resize( 2, 0 );
504  // must initialize this for comparison below
505  intersections[0] = -std::numeric_limits< double >::max();
506  }
507 
508  // negative case
509  if( 0.0 > dist )
510  {
511  set_intersection( 0, set, facet, dist );
512  search_win.second = &intersections[0];
513  // nonnegative case
514  }
515  else
516  {
517  set_intersection( 1, set, facet, dist );
518  search_win.first = &intersections[1];
519  // if the intersection is closer than the negative one, remove the negative one
520  if( dist < -*( search_win.second ) )
521  {
522  set_intersection( 0, 0, 0, -intersections[1] );
523  search_win.second = &intersections[0];
524  }
525  }
526  // std::cout << "add_intersection: dist = " << dist << " search_win.second=" <<
527  // *search_win.second
528  // << " search_win.first=" << *search_win.first << std::endl;
529  return;
530 }
531 
533  EntityHandle facet,
534  double dist,
536 {
537 
538  // Mode 1, detected by non-null neg_ray_len pointer
539  // keep the closest nonneg intersection and one negative intersection, if closer
540  if( search_win.second && search_win.first )
541  {
542  return add_mode1_intersection( set, facet, dist, search_win );
543  }
544 
545  // ---------------------------------------------------------------------------
546  /* Mode 2: Used if neg_ray_len is not specified
547  variables used: min_tol_int, tol, search_win.first
548  variables not used: neg_ray_len
549  1) if(min_tol_int<0) return all intersections
550  2) otherwise return all inside tolerance and unless there are >min_tol_int
551  inside of tolerance, return the closest outside of tolerance */
552  // Mode 2
553  // If minTolInt is less than zero, return all intersections
554  if( minTolInt < 0 && dist > -tol )
555  {
556  append_intersection( set, facet, dist );
557  neighborhoods.push_back( neighborhood );
558  return;
559  }
560 
561  // Check if the 'len' pointer is pointing into the intersection
562  // list. If this is the case, then the list contains, at that
563  // location, an intersection greater than the tolerance away from
564  // the base point of the ray.
565  int len_idx = -1;
566  if( search_win.first && search_win.first >= &intersections[0] &&
567  search_win.first < &intersections[0] + intersections.size() )
568  len_idx = search_win.first - &intersections[0];
569 
570  // If the intersection is within tol of the ray base point, we
571  // always add it to the list.
572  if( dist <= tol )
573  {
574  // If the list contains an intersection outside the tolerance...
575  if( len_idx >= 0 )
576  {
577  // If we no longer want an intersection outside the tolerance,
578  // remove it.
579  if( (int)intersections.size() >= minTolInt )
580  {
581  set_intersection( len_idx, set, facet, dist );
582  // From now on, we want only intersections within the tolerance,
583  // so update length accordingly
584  search_win.first = &tol;
585  }
586  // Otherwise appended to the list and update pointer
587  else
588  {
589  append_intersection( set, facet, dist );
590  search_win.first = &intersections[len_idx];
591  }
592  }
593  // Otherwise just append it
594  else
595  {
596  append_intersection( set, facet, dist );
597  // If we have all the intersections we want, set
598  // length such that we will only find further intersections
599  // within the tolerance
600  if( (int)intersections.size() >= minTolInt ) search_win.first = &tol;
601  }
602  }
603  // Otherwise the intersection is outside the tolerance
604  // If we already have an intersection outside the tolerance and
605  // this one is closer, replace the existing one with this one.
606  else if( len_idx >= 0 )
607  {
608  if( dist <= *search_win.first )
609  {
610  set_intersection( len_idx, set, facet, dist );
611  }
612  }
613  // Otherwise if we want an intersection outside the tolerance
614  // and don'thave one yet, add it.
615  else if( (int)intersections.size() < minTolInt )
616  {
617  append_intersection( set, facet, dist );
618  // update length. this is currently the closest intersection, so
619  // only want further intersections that are closer than this one.
620  search_win.first = &intersections.back();
621  }
622 }
623 
625  bool find_geomsets,
626  EntityHandle modelRootSet,
627  bool p_rootSets_vector,
628  bool restore_rootSets,
629  bool trace_counting,
630  double overlap_thickness,
631  double numerical_precision )
632  : verbose( false ), owns_gtt( true )
633 {
634  geomTopoTool = new GeomTopoTool( impl, find_geomsets, modelRootSet, p_rootSets_vector, restore_rootSets );
635 
637 
640 
641  counting = trace_counting;
642  overlapThickness = overlap_thickness;
643  numericalPrecision = numerical_precision;
644 
645  // reset query counters
646  n_pt_in_vol_calls = 0;
647  n_ray_fire_calls = 0;
648 }
649 
651  bool trace_counting,
652  double overlap_thickness,
653  double numerical_precision )
654  : verbose( false ), owns_gtt( false )
655 {
656 
657  geomTopoTool = geomtopotool;
658 
660 
663 
664  counting = trace_counting;
665  overlapThickness = overlap_thickness;
666  numericalPrecision = numerical_precision;
667 
668  // reset query counters
669  n_pt_in_vol_calls = 0;
670  n_ray_fire_calls = 0;
671 }
672 
674 {
675  if( owns_gtt )
676  {
677  delete geomTopoTool;
678  }
679 }
680 
682 {
683 
684  ErrorCode rval;
685 
686  rval = geomTopoTool->find_geomsets();MB_CHK_SET_ERR( rval, "Failed to find geometry sets" );
687 
688  rval = geomTopoTool->setup_implicit_complement();MB_CHK_SET_ERR( rval, "Couldn't setup the implicit complement" );
689 
690  rval = geomTopoTool->construct_obb_trees();MB_CHK_SET_ERR( rval, "Failed to construct OBB trees" );
691 
692  return MB_SUCCESS;
693 }
694 
696 {
697  prev_facets.clear();
698 }
699 
701 {
702 
703  if( prev_facets.size() > 1 )
704  {
705  prev_facets[0] = prev_facets.back();
706  prev_facets.resize( 1 );
707  }
708 }
709 
711 {
712  if( prev_facets.size() ) prev_facets.pop_back();
713 }
714 
716 {
717  if( prev_facets.size() > 0 )
718  {
719  last_facet_hit = prev_facets.back();
720  return MB_SUCCESS;
721  }
722  else
723  {
724  return MB_ENTITY_NOT_FOUND;
725  }
726 }
727 
729 {
730  return std::find( prev_facets.begin(), prev_facets.end(), ent ) != prev_facets.end();
731 }
732 
734 {
735  prev_facets.push_back( ent );
736 }
737 
739  const double point[3],
740  const double dir[3],
741  EntityHandle& next_surf,
742  double& next_surf_dist,
743  RayHistory* history,
744  double user_dist_limit,
745  int ray_orientation,
747 {
748 
749  // take some stats that are independent of nps
750  if( counting )
751  {
753  if( 0 == n_ray_fire_calls % 10000000 )
754  {
755  std::cout << "n_ray_fires=" << n_ray_fire_calls << " n_pt_in_vols=" << n_pt_in_vol_calls << std::endl;
756  }
757  }
758 
759 #ifdef MB_GQT_DEBUG
760  std::cout << "ray_fire:"
761  << " xyz=" << point[0] << " " << point[1] << " " << point[2] << " uvw=" << dir[0] << " " << dir[1] << " "
762  << dir[2] << " entity_handle=" << volume << std::endl;
763 #endif
764 
765  const double huge_val = std::numeric_limits< double >::max();
766  double dist_limit = huge_val;
767  if( user_dist_limit > 0 ) dist_limit = user_dist_limit;
768 
769  // don't recreate these every call
770  std::vector< double > dists;
771  std::vector< EntityHandle > surfs;
772  std::vector< EntityHandle > facets;
773 
774  EntityHandle root;
775  ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the obb tree root of the volume" );
776 
777  // check behind the ray origin for intersections
778  double neg_ray_len;
779  if( 0 == overlapThickness )
780  {
781  neg_ray_len = -numericalPrecision;
782  }
783  else
784  {
785  neg_ray_len = -overlapThickness;
786  }
787 
788  // optionally, limit the nonneg_ray_len with the distance to next collision.
789  double nonneg_ray_len = dist_limit;
790 
791  // the nonneg_ray_len should not be less than -neg_ray_len, or an overlap
792  // may be missed due to optimization within ray_intersect_sets
793  if( nonneg_ray_len < -neg_ray_len ) nonneg_ray_len = -neg_ray_len;
794  if( 0 > nonneg_ray_len || 0 <= neg_ray_len )
795  {
796  MB_SET_ERR( MB_FAILURE, "Incorrect ray length provided" );
797  }
798 
799  // min_tolerance_intersections is passed but not used in this call
800  const int min_tolerance_intersections = 0;
801 
802  // numericalPrecision is used for box.intersect_ray and find triangles in the
803  // neighborhood of edge/node intersections.
804  GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), point, dir, numericalPrecision, min_tolerance_intersections,
805  &root, &volume, &senseTag, &ray_orientation,
806  history ? &( history->prev_facets ) : NULL );
807 
808  OrientedBoxTreeTool::IntersectSearchWindow search_win( &nonneg_ray_len, &neg_ray_len );
809  rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, point, dir,
810  search_win, int_reg_ctxt, stats );
811 
812  MB_CHK_SET_ERR( rval, "Ray query failed" );
813 
814  // If no distances are returned, the particle is lost unless the physics limit
815  // is being used. If the physics limit is being used, there is no way to tell
816  // if the particle is lost. To avoid ambiguity, DO NOT use the distance limit
817  // unless you know lost particles do not occur.
818  if( dists.empty() )
819  {
820  next_surf = 0;
821 #ifdef MB_GQT_DEBUG
822  std::cout << " next_surf=0 dist=(undef)" << std::endl;
823 #endif
824  return MB_SUCCESS;
825  }
826 
827  // Assume that a (neg, nonneg) pair of RTIs could be returned,
828  // however, only one or the other may exist. dists[] may be populated, but
829  // intersections are ONLY indicated by nonzero surfs[] and facets[].
830  if( 2 != dists.size() || 2 != facets.size() )
831  {
832  MB_SET_ERR( MB_FAILURE, "Incorrect number of facets/distances" );
833  }
834  if( 0.0 < dists[0] || 0.0 > dists[1] )
835  {
836  MB_SET_ERR( MB_FAILURE, "Invalid intersection distance signs" );
837  }
838 
839  // If both negative and nonnegative RTIs are returned, the negative RTI must
840  // closer to the origin.
841  if( ( 0 != facets[0] && 0 != facets[1] ) && ( -dists[0] > dists[1] ) )
842  {
843  MB_SET_ERR( MB_FAILURE, "Invalid intersection distance values" );
844  }
845 
846  // If an RTI is found at negative distance, perform a PMT to see if the
847  // particle is inside an overlap.
848  int exit_idx = -1;
849  if( 0 != facets[0] )
850  {
851  // get the next volume
852  std::vector< EntityHandle > vols;
853  EntityHandle nx_vol;
854  rval = MBI->get_parent_meshsets( surfs[0], vols );MB_CHK_SET_ERR( rval, "Failed to get the parent meshsets" );
855  if( 2 != vols.size() )
856  {
857  MB_SET_ERR( MB_FAILURE, "Invaid number of parent volumes found" );
858  }
859  if( vols.front() == volume )
860  {
861  nx_vol = vols.back();
862  }
863  else
864  {
865  nx_vol = vols.front();
866  }
867  // Check to see if the point is actually in the next volume.
868  // The list of previous facets is used to topologically identify the
869  // "on_boundary" result of the PMT. This avoids a test that uses proximity
870  // (a tolerance).
871  int result;
872  rval = point_in_volume( nx_vol, point, result, dir, history );MB_CHK_SET_ERR( rval, "Point in volume query failed" );
873  if( 1 == result ) exit_idx = 0;
874  }
875 
876  // if the negative distance is not the exit, try the nonnegative distance
877  if( -1 == exit_idx && 0 != facets[1] ) exit_idx = 1;
878 
879  // if the exit index is still unknown, the particle is lost
880  if( -1 == exit_idx )
881  {
882  next_surf = 0;
883 #ifdef MB_GQT_DEBUG
884  std::cout << "next surf hit = 0, dist = (undef)" << std::endl;
885 #endif
886  return MB_SUCCESS;
887  }
888 
889  // return the intersection
890  next_surf = surfs[exit_idx];
891  next_surf_dist = ( 0 > dists[exit_idx] ? 0 : dists[exit_idx] );
892 
893  if( history )
894  {
895  history->prev_facets.push_back( facets[exit_idx] );
896  }
897 
898 #ifdef MB_GQT_DEBUG
899  if( 0 > dists[exit_idx] )
900  {
901  std::cout << " OVERLAP track length=" << dists[exit_idx] << std::endl;
902  }
903  std::cout << " next_surf = " << next_surf // todo: use geomtopotool to get id by entity handle
904  << ", dist = " << next_surf_dist << " new_pt=";
905  for( int i = 0; i < 3; ++i )
906  {
907  std::cout << point[i] + dir[i] * next_surf_dist << " ";
908  }
909  std::cout << std::endl;
910 #endif
911 
912  return MB_SUCCESS;
913 }
914 
916  const double xyz[3],
917  int& result,
918  const double* uvw,
919  const RayHistory* history )
920 {
921  // take some stats that are independent of nps
922  if( counting ) ++n_pt_in_vol_calls;
923 
924  // early fail for piv - see if point inside the root level obb
925  // if its not even in the box dont bother doing anything else
926  ErrorCode rval = point_in_box( volume, xyz, result );
927  if( !result )
928  {
929  result = 0;
930  return MB_SUCCESS;
931  }
932 
933  // get OBB Tree for volume
934  EntityHandle root;
935  rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to find the volume's obb tree root" );
936 
937  // Don't recreate these every call. These cannot be the same as the ray_fire
938  // vectors because both are used simultaneously.
939  std::vector< double > dists;
940  std::vector< EntityHandle > surfs;
941  std::vector< EntityHandle > facets;
942  std::vector< int > dirs;
943 
944  // if uvw is not given or is full of zeros, use a random direction
945  double u = 0, v = 0, w = 0;
946 
947  if( uvw )
948  {
949  u = uvw[0];
950  v = uvw[1], w = uvw[2];
951  }
952 
953  if( u == 0 && v == 0 && w == 0 )
954  {
955  u = rand();
956  v = rand();
957  w = rand();
958  const double magnitude = sqrt( u * u + v * v + w * w );
959  u /= magnitude;
960  v /= magnitude;
961  w /= magnitude;
962  }
963 
964  const double ray_direction[] = { u, v, w };
965 
966  // if overlaps, ray must be cast to infinity and all RTIs must be returned
967  const double large = 1e15;
968  double ray_length = large;
969 
970  // If overlaps occur, the pt is inside if traveling along the ray from the
971  // origin, there are ever more exits than entrances. In lieu of implementing
972  // that, all intersections to infinity are required if overlaps occur (expensive)
973  int min_tolerance_intersections;
974  if( 0 != overlapThickness )
975  {
976  min_tolerance_intersections = -1;
977  // only the first intersection is needed if overlaps do not occur (cheap)
978  }
979  else
980  {
981  min_tolerance_intersections = 1;
982  }
983 
984  // Get intersection(s) of forward and reverse orientation. Do not return
985  // glancing intersections or previous facets.
986  GQT_IntRegCtxt int_reg_ctxt( geomTopoTool->obb_tree(), xyz, ray_direction, numericalPrecision,
987  min_tolerance_intersections, &root, &volume, &senseTag, NULL,
988  history ? &( history->prev_facets ) : NULL );
989 
990  OrientedBoxTreeTool::IntersectSearchWindow search_win( &ray_length, (double*)NULL );
991  rval = geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, root, numericalPrecision, xyz,
992  ray_direction, search_win, int_reg_ctxt );MB_CHK_SET_ERR( rval, "Ray fire query failed" );
993 
994  // determine orientation of all intersections
995  // 1 for entering, 0 for leaving, -1 for tangent
996  // Tangent intersections are not returned from ray_tri_intersect.
997  dirs.resize( dists.size() );
998  for( unsigned i = 0; i < dists.size(); ++i )
999  {
1000  rval = boundary_case( volume, dirs[i], u, v, w, facets[i], surfs[i] );MB_CHK_SET_ERR( rval, "Failed to resolve boundary case" );
1001  }
1002 
1003  // count all crossings
1004  if( 0 != overlapThickness )
1005  {
1006  int sum = 0;
1007  for( unsigned i = 0; i < dirs.size(); ++i )
1008  {
1009  if( 1 == dirs[i] )
1010  sum += 1; // +1 for entering
1011  else if( 0 == dirs[i] )
1012  sum -= 1; // -1 for leaving
1013  else if( -1 == dirs[i] )
1014  { // 0 for tangent
1015  std::cout << "direction==tangent" << std::endl;
1016  sum += 0;
1017  }
1018  else
1019  {
1020  MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
1021  }
1022  }
1023 
1024  // inside/outside depends on the sum
1025  if( 0 < sum )
1026  result = 0; // pt is outside (for all vols)
1027  else if( 0 > sum )
1028  result = 1; // pt is inside (for all vols)
1029  else if( geomTopoTool->is_implicit_complement( volume ) )
1030  result = 1; // pt is inside (for impl_compl_vol)
1031  else
1032  result = 0; // pt is outside (for all other vols)
1033 
1034  // Only use the first crossing
1035  }
1036  else
1037  {
1038  if( dirs.empty() )
1039  {
1040  result = 0; // pt is outside
1041  }
1042  else
1043  {
1044  int smallest = std::min_element( dists.begin(), dists.end() ) - dists.begin();
1045  if( 1 == dirs[smallest] )
1046  result = 0; // pt is outside
1047  else if( 0 == dirs[smallest] )
1048  result = 1; // pt is inside
1049  else if( -1 == dirs[smallest] )
1050  {
1051  // Should not be here because Plucker ray-triangle test does not
1052  // return coplanar rays as intersections.
1053  std::cerr << "direction==tangent" << std::endl;
1054  result = -1;
1055  }
1056  else
1057  {
1058  MB_SET_ERR( MB_FAILURE, "Error: unknown direction" );
1059  }
1060  }
1061  }
1062 
1063 #ifdef MB_GQT_DEBUG
1064  std::cout << "pt_in_vol: result=" << result << " xyz=" << xyz[0] << " " << xyz[1] << " " << xyz[2] << " uvw=" << u
1065  << " " << v << " " << w << " vol_id=" << volume
1066  << std::endl; // todo: use geomtopotool to get id by entity handle
1067 #endif
1068 
1069  return MB_SUCCESS;
1070 }
1071 
1072 /**
1073  * \brief For the volume pointed to and the point wished to be tested, returns
1074  * whether the point is inside or outside the bounding box of the volume.
1075  * inside = 0, not inside, inside = 1, inside
1076  */
1077 ErrorCode GeomQueryTool::point_in_box( EntityHandle volume, const double point[3], int& inside )
1078 {
1079  double minpt[3];
1080  double maxpt[3];
1081  ErrorCode rval = geomTopoTool->get_bounding_coords( volume, minpt, maxpt );MB_CHK_SET_ERR( rval, "Failed to get the bounding coordinates of the volume" );
1082 
1083  // early exits
1084  if( point[0] > maxpt[0] || point[0] < minpt[0] )
1085  {
1086  inside = 0;
1087  return rval;
1088  }
1089  if( point[1] > maxpt[1] || point[1] < minpt[1] )
1090  {
1091  inside = 0;
1092  return rval;
1093  }
1094  if( point[2] > maxpt[2] || point[2] < minpt[2] )
1095  {
1096  inside = 0;
1097  return rval;
1098  }
1099  inside = 1;
1100  return rval;
1101 }
1102 
1104  const EntityHandle surface,
1105  const double xyz[3],
1106  const double uvw[3],
1107  int& result,
1108  const RayHistory* history )
1109 {
1110  ErrorCode rval;
1111  int dir;
1112 
1113  if( history && history->prev_facets.size() )
1114  {
1115  // the current facet is already available
1116  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" );
1117  }
1118  else
1119  {
1120  // look up nearest facet
1121 
1122  // Get OBB Tree for surface
1123  EntityHandle root;
1124  rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's OBB tree root" );
1125 
1126  // Get closest triangle on surface
1127  const CartVect point( xyz );
1128  CartVect nearest;
1129  EntityHandle facet_out;
1130  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" );
1131 
1132  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" );
1133  }
1134 
1135  result = dir;
1136 
1137  return MB_SUCCESS;
1138 }
1139 
1140 // use spherical area test to determine inside/outside of a polyhedron.
1141 ErrorCode GeomQueryTool::point_in_volume_slow( EntityHandle volume, const double xyz[3], int& result )
1142 {
1143  ErrorCode rval;
1144  Range faces;
1145  std::vector< EntityHandle > surfs;
1146  std::vector< int > senses;
1147  double sum = 0.0;
1148  const CartVect point( xyz );
1149 
1150  rval = MBI->get_child_meshsets( volume, surfs );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
1151 
1152  senses.resize( surfs.size() );
1153  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" );
1154 
1155  for( unsigned i = 0; i < surfs.size(); ++i )
1156  {
1157  if( !senses[i] ) // skip non-manifold surfaces
1158  continue;
1159 
1160  double surf_area = 0.0, face_area;
1161  faces.clear();
1162  rval = MBI->get_entities_by_dimension( surfs[i], 2, faces );MB_CHK_SET_ERR( rval, "Failed to get the surface entities by dimension" );
1163 
1164  for( Range::iterator j = faces.begin(); j != faces.end(); ++j )
1165  {
1166  rval = poly_solid_angle( *j, point, face_area );MB_CHK_SET_ERR( rval, "Failed to determin the polygon's solid angle" );
1167 
1168  surf_area += face_area;
1169  }
1170 
1171  sum += senses[i] * surf_area;
1172  }
1173 
1174  result = fabs( sum ) > 2.0 * M_PI;
1175  return MB_SUCCESS;
1176 }
1177 
1178 ErrorCode GeomQueryTool::find_volume( const double xyz[3], EntityHandle& volume, const double* dir )
1179 {
1180  ErrorCode rval;
1181  volume = 0;
1182 
1183  EntityHandle global_surf_tree_root = geomTopoTool->get_one_vol_root();
1184 
1185  // fast check - make sure point is in the implicit complement bounding box
1186  int ic_result;
1187  EntityHandle ic_handle;
1188  rval = geomTopoTool->get_implicit_complement( ic_handle );MB_CHK_SET_ERR( rval, "Failed to get the implicit complement handle" );
1189 
1190  rval = point_in_box( ic_handle, xyz, ic_result );MB_CHK_SET_ERR( rval, "Failed to check implicit complement for containment" );
1191  if( ic_result == 0 )
1192  {
1193  volume = 0;
1194  return MB_ENTITY_NOT_FOUND;
1195  }
1196 
1197  // if geomTopoTool doesn't have a global tree, use a loop over vols (slow)
1198  if( !global_surf_tree_root )
1199  {
1200  rval = find_volume_slow( xyz, volume, dir );
1201  return rval;
1202  }
1203 
1204  moab::CartVect uvw( 0.0 );
1205 
1206  if( dir )
1207  {
1208  uvw[0] = dir[0];
1209  uvw[1] = dir[1];
1210  uvw[2] = dir[2];
1211  }
1212 
1213  if( uvw == 0.0 )
1214  {
1215  uvw[0] = rand();
1216  uvw[1] = rand();
1217  uvw[2] = rand();
1218  }
1219 
1220  // always normalize direction
1221  uvw.normalize();
1222 
1223  // fire a ray along dir and get surface
1224  const double huge_val = std::numeric_limits< double >::max();
1225  double pos_ray_len = huge_val;
1226  double neg_ray_len = -huge_val;
1227 
1228  // RIS output data
1229  std::vector< double > dists;
1230  std::vector< EntityHandle > surfs;
1231  std::vector< EntityHandle > facets;
1232 
1233  FindVolumeIntRegCtxt find_vol_reg_ctxt;
1234  OrientedBoxTreeTool::IntersectSearchWindow search_win( &pos_ray_len, &neg_ray_len );
1235  rval =
1236  geomTopoTool->obb_tree()->ray_intersect_sets( dists, surfs, facets, global_surf_tree_root, numericalPrecision,
1237  xyz, uvw.array(), search_win, find_vol_reg_ctxt );MB_CHK_SET_ERR( rval, "Failed in global tree ray fire" );
1238 
1239  // if there was no intersection, no volume is found
1240  if( surfs.size() == 0 || surfs[0] == 0 )
1241  {
1242  volume = 0;
1243  return MB_ENTITY_NOT_FOUND;
1244  }
1245 
1246  // get the positive distance facet, surface hit
1247  EntityHandle facet = facets[0];
1248  EntityHandle surf = surfs[0];
1249 
1250  // get these now, we're going to use them no matter what
1251  EntityHandle fwd_vol, bwd_vol;
1252  rval = geomTopoTool->get_surface_senses( surf, fwd_vol, bwd_vol );MB_CHK_SET_ERR( rval, "Failed to get sense data" );
1253  EntityHandle parent_vols[2];
1254  parent_vols[0] = fwd_vol;
1255  parent_vols[1] = bwd_vol;
1256 
1257  // get triangle normal
1258  std::vector< EntityHandle > conn;
1259  CartVect coords[3];
1260  rval = MBI->get_connectivity( &facet, 1, conn );MB_CHK_SET_ERR( rval, "Failed to get triangle connectivity" );
1261 
1262  rval = MBI->get_coords( &conn[0], 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get triangle coordinates" );
1263 
1264  CartVect normal = ( coords[1] - coords[0] ) * ( coords[2] - coords[0] );
1265  normal.normalize();
1266 
1267  // reverse direction if a hit in the negative direction is found
1268  if( dists[0] < 0 )
1269  {
1270  uvw *= -1;
1271  }
1272 
1273  // if this is a "forward" intersection return the first sense entity
1274  // otherwise return the second, "reverse" sense entity
1275  double dot_prod = uvw % normal;
1276  int idx = dot_prod > 0.0 ? 0 : 1;
1277 
1278  if( dot_prod == 0.0 )
1279  {
1280  std::cerr << "Tangent dot product in find_volume. Shouldn't be here." << std::endl;
1281  volume = 0;
1282  return MB_FAILURE;
1283  }
1284 
1285  volume = parent_vols[idx];
1286 
1287  return MB_SUCCESS;
1288 }
1289 
1290 ErrorCode GeomQueryTool::find_volume_slow( const double xyz[3], EntityHandle& volume, const double* dir )
1291 {
1292  ErrorCode rval;
1293  volume = 0;
1294  // get all volumes
1295  Range all_vols;
1296  rval = geomTopoTool->get_gsets_by_dimension( 3, all_vols );MB_CHK_SET_ERR( rval, "Failed to get all volumes in the model" );
1297 
1298  Range::iterator it;
1299  int result = 0;
1300  for( it = all_vols.begin(); it != all_vols.end(); it++ )
1301  {
1302  rval = point_in_volume( *it, xyz, result, dir );MB_CHK_SET_ERR( rval, "Failed in point in volume loop" );
1303  if( result )
1304  {
1305  volume = *it;
1306  break;
1307  }
1308  }
1309  return volume ? MB_SUCCESS : MB_ENTITY_NOT_FOUND;
1310 }
1311 
1312 // detemine distance to nearest surface
1314  const double coords[3],
1315  double& result,
1316  EntityHandle* closest_surface )
1317 {
1318  // Get OBB Tree for volume
1319  EntityHandle root;
1320  ErrorCode rval = geomTopoTool->get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get the volume's obb tree root" );
1321 
1322  // Get closest triangles in volume
1323  const CartVect point( coords );
1324  CartVect nearest;
1325  EntityHandle facet_out;
1326 
1327  rval = geomTopoTool->obb_tree()->closest_to_location( point.array(), root, nearest.array(), facet_out,
1328  closest_surface );MB_CHK_SET_ERR( rval, "Failed to get the closest intersection to location" );
1329  // calculate distance between point and nearest facet
1330  result = ( point - nearest ).length();
1331 
1332  return MB_SUCCESS;
1333 }
1334 
1335 // calculate volume of polyhedron
1337 {
1338  ErrorCode rval;
1339  std::vector< EntityHandle > surfaces;
1340  result = 0.0;
1341 
1342  // don't try to calculate volume of implicit complement
1343  if( geomTopoTool->is_implicit_complement( volume ) )
1344  {
1345  result = 1.0;
1346  return MB_SUCCESS;
1347  }
1348 
1349  // get surfaces from volume
1350  rval = MBI->get_child_meshsets( volume, surfaces );MB_CHK_SET_ERR( rval, "Failed to get the volume's child surfaces" );
1351 
1352  // get surface senses
1353  std::vector< int > senses( surfaces.size() );
1354  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" );
1355 
1356  for( unsigned i = 0; i < surfaces.size(); ++i )
1357  {
1358  // skip non-manifold surfaces
1359  if( !senses[i] ) continue;
1360 
1361  // get triangles in surface
1362  Range triangles;
1363  rval = MBI->get_entities_by_dimension( surfaces[i], 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
1364 
1365  if( !triangles.all_of_type( MBTRI ) )
1366  {
1367  std::cout << "WARNING: Surface " << surfaces[i] // todo: use geomtopotool to get id by entity handle
1368  << " contains non-triangle elements. Volume calculation may be incorrect." << std::endl;
1369  triangles.clear();
1370  rval = MBI->get_entities_by_type( surfaces[i], MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface triangles" );
1371  }
1372 
1373  // calculate signed volume beneath surface (x 6.0)
1374  double surf_sum = 0.0;
1375  const EntityHandle* conn;
1376  int len;
1377  CartVect coords[3];
1378  for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
1379  {
1380  rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current triangle" );
1381  if( 3 != len )
1382  {
1383  MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1384  }
1385  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" );
1386 
1387  coords[1] -= coords[0];
1388  coords[2] -= coords[0];
1389  surf_sum += ( coords[0] % ( coords[1] * coords[2] ) );
1390  }
1391  result += senses[i] * surf_sum;
1392  }
1393 
1394  result /= 6.0;
1395  return MB_SUCCESS;
1396 }
1397 
1398 // sum area of elements in surface
1400 {
1401  // get triangles in surface
1402  Range triangles;
1403  ErrorCode rval = MBI->get_entities_by_dimension( surface, 2, triangles );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" );
1404  if( !triangles.all_of_type( MBTRI ) )
1405  {
1406  std::cout << "WARNING: Surface " << surface // todo: use geomtopotool to get id by entity handle
1407  << " contains non-triangle elements. Area calculation may be incorrect." << std::endl;
1408  triangles.clear();
1409  rval = MBI->get_entities_by_type( surface, MBTRI, triangles );MB_CHK_SET_ERR( rval, "Failed to the surface's triangle entities" );
1410  }
1411 
1412  // calculate sum of area of triangles
1413  result = 0.0;
1414  const EntityHandle* conn;
1415  int len;
1416  CartVect coords[3];
1417  for( Range::iterator j = triangles.begin(); j != triangles.end(); ++j )
1418  {
1419  rval = MBI->get_connectivity( *j, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's connectivity" );
1420  if( 3 != len )
1421  {
1422  MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1423  }
1424  rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get the current triangle's vertex coordinates" );
1425 
1426  // calculated area using cross product of triangle edges
1427  CartVect v1 = coords[1] - coords[0];
1428  CartVect v2 = coords[2] - coords[0];
1429  CartVect xp = v1 * v2;
1430  result += xp.length();
1431  }
1432  result *= 0.5;
1433  return MB_SUCCESS;
1434 }
1435 
1437  const double in_pt[3],
1438  double angle[3],
1439  const RayHistory* history )
1440 {
1441  EntityHandle root;
1442  ErrorCode rval = geomTopoTool->get_root( surf, root );MB_CHK_SET_ERR( rval, "Failed to get the surface's obb tree root" );
1443 
1444  std::vector< EntityHandle > facets;
1445 
1446  // if no history or history empty, use nearby facets
1447  if( !history || ( history->prev_facets.size() == 0 ) )
1448  {
1449  rval = geomTopoTool->obb_tree()->closest_to_location( in_pt, root, numericalPrecision, facets );MB_CHK_SET_ERR( rval, "Failed to get closest intersection to location" );
1450  }
1451  // otherwise use most recent facet in history
1452  else
1453  {
1454  facets.push_back( history->prev_facets.back() );
1455  }
1456 
1457  CartVect coords[3], normal( 0.0 );
1458  const EntityHandle* conn;
1459  int len;
1460  for( unsigned i = 0; i < facets.size(); ++i )
1461  {
1462  rval = MBI->get_connectivity( facets[i], conn, len );MB_CHK_SET_ERR( rval, "Failed to get facet connectivity" );
1463  if( 3 != len )
1464  {
1465  MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1466  }
1467 
1468  rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
1469 
1470  coords[1] -= coords[0];
1471  coords[2] -= coords[0];
1472  normal += coords[1] * coords[2];
1473  }
1474 
1475  normal.normalize();
1476  normal.get( angle );
1477 
1478  return MB_SUCCESS;
1479 }
1480 
1481 /* SECTION II (private) */
1482 
1483 // If point is on boundary, then this function is called to
1484 // discriminate cases in which the ray is entering or leaving.
1485 // result= 1 -> inside volume or entering volume
1486 // result= 0 -> outside volume or leaving volume
1487 // result=-1 -> on boundary with null or tangent uvw
1489  int& result,
1490  double u,
1491  double v,
1492  double w,
1493  EntityHandle facet,
1494  EntityHandle surface )
1495 {
1496  ErrorCode rval;
1497 
1498  // test to see if uvw is provided
1499  if( u <= 1.0 && v <= 1.0 && w <= 1.0 )
1500  {
1501 
1502  const CartVect ray_vector( u, v, w );
1503  CartVect coords[3], normal( 0.0 );
1504  const EntityHandle* conn;
1505  int len, sense_out;
1506 
1507  rval = MBI->get_connectivity( facet, conn, len );MB_CHK_SET_ERR( rval, "Failed to get the triangle's connectivity" );
1508  if( 3 != len )
1509  {
1510  MB_SET_ERR( MB_FAILURE, "Incorrect connectivity length for triangle" );
1511  }
1512 
1513  rval = MBI->get_coords( conn, 3, coords[0].array() );MB_CHK_SET_ERR( rval, "Failed to get vertex coordinates" );
1514 
1515  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" );
1516 
1517  coords[1] -= coords[0];
1518  coords[2] -= coords[0];
1519  normal = sense_out * ( coords[1] * coords[2] );
1520 
1521  double sense = ray_vector % normal;
1522 
1523  if( sense < 0.0 )
1524  {
1525  result = 1; // inside or entering
1526  }
1527  else if( sense > 0.0 )
1528  {
1529  result = 0; // outside or leaving
1530  }
1531  else if( sense == 0.0 )
1532  {
1533  result = -1; // tangent, therefore on boundary
1534  }
1535  else
1536  {
1537  result = -1; // failure
1538  MB_SET_ERR( MB_FAILURE, "Failed to resolve boundary case" );
1539  }
1540 
1541  // if uvw not provided, return on_boundary.
1542  }
1543  else
1544  {
1545  result = -1; // on boundary
1546  return MB_SUCCESS;
1547  }
1548 
1549  return MB_SUCCESS;
1550 }
1551 
1552 // point_in_volume_slow, including poly_solid_angle helper subroutine
1553 // are adapted from "Point in Polyhedron Testing Using Spherical Polygons", Paulo Cezar
1554 // Pinto Carvalho and Paulo Roma Cavalcanti, _Graphics Gems V_, pg. 42. Original algorithm
1555 // was described in "An Efficient Point In Polyhedron Algorithm", Jeff Lane, Bob Magedson,
1556 // and Mike Rarick, _Computer Vision, Graphics, and Image Processing 26_, pg. 118-225, 1984.
1557 
1558 // helper function for point_in_volume_slow. calculate area of a polygon
1559 // projected into a unit-sphere space
1561 {
1562  ErrorCode rval;
1563 
1564  // Get connectivity
1565  const EntityHandle* conn;
1566  int len;
1567  rval = MBI->get_connectivity( face, conn, len, true );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the polygon" );
1568 
1569  // Allocate space to store vertices
1570  CartVect coords_static[4];
1571  std::vector< CartVect > coords_dynamic;
1572  CartVect* coords = coords_static;
1573  if( (unsigned)len > ( sizeof( coords_static ) / sizeof( coords_static[0] ) ) )
1574  {
1575  coords_dynamic.resize( len );
1576  coords = &coords_dynamic[0];
1577  }
1578 
1579  // get coordinates
1580  rval = MBI->get_coords( conn, len, coords->array() );MB_CHK_SET_ERR( rval, "Failed to get the coordinates of the polygon vertices" );
1581 
1582  // calculate normal
1583  CartVect norm( 0.0 ), v1, v0 = coords[1] - coords[0];
1584  for( int i = 2; i < len; ++i )
1585  {
1586  v1 = coords[i] - coords[0];
1587  norm += v0 * v1;
1588  v0 = v1;
1589  }
1590 
1591  // calculate area
1592  double s, ang;
1593  area = 0.0;
1594  CartVect r, n1, n2, b, a = coords[len - 1] - coords[0];
1595  for( int i = 0; i < len; ++i )
1596  {
1597  r = coords[i] - point;
1598  b = coords[( i + 1 ) % len] - coords[i];
1599  n1 = a * r; // = norm1 (magnitude is important)
1600  n2 = r * b; // = norm2 (magnitude is important)
1601  s = ( n1 % n2 ) / ( n1.length() * n2.length() ); // = cos(angle between norm1,norm2)
1602  ang = s <= -1.0 ? M_PI : s >= 1.0 ? 0.0 : acos( s ); // = acos(s)
1603  s = ( b * a ) % norm; // =orientation of triangle wrt point
1604  area += s > 0.0 ? M_PI - ang : M_PI + ang;
1605  a = -b;
1606  }
1607 
1608  area -= M_PI * ( len - 2 );
1609  if( ( norm % r ) > 0 ) area = -area;
1610  return MB_SUCCESS;
1611 }
1612 
1613 void GeomQueryTool::set_overlap_thickness( double new_thickness )
1614 {
1615 
1616  if( new_thickness < 0 || new_thickness > 100 )
1617  {
1618  std::cerr << "Invalid overlap_thickness = " << new_thickness << std::endl;
1619  }
1620  else
1621  {
1622  overlapThickness = new_thickness;
1623  }
1624  if( verbose ) std::cout << "Set overlap thickness = " << overlapThickness << std::endl;
1625 }
1626 
1627 void GeomQueryTool::set_numerical_precision( double new_precision )
1628 {
1629 
1630  if( new_precision <= 0 || new_precision > 1 )
1631  {
1632  std::cerr << "Invalid numerical_precision = " << numericalPrecision << std::endl;
1633  }
1634  else
1635  {
1636  numericalPrecision = new_precision;
1637  }
1638  if( verbose ) std::cout << "Set numerical precision = " << numericalPrecision << std::endl;
1639 }
1640 
1641 } // namespace moab