 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
19 #include "moab/OrientedBoxTreeTool.hpp"
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
45 class FindVolumeIntRegCtxt : public OrientedBoxTreeTool::IntRegCtxt
46 {
47
48 public:
49 // Constructor
50 FindVolumeIntRegCtxt()
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
59 ErrorCode register_intersection( EntityHandle set,
60 EntityHandle tri,
61 double dist,
62 OrientedBoxTreeTool::IntersectSearchWindow& search_win,
63 GeomUtil::intersection_type /*it*/ )
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
123 class GQT_IntRegCtxt : public OrientedBoxTreeTool::IntRegCtxt
124 {
125
126 private:
127 // Input
128 OrientedBoxTreeTool* tool;
129 const CartVect ray_origin;
130 const CartVect ray_direction;
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,
157 OrientedBoxTreeTool::IntersectSearchWindow& search_win );
158 void append_intersection( EntityHandle set, EntityHandle facet, double dist );
159 void set_intersection( int len_idx, EntityHandle set, EntityHandle facet, double dist );
160 void add_mode1_intersection( EntityHandle set,
161 EntityHandle facet,
162 double dist,
163 OrientedBoxTreeTool::IntersectSearchWindow& search_win );
164 bool edge_node_piercing_intersect( const EntityHandle tri,
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:
176 GQT_IntRegCtxt( OrientedBoxTreeTool* obbtool,
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
192 virtual ErrorCode register_intersection( EntityHandle set,
193 EntityHandle triangle,
194 double distance,
195 OrientedBoxTreeTool::IntersectSearchWindow&,
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
205 ErrorCode GQT_IntRegCtxt::update_orient( EntityHandle set, int* surfTriOrient )
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
247 bool GQT_IntRegCtxt::in_prevFacets( const EntityHandle tri )
248 {
249 return ( prevFacets && ( ( *prevFacets ).end() != find( ( *prevFacets ).begin(), ( *prevFacets ).end(), tri ) ) );
250 }
251
252 bool GQT_IntRegCtxt::in_neighborhoods( const EntityHandle tri )
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 */
278 bool GQT_IntRegCtxt::edge_node_piercing_intersect( const EntityHandle tri,
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
422 ErrorCode GQT_IntRegCtxt::register_intersection( EntityHandle set,
423 EntityHandle t,
424 double int_dist,
425 OrientedBoxTreeTool::IntersectSearchWindow& search_win,
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,
456 tool->get_moab_instance(), &neighborhood ) )
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
472 void GQT_IntRegCtxt::append_intersection( EntityHandle set, EntityHandle facet, double dist )
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 */
494 void GQT_IntRegCtxt::add_mode1_intersection( EntityHandle set,
495 EntityHandle facet,
496 double dist,
497 OrientedBoxTreeTool::IntersectSearchWindow& search_win )
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
532 void GQT_IntRegCtxt::add_intersection( EntityHandle set,
533 EntityHandle facet,
534 double dist,
535 OrientedBoxTreeTool::IntersectSearchWindow& search_win )
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
624 GeomQueryTool::GeomQueryTool( Interface* impl,
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
636 senseTag = geomTopoTool->get_sense_tag();
637
638 obbTreeTool = geomTopoTool->obb_tree();
639 MBI = geomTopoTool->get_moab_instance();
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
650 GeomQueryTool::GeomQueryTool( GeomTopoTool* geomtopotool,
651 bool trace_counting,
652 double overlap_thickness,
653 double numerical_precision )
654 : verbose( false ), owns_gtt( false )
655 {
656
657 geomTopoTool = geomtopotool;
658
659 senseTag = geomTopoTool->get_sense_tag();
660
661 obbTreeTool = geomTopoTool->obb_tree();
662 MBI = geomTopoTool->get_moab_instance();
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
673 GeomQueryTool::~GeomQueryTool()
674 {
675 if( owns_gtt )
676 {
677 delete geomTopoTool;
678 }
679 }
680
681 ErrorCode GeomQueryTool::initialize()
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
695 void GeomQueryTool::RayHistory::reset()
696 {
697 prev_facets.clear();
698 }
699
700 void GeomQueryTool::RayHistory::reset_to_last_intersection()
701 {
702
703 if( prev_facets.size() > 1 )
704 {
705 prev_facets[0] = prev_facets.back();
706 prev_facets.resize( 1 );
707 }
708 }
709
710 void GeomQueryTool::RayHistory::rollback_last_intersection()
711 {
712 if( prev_facets.size() ) prev_facets.pop_back();
713 }
714
715 ErrorCode GeomQueryTool::RayHistory::get_last_intersection( EntityHandle& last_facet_hit ) const
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
728 bool GeomQueryTool::RayHistory::in_history( EntityHandle ent ) const
729 {
730 return std::find( prev_facets.begin(), prev_facets.end(), ent ) != prev_facets.end();
731 }
732
733 void GeomQueryTool::RayHistory::add_entity( EntityHandle ent )
734 {
735 prev_facets.push_back( ent );
736 }
737
738 ErrorCode GeomQueryTool::ray_fire( const EntityHandle volume,
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,
746 OrientedBoxTreeTool::TrvStats* stats )
747 {
748
749 // take some stats that are independent of nps
750 if( counting )
751 {
752 ++n_ray_fire_calls;
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
915 ErrorCode GeomQueryTool::point_in_volume( const EntityHandle volume,
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
1103 ErrorCode GeomQueryTool::test_volume_boundary( const EntityHandle volume,
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
1313 ErrorCode GeomQueryTool::closest_to_location( EntityHandle volume,
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
1336 ErrorCode GeomQueryTool::measure_volume( EntityHandle volume, double& result )
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
1399 ErrorCode GeomQueryTool::measure_area( EntityHandle surface, double& result )
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
1436 ErrorCode GeomQueryTool::get_normal( EntityHandle surf,
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
1488 ErrorCode GeomQueryTool::boundary_case( EntityHandle volume,
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
1560 ErrorCode GeomQueryTool::poly_solid_angle( EntityHandle face, const CartVect& point, double& area )
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
1.9.1.