 |
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 /*
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2008 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 /**\file BSPTree.cpp
17 *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
18 *\date 2008-05-13
19 */
20
21 #include "moab/BSPTree.hpp"
22 #include "moab/GeomUtil.hpp"
23 #include "moab/Range.hpp"
24 #include "Internals.hpp"
25 #include "moab/BSPTreePoly.hpp"
26 #include "moab/Util.hpp"
27
28 #include <cassert>
29 #include <cstring>
30 #include <algorithm>
31 #include <limits>
32
33 #if defined( MOAB_HAVE_IEEEFP_H )
34 #include <ieeefp.h>
35 #endif
36
37 #define MB_BSP_TREE_DEFAULT_TAG_NAME "BSPTree"
38
39 namespace moab
40 {
41
42 static void corners_from_box( const double box_min[3], const double box_max[3], double corners[8][3] )
43 {
44 const double* ranges[] = { box_min, box_max };
45 for( int z = 0; z < 2; ++z )
46 {
47 corners[4 * z][0] = box_min[0];
48 corners[4 * z][1] = box_min[1];
49 corners[4 * z][2] = ranges[z][2];
50
51 corners[4 * z + 1][0] = box_max[0];
52 corners[4 * z + 1][1] = box_min[1];
53 corners[4 * z + 1][2] = ranges[z][2];
54
55 corners[4 * z + 2][0] = box_max[0];
56 corners[4 * z + 2][1] = box_max[1];
57 corners[4 * z + 2][2] = ranges[z][2];
58
59 corners[4 * z + 3][0] = box_min[0];
60 corners[4 * z + 3][1] = box_max[1];
61 corners[4 * z + 3][2] = ranges[z][2];
62 }
63 }
64
65 // assume box has planar sides
66 // test if point is contained in box
67 static bool point_in_box( const double corners[8][3], const double point[3] )
68 {
69 const unsigned side_verts[6][3] = { { 0, 3, 1 }, { 4, 5, 7 }, { 0, 1, 4 }, { 1, 2, 5 }, { 2, 3, 6 }, { 3, 0, 7 } };
70 // If we assume planar sides, then the box is the intersection
71 // of 6 half-spaces defined by the planes of the sides.
72 const CartVect pt( point );
73 for( unsigned s = 0; s < 6; ++s )
74 {
75 CartVect v0( corners[side_verts[s][0]] );
76 CartVect v1( corners[side_verts[s][1]] );
77 CartVect v2( corners[side_verts[s][2]] );
78 CartVect N = ( v1 - v0 ) * ( v2 - v0 );
79 if( ( v0 - pt ) % N < 0 ) return false;
80 }
81 return true;
82 }
83
84 void BSPTree::Plane::set( const double pt1[3], const double pt2[3], const double pt3[3] )
85 {
86 const double v1[] = { pt2[0] - pt1[0], pt2[1] - pt1[1], pt2[2] - pt1[2] };
87 const double v2[] = { pt3[0] - pt1[0], pt3[1] - pt1[1], pt3[2] - pt1[2] };
88 const double nrm[] = { v1[1] * v2[2] - v1[2] * v2[1], v1[2] * v2[0] - v1[0] * v2[2],
89 v1[0] * v2[1] - v1[1] * v2[0] };
90 set( nrm, pt1 );
91 }
92
93 ErrorCode BSPTree::init_tags( const char* tagname )
94 {
95 if( !tagname ) tagname = MB_BSP_TREE_DEFAULT_TAG_NAME;
96
97 std::string rootname( tagname );
98 rootname += "_box";
99
100 ErrorCode rval = moab()->tag_get_handle( tagname, 4, MB_TYPE_DOUBLE, planeTag, MB_TAG_CREAT | MB_TAG_DENSE );
101 if( MB_SUCCESS != rval )
102 planeTag = 0;
103 else
104 rval = moab()->tag_get_handle( rootname.c_str(), 24, MB_TYPE_DOUBLE, rootTag, MB_TAG_CREAT | MB_TAG_SPARSE );
105 if( MB_SUCCESS != rval ) rootTag = 0;
106 return rval;
107 }
108
109 BSPTree::BSPTree( Interface* mb, const char* tagname, unsigned set_flags )
110 : mbInstance( mb ), meshSetFlags( set_flags ), cleanUpTrees( false )
111 {
112 init_tags( tagname );
113 }
114
115 BSPTree::BSPTree( Interface* mb, bool destroy_created_trees, const char* tagname, unsigned set_flags )
116 : mbInstance( mb ), meshSetFlags( set_flags ), cleanUpTrees( destroy_created_trees )
117 {
118 init_tags( tagname );
119 }
120
121 BSPTree::~BSPTree()
122 {
123 if( !cleanUpTrees ) return;
124
125 while( !createdTrees.empty() )
126 {
127 EntityHandle tree = createdTrees.back();
128 // make sure this is a tree (rather than some other, stale handle)
129 const void* data_ptr = 0;
130 ErrorCode rval = moab()->tag_get_by_ptr( rootTag, &tree, 1, &data_ptr );
131 if( MB_SUCCESS == rval ) rval = delete_tree( tree );
132 if( MB_SUCCESS != rval ) createdTrees.pop_back();
133 }
134 }
135
136 ErrorCode BSPTree::set_split_plane( EntityHandle node, const Plane& p )
137 {
138 // check for unit-length normal
139 const double lensqr = p.norm[0] * p.norm[0] + p.norm[1] * p.norm[1] + p.norm[2] * p.norm[2];
140 if( fabs( lensqr - 1.0 ) < std::numeric_limits< double >::epsilon() )
141 return moab()->tag_set_data( planeTag, &node, 1, &p );
142
143 const double inv_len = 1.0 / sqrt( lensqr );
144 Plane p2( p );
145 p2.norm[0] *= inv_len;
146 p2.norm[1] *= inv_len;
147 p2.norm[2] *= inv_len;
148 p2.coeff *= inv_len;
149
150 // check for zero-length normal
151 if( !Util::is_finite( p2.norm[0] + p2.norm[1] + p2.norm[2] + p2.coeff ) ) return MB_FAILURE;
152
153 // store plane
154 return moab()->tag_set_data( planeTag, &node, 1, &p2 );
155 }
156
157 ErrorCode BSPTree::set_tree_box( EntityHandle root_handle, const double box_min[3], const double box_max[3] )
158 {
159 double corners[8][3];
160 corners_from_box( box_min, box_max, corners );
161 return set_tree_box( root_handle, corners );
162 }
163
164 ErrorCode BSPTree::set_tree_box( EntityHandle root_handle, const double corners[8][3] )
165 {
166 return moab()->tag_set_data( rootTag, &root_handle, 1, corners );
167 }
168
169 ErrorCode BSPTree::get_tree_box( EntityHandle root_handle, double corners[8][3] )
170 {
171 return moab()->tag_get_data( rootTag, &root_handle, 1, corners );
172 }
173
174 ErrorCode BSPTree::get_tree_box( EntityHandle root_handle, double corners[24] )
175 {
176 return moab()->tag_get_data( rootTag, &root_handle, 1, corners );
177 }
178
179 ErrorCode BSPTree::create_tree( EntityHandle& root_handle )
180 {
181 const double min[3] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL };
182 const double max[3] = { HUGE_VAL, HUGE_VAL, HUGE_VAL };
183 return create_tree( min, max, root_handle );
184 }
185
186 ErrorCode BSPTree::create_tree( const double corners[8][3], EntityHandle& root_handle )
187 {
188 ErrorCode rval = moab()->create_meshset( meshSetFlags, root_handle );
189 if( MB_SUCCESS != rval ) return rval;
190
191 rval = set_tree_box( root_handle, corners );
192 if( MB_SUCCESS != rval )
193 {
194 moab()->delete_entities( &root_handle, 1 );
195 root_handle = 0;
196 return rval;
197 }
198
199 createdTrees.push_back( root_handle );
200 return MB_SUCCESS;
201 }
202
203 ErrorCode BSPTree::create_tree( const double box_min[3], const double box_max[3], EntityHandle& root_handle )
204 {
205 double corners[8][3];
206 corners_from_box( box_min, box_max, corners );
207 return create_tree( corners, root_handle );
208 }
209
210 ErrorCode BSPTree::delete_tree( EntityHandle root_handle )
211 {
212 ErrorCode rval;
213
214 std::vector< EntityHandle > children, dead_sets, current_sets;
215 current_sets.push_back( root_handle );
216 while( !current_sets.empty() )
217 {
218 EntityHandle set = current_sets.back();
219 current_sets.pop_back();
220 dead_sets.push_back( set );
221 rval = moab()->get_child_meshsets( set, children );
222 if( MB_SUCCESS != rval ) return rval;
223 std::copy( children.begin(), children.end(), std::back_inserter( current_sets ) );
224 children.clear();
225 }
226
227 rval = moab()->tag_delete_data( rootTag, &root_handle, 1 );
228 if( MB_SUCCESS != rval ) return rval;
229
230 createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), root_handle ), createdTrees.end() );
231 return moab()->delete_entities( &dead_sets[0], dead_sets.size() );
232 }
233
234 ErrorCode BSPTree::find_all_trees( Range& results )
235 {
236 return moab()->get_entities_by_type_and_tag( 0, MBENTITYSET, &rootTag, 0, 1, results );
237 }
238
239 ErrorCode BSPTree::get_tree_iterator( EntityHandle root, BSPTreeIter& iter )
240 {
241 ErrorCode rval = iter.initialize( this, root );
242 if( MB_SUCCESS != rval ) return rval;
243 return iter.step_to_first_leaf( BSPTreeIter::LEFT );
244 }
245
246 ErrorCode BSPTree::get_tree_end_iterator( EntityHandle root, BSPTreeIter& iter )
247 {
248 ErrorCode rval = iter.initialize( this, root );
249 if( MB_SUCCESS != rval ) return rval;
250 return iter.step_to_first_leaf( BSPTreeIter::RIGHT );
251 }
252
253 ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane, EntityHandle& left, EntityHandle& right )
254 {
255 ErrorCode rval;
256
257 rval = moab()->create_meshset( meshSetFlags, left );
258 if( MB_SUCCESS != rval ) return rval;
259
260 rval = moab()->create_meshset( meshSetFlags, right );
261 if( MB_SUCCESS != rval )
262 {
263 moab()->delete_entities( &left, 1 );
264 return rval;
265 }
266
267 if( MB_SUCCESS != set_split_plane( leaf.handle(), plane ) ||
268 MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), left ) ||
269 MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), right ) ||
270 MB_SUCCESS != leaf.step_to_first_leaf( BSPTreeIter::LEFT ) )
271 {
272 EntityHandle children[] = { left, right };
273 moab()->delete_entities( children, 2 );
274 return MB_FAILURE;
275 }
276
277 return MB_SUCCESS;
278 }
279
280 ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane )
281 {
282 EntityHandle left, right;
283 return split_leaf( leaf, plane, left, right );
284 }
285
286 ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane, const Range& left_entities, const Range& right_entities )
287 {
288 EntityHandle left, right, parent = leaf.handle();
289 ErrorCode rval = split_leaf( leaf, plane, left, right );
290 if( MB_SUCCESS != rval ) return rval;
291
292 if( MB_SUCCESS == moab()->add_entities( left, left_entities ) &&
293 MB_SUCCESS == moab()->add_entities( right, right_entities ) &&
294 MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
295 return MB_SUCCESS;
296
297 moab()->remove_child_meshset( parent, left );
298 moab()->remove_child_meshset( parent, right );
299 EntityHandle children[] = { left, right };
300 moab()->delete_entities( children, 2 );
301 return MB_FAILURE;
302 }
303
304 ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf,
305 Plane plane,
306 const std::vector< EntityHandle >& left_entities,
307 const std::vector< EntityHandle >& right_entities )
308 {
309 EntityHandle left, right, parent = leaf.handle();
310 ErrorCode rval = split_leaf( leaf, plane, left, right );
311 if( MB_SUCCESS != rval ) return rval;
312
313 if( MB_SUCCESS == moab()->add_entities( left, &left_entities[0], left_entities.size() ) &&
314 MB_SUCCESS == moab()->add_entities( right, &right_entities[0], right_entities.size() ) &&
315 MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
316 return MB_SUCCESS;
317
318 moab()->remove_child_meshset( parent, left );
319 moab()->remove_child_meshset( parent, right );
320 EntityHandle children[] = { left, right };
321 moab()->delete_entities( children, 2 );
322 return MB_FAILURE;
323 }
324
325 ErrorCode BSPTree::merge_leaf( BSPTreeIter& iter )
326 {
327 ErrorCode rval;
328 if( iter.depth() == 1 ) // at root
329 return MB_FAILURE;
330
331 // Move iter to parent
332 iter.up();
333
334 // Get sets to merge
335 EntityHandle parent = iter.handle();
336 iter.childVect.clear();
337 rval = moab()->get_child_meshsets( parent, iter.childVect );
338 if( MB_SUCCESS != rval ) return rval;
339
340 // Remove child links
341 moab()->remove_child_meshset( parent, iter.childVect[0] );
342 moab()->remove_child_meshset( parent, iter.childVect[1] );
343 std::vector< EntityHandle > stack( iter.childVect );
344
345 // Get all entities from children and put them in parent
346 Range range;
347 while( !stack.empty() )
348 {
349 EntityHandle h = stack.back();
350 stack.pop_back();
351 range.clear();
352 rval = moab()->get_entities_by_handle( h, range );
353 if( MB_SUCCESS != rval ) return rval;
354 rval = moab()->add_entities( parent, range );
355 if( MB_SUCCESS != rval ) return rval;
356
357 iter.childVect.clear();
358 rval = moab()->get_child_meshsets( h, iter.childVect );MB_CHK_ERR( rval );
359 if( !iter.childVect.empty() )
360 {
361 moab()->remove_child_meshset( h, iter.childVect[0] );
362 moab()->remove_child_meshset( h, iter.childVect[1] );
363 stack.push_back( iter.childVect[0] );
364 stack.push_back( iter.childVect[1] );
365 }
366
367 rval = moab()->delete_entities( &h, 1 );
368 if( MB_SUCCESS != rval ) return rval;
369 }
370
371 return MB_SUCCESS;
372 }
373
374 ErrorCode BSPTreeIter::initialize( BSPTree* btool, EntityHandle root, const double* /*point*/ )
375 {
376 treeTool = btool;
377 mStack.clear();
378 mStack.push_back( root );
379 return MB_SUCCESS;
380 }
381
382 ErrorCode BSPTreeIter::step_to_first_leaf( Direction direction )
383 {
384 ErrorCode rval;
385 for( ;; )
386 {
387 childVect.clear();
388 rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
389 if( MB_SUCCESS != rval ) return rval;
390 if( childVect.empty() ) // leaf
391 break;
392
393 mStack.push_back( childVect[direction] );
394 }
395 return MB_SUCCESS;
396 }
397
398 ErrorCode BSPTreeIter::step( Direction direction )
399 {
400 EntityHandle node, parent;
401 ErrorCode rval;
402 const Direction opposite = static_cast< Direction >( 1 - direction );
403
404 // If stack is empty, then either this iterator is uninitialized
405 // or we reached the end of the iteration (and return
406 // MB_ENTITY_NOT_FOUND) already.
407 if( mStack.empty() ) return MB_FAILURE;
408
409 // Pop the current node from the stack.
410 // The stack should then contain the parent of the current node.
411 // If the stack is empty after this pop, then we've reached the end.
412 node = mStack.back();
413 mStack.pop_back();
414
415 while( !mStack.empty() )
416 {
417 // Get data for parent entity
418 parent = mStack.back();
419 childVect.clear();
420 rval = tool()->moab()->get_child_meshsets( parent, childVect );
421 if( MB_SUCCESS != rval ) return rval;
422
423 // If we're at the left child
424 if( childVect[opposite] == node )
425 {
426 // push right child on stack
427 mStack.push_back( childVect[direction] );
428 // descend to left-most leaf of the right child
429 return step_to_first_leaf( opposite );
430 }
431
432 // The current node is the right child of the parent,
433 // continue up the tree.
434 assert( childVect[direction] == node );
435 node = parent;
436 mStack.pop_back();
437 }
438
439 return MB_ENTITY_NOT_FOUND;
440 }
441
442 ErrorCode BSPTreeIter::up()
443 {
444 if( mStack.size() < 2 ) return MB_ENTITY_NOT_FOUND;
445 mStack.pop_back();
446 return MB_SUCCESS;
447 }
448
449 ErrorCode BSPTreeIter::down( const BSPTree::Plane& /*plane*/, Direction dir )
450 {
451 childVect.clear();
452 ErrorCode rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
453 if( MB_SUCCESS != rval ) return rval;
454 if( childVect.empty() ) return MB_ENTITY_NOT_FOUND;
455
456 mStack.push_back( childVect[dir] );
457 return MB_SUCCESS;
458 }
459
460 ErrorCode BSPTreeIter::get_parent_split_plane( BSPTree::Plane& plane ) const
461 {
462 if( mStack.size() < 2 ) // at tree root
463 return MB_ENTITY_NOT_FOUND;
464
465 EntityHandle parent = mStack[mStack.size() - 2];
466 return tool()->get_split_plane( parent, plane );
467 }
468
469 double BSPTreeIter::volume() const
470 {
471 BSPTreePoly polyhedron;
472 ErrorCode rval = calculate_polyhedron( polyhedron );
473 return MB_SUCCESS == rval ? polyhedron.volume() : -1.0;
474 }
475
476 bool BSPTreeIter::is_sibling( const BSPTreeIter& other_leaf ) const
477 {
478 const size_t s = mStack.size();
479 return ( s > 1 ) && ( s == other_leaf.mStack.size() ) && ( other_leaf.mStack[s - 2] == mStack[s - 2] ) &&
480 other_leaf.handle() != handle();
481 }
482
483 bool BSPTreeIter::is_sibling( EntityHandle other_leaf ) const
484 {
485 if( mStack.size() < 2 || other_leaf == handle() ) return false;
486 EntityHandle parent = mStack[mStack.size() - 2];
487 childVect.clear();
488 ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
489 if( MB_SUCCESS != rval || childVect.size() != 2 )
490 {
491 assert( false );
492 return false;
493 }
494 return childVect[0] == other_leaf || childVect[1] == other_leaf;
495 }
496
497 bool BSPTreeIter::sibling_is_forward() const
498 {
499 if( mStack.size() < 2 ) // if root
500 return false;
501 EntityHandle parent = mStack[mStack.size() - 2];
502 childVect.clear();
503 ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
504 if( MB_SUCCESS != rval || childVect.size() != 2 )
505 {
506 assert( false );
507 return false;
508 }
509 return childVect[0] == handle();
510 }
511
512 ErrorCode BSPTreeIter::calculate_polyhedron( BSPTreePoly& poly_out ) const
513 {
514 ErrorCode rval;
515
516 assert( sizeof( CartVect ) == 3 * sizeof( double ) );
517 CartVect corners[8];
518 rval = treeTool->get_tree_box( mStack.front(), corners[0].array() );
519 if( MB_SUCCESS != rval ) return rval;
520
521 rval = poly_out.set( corners );
522 if( MB_SUCCESS != rval ) return rval;
523
524 BSPTree::Plane plane;
525 std::vector< EntityHandle >::const_iterator i = mStack.begin();
526 std::vector< EntityHandle >::const_iterator here = mStack.end() - 1;
527 while( i != here )
528 {
529 rval = treeTool->get_split_plane( *i, plane );
530 if( MB_SUCCESS != rval ) return rval;
531
532 childVect.clear();
533 rval = treeTool->moab()->get_child_meshsets( *i, childVect );
534 if( MB_SUCCESS != rval ) return rval;
535 if( childVect.size() != 2 ) return MB_FAILURE;
536
537 ++i;
538 if( childVect[1] == *i ) plane.flip();
539
540 CartVect norm( plane.norm );
541 poly_out.cut_polyhedron( norm, plane.coeff );
542 }
543
544 return MB_SUCCESS;
545 }
546
547 ErrorCode BSPTreeBoxIter::initialize( BSPTree* tool_ptr, EntityHandle root, const double* point )
548 {
549 ErrorCode rval = BSPTreeIter::initialize( tool_ptr, root );
550 if( MB_SUCCESS != rval ) return rval;
551
552 rval = tool()->get_tree_box( root, leafCoords );
553 if( MB_SUCCESS != rval ) return rval;
554
555 if( point && !point_in_box( leafCoords, point ) ) return MB_ENTITY_NOT_FOUND;
556
557 stackData.resize( 1 );
558 return MB_SUCCESS;
559 }
560
561 BSPTreeBoxIter::SideBits BSPTreeBoxIter::side_above_plane( const double hex_coords[8][3], const BSPTree::Plane& plane )
562 {
563 unsigned result = 0;
564 for( unsigned i = 0; i < 8u; ++i )
565 result |= plane.above( hex_coords[i] ) << i;
566 return (BSPTreeBoxIter::SideBits)result;
567 }
568
569 BSPTreeBoxIter::SideBits BSPTreeBoxIter::side_on_plane( const double hex_coords[8][3], const BSPTree::Plane& plane )
570 {
571 unsigned result = 0;
572 for( unsigned i = 0; i < 8u; ++i )
573 {
574 bool on = plane.distance( hex_coords[i] ) <= BSPTree::epsilon();
575 result |= on << i;
576 }
577 return (BSPTreeBoxIter::SideBits)result;
578 }
579
580 static inline void copy_coords( const double src[3], double dest[3] )
581 {
582 dest[0] = src[0];
583 dest[1] = src[1];
584 dest[2] = src[2];
585 }
586
587 ErrorCode BSPTreeBoxIter::face_corners( const SideBits face, const double hex_corners[8][3], double face_corners[4][3] )
588 {
589 switch( face )
590 {
591 case BSPTreeBoxIter::B0154:
592 copy_coords( hex_corners[0], face_corners[0] );
593 copy_coords( hex_corners[1], face_corners[1] );
594 copy_coords( hex_corners[5], face_corners[2] );
595 copy_coords( hex_corners[4], face_corners[3] );
596 break;
597 case BSPTreeBoxIter::B1265:
598 copy_coords( hex_corners[1], face_corners[0] );
599 copy_coords( hex_corners[2], face_corners[1] );
600 copy_coords( hex_corners[6], face_corners[2] );
601 copy_coords( hex_corners[5], face_corners[3] );
602 break;
603 case BSPTreeBoxIter::B2376:
604 copy_coords( hex_corners[2], face_corners[0] );
605 copy_coords( hex_corners[3], face_corners[1] );
606 copy_coords( hex_corners[7], face_corners[2] );
607 copy_coords( hex_corners[6], face_corners[3] );
608 break;
609 case BSPTreeBoxIter::B3047:
610 copy_coords( hex_corners[3], face_corners[0] );
611 copy_coords( hex_corners[0], face_corners[1] );
612 copy_coords( hex_corners[4], face_corners[2] );
613 copy_coords( hex_corners[7], face_corners[3] );
614 break;
615 case BSPTreeBoxIter::B3210:
616 copy_coords( hex_corners[3], face_corners[0] );
617 copy_coords( hex_corners[2], face_corners[1] );
618 copy_coords( hex_corners[1], face_corners[2] );
619 copy_coords( hex_corners[0], face_corners[3] );
620 break;
621 case BSPTreeBoxIter::B4567:
622 copy_coords( hex_corners[4], face_corners[0] );
623 copy_coords( hex_corners[5], face_corners[1] );
624 copy_coords( hex_corners[6], face_corners[2] );
625 copy_coords( hex_corners[7], face_corners[3] );
626 break;
627 default:
628 return MB_FAILURE; // child is not a box
629 }
630
631 return MB_SUCCESS;
632 }
633
634 /** \brief Clip an edge using a plane
635 *
636 * Given an edge from keep_end_coords to cut_end_coords,
637 * cut the edge using the passed plane, such that cut_end_coords
638 * is updated with a new location on the plane, and old_coords_out
639 * contains the original value of cut_end_coords.
640 */
641 static inline void plane_cut_edge( double old_coords_out[3],
642 const double keep_end_coords[3],
643 double cut_end_coords[3],
644 const BSPTree::Plane& plane )
645 {
646 const CartVect start( keep_end_coords ), end( cut_end_coords );
647 const CartVect norm( plane.norm );
648 CartVect xsect_point;
649
650 const CartVect m = end - start;
651 const double t = -( norm % start + plane.coeff ) / ( norm % m );
652 assert( t > 0.0 && t < 1.0 );
653 xsect_point = start + t * m;
654
655 end.get( old_coords_out );
656 xsect_point.get( cut_end_coords );
657 }
658
659 /** Given the corners of a hexahedron in corners_input and a
660 * plane, cut the hex with the plane, updating corners_input
661 * and storing the original,cut-off side of the hex in cut_face_out.
662 *
663 * The portion of the hex below the plane is retained. cut_face_out
664 * will contain the side of the hex that is entirely above the plane.
665 *\return MB_FAILURE if plane/hex intersection is not a quadrilateral.
666 */
667 static ErrorCode plane_cut_box( double cut_face_out[4][3], double corners_inout[8][3], const BSPTree::Plane& plane )
668 {
669 switch( BSPTreeBoxIter::side_above_plane( corners_inout, plane ) )
670 {
671 case BSPTreeBoxIter::B0154:
672 plane_cut_edge( cut_face_out[0], corners_inout[3], corners_inout[0], plane );
673 plane_cut_edge( cut_face_out[1], corners_inout[2], corners_inout[1], plane );
674 plane_cut_edge( cut_face_out[2], corners_inout[6], corners_inout[5], plane );
675 plane_cut_edge( cut_face_out[3], corners_inout[7], corners_inout[4], plane );
676 break;
677 case BSPTreeBoxIter::B1265:
678 plane_cut_edge( cut_face_out[0], corners_inout[0], corners_inout[1], plane );
679 plane_cut_edge( cut_face_out[1], corners_inout[3], corners_inout[2], plane );
680 plane_cut_edge( cut_face_out[2], corners_inout[7], corners_inout[6], plane );
681 plane_cut_edge( cut_face_out[3], corners_inout[4], corners_inout[5], plane );
682 break;
683 case BSPTreeBoxIter::B2376:
684 plane_cut_edge( cut_face_out[0], corners_inout[1], corners_inout[2], plane );
685 plane_cut_edge( cut_face_out[1], corners_inout[0], corners_inout[3], plane );
686 plane_cut_edge( cut_face_out[2], corners_inout[4], corners_inout[7], plane );
687 plane_cut_edge( cut_face_out[3], corners_inout[5], corners_inout[6], plane );
688 break;
689 case BSPTreeBoxIter::B3047:
690 plane_cut_edge( cut_face_out[0], corners_inout[2], corners_inout[3], plane );
691 plane_cut_edge( cut_face_out[1], corners_inout[1], corners_inout[0], plane );
692 plane_cut_edge( cut_face_out[2], corners_inout[5], corners_inout[4], plane );
693 plane_cut_edge( cut_face_out[3], corners_inout[6], corners_inout[7], plane );
694 break;
695 case BSPTreeBoxIter::B3210:
696 plane_cut_edge( cut_face_out[0], corners_inout[7], corners_inout[3], plane );
697 plane_cut_edge( cut_face_out[1], corners_inout[6], corners_inout[2], plane );
698 plane_cut_edge( cut_face_out[2], corners_inout[5], corners_inout[1], plane );
699 plane_cut_edge( cut_face_out[3], corners_inout[4], corners_inout[0], plane );
700 break;
701 case BSPTreeBoxIter::B4567:
702 plane_cut_edge( cut_face_out[0], corners_inout[0], corners_inout[4], plane );
703 plane_cut_edge( cut_face_out[1], corners_inout[1], corners_inout[5], plane );
704 plane_cut_edge( cut_face_out[2], corners_inout[2], corners_inout[6], plane );
705 plane_cut_edge( cut_face_out[3], corners_inout[3], corners_inout[7], plane );
706 break;
707 default:
708 return MB_FAILURE; // child is not a box
709 }
710
711 return MB_SUCCESS;
712 }
713
714 static inline void copy_coords( double dest[3], const double source[3] )
715 {
716 dest[0] = source[0];
717 dest[1] = source[1];
718 dest[2] = source[2];
719 }
720
721 /** reverse of plane_cut_box */
722 static inline ErrorCode plane_uncut_box( const double cut_face_in[4][3],
723 double corners_inout[8][3],
724 const BSPTree::Plane& plane )
725 {
726 switch( BSPTreeBoxIter::side_on_plane( corners_inout, plane ) )
727 {
728 case BSPTreeBoxIter::B0154:
729 copy_coords( corners_inout[0], cut_face_in[0] );
730 copy_coords( corners_inout[1], cut_face_in[1] );
731 copy_coords( corners_inout[5], cut_face_in[2] );
732 copy_coords( corners_inout[4], cut_face_in[3] );
733 break;
734 case BSPTreeBoxIter::B1265:
735 copy_coords( corners_inout[1], cut_face_in[0] );
736 copy_coords( corners_inout[2], cut_face_in[1] );
737 copy_coords( corners_inout[6], cut_face_in[2] );
738 copy_coords( corners_inout[5], cut_face_in[3] );
739 break;
740 case BSPTreeBoxIter::B2376:
741 copy_coords( corners_inout[2], cut_face_in[0] );
742 copy_coords( corners_inout[3], cut_face_in[1] );
743 copy_coords( corners_inout[7], cut_face_in[2] );
744 copy_coords( corners_inout[6], cut_face_in[3] );
745 break;
746 case BSPTreeBoxIter::B3047:
747 copy_coords( corners_inout[3], cut_face_in[0] );
748 copy_coords( corners_inout[0], cut_face_in[1] );
749 copy_coords( corners_inout[4], cut_face_in[2] );
750 copy_coords( corners_inout[7], cut_face_in[3] );
751 break;
752 case BSPTreeBoxIter::B3210:
753 copy_coords( corners_inout[3], cut_face_in[0] );
754 copy_coords( corners_inout[2], cut_face_in[1] );
755 copy_coords( corners_inout[1], cut_face_in[2] );
756 copy_coords( corners_inout[0], cut_face_in[3] );
757 break;
758 case BSPTreeBoxIter::B4567:
759 copy_coords( corners_inout[4], cut_face_in[0] );
760 copy_coords( corners_inout[5], cut_face_in[1] );
761 copy_coords( corners_inout[6], cut_face_in[2] );
762 copy_coords( corners_inout[7], cut_face_in[3] );
763 break;
764 default:
765 return MB_FAILURE; // child is not a box
766 }
767
768 return MB_SUCCESS;
769 }
770
771 ErrorCode BSPTreeBoxIter::step_to_first_leaf( Direction direction )
772 {
773 ErrorCode rval;
774 BSPTree::Plane plane;
775 Corners clipped_corners;
776
777 for( ;; )
778 {
779 childVect.clear();
780 rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
781 if( MB_SUCCESS != rval ) return rval;
782 if( childVect.empty() ) // leaf
783 break;
784
785 rval = tool()->get_split_plane( mStack.back(), plane );
786 if( MB_SUCCESS != rval ) return rval;
787
788 if( direction == RIGHT ) plane.flip();
789 rval = plane_cut_box( clipped_corners.coords, leafCoords, plane );
790 if( MB_SUCCESS != rval ) return rval;
791 mStack.push_back( childVect[direction] );
792 stackData.push_back( clipped_corners );
793 }
794 return MB_SUCCESS;
795 }
796
797 ErrorCode BSPTreeBoxIter::up()
798 {
799 ErrorCode rval;
800 if( mStack.size() == 1 ) return MB_ENTITY_NOT_FOUND;
801
802 EntityHandle node = mStack.back();
803 Corners clipped_face = stackData.back();
804 mStack.pop_back();
805 stackData.pop_back();
806
807 BSPTree::Plane plane;
808 rval = tool()->get_split_plane( mStack.back(), plane );
809 if( MB_SUCCESS != rval )
810 {
811 mStack.push_back( node );
812 stackData.push_back( clipped_face );
813 return rval;
814 }
815
816 rval = plane_uncut_box( clipped_face.coords, leafCoords, plane );
817 if( MB_SUCCESS != rval )
818 {
819 mStack.push_back( node );
820 stackData.push_back( clipped_face );
821 return rval;
822 }
823
824 return MB_SUCCESS;
825 }
826
827 ErrorCode BSPTreeBoxIter::down( const BSPTree::Plane& plane_ref, Direction direction )
828 {
829 childVect.clear();
830 ErrorCode rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
831 if( MB_SUCCESS != rval ) return rval;
832 if( childVect.empty() ) return MB_ENTITY_NOT_FOUND;
833
834 BSPTree::Plane plane( plane_ref );
835 if( direction == RIGHT ) plane.flip();
836
837 Corners clipped_face;
838 rval = plane_cut_box( clipped_face.coords, leafCoords, plane );
839 if( MB_SUCCESS != rval ) return rval;
840
841 mStack.push_back( childVect[direction] );
842 stackData.push_back( clipped_face );
843 return MB_SUCCESS;
844 }
845
846 ErrorCode BSPTreeBoxIter::step( Direction direction )
847 {
848 EntityHandle node, parent;
849 Corners clipped_face;
850 ErrorCode rval;
851 BSPTree::Plane plane;
852 const Direction opposite = static_cast< Direction >( 1 - direction );
853
854 // If stack is empty, then either this iterator is uninitialized
855 // or we reached the end of the iteration (and return
856 // MB_ENTITY_NOT_FOUND) already.
857 if( mStack.empty() ) return MB_FAILURE;
858
859 // Pop the current node from the stack.
860 // The stack should then contain the parent of the current node.
861 // If the stack is empty after this pop, then we've reached the end.
862 node = mStack.back();
863 mStack.pop_back();
864 clipped_face = stackData.back();
865 stackData.pop_back();
866
867 while( !mStack.empty() )
868 {
869 // Get data for parent entity
870 parent = mStack.back();
871 childVect.clear();
872 rval = tool()->moab()->get_child_meshsets( parent, childVect );
873 if( MB_SUCCESS != rval ) return rval;
874 rval = tool()->get_split_plane( parent, plane );
875 if( MB_SUCCESS != rval ) return rval;
876 if( direction == LEFT ) plane.flip();
877
878 // If we're at the left child
879 if( childVect[opposite] == node )
880 {
881 // change from box of left child to box of parent
882 plane_uncut_box( clipped_face.coords, leafCoords, plane );
883 // change from box of parent to box of right child
884 plane.flip();
885 plane_cut_box( clipped_face.coords, leafCoords, plane );
886 // push right child on stack
887 mStack.push_back( childVect[direction] );
888 stackData.push_back( clipped_face );
889 // descend to left-most leaf of the right child
890 return step_to_first_leaf( opposite );
891 }
892
893 // The current node is the right child of the parent,
894 // continue up the tree.
895 assert( childVect[direction] == node );
896 plane.flip();
897 plane_uncut_box( clipped_face.coords, leafCoords, plane );
898 node = parent;
899 clipped_face = stackData.back();
900 mStack.pop_back();
901 stackData.pop_back();
902 }
903
904 return MB_ENTITY_NOT_FOUND;
905 }
906
907 ErrorCode BSPTreeBoxIter::get_box_corners( double coords[8][3] ) const
908 {
909 memcpy( coords, leafCoords, 24 * sizeof( double ) );
910 return MB_SUCCESS;
911 }
912
913 // result = a - b
914 static void subtr( double result[3], const double a[3], const double b[3] )
915 {
916 result[0] = a[0] - b[0];
917 result[1] = a[1] - b[1];
918 result[2] = a[2] - b[2];
919 }
920
921 // result = a + b + c + d
922 static void sum( double result[3], const double a[3], const double b[3], const double c[3], const double d[3] )
923 {
924 result[0] = a[0] + b[0] + c[0] + d[0];
925 result[1] = a[1] + b[1] + c[1] + d[1];
926 result[2] = a[2] + b[2] + c[2] + d[2];
927 }
928
929 // result = a cross b
930 static void cross( double result[3], const double a[3], const double b[3] )
931 {
932 result[0] = a[1] * b[2] - a[2] * b[1];
933 result[1] = a[2] * b[0] - a[0] * b[2];
934 result[2] = a[0] * b[1] - a[1] * b[0];
935 }
936
937 static double dot( const double a[3], const double b[3] )
938 {
939 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
940 }
941
942 double BSPTreeBoxIter::volume() const
943 {
944 // have planar sides, so use mid-face tripple product
945 double f1[3], f2[3], f3[3], f4[3], f5[3], f6[3];
946 sum( f1, leafCoords[0], leafCoords[1], leafCoords[4], leafCoords[5] );
947 sum( f2, leafCoords[1], leafCoords[2], leafCoords[5], leafCoords[6] );
948 sum( f3, leafCoords[2], leafCoords[3], leafCoords[6], leafCoords[7] );
949 sum( f4, leafCoords[0], leafCoords[3], leafCoords[4], leafCoords[7] );
950 sum( f5, leafCoords[0], leafCoords[1], leafCoords[2], leafCoords[3] );
951 sum( f6, leafCoords[4], leafCoords[5], leafCoords[6], leafCoords[7] );
952 double v13[3], v24[3], v65[3];
953 subtr( v13, f1, f3 );
954 subtr( v24, f2, f4 );
955 subtr( v65, f6, f5 );
956 double cr[3];
957 cross( cr, v13, v24 );
958 return ( 1. / 64 ) * dot( cr, v65 );
959 }
960
961 BSPTreeBoxIter::XSect BSPTreeBoxIter::splits( const BSPTree::Plane& plane ) const
962 {
963 // test each corner relative to the plane
964 unsigned result = 0;
965 for( unsigned i = 0; i < 8u; ++i )
966 {
967 double d = plane.signed_distance( leafCoords[i] );
968 // if corner is on plane, than intersection
969 // will result in a degenerate hex
970 if( fabs( d ) < BSPTree::epsilon() ) return NONHEX;
971 // if mark vertices above plane
972 if( d > 0.0 ) result |= 1 << i;
973 }
974
975 switch( result )
976 {
977 // if all vertices or no vertices above plane,
978 // then plane doesn't intersect
979 case 0:
980 case 0xFF:
981 return MISS;
982
983 // if there are four vertices above the plane
984 // and they compose a single face of the hex,
985 // then the cut will result in two hexes
986 case B0154:
987 case B1265:
988 case B2376:
989 case B3047:
990 case B3210:
991 case B4567:
992 return SPLIT;
993
994 // otherwise intersects, but split would not result
995 // in two hexahedrons
996 default:
997 return NONHEX;
998 }
999 }
1000
1001 bool BSPTreeBoxIter::intersects( const BSPTree::Plane& plane ) const
1002 {
1003 // test each corner relative to the plane
1004 unsigned count = 0;
1005 for( unsigned i = 0; i < 8u; ++i )
1006 count += plane.above( leafCoords[i] );
1007 return count > 0 && count < 8u;
1008 }
1009
1010 ErrorCode BSPTreeBoxIter::sibling_side( SideBits& side_out ) const
1011 {
1012 if( mStack.size() < 2 ) // at tree root
1013 return MB_ENTITY_NOT_FOUND;
1014
1015 EntityHandle parent = mStack[mStack.size() - 2];
1016 BSPTree::Plane plane;
1017 ErrorCode rval = tool()->get_split_plane( parent, plane );
1018 if( MB_SUCCESS != rval ) return MB_FAILURE;
1019
1020 side_out = side_on_plane( leafCoords, plane );
1021 return MB_SUCCESS;
1022 }
1023
1024 ErrorCode BSPTreeBoxIter::get_neighbors( SideBits side, std::vector< BSPTreeBoxIter >& results, double epsilon ) const
1025 {
1026 EntityHandle tmp_handle;
1027 BSPTree::Plane plane;
1028 ErrorCode rval;
1029 int n;
1030
1031 Corners face;
1032 rval = face_corners( side, leafCoords, face.coords );
1033 if( MB_SUCCESS != rval ) return rval;
1034
1035 // Move up tree until we find the split that created the specified side.
1036 // Push the sibling at that level onto the iterator stack as
1037 // all neighbors will be rooted at that node.
1038 BSPTreeBoxIter iter( *this ); // temporary iterator (don't modifiy *this)
1039 for( ;; )
1040 {
1041 tmp_handle = iter.handle();
1042
1043 rval = iter.up();
1044 if( MB_SUCCESS != rval ) // reached root - no neighbors on that side
1045 return ( rval == MB_ENTITY_NOT_FOUND ) ? MB_SUCCESS : rval;
1046
1047 iter.childVect.clear();
1048 rval = tool()->moab()->get_child_meshsets( iter.handle(), iter.childVect );
1049 if( MB_SUCCESS != rval ) return rval;
1050
1051 rval = tool()->get_split_plane( iter.handle(), plane );
1052 if( MB_SUCCESS != rval ) return rval;
1053 SideBits s = side_above_plane( iter.leafCoords, plane );
1054
1055 if( tmp_handle == iter.childVect[0] && s == side )
1056 {
1057 rval = iter.down( plane, RIGHT );
1058 if( MB_SUCCESS != rval ) return rval;
1059 break;
1060 }
1061 else if( tmp_handle == iter.childVect[1] && opposite_face( s ) == side )
1062 {
1063 rval = iter.down( plane, LEFT );
1064 if( MB_SUCCESS != rval ) return rval;
1065 break;
1066 }
1067 }
1068
1069 // now move down tree, searching for adjacent boxes
1070 std::vector< BSPTreeBoxIter > list;
1071 // loop over all potential paths to neighbors (until list is empty)
1072 for( ;; )
1073 {
1074 // follow a single path to a leaf, append any other potential
1075 // paths to neighbors to 'list'
1076 for( ;; )
1077 {
1078 rval = tool()->moab()->num_child_meshsets( iter.handle(), &n );
1079 if( MB_SUCCESS != rval ) return rval;
1080
1081 // if leaf
1082 if( !n )
1083 {
1084 results.push_back( iter );
1085 break;
1086 }
1087
1088 rval = tool()->get_split_plane( iter.handle(), plane );
1089 if( MB_SUCCESS != rval ) return rval;
1090
1091 bool some_above = false, some_below = false;
1092 for( int i = 0; i < 4; ++i )
1093 {
1094 double signed_d = plane.signed_distance( face.coords[i] );
1095 if( signed_d > -epsilon ) some_above = true;
1096 if( signed_d < epsilon ) some_below = true;
1097 }
1098
1099 if( some_above && some_below )
1100 {
1101 list.push_back( iter );
1102 list.back().down( plane, RIGHT );
1103 iter.down( plane, LEFT );
1104 }
1105 else if( some_above )
1106 {
1107 iter.down( plane, RIGHT );
1108 }
1109 else if( some_below )
1110 {
1111 iter.down( plane, LEFT );
1112 }
1113 else
1114 {
1115 // tolerance issue -- epsilon to small? 2D box?
1116 return MB_FAILURE;
1117 }
1118 }
1119
1120 if( list.empty() ) break;
1121
1122 iter = list.back();
1123 list.pop_back();
1124 }
1125
1126 return MB_SUCCESS;
1127 }
1128
1129 ErrorCode BSPTreeBoxIter::calculate_polyhedron( BSPTreePoly& poly_out ) const
1130 {
1131 const CartVect* ptr = reinterpret_cast< const CartVect* >( leafCoords );
1132 return poly_out.set( ptr );
1133 }
1134
1135 ErrorCode BSPTree::leaf_containing_point( EntityHandle tree_root, const double point[3], EntityHandle& leaf_out )
1136 {
1137 std::vector< EntityHandle > children;
1138 Plane plane;
1139 EntityHandle node = tree_root;
1140 ErrorCode rval = moab()->get_child_meshsets( node, children );
1141 if( MB_SUCCESS != rval ) return rval;
1142 while( !children.empty() )
1143 {
1144 rval = get_split_plane( node, plane );
1145 if( MB_SUCCESS != rval ) return rval;
1146
1147 node = children[plane.above( point )];
1148 children.clear();
1149 rval = moab()->get_child_meshsets( node, children );
1150 if( MB_SUCCESS != rval ) return rval;
1151 }
1152 leaf_out = node;
1153 return MB_SUCCESS;
1154 }
1155
1156 ErrorCode BSPTree::leaf_containing_point( EntityHandle root, const double point[3], BSPTreeIter& iter )
1157 {
1158 ErrorCode rval;
1159
1160 rval = iter.initialize( this, root, point );
1161 if( MB_SUCCESS != rval ) return rval;
1162
1163 for( ;; )
1164 {
1165 iter.childVect.clear();
1166 rval = moab()->get_child_meshsets( iter.handle(), iter.childVect );
1167 if( MB_SUCCESS != rval || iter.childVect.empty() ) return rval;
1168
1169 Plane plane;
1170 rval = get_split_plane( iter.handle(), plane );
1171 if( MB_SUCCESS != rval ) return rval;
1172
1173 rval = iter.down( plane, ( BSPTreeIter::Direction )( plane.above( point ) ) );
1174 if( MB_SUCCESS != rval ) return rval;
1175 }
1176 }
1177
1178 template < typename PlaneIter >
1179 static inline bool ray_intersect_halfspaces( const CartVect& ray_pt,
1180 const CartVect& ray_dir,
1181 const PlaneIter& begin,
1182 const PlaneIter& end,
1183 double& t_enter,
1184 double& t_exit )
1185 {
1186 const double epsilon = 1e-12;
1187
1188 // begin with inifinite ray
1189 t_enter = 0.0;
1190 t_exit = std::numeric_limits< double >::infinity();
1191
1192 // cut ray such that we keep only the portion contained
1193 // in each halfspace
1194 for( PlaneIter i = begin; i != end; ++i )
1195 {
1196 CartVect norm( i->norm );
1197 double coeff = i->coeff;
1198 double den = norm % ray_dir;
1199 if( fabs( den ) < epsilon )
1200 { // ray is parallel to plane
1201 if( i->above( ray_pt.array() ) ) return false; // ray entirely outside half-space
1202 }
1203 else
1204 {
1205 double t_xsect = ( -coeff - ( norm % ray_pt ) ) / den;
1206 // keep portion of ray/segment below plane
1207 if( den > 0 )
1208 {
1209 if( t_xsect < t_exit ) t_exit = t_xsect;
1210 }
1211 else
1212 {
1213 if( t_xsect > t_enter ) t_enter = t_xsect;
1214 }
1215 }
1216 }
1217
1218 return t_exit >= t_enter;
1219 }
1220
1221 class BoxPlaneIter
1222 {
1223 int faceNum;
1224 BSPTree::Plane facePlanes[6];
1225
1226 public:
1227 BoxPlaneIter( const double coords[8][3] );
1228 BoxPlaneIter() : faceNum( 6 ) {} // initialize to 'end'
1229 const BSPTree::Plane* operator->() const
1230 {
1231 return facePlanes + faceNum;
1232 }
1233 bool operator==( const BoxPlaneIter& other ) const
1234 {
1235 return faceNum == other.faceNum;
1236 }
1237 bool operator!=( const BoxPlaneIter& other ) const
1238 {
1239 return faceNum != other.faceNum;
1240 }
1241 BoxPlaneIter& operator++()
1242 {
1243 ++faceNum;
1244 return *this;
1245 }
1246 };
1247
1248 static const int box_face_corners[6][4] = { { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 3, 7, 6 },
1249 { 3, 0, 4, 7 }, { 3, 2, 1, 0 }, { 4, 5, 6, 7 } };
1250
1251 BoxPlaneIter::BoxPlaneIter( const double coords[8][3] ) : faceNum( 0 )
1252 {
1253 // NOTE: In the case of a BSP tree, all sides of the
1254 // leaf will planar.
1255 assert( sizeof( CartVect ) == sizeof( coords[0] ) );
1256 const CartVect* corners = reinterpret_cast< const CartVect* >( coords );
1257 for( int i = 0; i < 6; ++i )
1258 {
1259 const int* indices = box_face_corners[i];
1260 CartVect v1 = corners[indices[1]] - corners[indices[0]];
1261 CartVect v2 = corners[indices[3]] - corners[indices[0]];
1262 CartVect n = v1 * v2;
1263 facePlanes[i] = BSPTree::Plane( n.array(), -( n % corners[indices[2]] ) );
1264 }
1265 }
1266
1267 bool BSPTreeBoxIter::intersect_ray( const double ray_point[3],
1268 const double ray_vect[3],
1269 double& t_enter,
1270 double& t_exit ) const
1271 {
1272 BoxPlaneIter iter( this->leafCoords ), end;
1273 return ray_intersect_halfspaces( CartVect( ray_point ), CartVect( ray_vect ), iter, end, t_enter, t_exit );
1274 }
1275
1276 class BSPTreePlaneIter
1277 {
1278 BSPTree* toolPtr;
1279 const EntityHandle* const pathToRoot;
1280 int pathPos;
1281 BSPTree::Plane tmpPlane;
1282 std::vector< EntityHandle > tmpChildren;
1283
1284 public:
1285 BSPTreePlaneIter( BSPTree* tool, const EntityHandle* path, int path_len )
1286 : toolPtr( tool ), pathToRoot( path ), pathPos( path_len - 1 )
1287 {
1288 operator++();
1289 }
1290 BSPTreePlaneIter() // initialize to 'end'
1291 : toolPtr( 0 ), pathToRoot( 0 ), pathPos( -1 )
1292 {
1293 }
1294
1295 const BSPTree::Plane* operator->() const
1296 {
1297 return &tmpPlane;
1298 }
1299 bool operator==( const BSPTreePlaneIter& other ) const
1300 {
1301 return pathPos == other.pathPos;
1302 }
1303 bool operator!=( const BSPTreePlaneIter& other ) const
1304 {
1305 return pathPos != other.pathPos;
1306 }
1307 BSPTreePlaneIter& operator++();
1308 };
1309
1310 BSPTreePlaneIter& BSPTreePlaneIter::operator++()
1311 {
1312 if( --pathPos < 0 ) return *this;
1313
1314 EntityHandle prev = pathToRoot[pathPos + 1];
1315 EntityHandle curr = pathToRoot[pathPos];
1316
1317 ErrorCode rval = toolPtr->get_split_plane( curr, tmpPlane );
1318 if( MB_SUCCESS != rval )
1319 {
1320 assert( false );
1321 pathPos = 0;
1322 return *this;
1323 }
1324
1325 tmpChildren.clear();
1326 rval = toolPtr->moab()->get_child_meshsets( curr, tmpChildren );
1327 if( MB_SUCCESS != rval || tmpChildren.size() != 2 )
1328 {
1329 assert( false );
1330 pathPos = 0;
1331 return *this;
1332 }
1333
1334 if( tmpChildren[1] == prev ) tmpPlane.flip();
1335 return *this;
1336 }
1337
1338 bool BSPTreeIter::intersect_ray( const double ray_point[3],
1339 const double ray_vect[3],
1340 double& t_enter,
1341 double& t_exit ) const
1342 {
1343 // intersect with half-spaces defining tree
1344 BSPTreePlaneIter iter1( tool(), &mStack[0], mStack.size() ), end1;
1345 if( !ray_intersect_halfspaces( CartVect( ray_point ), CartVect( ray_vect ), iter1, end1, t_enter, t_exit ) )
1346 return false;
1347
1348 // itersect with box bounding entire tree
1349 double corners[8][3];
1350 ErrorCode rval = tool()->get_tree_box( mStack.front(), corners );
1351 if( MB_SUCCESS != rval )
1352 {
1353 assert( false );
1354 return false;
1355 }
1356
1357 BoxPlaneIter iter2( corners ), end2;
1358 double t2_enter, t2_exit;
1359 if( !ray_intersect_halfspaces( CartVect( ray_point ), CartVect( ray_vect ), iter2, end2, t2_enter, t2_exit ) )
1360 return false;
1361
1362 // if intersect both box and halfspaces, check that
1363 // two intersections overlap
1364 if( t_enter < t2_enter ) t_enter = t2_enter;
1365 if( t_exit > t2_exit ) t_exit = t2_exit;
1366 return t_enter <= t_exit;
1367 }
1368
1369 } // namespace moab
1.9.1.