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