 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "moab/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
12 struct BSPTreePoly::Vertex : public CartVect
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 }
22 ~Vertex()
23 {
24 assert( !usePtr );
25 }
26 BSPTreePoly::VertexUse* usePtr;
27 int markVal;
28 #ifdef DEBUG_IDS
29 int id;
30 static int nextID;
31 #endif
32 };
33
34 struct BSPTreePoly::VertexUse
35 {
36 VertexUse( Edge* edge, Vertex* vtx );
37 ~VertexUse();
38
39 void set_vertex( BSPTreePoly::Vertex*& vtx_ptr );
40
41 BSPTreePoly::VertexUse *nextPtr, *prevPtr;
42 BSPTreePoly::Vertex* vtxPtr;
43 BSPTreePoly::Edge* edgePtr;
44 };
45
46 struct BSPTreePoly::EdgeUse
47 {
48 EdgeUse( Edge* edge );
49 EdgeUse( Edge* edge, Face* face );
50 ~EdgeUse();
51
52 BSPTreePoly::EdgeUse *prevPtr, *nextPtr;
53 BSPTreePoly::Edge* edgePtr;
54 BSPTreePoly::Face* facePtr;
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
64 struct BSPTreePoly::Edge
65 {
66 BSPTreePoly::VertexUse *startPtr, *endPtr;
67 BSPTreePoly::EdgeUse *forwardPtr, *reversePtr;
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
86 BSPTreePoly::Vertex* start() const
87 {
88 return startPtr->vtxPtr;
89 }
90 BSPTreePoly::Vertex* end() const
91 {
92 return endPtr->vtxPtr;
93 }
94
95 BSPTreePoly::Face* forward() const
96 {
97 return forwardPtr ? forwardPtr->facePtr : 0;
98 }
99 BSPTreePoly::Face* reverse() const
100 {
101 return reversePtr ? reversePtr->facePtr : 0;
102 }
103
104 BSPTreePoly::VertexUse* use( BSPTreePoly::Vertex* vtx ) const
105 {
106 return ( vtx == startPtr->vtxPtr ) ? startPtr : ( vtx == endPtr->vtxPtr ) ? endPtr : 0;
107 }
108 BSPTreePoly::Edge* next( BSPTreePoly::Vertex* about ) const
109 {
110 return use( about )->nextPtr->edgePtr;
111 }
112 BSPTreePoly::Edge* prev( BSPTreePoly::Vertex* about ) const
113 {
114 return use( about )->prevPtr->edgePtr;
115 }
116
117 BSPTreePoly::EdgeUse* use( BSPTreePoly::Face* face ) const
118 {
119 return ( face == forwardPtr->facePtr ) ? forwardPtr : ( face == reversePtr->facePtr ) ? reversePtr : 0;
120 }
121 BSPTreePoly::Edge* next( BSPTreePoly::Face* about ) const
122 {
123 return use( about )->nextPtr->edgePtr;
124 }
125 BSPTreePoly::Edge* prev( BSPTreePoly::Face* about ) const
126 {
127 return use( about )->prevPtr->edgePtr;
128 }
129
130 BSPTreePoly::VertexUse* other( BSPTreePoly::VertexUse* vuse ) const
131 {
132 return vuse == startPtr ? endPtr : vuse == endPtr ? startPtr : 0;
133 }
134 BSPTreePoly::EdgeUse* other( BSPTreePoly::EdgeUse* vuse ) const
135 {
136 return vuse == forwardPtr ? reversePtr : vuse == reversePtr ? forwardPtr : 0;
137 }
138 BSPTreePoly::Vertex* other( BSPTreePoly::Vertex* vtx ) const
139 {
140 return vtx == startPtr->vtxPtr ? endPtr->vtxPtr : vtx == endPtr->vtxPtr ? startPtr->vtxPtr : 0;
141 }
142 BSPTreePoly::Vertex* common( BSPTreePoly::Edge* eother ) const
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
151 void remove_from_vertex( BSPTreePoly::Vertex*& vtx_ptr );
152 void remove_from_face( BSPTreePoly::Face*& face_ptr );
153 void add_to_vertex( BSPTreePoly::Vertex* vtx_ptr );
154 };
155
156 struct BSPTreePoly::Face
157 {
158 Face( Face* next )
159 : usePtr( 0 ), nextPtr( next )
160 #ifdef DEBUG_IDS
161 ,
162 id( nextID++ )
163 #endif
164 {
165 }
166 Face()
167 : usePtr( 0 ), nextPtr( 0 )
168 #ifdef DEBUG_IDS
169 ,
170 id( nextID++ )
171 #endif
172 {
173 }
174 ~Face();
175 BSPTreePoly::EdgeUse* usePtr;
176 BSPTreePoly::Face* nextPtr;
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
189 void BSPTreePoly::reset_debug_ids()
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
203 BSPTreePoly::VertexUse::VertexUse( BSPTreePoly::Edge* edge, BSPTreePoly::Vertex* vtx ) : vtxPtr( vtx ), edgePtr( edge )
204 {
205 if( !vtx->usePtr )
206 {
207 vtx->usePtr = prevPtr = nextPtr = this;
208 return;
209 }
210
211 nextPtr = vtx->usePtr;
212 prevPtr = nextPtr->prevPtr;
213 assert( prevPtr->nextPtr == nextPtr );
214 nextPtr->prevPtr = this;
215 prevPtr->nextPtr = this;
216 }
217
218 BSPTreePoly::VertexUse::~VertexUse()
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
235 void BSPTreePoly::VertexUse::set_vertex( BSPTreePoly::Vertex*& vtx )
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
266 BSPTreePoly::EdgeUse::EdgeUse( BSPTreePoly::Edge* edge, BSPTreePoly::Face* face ) : edgePtr( edge ), facePtr( face )
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;
278 prevPtr = nextPtr->prevPtr;
279 assert( prevPtr->nextPtr == nextPtr );
280 nextPtr->prevPtr = this;
281 prevPtr->nextPtr = this;
282 }
283
284 void BSPTreePoly::EdgeUse::insert_after( BSPTreePoly::EdgeUse* prev )
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
298 void BSPTreePoly::EdgeUse::insert_before( BSPTreePoly::EdgeUse* next )
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
312 BSPTreePoly::EdgeUse::~EdgeUse()
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
326 int BSPTreePoly::EdgeUse::sense() const
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
336 BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::start() const
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
346 BSPTreePoly::Vertex* BSPTreePoly::EdgeUse::end() const
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
356 BSPTreePoly::Edge::~Edge()
357 {
358 delete startPtr;
359 delete endPtr;
360 delete forwardPtr;
361 delete reversePtr;
362 }
363
364 int BSPTreePoly::Edge::sense( BSPTreePoly::Face* face ) const
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
374 BSPTreePoly::Face::~Face()
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
388 void BSPTreePoly::clear()
389 {
390 while( faceList )
391 {
392 Face* face = faceList;
393 faceList = faceList->nextPtr;
394 delete face;
395 }
396 }
397
398 ErrorCode BSPTreePoly::set( const CartVect hex_corners[8] )
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
486 double BSPTreePoly::Face::signed_volume() const
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
508 void BSPTreePoly::set_vertex_marks( int value )
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 */
546 static BSPTreePoly::Edge* split_edge( BSPTreePoly::Vertex*& new_vtx, BSPTreePoly::Edge* into_edge )
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
567 static BSPTreePoly::Face* split_face( BSPTreePoly::EdgeUse* start, BSPTreePoly::EdgeUse* end )
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
805 bool BSPTreePoly::is_valid() const
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
1.9.1.