Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
BSPTree.cpp
Go to the documentation of this file.
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 ([email protected])
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 
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 
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 
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 
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 
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 
235 {
236  return moab()->get_entities_by_type_and_tag( 0, MBENTITYSET, &rootTag, 0, 1, results );
237 }
238 
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 
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 
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 ) ||
271  {
272  EntityHandle children[] = { left, right };
273  moab()->delete_entities( children, 2 );
274  return MB_FAILURE;
275  }
276 
277  return MB_SUCCESS;
278 }
279 
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 
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 
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 
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 
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 
443 {
444  if( mStack.size() < 2 ) return MB_ENTITY_NOT_FOUND;
445  mStack.pop_back();
446  return MB_SUCCESS;
447 }
448 
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 
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 
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 
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  {
596  break;
602  break;
608  break;
614  break;
620  break;
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  {
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;
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;
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;
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;
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;
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  {
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;
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;
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;
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;
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;
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
1222 {
1223  int faceNum;
1225 
1226  public:
1227  BoxPlaneIter( const double coords[8][3] );
1228  BoxPlaneIter() : faceNum( 6 ) {} // initialize to 'end'
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  }
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 
1277 {
1279  const EntityHandle* const pathToRoot;
1280  int pathPos;
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 
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  }
1308 };
1309 
1311 {
1312  if( --pathPos < 0 ) return *this;
1313 
1314  EntityHandle prev = pathToRoot[pathPos + 1];
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