Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
BSPTreePoly.cpp
Go to the documentation of this file.
1 #include "moab/CartVect.hpp"
2 #include "moab/BSPTreePoly.hpp"
3 #include <cassert>
4 #include <cstdlib>
5 #include <set>
6 
7 #undef DEBUG_IDS
8 
9 namespace moab
10 {
11 
13 {
14  Vertex( const CartVect& v )
15  : CartVect( v ), usePtr( 0 ), markVal( 0 )
16 #ifdef DEBUG_IDS
17  ,
18  id( nextID++ )
19 #endif
20  {
21  }
23  {
24  assert( !usePtr );
25  }
27  int markVal;
28 #ifdef DEBUG_IDS
29  int id;
30  static int nextID;
31 #endif
32 };
33 
35 {
36  VertexUse( Edge* edge, Vertex* vtx );
37  ~VertexUse();
38 
39  void set_vertex( BSPTreePoly::Vertex*& vtx_ptr );
40 
44 };
45 
47 {
48  EdgeUse( Edge* edge );
49  EdgeUse( Edge* edge, Face* face );
50  ~EdgeUse();
51 
55 
56  inline BSPTreePoly::Vertex* start() const;
57  inline BSPTreePoly::Vertex* end() const;
58  int sense() const;
59 
60  void insert_after( BSPTreePoly::EdgeUse* prev );
61  void insert_before( BSPTreePoly::EdgeUse* next );
62 };
63 
65 {
68 #ifdef DEBUG_IDS
69  int id;
70  static int nextID;
71 #endif
72 
73  Edge( Vertex* vstart, Vertex* vend )
74  : forwardPtr( 0 ), reversePtr( 0 )
75 #ifdef DEBUG_IDS
76  ,
77  id( nextID++ )
78 #endif
79  {
80  startPtr = new VertexUse( this, vstart );
81  endPtr = new VertexUse( this, vend );
82  }
83 
84  ~Edge();
85 
87  {
88  return startPtr->vtxPtr;
89  }
91  {
92  return endPtr->vtxPtr;
93  }
94 
96  {
97  return forwardPtr ? forwardPtr->facePtr : 0;
98  }
100  {
101  return reversePtr ? reversePtr->facePtr : 0;
102  }
103 
105  {
106  return ( vtx == startPtr->vtxPtr ) ? startPtr : ( vtx == endPtr->vtxPtr ) ? endPtr : 0;
107  }
109  {
110  return use( about )->nextPtr->edgePtr;
111  }
113  {
114  return use( about )->prevPtr->edgePtr;
115  }
116 
118  {
119  return ( face == forwardPtr->facePtr ) ? forwardPtr : ( face == reversePtr->facePtr ) ? reversePtr : 0;
120  }
122  {
123  return use( about )->nextPtr->edgePtr;
124  }
126  {
127  return use( about )->prevPtr->edgePtr;
128  }
129 
131  {
132  return vuse == startPtr ? endPtr : vuse == endPtr ? startPtr : 0;
133  }
135  {
136  return vuse == forwardPtr ? reversePtr : vuse == reversePtr ? forwardPtr : 0;
137  }
139  {
140  return vtx == startPtr->vtxPtr ? endPtr->vtxPtr : vtx == endPtr->vtxPtr ? startPtr->vtxPtr : 0;
141  }
143  {
144  return start() == eother->start() || start() == eother->end() ? start()
145  : end() == eother->start() || end() == eother->end() ? end()
146  : 0;
147  }
148 
149  int sense( BSPTreePoly::Face* face ) const;
150 
154 };
155 
157 {
158  Face( Face* next )
159  : usePtr( 0 ), nextPtr( next )
160 #ifdef DEBUG_IDS
161  ,
162  id( nextID++ )
163 #endif
164  {
165  }
167  : usePtr( 0 ), nextPtr( 0 )
168 #ifdef DEBUG_IDS
169  ,
170  id( nextID++ )
171 #endif
172  {
173  }
174  ~Face();
177 #ifdef DEBUG_IDS
178  int id;
179  static int nextID;
180 #endif
181  double signed_volume() const;
182 };
183 
184 #ifdef DEBUG_IDS
185 int BSPTreePoly::Vertex::nextID = 1;
186 int BSPTreePoly::Edge::nextID = 1;
187 int BSPTreePoly::Face::nextID = 1;
188 #endif
190 {
191 #ifdef DEBUG_IDS
192  BSPTreePoly::Vertex::nextID = 1;
193  BSPTreePoly::Edge::nextID = 1;
194  BSPTreePoly::Face::nextID = 1;
195 #endif
196 }
197 
198 // static void merge_edges( BSPTreePoly::Edge* keep_edge,
199 // BSPTreePoly::Edge* dead_edge );
200 
201 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge );
202 
204 {
205  if( !vtx->usePtr )
206  {
207  vtx->usePtr = prevPtr = nextPtr = this;
208  return;
209  }
210 
211  nextPtr = vtx->usePtr;
213  assert( prevPtr->nextPtr == nextPtr );
214  nextPtr->prevPtr = this;
215  prevPtr->nextPtr = this;
216 }
217 
219 {
220  if( nextPtr == this )
221  {
222  assert( prevPtr == this );
223  assert( vtxPtr->usePtr == this );
224  vtxPtr->usePtr = 0;
225  delete vtxPtr;
226  }
227  else if( vtxPtr->usePtr == this )
228  vtxPtr->usePtr = nextPtr;
229 
230  nextPtr->prevPtr = prevPtr;
231  prevPtr->nextPtr = nextPtr;
232  nextPtr = prevPtr = 0;
233 }
234 
236 {
237  if( vtxPtr )
238  {
239  if( nextPtr == prevPtr )
240  {
241  assert( nextPtr == this );
242  vtxPtr->usePtr = 0;
243  delete vtx;
244  vtx = 0;
245  }
246  else
247  {
248  nextPtr->prevPtr = prevPtr;
249  prevPtr->nextPtr = nextPtr;
250  if( vtxPtr->usePtr == this ) vtxPtr->usePtr = nextPtr;
251  }
252  }
253 
254  if( vtx )
255  {
256  vtxPtr = vtx;
257  nextPtr = vtxPtr->usePtr->nextPtr;
258  prevPtr = vtxPtr->usePtr;
259  nextPtr->prevPtr = this;
260  vtxPtr->usePtr->nextPtr = this;
261  }
262 }
263 
264 BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge ) : prevPtr( 0 ), nextPtr( 0 ), edgePtr( edge ), facePtr( 0 ) {}
265 
267 {
268  assert( !face->usePtr );
269  face->usePtr = prevPtr = nextPtr = this;
270 
271  if( !face->usePtr )
272  {
273  face->usePtr = prevPtr = nextPtr = this;
274  return;
275  }
276 
277  nextPtr = face->usePtr;
279  assert( prevPtr->nextPtr == nextPtr );
280  nextPtr->prevPtr = this;
281  prevPtr->nextPtr = this;
282 }
283 
285 {
286  // shouldn't already be in a face
287  assert( !facePtr );
288  // adjacent edges should share vertices
289  assert( start() == prev->end() );
290 
291  facePtr = prev->facePtr;
292  nextPtr = prev->nextPtr;
293  prevPtr = prev;
294  nextPtr->prevPtr = this;
295  prevPtr->nextPtr = this;
296 }
297 
299 {
300  // shouldn't already be in a face
301  assert( !facePtr );
302  // adjacent edges should share vertices
303  assert( end() == next->start() );
304 
305  facePtr = next->facePtr;
306  prevPtr = next->prevPtr;
307  nextPtr = next;
308  nextPtr->prevPtr = this;
309  prevPtr->nextPtr = this;
310 }
311 
313 {
314  if( facePtr->usePtr == this ) facePtr->usePtr = ( nextPtr == this ) ? 0 : nextPtr;
315 
316  if( edgePtr->forwardPtr == this ) edgePtr->forwardPtr = 0;
317  if( edgePtr->reversePtr == this ) edgePtr->reversePtr = 0;
318 
319  if( !edgePtr->forwardPtr && !edgePtr->reversePtr ) delete edgePtr;
320 
321  nextPtr->prevPtr = prevPtr;
322  prevPtr->nextPtr = nextPtr;
323  nextPtr = prevPtr = 0;
324 }
325 
327 {
328  if( edgePtr->forwardPtr == this )
329  return 1;
330  else if( edgePtr->reversePtr == this )
331  return -1;
332  else
333  return 0;
334 }
335 
337 {
338  if( edgePtr->forwardPtr == this )
339  return edgePtr->start();
340  else if( edgePtr->reversePtr == this )
341  return edgePtr->end();
342  else
343  return 0;
344 }
345 
347 {
348  if( edgePtr->forwardPtr == this )
349  return edgePtr->end();
350  else if( edgePtr->reversePtr == this )
351  return edgePtr->start();
352  else
353  return 0;
354 }
355 
357 {
358  delete startPtr;
359  delete endPtr;
360  delete forwardPtr;
361  delete reversePtr;
362 }
363 
365 {
366  if( forwardPtr && forwardPtr->facePtr == face )
367  return 1;
368  else if( reversePtr && reversePtr->facePtr == face )
369  return -1;
370  else
371  return 0;
372 }
373 
375 {
376  BSPTreePoly::EdgeUse* nextEdgeUsePtr = usePtr;
377  while( nextEdgeUsePtr )
378  {
379  delete nextEdgeUsePtr; // This is tricky: ~EdgeUse() might change the value of usePtr
380  if( usePtr && usePtr != nextEdgeUsePtr )
381  nextEdgeUsePtr = usePtr;
382  else
383  nextEdgeUsePtr = 0;
384  }
385  usePtr = 0;
386 }
387 
389 {
390  while( faceList )
391  {
392  Face* face = faceList;
394  delete face;
395  }
396 }
397 
399 {
400  clear();
401 
402  Vertex* vertices[8];
403  for( int i = 0; i < 8; ++i )
404  vertices[i] = new Vertex( hex_corners[i] );
405 
406  Edge* edges[12];
407 #ifdef DEBUG_IDS
408  int start_id = Edge::nextID;
409 #endif
410  for( int i = 0; i < 4; ++i )
411  {
412  int j = ( i + 1 ) % 4;
413  edges[i] = new Edge( vertices[i], vertices[j] );
414  edges[i + 4] = new Edge( vertices[i], vertices[i + 4] );
415  edges[i + 8] = new Edge( vertices[i + 4], vertices[j + 4] );
416  }
417 #ifdef DEBUG_IDS
418  for( int i = 0; i < 12; ++i )
419  edges[i]->id = start_id++;
420 #endif
421 
422  static const int face_conn[6][4] = { { 0, 5, -8, -4 }, { 1, 6, -9, -5 }, { 2, 7, -10, -6 },
423  { 3, 4, -11, -7 }, { -3, -2, -1, -12 }, { 8, 9, 10, 11 } };
424  for( int i = 0; i < 6; ++i )
425  {
426  faceList = new Face( faceList );
427  EdgeUse* prev = 0;
428  for( int j = 0; j < 4; ++j )
429  {
430  int e = face_conn[i][j];
431  if( e < 0 )
432  {
433  e = ( -e ) % 12;
434  assert( !edges[e]->reversePtr );
435  if( !prev )
436  {
437  edges[e]->reversePtr = new EdgeUse( edges[e], faceList );
438  }
439  else
440  {
441  edges[e]->reversePtr = new EdgeUse( edges[e] );
442  edges[e]->reversePtr->insert_after( prev );
443  }
444  prev = edges[e]->reversePtr;
445  }
446  else
447  {
448  assert( !edges[e]->forwardPtr );
449  if( !prev )
450  {
451  edges[e]->forwardPtr = new EdgeUse( edges[e], faceList );
452  }
453  else
454  {
455  edges[e]->forwardPtr = new EdgeUse( edges[e] );
456  edges[e]->forwardPtr->insert_after( prev );
457  }
458  prev = edges[e]->forwardPtr;
459  }
460  }
461  }
462 
463  return MB_SUCCESS;
464 }
465 
466 void BSPTreePoly::get_faces( std::vector< const Face* >& face_list ) const
467 {
468  face_list.clear();
469  for( Face* face = faceList; face; face = face->nextPtr )
470  face_list.push_back( face );
471 }
472 
473 void BSPTreePoly::get_vertices( const Face* face, std::vector< CartVect >& vertices ) const
474 {
475  vertices.clear();
476  if( !face || !face->usePtr ) return;
477 
478  EdgeUse* coedge = face->usePtr;
479  do
480  {
481  vertices.push_back( *coedge->end() );
482  coedge = coedge->nextPtr;
483  } while( coedge != face->usePtr );
484 }
485 
487 {
488  CartVect sum( 0.0 );
489  const CartVect* base = usePtr->start();
490  CartVect d1 = ( *usePtr->end() - *base );
491  for( EdgeUse* coedge = usePtr->nextPtr; coedge != usePtr; coedge = coedge->nextPtr )
492  {
493  CartVect d2 = ( *coedge->end() - *base );
494  sum += d1 * d2;
495  d1 = d2;
496  }
497  return ( 1.0 / 6.0 ) * ( sum % *base );
498 }
499 
500 double BSPTreePoly::volume() const
501 {
502  double result = 0;
503  for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
504  result += ptr->signed_volume();
505  return result;
506 }
507 
509 {
510  for( Face* face = faceList; face; face = face->nextPtr )
511  {
512  EdgeUse* edge = face->usePtr;
513  do
514  {
515  edge->edgePtr->start()->markVal = value;
516  edge->edgePtr->end()->markVal = value;
517  edge = edge->nextPtr;
518  } while( edge && edge != face->usePtr );
519  }
520 }
521 /*
522 static void merge_edges( BSPTreePoly::Edge* keep_edge,
523  BSPTreePoly::Edge* dead_edge )
524 {
525  // edges must share a vertex
526  BSPTreePoly::Vertex* dead_vtx = keep_edge->common(dead_edge);
527  assert(dead_vtx);
528  // vertex may have only two adjacent edges
529  BSPTreePoly::VertexUse* dead_vtxuse = dead_edge->use(dead_vtx);
530  assert(dead_vtxuse);
531  BSPTreePoly::VertexUse* keep_vtxuse = dead_vtxuse->nextPtr;
532  assert(keep_vtxuse);
533  assert(keep_vtxuse->edgePtr == keep_edge);
534  assert(keep_vtxuse->nextPtr == dead_vtxuse);
535  assert(keep_vtxuse->prevPtr == dead_vtxuse);
536  assert(dead_vtxuse->prevPtr == keep_vtxuse);
537 
538  // kept edge now ends with the kept vertex on the dead edge
539  keep_vtxuse->set_vertex( dead_edge->other(dead_vtx) );
540 
541  // destructors should take care of everything else
542  // (including removing dead edge from face loops)
543  delete dead_edge;
544 }
545 */
547 {
548  // split edge, creating new edge
549  BSPTreePoly::Edge* new_edge = new BSPTreePoly::Edge( new_vtx, into_edge->end() );
550  into_edge->endPtr->set_vertex( new_vtx ); // This call might delete new_vtx
551 
552  // update coedge loops in faces
553  if( into_edge->forwardPtr )
554  {
555  new_edge->forwardPtr = new BSPTreePoly::EdgeUse( new_edge );
556  new_edge->forwardPtr->insert_after( into_edge->forwardPtr );
557  }
558  if( into_edge->reversePtr )
559  {
560  new_edge->reversePtr = new BSPTreePoly::EdgeUse( new_edge );
561  new_edge->reversePtr->insert_before( into_edge->reversePtr );
562  }
563 
564  return new_edge;
565 }
566 
568 {
569  BSPTreePoly::Face* face = start->facePtr;
570  assert( face == end->facePtr );
571  BSPTreePoly::Face* new_face = new BSPTreePoly::Face;
572  BSPTreePoly::EdgeUse* keep_start = start->prevPtr;
573  BSPTreePoly::EdgeUse* keep_end = end->nextPtr;
574  for( BSPTreePoly::EdgeUse* ptr = start; ptr != keep_end; ptr = ptr->nextPtr )
575  {
576  if( face->usePtr == ptr ) face->usePtr = keep_start;
577  ptr->facePtr = new_face;
578  }
579  new_face->usePtr = start;
580  BSPTreePoly::Edge* edge = new BSPTreePoly::Edge( start->start(), end->end() );
581  edge->forwardPtr = new BSPTreePoly::EdgeUse( edge );
582  edge->reversePtr = new BSPTreePoly::EdgeUse( edge );
583 
584  edge->forwardPtr->facePtr = face;
585  edge->forwardPtr->prevPtr = keep_start;
586  keep_start->nextPtr = edge->forwardPtr;
587  edge->forwardPtr->nextPtr = keep_end;
588  keep_end->prevPtr = edge->forwardPtr;
589 
590  edge->reversePtr->facePtr = new_face;
591  edge->reversePtr->nextPtr = start;
592  start->prevPtr = edge->reversePtr;
593  edge->reversePtr->prevPtr = end;
594  end->nextPtr = edge->reversePtr;
595 
596  return new_face;
597 }
598 
599 bool BSPTreePoly::cut_polyhedron( const CartVect& plane_normal, double plane_coeff )
600 {
601  const double EPSILON = 1e-6; // points this close are considered coincident
602 
603  // scale epsilon rather than normalizing normal vector
604  const double epsilon = EPSILON * ( plane_normal % plane_normal );
605 
606  // Classify all points above/below plane and destroy any faces
607  // that have no vertices below the plane.
608  const int UNKNOWN = 0;
609  const int ABOVE = 1;
610  const int ON = 2;
611  const int BELOW = 3;
612  int num_above = 0;
613  set_vertex_marks( UNKNOWN );
614 
615  // Classify all points above/below plane and
616  // split any edge that intersect the plane.
617  for( Face* face = faceList; face; face = face->nextPtr )
618  {
619  EdgeUse* edge = face->usePtr;
620 
621  do
622  {
623  Vertex* start = edge->edgePtr->start();
624  Vertex* end = edge->edgePtr->end();
625 
626  if( !start->markVal )
627  {
628  double d = plane_normal % *start + plane_coeff;
629  if( d * d <= epsilon )
630  start->markVal = ON;
631  else if( d < 0.0 )
632  start->markVal = BELOW;
633  else
634  {
635  start->markVal = ABOVE;
636  ++num_above;
637  }
638  }
639 
640  if( !end->markVal )
641  {
642  double d = plane_normal % *end + plane_coeff;
643  if( d * d <= epsilon )
644  end->markVal = ON;
645  else if( d < 0.0 )
646  end->markVal = BELOW;
647  else
648  {
649  end->markVal = ABOVE;
650  ++num_above;
651  }
652  }
653 
654  if( ( end->markVal == ABOVE && start->markVal == BELOW ) ||
655  ( end->markVal == BELOW && start->markVal == ABOVE ) )
656  {
657  CartVect dir = *end - *start;
658  double t = -( plane_normal % *start + plane_coeff ) / ( dir % plane_normal );
659  Vertex* new_vtx = new Vertex( *start + t * dir );
660  new_vtx->markVal = ON;
661  split_edge( new_vtx, edge->edgePtr ); // This call might delete new_vtx
662  end = new_vtx;
663  }
664 
665  edge = edge->nextPtr;
666  } while( edge && edge != face->usePtr );
667  }
668 
669  if( !num_above ) return false;
670 
671  // Split faces
672  for( Face* face = faceList; face; face = face->nextPtr )
673  {
674  EdgeUse* edge = face->usePtr;
675 
676  EdgeUse *split_start = 0, *split_end = 0, *other_split = 0;
677  do
678  {
679  if( edge->end()->markVal == ON && edge->start()->markVal != ON )
680  {
681  if( !split_start )
682  split_start = edge->nextPtr;
683  else if( !split_end )
684  split_end = edge;
685  else
686  other_split = edge;
687  }
688 
689  edge = edge->nextPtr;
690  } while( edge && edge != face->usePtr );
691 
692  // If two vertices are on plane (but not every vertex)
693  // then split the face
694  if( split_end && !other_split )
695  {
696  assert( split_start );
697  Face* new_face = split_face( split_start, split_end );
698  new_face->nextPtr = faceList;
699  faceList = new_face;
700  }
701  }
702 
703  // Destroy all faces that are above the plane
704  Face** lptr = &faceList;
705  while( *lptr )
706  {
707  EdgeUse* edge = ( *lptr )->usePtr;
708  bool some_above = false;
709  do
710  {
711  if( edge->start()->markVal == ABOVE )
712  {
713  some_above = true;
714  break;
715  }
716  edge = edge->nextPtr;
717  } while( edge && edge != ( *lptr )->usePtr );
718 
719  if( some_above )
720  {
721  Face* dead = *lptr;
722  *lptr = ( *lptr )->nextPtr;
723  delete dead;
724  }
725  else
726  {
727  lptr = &( ( *lptr )->nextPtr );
728  }
729  }
730 
731  // Construct a new face in the cut plane
732 
733  // First find an edge to start at
734  Edge* edge_ptr = 0;
735  for( Face* face = faceList; face && !edge_ptr; face = face->nextPtr )
736  {
737  EdgeUse* co_edge = face->usePtr;
738  do
739  {
740  if( 0 == co_edge->edgePtr->other( co_edge ) )
741  {
742  edge_ptr = co_edge->edgePtr;
743  break;
744  }
745  co_edge = co_edge->nextPtr;
746  } while( co_edge && co_edge != face->usePtr );
747  }
748  if( !edge_ptr ) return false;
749 
750  // Constuct new face and first CoEdge
751  faceList = new Face( faceList );
752  Vertex *next_vtx, *start_vtx;
753  EdgeUse* prev_coedge;
754  if( edge_ptr->forwardPtr )
755  {
756  next_vtx = edge_ptr->start();
757  start_vtx = edge_ptr->end();
758  assert( !edge_ptr->reversePtr );
759  prev_coedge = edge_ptr->reversePtr = new EdgeUse( edge_ptr, faceList );
760  }
761  else
762  {
763  next_vtx = edge_ptr->end();
764  start_vtx = edge_ptr->start();
765  prev_coedge = edge_ptr->forwardPtr = new EdgeUse( edge_ptr, faceList );
766  }
767 
768  // Construct coedges until loop is closed
769  while( next_vtx != start_vtx )
770  {
771  // find next edge adjacent to vertex with only one adjacent face
772  VertexUse* this_use = edge_ptr->use( next_vtx );
773  VertexUse* use = this_use->nextPtr;
774  while( use != this_use )
775  {
776  if( use->edgePtr->forwardPtr == 0 )
777  {
778  edge_ptr = use->edgePtr;
779  assert( edge_ptr->start() == next_vtx );
780  next_vtx = edge_ptr->end();
781  edge_ptr->forwardPtr = new EdgeUse( edge_ptr );
782  edge_ptr->forwardPtr->insert_after( prev_coedge );
783  prev_coedge = edge_ptr->forwardPtr;
784  break;
785  }
786  else if( use->edgePtr->reversePtr == 0 )
787  {
788  edge_ptr = use->edgePtr;
789  assert( edge_ptr->end() == next_vtx );
790  next_vtx = edge_ptr->start();
791  edge_ptr->reversePtr = new EdgeUse( edge_ptr );
792  edge_ptr->reversePtr->insert_after( prev_coedge );
793  prev_coedge = edge_ptr->reversePtr;
794  break;
795  }
796 
797  use = use->nextPtr;
798  assert( use != this_use ); // failed to close loop!
799  }
800  }
801 
802  return true;
803 }
804 
806 {
807  std::set< Face* > list_faces;
808 
809  int i = 0;
810  for( Face* ptr = faceList; ptr; ptr = ptr->nextPtr )
811  {
812  if( ++i > 10000 ) return false;
813  if( !list_faces.insert( ptr ).second ) return false;
814  }
815 
816  std::set< Vertex* > vertices;
817  for( Face* face = faceList; face; face = face->nextPtr )
818  {
819  i = 0;
820  EdgeUse* coedge = face->usePtr;
821  do
822  {
823  if( ++i > 10000 ) return false;
824 
825  if( coedge->facePtr != face ) return false;
826 
827  Edge* edge = coedge->edgePtr;
828  if( !edge->startPtr || !edge->endPtr ) return false;
829 
830  vertices.insert( edge->start() );
831  vertices.insert( edge->end() );
832 
833  EdgeUse* other;
834  if( edge->forwardPtr == coedge )
835  other = edge->reversePtr;
836  else if( edge->reversePtr != coedge )
837  return false;
838  else
839  other = edge->forwardPtr;
840  if( !other ) return false;
841  if( list_faces.find( other->facePtr ) == list_faces.end() ) return false;
842 
843  EdgeUse* next = coedge->nextPtr;
844  if( next->prevPtr != coedge ) return false;
845  if( coedge->end() != next->start() ) return false;
846 
847  coedge = next;
848  } while( coedge != face->usePtr );
849  }
850 
851  for( std::set< Vertex* >::iterator j = vertices.begin(); j != vertices.end(); ++j )
852  {
853  Vertex* vtx = *j;
854 
855  i = 0;
856  VertexUse* use = vtx->usePtr;
857  do
858  {
859  if( ++i > 10000 ) return false;
860 
861  if( use->vtxPtr != vtx ) return false;
862 
863  Edge* edge = use->edgePtr;
864  if( !edge ) return false;
865  if( edge->startPtr != use && edge->endPtr != use ) return false;
866 
867  VertexUse* next = use->nextPtr;
868  if( next->prevPtr != use ) return false;
869 
870  use = next;
871  } while( use != vtx->usePtr );
872  }
873 
874  return true;
875 }
876 
877 bool BSPTreePoly::is_point_contained( const CartVect& point ) const
878 {
879  if( !faceList ) // empty (zero-dimension) polyhedron
880  return false;
881 
882  const double EPSILON = 1e-6;
883  // Test that point is below the plane of each face
884  // NOTE: This will NOT work for polyhedra w/ concavities
885  for( Face* face = faceList; face; face = face->nextPtr )
886  {
887  Vertex *pt1, *pt2, *pt3;
888  pt1 = face->usePtr->start();
889  pt2 = face->usePtr->end();
890  pt3 = face->usePtr->nextPtr->end();
891 
892  if( pt3 == pt1 ) // degenerate
893  continue;
894 
895  CartVect norm = ( *pt3 - *pt2 ) * ( *pt1 - *pt2 );
896  double coeff = -( norm % *pt2 );
897  if( ( norm % point + coeff ) > EPSILON ) // if above plane, with some -epsilon
898  return false;
899  }
900 
901  return true;
902 }
903 
904 } // namespace moab