cgma
FacetProjectTool.cpp
Go to the documentation of this file.
00001 #include "FacetProjectTool.hpp"
00002 #include "PST_Data.hpp"
00003 #include "GeometryDefines.h"
00004 #include "CubitMessage.hpp"
00005 #include "CubitPlane.hpp"
00006 
00007 #include "GfxDebug.hpp" /* for debugging output */
00008 
00009 const int VG_FACET_DEBUG = 145;
00010 const int VG_FACET_BOUNDARY_COLOR = CUBIT_RED_INDEX;
00011 const int VG_FACET_FACET_COLOR = CUBIT_GREEN_INDEX;
00012 const int VG_FACET_IMPRINT_COLOR = CUBIT_WHITE_INDEX;
00013 const int VG_FACET_SEGMENT_COLOR = CUBIT_CYAN_INDEX;
00014 #define VG_FACET_PRINT PRINT_DEBUG(VG_FACET_DEBUG)
00015 
00016 
00017 //-------------------------------------------------------------------------
00018 // Purpose       : Constructor
00019 //
00020 // Special Notes :  
00021 //
00022 // Creator       : Jason Kraftcheck
00023 //
00024 // Creation Date : 05/11/02
00025 //-------------------------------------------------------------------------
00026 FacetProjectTool::FacetProjectTool( ) 
00027 { }
00028 
00029 
00030 static int parent_compare_facets( PST_Face*& f1, PST_Face*& f2 )
00031 {
00032   return (f1->parent < f2->parent) ? -1 :
00033          (f1->parent > f2->parent) ?  1 : 0;
00034 }
00035 static int sequence_compare_pts( PST_Point*& p1, PST_Point*& p2 )
00036 {
00037   return (p1->sequence < p2->sequence) ? -1 :
00038          (p1->sequence > p2->sequence) ?  1 : 0;
00039 }
00040 
00041 //-------------------------------------------------------------------------
00042 // Purpose       : takes a list of edge segments and a surface and returns
00043 //                 the new facets, dudded facets, new points, and new edges
00044 // Special Notes :  
00045 //                 
00046 //
00047 // Creator       : John Fowler
00048 //
00049 // Creation Date : 12/03/02
00050 //-------------------------------------------------------------------------
00051 CubitStatus FacetProjectTool::project(
00052     DLIList<CubitVector*>& segments,                // in
00053     const std::vector<double>& coordinates,     // in
00054     const std::vector<int>& connections,        // in
00055     std::vector<int>& duddedFacets,             // out
00056     std::vector<int>& newFacets,            // out
00057     std::vector<int>& newFacetsIndex,       // out
00058     std::vector<double>& newPoints,             // out
00059     std::vector<int>& edges,                // out
00060     std::vector<int>& segmentPoints,
00061     const double *tolerance_length
00062  )
00063 {
00064   int i;
00065 //   const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS; 
00066                      
00067   assert(connections.size()%3 == 0);  //  For triangles these must be triples
00068   assert(coordinates.size()%3 == 0);  //  Coordinates must come in triples   
00069   DLIList<PST_Edge*> facet_edges;  
00070 
00071     // make the PST edges and faces from the coordinates and connections.
00072   PST_Edge::make_facets( coordinates, connections, GEOMETRY_RESABS, facet_edges );
00073   DLIList<PST_Face*> faces;
00074 
00075   CubitStatus success = CUBIT_SUCCESS;
00076   if(facet_edges.size() > 0)
00077   {
00078 
00079     PST_Edge::faces( facet_edges, faces );
00080     populate_data(faces);
00081 
00082       // Order the orignal points by sequence number.  New points
00083       // will be appended in order.
00084     pointList.sort(&sequence_compare_pts);
00085         
00086       // Do the work
00087     CubitBoolean point_changed;
00088     success = do_projection(segments, point_changed, tolerance_length);
00089     
00090     bool debug = false;
00091     if( debug )
00092     {
00093       GfxDebug::clear();
00094       for( int j=0; j<facetList.size(); j++ )
00095         facetList.get_and_step()->debug_draw( CUBIT_CYAN_INDEX );
00096       GfxDebug::mouse_xforms();
00097     }
00098 
00099       // fill in segmentPoints
00100 //    assert (segPoints.size() == segments.size());
00101     segmentPoints.resize( segPoints.size() );
00102     segPoints.reset();
00103     for ( i = 0; i < segPoints.size(); i++ )
00104     {
00105       PST_Point* pt = segPoints.get_and_step();     
00106       segmentPoints[i] = pt ? pt->sequence : -1;
00107       if (point_changed)
00108         success = CUBIT_FAILURE;
00109     }
00110   }
00111   else
00112     success = CUBIT_FAILURE;
00113   
00114   if (!success)
00115   {
00116      cleanup();
00117      return success;
00118   }
00119                                                                                 
00120     //  Now put the points with sequence > coordinates.size()/3 in newPoints.
00121   int orig_num_points = coordinates.size()/3;
00122   int num_new_points = pointList.size() - orig_num_points;
00123   pointList.reset();
00124   pointList.step(orig_num_points);
00125   newPoints.resize(num_new_points*3);
00126   std::vector<double>::iterator ditor = newPoints.begin();
00127   while( num_new_points-- ) {
00128     PST_Point* pt = pointList.get_and_step();
00129     *ditor++ = pt->x();
00130     *ditor++ = pt->y();
00131     *ditor++ = pt->z();
00132   }
00133 
00134   if ( facetList.size() > 1 ) {  //  If only one facet, skip finding dudded facets
00135     // Sort (group) facets by parent.
00136     facetList.sort(&parent_compare_facets);
00137 
00138     // Fill in the duddedFacets, newFacets, and newFacetsIndex vectors.  
00139 
00140 #ifndef NDEBUG
00141     int orig_num_facets = connections.size()/3;
00142 #endif
00143     int new_facet_index = 0;
00144     DLIList<PST_Point*> facet_pts(3);
00145 
00146     facetList.reset();
00147     duddedFacets.clear();
00148     newFacetsIndex.clear();
00149     newFacets.clear();
00150 
00151     for ( i = facetList.size(); i > 0; ) {
00152       PST_Face* dudded_facet = 0;
00153       if( facetList.get()->parent != facetList.next()->parent ) {
00154 
00155           // unmodified original facet - skip it
00156 #ifndef NDEBUG
00157         PST_Face* facet = facetList.get();
00158         assert(facet->sequence < orig_num_facets &&
00159                facet->sequence == facet->parent);
00160 #endif
00161         facetList.step();
00162         i--;
00163       } 
00164       else {
00165           // new facets
00166         
00167           // put all new facets split from the
00168           // same original facet (having same parent),
00169           // including the orignal facet which was 
00170           // recycled, into newFacets
00171         PST_Face* facet = 0;
00172         do {
00173           facet = facetList.get_and_step();
00174           i--;
00175         
00176           if ( facet->parent == facet->sequence ) {
00177             assert(!dudded_facet);
00178             dudded_facet = facet;
00179           }
00180       
00181           facet_pts.clean_out();
00182           facet->append_points(facet_pts);
00183           facet_pts.reset();
00184           newFacets.push_back(facet_pts.get_and_step()->sequence);
00185           newFacets.push_back(facet_pts.get_and_step()->sequence);
00186           newFacets.push_back(facet_pts.get_and_step()->sequence);
00187           new_facet_index += 3;
00188         } while ( (facet->parent == facetList.get()->parent) && (i > 0) );
00189       
00190           // add replaced facet to duddedFacets
00191         assert(dudded_facet && dudded_facet->sequence < orig_num_facets);
00192         duddedFacets.push_back(dudded_facet->sequence);
00193           // add end position in new_facets for the
00194           // set of replacement facets to newFacetsIndex
00195         newFacetsIndex.push_back(new_facet_index);
00196       }
00197     }  
00198   }
00199   
00200   DLIList<PST_CoEdge*> coedge_list;
00201   edges.clear();
00202   finalize_results();
00203   while( get_result_set(coedge_list) ) {
00204 
00205       // add count and first point
00206     coedge_list.reset();
00207     edges.push_back(coedge_list.size() + 1);
00208     PST_Point* pt = coedge_list.get()->start_point();
00209     edges.push_back(pt->sequence);
00210     
00211       // add all other points
00212     for( i = coedge_list.size(); i--; ) {
00213       pt = coedge_list.get_and_step()->end_point();
00214       edges.push_back(pt->sequence);
00215     }
00216   
00217       // clean up for next iteration
00218     coedge_list.clean_out();
00219   }
00220   edges.push_back(0);  //  terminating zero for next segment length      
00221   
00222   cleanup();
00223   return CUBIT_SUCCESS;
00224 }
00225 
00226  
00227 //-------------------------------------------------------------------------
00228 // Purpose       : fills in the data
00229 //
00230 // Special Notes : populate private data lists and set marks on facet 
00231 //                 entities
00232 //
00233 // Creator       : Jason Kraftcheck, modified by John Fowler --
00234 //         used to be the class constructor
00235 //
00236 // Creation Date : 05/11/02, 12/03/02
00237 //-------------------------------------------------------------------------
00238 CubitStatus FacetProjectTool::populate_data( const DLIList<PST_Face*>& facets )
00239 {
00240   facetList = facets;
00241   int i;
00242   
00243     // Mark all connected edges with a 3, and initialize 
00244     // parent number to sequence number.
00245   facetList.last();
00246   for( i = facetList.size(); i--; )
00247   {
00248     PST_Face* face = facetList.step_and_get();
00249     face->parent = face->sequence;
00250     PST_CoEdge* coe = face->first_coedge();
00251     do
00252     {
00253       coe->edge()->mark = 3;
00254       coe = coe->next();
00255     } while( coe != face->first_coedge() );
00256   }
00257   
00258     // Decrement the edge mark by one for each
00259     // face that contains the edge.
00260   for( i = facetList.size(); i--; )
00261   {
00262     PST_Face* face = facetList.step_and_get();
00263     PST_CoEdge* coe = face->first_coedge();
00264     do
00265     {
00266       coe->edge()->mark--;
00267       coe = coe->next();
00268     } while( coe != face->first_coedge() );
00269   }
00270   
00271     // Populate edgeList and boundaryList with edges
00272   for( i = facetList.size(); i--; )
00273   {
00274     PST_Face* face = facetList.step_and_get();
00275     PST_CoEdge* coe = face->first_coedge();
00276     do
00277     {
00278       PST_Edge* edge = coe->edge();
00279         // If mark is 2, edge was only in one face in the
00280         // passed list, and is thus on the boundary of the
00281         // passed patch of faces.
00282       if( edge->mark == 2 )
00283       {
00284         boundaryList.append( edge );
00285         edgeList.append( edge );
00286         edge->mark = 0;
00287       }
00288         // If mark is 1, the edge is an interior edge in the
00289         // patch of faces
00290       else if( edge->mark == 1 )
00291       {
00292         edgeList.append( edge );
00293         edge->mark = 0;
00294       }
00295         // If the mark on the edge is zero, we have already
00296         // added it to the list(s).
00297       else
00298       {
00299         assert( edge->mark == 0 );  
00300       }
00301       
00302       coe = coe->next();
00303     } while( coe != face->first_coedge() );
00304   }
00305   
00306     // Get all the points
00307   PST_Edge::points( edgeList, pointList );
00308   
00309     // Clear marks on all faces connected to our edge
00310     // list (the passed faces and any faces connected
00311     // at the boundary edges.)
00312   for( i = edgeList.size(); i--; )
00313   {
00314     PST_Edge* edge = edgeList.step_and_get();
00315     edge->mark = 0;
00316     if( edge->forward()->face() )
00317       edge->forward()->face()->mark = 0;
00318     if( edge->reverse()->face() )
00319       edge->reverse()->face()->mark = 0;
00320   }
00321 
00322     // clear marks on faces
00323   for( i = facetList.size(); i--; )
00324     facetList.step_and_get()->mark = 0;
00325     // Set mark to 1 on boundary edges
00326   for( i = boundaryList.size(); i--; )
00327     boundaryList.step_and_get()->mark = 1;
00328     
00329   if( DEBUG_FLAG( VG_FACET_DEBUG ) )
00330   {
00331 //    PST_Face::draw_faces( facetList, VG_FACET_FACET_COLOR, false );
00332     PST_Edge::debug_draw_edges( boundaryList, VG_FACET_BOUNDARY_COLOR );
00333   }
00334 
00335   return CUBIT_SUCCESS;
00336 } 
00337  
00338 //-------------------------------------------------------------------------
00339 // Purpose       : Find the face for which the projection of the passed
00340 //                 position onto the face is closest to the passed
00341 //                 position. 
00342 //
00343 // Special Notes : 
00344 //
00345 // Creator       : Jason Kraftcheck
00346 //
00347 // Creation Date : 05/11/02
00348 //-------------------------------------------------------------------------
00349 PST_Face* FacetProjectTool::closest_face( const CubitVector& position )
00350 {
00351   PST_Face* result = 0;
00352   double dist_sqr = 0;
00353   CubitVector projected;
00354   
00355   for( int i = facetList.size(); i--; )
00356   {
00357     PST_Face* facet = facetList.step_and_get();
00358     if( project_to_face( facet, position, projected ) )
00359     {
00360       double d = (projected - position).length_squared();
00361       if( !result || d < dist_sqr )
00362       {
00363         result = facet;
00364         dist_sqr = d;
00365       }
00366     }
00367   }
00368   
00369   return result;
00370 }
00371 
00372 
00373 //-------------------------------------------------------------------------
00374 // Purpose       : Find the edge for which the projection of the passed
00375 //                 position onto the edge is closest to the passed
00376 //                 position.
00377 //
00378 // Special Notes : 
00379 //
00380 // Creator       : Jason Kraftcheck
00381 //
00382 // Creation Date : 05/11/02
00383 //-------------------------------------------------------------------------
00384 PST_Edge* FacetProjectTool::closest_edge( const CubitVector& position )
00385 {
00386   PST_Edge* result = 0;
00387   double dist_sqr = 0;
00388   
00389   for( int i = edgeList.size(); i--; )
00390   {
00391     PST_Edge* edge = edgeList.step_and_get();
00392     
00393     double t = edge->closest_on_line( position );
00394     if( (t > -CUBIT_RESABS) && (t < (1.0 + CUBIT_RESABS)) )
00395     {
00396       double d = (edge->position(t) - position).length_squared();
00397       if( !result || d < dist_sqr )
00398       {
00399         result = edge;
00400         dist_sqr = d;
00401       }
00402     }
00403   }
00404 
00405   return result;
00406 }
00407 
00408 //-------------------------------------------------------------------------
00409 // Purpose       : Find the point closest to the passed position.
00410 //
00411 // Special Notes : 
00412 //
00413 // Creator       : Jason Kraftcheck
00414 //
00415 // Creation Date : 05/11/02
00416 //-------------------------------------------------------------------------
00417 PST_Point* FacetProjectTool::closest_point( const CubitVector& position )
00418 {
00419   PST_Point* result = 0;
00420   double dist_sqr = 0;
00421   
00422   for( int i = pointList.size(); i--; )
00423   {
00424     PST_Point* pt = pointList.step_and_get();
00425     double d = (position - *pt).length_squared();
00426     if( !result || d < dist_sqr )
00427     {
00428       dist_sqr = d;
00429       result = pt;
00430     }
00431   }
00432   
00433   return result;
00434 }
00435 
00436 
00437 //-------------------------------------------------------------------------
00438 // Purpose       : Project the passed position onto the plane of the 
00439 //                 passed facet.  Returns false if projection is outside
00440 //                 the facet.
00441 //
00442 // Special Notes : "result" is unchanged if false is returned.
00443 //
00444 // Creator       : Jason Kraftcheck
00445 //
00446 // Creation Date : 05/11/02
00447 //-------------------------------------------------------------------------
00448 bool FacetProjectTool::project_to_face( PST_Face* triangle,
00449                                         const CubitVector& position,
00450                                         CubitVector& result )
00451 {
00452     // First check if projection of the point will be within
00453     // the bounds of the triangle.  
00454   if( ! inside_face( triangle, position ) )
00455     return false;
00456   
00457     // Calculate the actual projection
00458   result = triangle->plane().project( position );
00459   return true;
00460 }
00461 
00462 //-------------------------------------------------------------------------
00463 // Purpose       : Destructor.  Clear marks on all facet entities.
00464 //
00465 // Special Notes : 
00466 //
00467 // Creator       : Jason Kraftcheck
00468 //
00469 // Creation Date : 05/11/02
00470 //-------------------------------------------------------------------------
00471 FacetProjectTool::~FacetProjectTool()
00472 {
00473   cleanup();
00474 }
00475 
00476 
00477 //-------------------------------------------------------------------------
00478 // Purpose       : clean out internal data
00479 //
00480 // Special Notes : 
00481 //
00482 // Creator       : Jason Kraftcheck
00483 //
00484 // Creation Date : 02/13/03
00485 //-------------------------------------------------------------------------
00486 void FacetProjectTool::cleanup()
00487 {
00488   while( pointList.size() )
00489     delete pointList.pop();
00490   facetList.clean_out();
00491   edgeList.clean_out();
00492   boundaryList.clean_out();
00493   imprintResult.clean_out();
00494   segPoints.clean_out();
00495 }
00496 
00497 
00498 //-------------------------------------------------------------------------
00499 // Purpose       : Inset a point in the interior of a facet.
00500 //
00501 // Special Notes : Assumes input facet is a triangle, and splits facet
00502 //                 as necessary to ensure that resulting facets are tris.
00503 //
00504 // Creator       : Jason Kraftcheck
00505 //
00506 // Creation Date : 05/11/02
00507 //-------------------------------------------------------------------------
00508 PST_Point* FacetProjectTool::insert_in_face( PST_Face* face, 
00509                                              const CubitVector& position )
00510 {
00511   PST_Point* new_point = new PST_Point( position );
00512   PST_CoEdge* ce1 = face->first_coedge();
00513   PST_CoEdge* ce2 = ce1->next();
00514   PST_CoEdge* ce3 = ce2->next();
00515 
00516     // insert non-mainifold edge in face after ce1
00517   PST_Edge* edge1 = PST_Edge::insert_in_face( new_point, ce1 );
00518   if( ! edge1 )
00519   {
00520     delete new_point;
00521     return 0;
00522   }
00523   new_point->sequence = pointList.size();
00524   pointList.append(new_point);
00525   edgeList.append( edge1 );
00526   if( DEBUG_FLAG( VG_FACET_DEBUG ) )
00527     edge1->debug_draw( VG_FACET_FACET_COLOR );
00528   
00529     // connect non-manifold edge to opposite side of the face
00530     // to make manifold topology.
00531   PST_Edge* edge2 = PST_Edge::split_face( ce2, edge1->start_point() == new_point ?
00532                                                edge1->reverse() : edge1->forward() );   
00533   if( ! edge2 ) return new_point;
00534   edgeList.append( edge2 );
00535   PST_Face* new_face = edge2->other(face);
00536   new_face->sequence = facetList.size();
00537   new_face->parent = face->parent;
00538   facetList.append( new_face );
00539   if( face->owner() )
00540     face->owner()->notify_split( face, new_face );
00541     
00542   if( DEBUG_FLAG( VG_FACET_DEBUG ) )
00543     edge2->debug_draw( VG_FACET_FACET_COLOR );
00544   
00545     // Split face so that all faces are triangles.
00546   PST_Face* split_face = ce3->face();
00547   PST_Edge* edge3 = PST_Edge::split_face( new_point, ce3->end_point(), split_face );
00548   if( edge3 )
00549   {
00550     edgeList.append( edge3 );
00551     new_face = edge3->other(split_face);
00552     new_face->sequence = facetList.size();
00553     new_face->parent = face->parent;
00554     facetList.append( new_face );
00555     if( split_face->owner() )
00556       split_face->owner()->notify_split( split_face, new_face );
00557   }
00558   
00559   return new_point;
00560 }
00561 
00562 
00563 //-------------------------------------------------------------------------
00564 // Purpose       : Add a point where the passed position is projected
00565 //                 into the facet patch.
00566 //
00567 // Special Notes : May return NULL if projection does not lie within 
00568 //                 facet patch.
00569 //
00570 // Creator       : Jason Kraftcheck
00571 //
00572 // Creation Date : 05/11/02
00573 //-------------------------------------------------------------------------
00574 PST_Point* FacetProjectTool::project( const CubitVector& pos,
00575                                       const double* tolerance_length )
00576 {
00577   PST_Point* point = closest_point( pos );
00578   PST_Edge *  edge = closest_edge ( pos );
00579   PST_Face *  face = closest_face ( pos );
00580   
00581     // find closest facet
00582   double face_dist = CUBIT_DBL_MAX;
00583   CubitVector face_pos;
00584   if( face )
00585   {
00586     project_to_face( face, pos, face_pos );
00587     face_dist = (pos - face_pos).length_squared();
00588   }
00589   
00590     // find closest edge
00591   double edge_dist = CUBIT_DBL_MAX;
00592   double edge_pos  = CUBIT_DBL_MAX;
00593   if( edge )
00594   {
00595     edge_pos = edge->closest_on_edge( pos );
00596     edge_dist = (pos - edge->position(edge_pos)).length_squared();
00597   }
00598   
00599     // find closest point
00600   double point_dist = CUBIT_DBL_MAX;
00601   if( point )
00602   {
00603     point_dist = (pos - *point).length_squared();
00604   }
00605   
00606     // if the position is closer to a facet than a point or edge...
00607   if( face && (face_dist < point_dist) && (face_dist < edge_dist) )
00608   {
00609       // get a tolerance based on a sample face size. I attempted
00610       // to do this with a representative facet size but still 
00611       // ended up with problems.  While this is slower is more
00612       // robust.  KGM
00613     double facet_length;
00614 
00615     if(tolerance_length)
00616     {
00617       facet_length = *tolerance_length;
00618     }
00619     else
00620     {
00621       facet_length = face->bounding_length();
00622     }
00623 
00624     // use the square of the facet length times an arbitrary factor to
00625     // determine the tolerance
00626     double tol_sqr = facet_length*facet_length*.00001;
00627 
00628       // check if projected position is within tolerance of
00629       // any bounding edge of the face.
00630     PST_CoEdge* coe = face->first_coedge();
00631     bool insert = true;
00632     do {
00633       PST_Edge* e = coe->edge();
00634       CubitVector p = e->position( e->closest_on_line( face_pos ) );
00635       double d = (p - face_pos).length_squared();
00636 
00637       // use a relative tolerance here - KGM
00638       if( d < tol_sqr ) 
00639       {
00640         edge = e;
00641         edge_pos = edge->closest_on_edge( pos );
00642         face_dist = d;
00643         edge_dist = d;
00644         insert = false;
00645       }
00646       coe = coe->next();
00647     } while( coe != face->first_coedge() );
00648     
00649       // If projected position was sufficiently far from the edges
00650       // of the facet, insert an interior point in the facet.  Otherwise
00651       // fall through and split edge.
00652     if( insert )
00653     {
00654       return insert_in_face( face, face_pos );
00655     }
00656   }
00657   
00658     // If the point was closer to an edge than a point
00659   if( edge && (edge_dist <= point_dist) )
00660   {
00661       // If within tolerance of an end point of the edge, return
00662       // the end point.
00663     CubitVector pt = edge->position( edge_pos );
00664     if( (pt - *edge->start_point()).length_squared() < GEOMETRY_RESABS*GEOMETRY_RESABS )
00665       point = edge->start_point();
00666     else if( (pt - *edge->end_point()).length_squared() < GEOMETRY_RESABS*GEOMETRY_RESABS )
00667       point = edge->end_point();
00668       // Otherwise split the edge at the projected position
00669     else point = split_edge( edge, edge_pos );
00670   }
00671   
00672   if( DEBUG_FLAG( VG_FACET_DEBUG ) && point )
00673     point->debug_draw( VG_FACET_IMPRINT_COLOR );
00674   
00675     // If this was not set in the above if-block, then it is either
00676     // the closest point in the facet patch or NULL if the position
00677     // projected outside the facet patch.
00678   return point;
00679 }    
00680   
00681 //-------------------------------------------------------------------------
00682 // Purpose       : Split an edge at the specified parameter value on the
00683 //                 edge.  Splits connected facets to maintain triangular
00684 //                 facets.
00685 //
00686 // Special Notes : 
00687 //
00688 // Creator       : Jason Kraftcheck
00689 //
00690 // Creation Date : 05/11/02
00691 //-------------------------------------------------------------------------
00692 PST_Point* FacetProjectTool::split_edge( PST_Edge* edge, double t )
00693 {
00694   assert( t > 0.0 && t < 1.0 );
00695   
00696     // create point
00697   PST_Point* point = new PST_Point( edge->position(t) );
00698   
00699     // get connected faces 
00700   PST_Face* face1 = edge->forward()->face();
00701   PST_Face* face2 = edge->reverse()->face();
00702   
00703     // point on each triangle opposite the edge to be split
00704   PST_Point* face1_pt = face1 ? edge->forward()->next()->end_point() : 0;
00705   PST_Point* face2_pt = face2 ? edge->reverse()->next()->end_point() : 0;
00706 
00707   PST_Edge* face1_edge = 0, *face2_edge = 0;
00708   
00709     // split edge edge
00710   PST_Edge* split_edge = edge->split( point );
00711   if( edge->owner() )
00712     edge->owner()->notify_split( edge, split_edge );
00713     // if we split a boundary edge, put the new edge in the
00714     // list of boundary edges too.
00715   if( edge->mark )
00716   {
00717     split_edge->mark = edge->mark;
00718     boundaryList.append( split_edge );
00719   }
00720     // add edge to our list of edges
00721   edgeList.append( split_edge );
00722     // If the first face was not a boundary face (different
00723     // than the boundary of our patch of facets), then split
00724     // it so that it remains a triangle.
00725   if( face1 ) 
00726   {
00727     face1_edge = PST_Edge::split_face( point, face1_pt, face1 );
00728     
00729     if( face1_edge )
00730     {
00731       if( DEBUG_FLAG( VG_FACET_DEBUG ) )
00732         face1_edge->debug_draw( VG_FACET_FACET_COLOR );
00733 
00734       edgeList.append( face1_edge );
00735         PST_Face* new_face = face1_edge->other(face1);
00736       new_face->sequence = facetList.size();
00737       new_face->parent = face1->parent;
00738       facetList.append(new_face );
00739       
00740       if( face1->owner() )
00741         face1->owner()->notify_split( face1, face1_edge->other(face1) );
00742     }
00743   }
00744     // If the second face is not a boundary face, split it
00745     // so that we still have triangles.
00746   if( face2 )
00747   {
00748     face2_edge = PST_Edge::split_face( point, face2_pt, face2 );
00749     if( face2_edge )
00750     {
00751       if( DEBUG_FLAG( VG_FACET_DEBUG ) )
00752         face2_edge->debug_draw( VG_FACET_FACET_COLOR );
00753 
00754       edgeList.append( face2_edge );
00755         PST_Face* new_face = face2_edge->other(face2);
00756       new_face->sequence = facetList.size();
00757       new_face->parent = face2->parent;
00758       facetList.append( new_face );
00759       
00760       if( face2->owner() )
00761         face2->owner()->notify_split( face2, face2_edge->other(face2) );
00762     }
00763   }
00764   point->sequence = pointList.size();
00765   pointList.append(point);
00766   
00767   return point;
00768 }
00769 
00770 //-------------------------------------------------------------------------
00771 // Purpose       : This is the main or outer-loop function for the
00772 //                 polyline projection.  It is called after the input data
00773 //                 is converted to the working representation.  The passed
00774 //                 list of CubitVectors is the polyline.
00775 //
00776 // Special Notes : 
00777 //
00778 // Creator       : Jason Kraftcheck
00779 //
00780 // Creation Date : ??
00781 //-------------------------------------------------------------------------
00782 
00783 CubitStatus 
00784 FacetProjectTool::do_projection( const DLIList<CubitVector*>& passed_list, 
00785                                  CubitBoolean & point_changed,
00786                                  const double *tolerance_length)
00787 {
00788   const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;
00789   point_changed = CUBIT_FALSE;
00790   DLIList<CubitVector*> segments( passed_list );  
00791   
00792   bool first_and_last_points_coincident = false;  
00793   if( segments[0]->distance_between_squared( *(segments[ segments.size()-1]) ) < TOL_SQR )
00794     first_and_last_points_coincident = true; 
00795 
00796   segments.reset();
00797   CubitVector end = *(segments.get_and_step());
00798   PST_Point* last_point = 0;
00799 
00800   bool did_any_projection = false;
00801   
00802   int where_is_last_pt;
00803   int where_is_end_pt;
00804   
00805   bool debug = false;  
00806   
00807   for( int i = segments.size(); i > 1; i-- )
00808   {
00809     CubitVector start(end);
00810     end = *(segments.get_and_step());
00811     
00812     if( debug )
00813     {
00814       GfxDebug::clear();       
00815       for( int j=facetList.size(); j--; )
00816         facetList.get_and_step()->debug_draw( CUBIT_CYAN_INDEX+j, false );        
00817       GfxDebug::draw_point( start, CUBIT_GREEN_INDEX );
00818       GfxDebug::draw_point( end, CUBIT_RED_INDEX );
00819       GfxDebug::flush();
00820       GfxDebug::mouse_xforms();      
00821     }
00822 
00823     if( DEBUG_FLAG( VG_FACET_DEBUG ) )
00824     {
00825       GfxDebug::draw_line( 
00826         (float)start.x(), (float)start.y(), (float)start.z(),
00827         (float)end.x()  , (float)end.y()  , (float)end.z(),
00828         VG_FACET_SEGMENT_COLOR );     
00829       GfxDebug::flush();
00830       GfxDebug::mouse_xforms();
00831     }
00832 
00833     double facet_dist_tol = start.distance_between_squared( end );
00834     facet_dist_tol *= 0.0001;
00835 
00836     //determine if start is on/off/on-boundary of facets
00837     if( !last_point )
00838     {
00839       where_is_last_pt = point_on_facets( &start, facet_dist_tol );
00840 
00841       if( 1 == where_is_last_pt || 
00842           2 == where_is_last_pt  )
00843       {
00844         last_point = project( start, tolerance_length );
00845         
00846         segPoints.append(last_point); 
00847         start = *last_point;
00848 
00849         if( first_and_last_points_coincident )
00850           *(segments[segments.size()-1]) = start;
00851       }
00852       else
00853         segPoints.append(NULL);  //it's not on the facets
00854     } 
00855 
00856     where_is_end_pt = point_on_facets( &end, facet_dist_tol );
00857 
00858     //both points off the facets
00859     //or one-on, one-off
00860     PST_Point *end_pt = NULL;
00861 
00862     //both points off
00863     //one-point-on && one-point-off
00864     if( (0 == where_is_last_pt && 0 == where_is_end_pt ) || 
00865         (0 == where_is_last_pt && 1 == where_is_end_pt ) ||
00866         (1 == where_is_last_pt && 0 == where_is_end_pt ) )
00867     {
00868       //see if segment intersects boundary
00869       CubitVector int_pt;
00870       CubitStatus my_status = find_closest_int_point_with_boundary( start, end, int_pt );
00871       if( CUBIT_FAILURE == my_status )      
00872       {
00873         segPoints.append( NULL );
00874         continue;
00875       }      
00876 
00877       if( 0 == where_is_last_pt )
00878       {
00879         last_point = project( int_pt, tolerance_length );
00880 
00881         if( 0 == where_is_end_pt )
00882         {
00883           my_status = find_closest_int_point_with_boundary( end, *last_point, int_pt );
00884           if( CUBIT_SUCCESS == my_status )
00885           {
00886             end_pt = project( int_pt, tolerance_length );
00887             i++;
00888             segments.back();
00889           }          
00890         } 
00891         else
00892         {
00893           end_pt = project( end, tolerance_length );
00894           segPoints.append( end_pt );
00895         }
00896       }
00897       else if( 0 == where_is_end_pt )
00898       {
00899         end_pt = project( int_pt, tolerance_length );
00900         i++;
00901         segments.back();
00902       }      
00903     }
00904     else if( (1 == where_is_last_pt || 2 == where_is_last_pt) && 
00905              (1 == where_is_end_pt || 2 == where_is_end_pt) )
00906     {
00907       end_pt = project( end, tolerance_length );
00908       segPoints.append( end_pt );
00909     }         
00910     else if( (0 == where_is_last_pt && 2 == where_is_end_pt ) || 
00911              (2 == where_is_last_pt && 0 == where_is_end_pt ) )
00912     {
00913       if( 0 == where_is_last_pt && 2 == where_is_end_pt )
00914       {
00915         //see if segment intersects boundary
00916         CubitVector int_pt;
00917         CubitStatus my_status = find_closest_int_point_with_boundary( start, end, int_pt );
00918         if( CUBIT_SUCCESS == my_status )
00919         {
00920           end_pt = project( end, tolerance_length );
00921           last_point = project( int_pt, tolerance_length );
00922           segPoints.append( end_pt );
00923         }
00924         else
00925         {
00926           segPoints.append(NULL);
00927           continue;
00928         }
00929       }
00930       else
00931       {
00932         //see if segment intersects boundary
00933         CubitVector int_pt;
00934         CubitStatus my_status = find_closest_int_point_with_boundary( end, start, int_pt );      
00935 
00936         if( CUBIT_SUCCESS == my_status )
00937         {
00938           end_pt = project( int_pt, tolerance_length );                        
00939           segPoints.append( NULL );
00940         }
00941         else
00942         {
00943           segPoints.append( NULL );
00944           continue;
00945         }
00946       }
00947     }
00948 
00949         
00950     PST_Point* start_point = last_point;
00951     PST_Point* end_point = end_pt;
00952 
00953     if( debug )
00954     {
00955       GfxDebug::draw_line( *last_point, *end_pt, CUBIT_MAGENTA_INDEX );
00956       GfxDebug::flush();
00957       GfxDebug::mouse_xforms();
00958     }
00959 
00960     if(project( last_point, end_pt ) == CUBIT_SUCCESS)
00961     {
00962       did_any_projection = true;
00963       if ((*start_point - *last_point).length_squared() > TOL_SQR)
00964       {
00965         segPoints.move_to(start_point);
00966         segPoints.change_to(last_point);
00967         point_changed = CUBIT_TRUE;
00968       }
00969       if ((*end_pt - *end_point).length_squared() > TOL_SQR)
00970       {
00971         segPoints.move_to(end_point);
00972         segPoints.change_to(end_pt);
00973         point_changed = CUBIT_TRUE;
00974       }    
00975     }
00976 
00977     last_point = end_pt;
00978     where_is_last_pt = where_is_end_pt;    
00979   
00980   } // end for( i = segments )
00981   
00982   if( !did_any_projection )
00983     return CUBIT_FAILURE;
00984 
00985   return CUBIT_SUCCESS;
00986 }
00987 
00988 //-------------------------------------------------------------------------
00989 // Purpose       : Determines if a position is off/on/on-boundary of the
00990 //                 facets in this tool.  For off/on/on-boundary returns
00991 //                 0, 1, and 2 respectively.  The passed in tolerance_sq is 
00992 //                 used as a coincidence tolerance.
00993 //
00994 // Special Notes : 
00995 //
00996 // Creator       : Corey Ernst
00997 //
00998 // Creation Date : 24-April-2013
00999 //-------------------------------------------------------------------------
01000 
01001 int FacetProjectTool::point_on_facets( 
01002   CubitVector *pos,  
01003   double tolerance_sq )
01004 {
01005   //Return Values:
01006   //0 is outside of facets
01007   //1 is inside of facets
01008   //2 is on boundary of facets
01009 
01010   PST_Point* point = closest_point( *pos );
01011   PST_Edge *  edge = closest_edge ( *pos );
01012   PST_Face *  face = closest_face ( *pos );
01013 
01014   // find closest facet
01015   double face_dist = CUBIT_DBL_MAX;
01016   CubitVector face_pos;
01017   if( face )
01018   {
01019     project_to_face( face, *pos, face_pos );
01020     face_dist = (*pos - face_pos).length_squared();
01021   }
01022 
01023   // find closest edge
01024   double edge_dist = CUBIT_DBL_MAX;
01025   double edge_pos  = CUBIT_DBL_MAX;
01026   if( edge )
01027   {
01028     edge_pos = edge->closest_on_edge( *pos );
01029     edge_dist = (*pos - edge->position(edge_pos)).length_squared();
01030   }
01031 
01032   // find closest point
01033   double point_dist = CUBIT_DBL_MAX;
01034   if( point )
01035   {
01036     point_dist = (*pos - *point).length_squared();
01037   }  
01038 
01039   // if the position is closer to a facet than a point or edge...
01040   if( face && (face_dist < point_dist) && (face_dist < edge_dist) )
01041   {
01042     // check if projected position is within tolerance of
01043     // any bounding edge of the face.
01044     PST_CoEdge* coe = face->first_coedge();
01045     bool edge_closer = false;
01046     do {
01047       PST_Edge* e = coe->edge();
01048       CubitVector p = e->position( e->closest_on_line( face_pos ) );
01049       double d = (p - face_pos).length_squared();
01050 
01051       if( d < tolerance_sq ) 
01052       {
01053         edge = e;
01054         edge_pos = edge->closest_on_edge( *pos );
01055         face_dist = d;
01056         edge_dist = d;
01057         edge_closer = true;
01058       }
01059       coe = coe->next();
01060     } while( coe != face->first_coedge() );
01061     
01062     if( edge_closer )
01063     {      
01064       if( boundaryList.is_in_list( edge ) )
01065         return 2;
01066       else
01067         return 1;    
01068     }      
01069     else
01070       return 1;    
01071     
01072     return 0;
01073   }
01074 
01075   // If the point was closer to an edge than a point
01076   if( edge && (edge_dist < point_dist) )
01077   {
01078     if( edge_dist < tolerance_sq )
01079     {
01080       if( boundaryList.is_in_list( edge ) )
01081         return 2;
01082       else
01083         return 1;
01084     }
01085     else
01086       return 0;
01087   }
01088 
01089   PST_Edge *first_edge = point->edge();
01090   PST_Edge *tmp_edge = first_edge;
01091   if( point_dist < tolerance_sq )
01092   {
01093     do
01094     {
01095       if( boundaryList.is_in_list( tmp_edge ) )
01096         return 2;
01097       tmp_edge = point->next( tmp_edge );
01098     }
01099     while( first_edge != tmp_edge );
01100 
01101     return 1;
01102   }
01103   return 0;
01104 }
01105 
01106 
01107 //-------------------------------------------------------------------------
01108 // Purpose       : Finds the closest intersection point to 'start_pt' 
01109 //                 between the boundary facet edges in this tool and a line
01110 //                 segement defined by 'start_pt' and 'end_pt'.  The 
01111 //                 intersection must be at or between the two end points
01112 //                 of the segment.
01113 //
01114 // Special Notes : 
01115 //
01116 // Creator       : Corey Ernst
01117 //
01118 // Creation Date : 24-April-2013
01119 //-------------------------------------------------------------------------
01120 CubitStatus FacetProjectTool::find_closest_int_point_with_boundary( 
01121   CubitVector &start_pt, CubitVector &end_pt, CubitVector &closest_pt  )
01122 {
01123   CubitVector seg_direction = end_pt - start_pt;   
01124   double seg_length_sq = seg_direction.length_squared();
01125   double shortest_length = CUBIT_DBL_MAX;    
01126   
01127   const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;
01128 
01129   bool debug = false;
01130   
01131   //find the farthest boundary edge that is on the line segment    
01132   CubitVector split_position;
01133   for( int k=boundaryList.size(); k--; )
01134   {
01135     PST_Edge *tmp_edge = boundaryList.get_and_step();
01136 
01137     //if lines are parallel
01138     CubitVector tmp_vec1 = end_pt - start_pt;
01139     CubitVector tmp_vec2 = tmp_edge->end_coord() - tmp_edge->start_coord();
01140     tmp_vec1 *= tmp_vec2;
01141     if( tmp_vec1.length_squared() < TOL_SQR )
01142       continue;
01143 
01144     double t1 = tmp_edge->closest_on_line( start_pt, end_pt - start_pt );
01145     
01146     //ignore points that are coincident with start_pt or end_pt
01147     CubitVector position_on_edge = tmp_edge->position( t1 );
01148     double dist_sq_start = position_on_edge.distance_between_squared( start_pt );
01149     double dist_sq_end = position_on_edge.distance_between_squared( end_pt );
01150     if( dist_sq_start < TOL_SQR || dist_sq_end < TOL_SQR )      
01151       continue;
01152 
01153     if( t1 >= 0 && t1 <= 1 )    
01154     {       
01155       if( debug )
01156       {
01157         GfxDebug::draw_line( start_pt, end_pt, CUBIT_MAGENTA_INDEX );      
01158         PRINT_INFO("t1 = %f\n", t1 );
01159         tmp_edge->debug_draw( CUBIT_YELLOW_INDEX ); 
01160       }
01161   
01162       CubitVector tmp_pos = tmp_edge->position( t1 );
01163 
01164       //This checks if the intersection between the edge and segment
01165       //is not within the segment.
01166       if( tmp_pos.distance_between_squared( start_pt ) > seg_length_sq || 
01167           tmp_pos.distance_between_squared( end_pt ) > seg_length_sq )      
01168         continue;              
01169       
01170       if( debug )
01171       {
01172         GfxDebug::draw_point( tmp_pos, CUBIT_RED_INDEX );        
01173         GfxDebug::flush();
01174         GfxDebug::mouse_xforms(); 
01175       }
01176 
01177       double tmp_dist = tmp_pos.distance_between( start_pt );
01178       
01179       if( tmp_dist < shortest_length )
01180       {
01181         shortest_length = tmp_dist;        
01182         closest_pt = tmp_pos;
01183       }    
01184     }
01185   }
01186 
01187   if( shortest_length == CUBIT_DBL_MAX )
01188     return CUBIT_FAILURE;
01189 
01190   return CUBIT_SUCCESS;
01191 }
01192 
01193 
01194 //-------------------------------------------------------------------------
01195 // Purpose       : Given two existing points in a triangulation, add edges
01196 //                 to the triangulation to connect those two points by the
01197 //                 most direct path that does not leave the triangulation.
01198 //
01199 // Special Notes : Thisis the second step of the facet-polyline projection
01200 //                 algorithm.  First each point of the polyline is 
01201 //                 projected onto the facetting.  Second this function is
01202 //                 called to find/create the set of facet edges connecting
01203 //                 those points.
01204 //
01205 // Creator       : Jason Kraftcheck
01206 //
01207 // Creation Date : ??
01208 //-------------------------------------------------------------------------
01209 CubitStatus FacetProjectTool::project( PST_Point* &start, PST_Point* &end )
01210 {
01211   PST_Edge* edge,* closest_edge, *last_closest_edge=NULL;
01212   PST_Face* closest_face;
01213   PST_Point* point;
01214   PST_Point* start_point = start;  
01215   CubitStatus stat;
01216   
01217   const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;
01218 
01219   while( start_point != end )
01220   {
01221       // If the points share an edge, that edge is the result.
01222       // Return it.
01223     CubitBoolean is_boundary_edge;
01224     if( (edge = start_point->common( end )) )
01225     {
01226       PST_CoEdge* coedge = edge->end_point() == end ? 
01227                            edge->forward() : edge->reverse();
01228       imprintResult.append( coedge );
01229 
01230       if( DEBUG_FLAG( VG_FACET_DEBUG ) )
01231       {
01232         edge->debug_draw( VG_FACET_IMPRINT_COLOR );        
01233       }
01234     
01235       return CUBIT_SUCCESS;
01236     }    
01237 
01238       // Draw the current segment
01239     if( DEBUG_FLAG( VG_FACET_DEBUG ) )
01240     {
01241       GfxDebug::draw_line( 
01242         (float)start_point->x(), (float)start_point->y(), 
01243         (float)start_point->z(),
01244         (float)end->x()  , (float)end->y()  , (float)end->z(),
01245         VG_FACET_SEGMENT_COLOR );
01246       GfxDebug::flush();
01247       GfxDebug::mouse_xforms();
01248     }
01249     
01250       // Get face or edge adjacent to start and closest to the polyline segment
01251     stat = next_around_point( start_point, *end, closest_face, closest_edge, 
01252                               is_boundary_edge, last_closest_edge );
01253     
01254     if(closest_edge)    
01255       last_closest_edge = closest_edge;    
01256 
01257     if (!stat || (closest_face && closest_edge))
01258     {
01259  //     assert(false);
01260       return stat;
01261     }
01262 
01263     const CubitVector seg_tan = *end - *start_point;
01264     double t_seg = 0.0;
01265     point = 0;
01266     
01267     if (closest_face)
01268     {
01269         // calculate intersection with triangle edge opposite
01270         // 'start_point'.
01271       double t_edge;
01272       edge = closest_face->opposite(start_point);
01273       PST_Point* edge_start = edge->start_point();
01274       PST_Point* edge_end   = edge->end_point();
01275       closest_on_lines( *edge_start, *edge_end - *edge_start,
01276                         *start_point,  seg_tan, t_edge, t_seg);
01277       
01278       const CubitVector  seg_pt = *start_point + t_seg  *  seg_tan;
01279       const CubitVector edge_pt = (1.0 - t_edge) * *edge_start
01280                                        + t_edge  * *edge_end;
01281       
01282         // if line segment intersects opposite edge of triangle...
01283       if ( (t_seg < 1.0) || ((seg_pt - edge_pt).length_squared() < TOL_SQR) )
01284       {
01285           // Check if intersection is near start of opposite edge
01286         if ((t_edge <= 0.0) || 
01287             ((*edge_start - edge_pt).length_squared() < TOL_SQR) ||
01288             ((*edge_start - seg_pt ).length_squared() < TOL_SQR) )
01289           point = edge_start;
01290           // Check if intersection is near end of opposite edge
01291         else if ((t_edge >= 1.0) || 
01292                  ((*edge_end - edge_pt).length_squared() < TOL_SQR) ||
01293                  ((*edge_end -  seg_pt).length_squared() < TOL_SQR))
01294           point = edge_end;
01295           // Otherwise intersection is in the interior of the edge
01296         else
01297           point = split_edge(edge, t_edge);
01298       }
01299         // otherwise segment end is inside triangle
01300       else
01301       {
01302         t_seg = 1.0;
01303         const CubitVector proj = closest_face->plane().project(*end);
01304         PST_Edge *edge1, *edge2;
01305         closest_face->two_edges (start_point, edge1, edge2);
01306           // Too close to start position - skip segment
01307         if ( (*start_point - proj).length_squared() < TOL_SQR )
01308           point = start_point;
01309           // If close to side of triangle, fall through to edge
01310           // handling code
01311         else if (within_tolerance(edge1, *end, GEOMETRY_RESABS) ||
01312                  within_tolerance(edge1, proj, GEOMETRY_RESABS))
01313           closest_edge = edge1;
01314           // If close to side of triangle, fall through to edge
01315           // handling code
01316         else if (within_tolerance(edge2, *end, GEOMETRY_RESABS) ||
01317                  within_tolerance(edge2, proj, GEOMETRY_RESABS))
01318           closest_edge = edge2;
01319           // Insert new point in triangle interior...
01320         else
01321           point = insert_in_face (closest_face, proj);
01322       }
01323     } // if(closest_face)
01324     
01325     if (closest_edge)
01326     {
01327       PST_Point *const edge_end = closest_edge->other(start_point);
01328         // If edge and segment ends are equal...
01329       if ( (*edge_end - *end).length_squared() < TOL_SQR )
01330       {
01331         point = edge_end;
01332         edge = closest_edge;
01333         if (is_boundary_edge)
01334           end = edge_end;
01335       }
01336         // Otherwise calculate closest point
01337       else
01338       {
01339           // edge and segment share a start point so just need
01340           // to calculate the projection of each on the other.
01341         const CubitVector edge_tan = *edge_end - *start_point;
01342         const double dot_prod = seg_tan % edge_tan;
01343         const double t_seg  = dot_prod /  seg_tan.length_squared();
01344         const double t_edge = dot_prod / edge_tan.length_squared();
01345         
01346         const CubitVector  seg_pt = *start_point + t_seg  *  seg_tan;
01347         const CubitVector edge_pt = *start_point + t_edge * edge_tan;
01348         
01349           // Too close to start point -- skip segment
01350         if (t_edge <= 0.0 || t_seg <= 0.0 ||
01351             (*start_point - edge_pt).length_squared() < TOL_SQR ||
01352             (*start_point -  seg_pt).length_squared() < TOL_SQR )
01353           point = start_point;
01354           // If segment end is near or passed edge end, 
01355           // then the edge end is the intersection with this triangle
01356         else if (t_edge > 1.0 || 
01357                  (*edge_end - edge_pt).length_squared() < TOL_SQR ||
01358                  (*edge_end -  seg_pt).length_squared() < TOL_SQR)
01359         {
01360           //make sure the facet edge is not a boundary edge.
01361           if (is_boundary_edge && start_point == start)
01362             start =  edge_end;
01363           point = edge_end;
01364         }
01365           // Otherwise closest point to segment end is in the interior
01366           // of the edge -- split the edge.
01367         else
01368         {
01369           if(start_point != closest_edge->start_point())
01370             point = split_edge( closest_edge, 1.0 - t_edge );
01371           else
01372             point = split_edge( closest_edge, t_edge );
01373         }
01374       }
01375     } // if(closest_edge)
01376 
01377     if (point == start_point) // skip perpendicular or very short segment
01378       continue;
01379 
01380     edge = start_point->common( point );
01381     if (!edge)
01382     {
01383       assert(false);
01384       continue;
01385     }
01386     
01387     if ( !is_boundary_edge || start_point != start )
01388     { 
01389        PST_CoEdge* coedge = edge->end_point() == point ?
01390                             edge->forward() : edge->reverse();
01391        imprintResult.append(coedge);
01392     }
01393 
01394     if( DEBUG_FLAG( VG_FACET_DEBUG ) )
01395     {
01396       edge->debug_draw( VG_FACET_IMPRINT_COLOR );
01397       point->debug_draw( VG_FACET_IMPRINT_COLOR );
01398     }    
01399     start_point = point;
01400   }
01401   
01402   return CUBIT_SUCCESS;
01403 }
01404 
01405 //-------------------------------------------------------------------------
01406 // Purpose       : Check if a point is within the passed tolerance
01407 //                 of a bounded edge.
01408 //
01409 // Special Notes : 
01410 //
01411 // Creator       : Jason Kraftcheck
01412 //
01413 // Creation Date : 09/04/03
01414 //-------------------------------------------------------------------------
01415 bool FacetProjectTool::within_tolerance( PST_Edge* edge,
01416                                          const CubitVector& point,
01417                                          double tolerance )
01418 {
01419   const double t_edge = edge->closest_on_edge( point );
01420   const CubitVector edge_pos = edge->position(t_edge);
01421   return (edge_pos - point).length_squared() < tolerance*tolerance;
01422 }
01423 
01424 //-------------------------------------------------------------------------
01425 // Purpose       : Given a line segment beginning at a point in the
01426 //                 triangulation, find the face or edge adjacent to the
01427 //                 point and closest to the line segment.
01428 //
01429 // Special Notes : 
01430 //
01431 // Creator       : Jason Kraftcheck
01432 //
01433 // Creation Date : 09/04/03
01434 //-------------------------------------------------------------------------
01435 CubitStatus FacetProjectTool::next_around_point ( PST_Point*& start,
01436                                         const CubitVector& seg_end,
01437                                         PST_Face*& closest_face,
01438                                         PST_Edge*& closest_edge,
01439                                         CubitBoolean & is_boundary_edge,
01440                                         PST_Edge *last_closest_edge)
01441 {
01442   DLIList<PST_Face*> face_list;
01443   DLIList<PST_Edge*> boundary_edges;
01444   is_boundary_edge = CUBIT_FALSE; 
01445   double shortest_dist_sqr = CUBIT_DBL_MAX;
01446   closest_face = 0;
01447   closest_edge = 0;
01448   
01449   //const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;
01450   
01451   const CubitVector tangent = seg_end - *start;
01452    
01453     // To traverse the edges in clockwise order about the point,
01454     // it is necessary to start with the appropriate boundary edge
01455     // if there is one.  If this point is not on the boundary of 
01456     // the facet patch, any edge is acceptable.
01457 
01458     // First mark all faces adjacent to the vertex.  This is to make
01459     // sure disjoint sets of faces (sets of faces with no edges in
01460     // common) are not missed when traversing edge-face adjacencies.
01461   start->faces( &face_list );
01462   for (int i = face_list.size(); i--; )  
01463     face_list.step_and_get()->mark = true;
01464   
01465     // Loop until all faces in face_list have been un-marked
01466   boundary_edges.clean_out();
01467   closest_edge = 0;
01468   closest_face = 0; 
01469   while (face_list.size())
01470   {
01471     PST_Face* face = face_list.pop();
01472     if (!face->mark)
01473       continue;
01474 
01475     PST_CoEdge* coedge = face->first_coedge();
01476     if (coedge->edge()->other(start) == 0)
01477       coedge = coedge->next();
01478 
01479       // search for the 'first' coedge in the clockwise list.
01480     PST_CoEdge* first = coedge;
01481     while (coedge->face())
01482     {
01483       if (coedge->end_point() == start)
01484         coedge = coedge->other();
01485       else
01486         coedge = coedge->previous();
01487 
01488       if (coedge == first) // no start, point is not on boundary
01489         break;
01490     }
01491 
01492       // keep boundary edges encountered for later
01493     if (!coedge->face())
01494       boundary_edges.append(coedge->edge());
01495 
01496       // Of the two edges on the face that are adjacent to the
01497       // point, begin with the one that is first in clockwise
01498       // order around the point.
01499     if (coedge->start_point() == start)
01500       coedge = coedge->other();
01501     first = coedge;
01502 
01503       // Traverse adjacent facets in clockwise order
01504       // around point.
01505     bool clockwise_of_prev = false;
01506     while ( (face = coedge->face()) )
01507     {
01508       face->mark = false;
01509 
01510         // get vectors for edges of face pointing outwards
01511         // from vertex.
01512       assert(coedge->end_point() == start && coedge->next()->start_point() == start);
01513       CubitVector prev = coedge->other()->direction();
01514       CubitVector next = coedge->next()->direction();
01515 
01516         // calculate if segment is to the 'inside' of each
01517         // adjacent edge.
01518       CubitVector normal = next * prev;  // ccw 
01519       bool inside_prev = ((tangent * prev) % normal) >= 0.0;
01520       bool inside_next = ((next * tangent) % normal) >= 0.0;
01521 
01522         // If edge is inside face, check distance from face
01523       if( (inside_prev && inside_next) )
01524       {
01525         // calculate distance from end point to plane of face
01526         const double plane_coeff = normal % *start;
01527         const double end_coeff = normal % seg_end;
01528         const double numerator = end_coeff - plane_coeff;
01529         const double dist_sqr = (numerator * numerator) / normal.length_squared();         
01530         
01531         if (dist_sqr < shortest_dist_sqr)
01532         {
01533           closest_face = face;
01534           closest_edge = 0;
01535           shortest_dist_sqr = dist_sqr;
01536         }
01537 
01538         clockwise_of_prev = false;
01539       }
01540         // If the edge is above a peak formed by the shared edge of two faces
01541       else if (inside_next && clockwise_of_prev)
01542       {
01543           // calculate distance from end of segment to line of facet edge
01544         const double dot_prod = tangent % prev;
01545         const double dist_sqr = 
01546           tangent.length_squared() - dot_prod * dot_prod / prev.length_squared();
01547 
01548         if (dist_sqr <= shortest_dist_sqr && coedge->edge() != last_closest_edge)
01549         {
01550           closest_face = 0;
01551           closest_edge = coedge->edge();
01552           shortest_dist_sqr = dist_sqr;
01553         }
01554         clockwise_of_prev = false;
01555       }
01556         // If clockwise of the face, note for next iteration
01557       else if (inside_prev)
01558       {
01559         clockwise_of_prev = true;
01560       }
01561       else
01562       {
01563         clockwise_of_prev = false;
01564       }  
01565 
01566       coedge = coedge->next()->other();
01567       if (coedge == first)
01568         break;
01569     } // while(coedge->face())
01570     
01571     if (!coedge->face())
01572       boundary_edges.append(coedge->edge());
01573   
01574   } // while(face_list.size())
01575   
01576     // If a closest entity was found, then done.
01577   if (closest_face || closest_edge)
01578     return CUBIT_SUCCESS;
01579     
01580 
01581   //Here we try to intersect the segment with the boundary 
01582   //facet edge closest to 'start' (farthest from 'end')
01583   bool debug = false;
01584   if( !closest_edge )
01585   {
01586     CubitVector seg_direction = seg_end - *start;   
01587     double seg_length = seg_direction.length();
01588     double longest_length = 0;            
01589     
01590     if( debug )
01591     {
01592       GfxDebug::clear(); 
01593       start->debug_draw(CUBIT_GREEN_INDEX);
01594       GfxDebug::draw_point( seg_end, CUBIT_RED_INDEX );
01595       for( int j=facetList.size(); j--; )
01596         facetList.get_and_step()->debug_draw( CUBIT_CYAN_INDEX );
01597       GfxDebug::mouse_xforms(); 
01598     }
01599     
01600     //find the farthest boundary edge that is on the line segment    
01601     CubitVector split_position;
01602     for( int k=boundaryList.size(); k--; )
01603     {
01604       PST_Edge *tmp_edge = boundaryList.get_and_step();
01605 
01606       if( tmp_edge->start_point() == start ||
01607           tmp_edge->end_point() == start )
01608           continue;      
01609 
01610       double t1 = tmp_edge->closest_on_line( *start, seg_end - *start );
01611       
01612       if( t1 > 1e-6 && t1 < 1 )
01613       { 
01614         CubitVector tmp_pos = tmp_edge->position( t1 );
01615 
01616         if( debug )
01617         {
01618           GfxDebug::draw_line( *start, seg_end, CUBIT_MAGENTA_INDEX );
01619           PRINT_INFO("t1 = %f\n", t1 );
01620           tmp_edge->debug_draw( CUBIT_YELLOW_INDEX );
01621           GfxDebug::draw_point( tmp_pos, CUBIT_RED_INDEX );        
01622           GfxDebug::flush();
01623           GfxDebug::mouse_xforms(); 
01624         }
01625 
01626         double tmp_dist = tmp_pos.distance_between( seg_end );
01627 
01628         if( tmp_dist > longest_length && tmp_dist < seg_length )
01629         {
01630           longest_length = tmp_dist;
01631           closest_edge = tmp_edge;
01632           split_position = tmp_pos;
01633         }
01634       }
01635     }
01636     
01637     if( closest_edge )
01638     {
01639       start = project( split_position );    
01640 
01641       CubitStatus my_stat = next_around_point( start, seg_end, closest_face, 
01642         closest_edge, is_boundary_edge, last_closest_edge );
01643 
01644       if( CUBIT_FAILURE == my_stat )
01645         return CUBIT_FAILURE;
01646      
01647       return closest_face ? CUBIT_SUCCESS : CUBIT_FAILURE;      
01648     }   
01649   } 
01650   
01651   return closest_edge ? CUBIT_SUCCESS : CUBIT_FAILURE;
01652 }
01653 
01654   
01655   
01656     
01657   
01658 
01659 //-------------------------------------------------------------------------
01660 // Purpose       : Remove duplciate edges from results if input polyline
01661 //                 was self-intersecting.
01662 //
01663 // Special Notes : 
01664 //
01665 // Creator       : Jason Kraftcheck
01666 //
01667 // Creation Date : ??
01668 //-------------------------------------------------------------------------
01669 void FacetProjectTool::finalize_results()
01670 {
01671   imprintResult.reset();
01672   for ( int i = imprintResult.size(); i--; )
01673   {
01674     PST_Edge* edge = imprintResult.step_and_get()->edge();
01675     if (edge->mark)
01676       imprintResult.change_to(0);
01677     else
01678       edge->mark = 1;
01679   }
01680   imprintResult.remove_all_with_value(0);
01681   imprintResult.reset();
01682 }
01683 
01684 //-------------------------------------------------------------------------
01685 // Purpose       : Return non-intersecting chains of result edges, in the
01686 //                 same sequence as the input polyline and with the same
01687 //                 orientation.
01688 //
01689 // Special Notes : This is called after the projection is done to get the
01690 //                 resulting facet edge chains.
01691 //
01692 // Creator       : Jason Kraftcheck
01693 //
01694 // Creation Date : 09/04/03
01695 //-------------------------------------------------------------------------
01696 CubitStatus FacetProjectTool::get_result_set( DLIList<PST_CoEdge*>& coedges )
01697 {
01698   if ( !imprintResult.size() )
01699     return CUBIT_FAILURE;
01700   
01701   imprintResult.reset();
01702   PST_CoEdge* prev = imprintResult.get();
01703   imprintResult.change_to(0);
01704   coedges.append(prev);
01705   
01706   for ( int i = imprintResult.size() - 1; i--; )
01707   {
01708     PST_CoEdge* coedge = imprintResult.step_and_get();
01709     PST_Point* pt = prev->end_point();
01710     if (pt != coedge->start_point())
01711       break;
01712     
01713     int count = 1;
01714     for (PST_Edge* pt_edge = pt->next(coedge->edge());
01715          pt_edge != coedge->edge();
01716          pt_edge = pt->next(pt_edge))
01717     {
01718       if (pt_edge->mark)
01719         count++;
01720     }
01721      
01722     if (count != 2)
01723       break;
01724     
01725     coedges.append(coedge);
01726     prev = coedge;
01727     imprintResult.change_to(0);
01728   }
01729   
01730   imprintResult.remove_all_with_value(0);
01731   return CUBIT_SUCCESS;
01732 }
01733 
01734 
01735 double FacetProjectTool::closest_on_line( const CubitVector& B,
01736                                       const CubitVector& M,
01737                                       const CubitVector& P ) 
01738 { 
01739   const double dist_tol_sqr = GEOMETRY_RESABS*GEOMETRY_RESABS;
01740   double D = M.length_squared();
01741   return (D < dist_tol_sqr ) ? 0.0 : (M % (P - B)) / D;
01742 }
01743 
01744 
01745 void FacetProjectTool::closest_on_lines( const CubitVector& B1,
01746                                      const CubitVector& M1,
01747                                      const CubitVector& B2,
01748                                      const CubitVector& M2,
01749                                      double& t1, double& t2 ) 
01750 {
01751   CubitVector N = M2 * (M2 * M1);
01752   double d = N % M1;
01753   if( fabs(d) < CUBIT_RESABS ) //parallel lines
01754     t1 = 0.0;
01755   else
01756     t1 = ((N % B2) - (N % B1)) / d;
01757 
01758   t2 = closest_on_line( B2, M2, B1 + t1 * M1 );
01759 }
01760   
01761       
01762 bool FacetProjectTool::inside_face( PST_Face* face, const CubitVector& pos )
01763 {      
01764     // This code checks if the position is in the interior of 
01765     // the prism formed by sweeping the facet infinetly in the
01766     // direction of, and opposite direction of the facet normal.
01767   PST_CoEdge* coe = face->first_coedge();
01768   do {
01769     CubitVector vect( pos );
01770     vect -= *coe->start_point();
01771     vect *= coe->direction();
01772     if( vect % face->normal() > -CUBIT_RESABS )
01773       return false;
01774     coe = coe->next();
01775   } while( coe != face->first_coedge() );
01776   return true;
01777 }
01778 
01779 /*
01780 static void draw_pt( const std::vector<double>& list, int index, int color )
01781 {
01782   GfxDebug::draw_point( 
01783     (float)list[index], (float)list[index+1], (float)list[index+2], color );
01784 }
01785 
01786 static void draw_edge( const std::vector<double>& list, int i1, int i2, int color )
01787 {
01788   i1 *= 3; i2 *= 3;
01789   GfxDebug::draw_line(
01790     (float)list[i1], (float)list[i1+1], (float)list[i1+2],
01791     (float)list[i2], (float)list[i2+1], (float)list[i2+2],
01792     color );
01793 }
01794 
01795 static void draw_tri( const std::vector<double>& pts,
01796                       const std::vector<int>& tris, int index, 
01797                       const std::vector<int>& seq, int color )
01798 {
01799   for ( int i = 0; i < 3; i++ )
01800   {
01801     int i1 = tris[index + i];
01802     int i2 = tris[index + (i + 1) % 3];
01803     if ( !seq[i1] || !seq[i2] || abs(seq[i1]-seq[i2]) != 1 )
01804       draw_edge( pts, i1, i2, color );
01805   }
01806 }
01807 */
01808 
01809 void FacetProjectTool::debug( DLIList<CubitVector*>& segments,
01810                               std::vector<double>& coordinates,
01811                               std::vector<int>& connections )
01812 {
01813 #if 1 /* test with old interface */  
01814   DLIList<PST_Edge*> edges;
01815   PST_Edge::make_facets(coordinates, connections, GEOMETRY_RESABS, edges);
01816   
01817   DLIList<PST_Face*> faces;
01818   PST_Edge::faces( edges, faces );
01819   populate_data(faces);
01820   
01821   CubitBoolean point_changed;
01822   do_projection( segments, point_changed );
01823   
01824   PST_Edge::debug_draw_edges( edgeList, CUBIT_BLUE_INDEX );
01825   
01826   DLIList<PST_CoEdge*> coedges;
01827   finalize_results();
01828   while( get_result_set( coedges ) ) {
01829     edges.clean_out();
01830     coedges.reverse();
01831     while (coedges.size()) edges.append(coedges.pop()->edge());
01832     PST_Edge::debug_draw_edges( edges, CUBIT_WHITE_INDEX );
01833     PST_Edge::debug_draw_points( edges, CUBIT_WHITE_INDEX );
01834   }
01835 
01836 #else /* test with new interface */
01837 
01838   int i, j;
01839 
01840     // do projection
01841   std::vector<double> new_pts;
01842   std::vector<int> dudded, new_facets, facet_index, polyline;
01843   if( !project( segments, coordinates, connections, dudded, 
01844                 new_facets, facet_index, new_pts, polyline ) )
01845   {
01846     PRINT_ERROR("FacetProjectTool::project returned CUBIT_FAILURE\n");
01847     return;
01848   } 
01849   
01850   assert( connections.size() % 3 == 0 );
01851   assert( coordinates.size() % 3 == 0 );
01852   assert( new_facets.size() % 3 == 0 );
01853   assert( new_pts.size() % 3 == 0 );
01854   assert( dudded.size() == facet_index.size() );
01855   int num_old_tri = connections.size() / 3;
01856   int num_new_tri = new_facets.size() / 3;
01857   int num_old_pts = coordinates.size() / 3;
01858   int num_new_pts = new_pts.size() / 3;
01859   int num_tot_tri = num_new_tri + num_old_tri;
01860   int num_tot_pts = num_new_pts + num_old_pts;
01861   
01862   PRINT_INFO("%d input triangles using %d points.\n", num_old_tri, num_old_pts );
01863   PRINT_INFO("%d new triangles and %d new poitns.\n", num_new_tri, num_new_pts );
01864   PRINT_INFO("%d total triangles using %d total points.\n", num_tot_tri, num_tot_pts );
01865   PRINT_INFO("connections.size() = %d\n", (int)connections.size());
01866   PRINT_INFO("coordinates.size() = %d\n", (int)coordinates.size());
01867   PRINT_INFO("dudded.size() = %d\n", (int)dudded.size());
01868   PRINT_INFO("new_facets.size() = %d\n", (int)new_facets.size());
01869   PRINT_INFO("facet_index.size() = %d\n", (int)facet_index.size());
01870   PRINT_INFO("new_pts.size() = %d\n", (int)new_pts.size());
01871   PRINT_INFO("polyline.size() = %d\n", (int)polyline.size());
01872   
01873   for ( i = 0; i < (int)new_facets.size(); i++ )
01874     if ( new_facets[i] >= num_tot_pts || new_facets[i] < 0 )
01875       PRINT_ERROR("Invalid value %d at new_facets[%d]\n", new_facets[i], i);
01876   for ( i = 0; i < (int)dudded.size(); i++ )
01877     if ( dudded[i] >= num_old_tri || dudded[i] < 0 )
01878       PRINT_ERROR("Invalid value %d at dudded[%d]\n", dudded[i], i);
01879   for ( i = 0; i < (int)polyline.size()-1; )
01880   {
01881     int count = polyline[i];
01882     if ( count < 2 )
01883       PRINT_ERROR("Invalid polyline length %d at polyline[%d]\n", count, i );
01884     i++;
01885     for ( j = 0; j < count; j++, i++ )
01886     {
01887       if( i >= (int)polyline.size() )
01888         PRINT_ERROR("Ran out of indices in polyline.\n");
01889       if ( polyline[i] >= num_tot_pts || polyline[i] < 0 )
01890         PRINT_ERROR("Invalid value %d at dudded[%d]\n", dudded[i], i);
01891     }
01892   }
01893   
01894     // draw original points
01895   for ( i = 0; i < (int)coordinates.size(); i += 3 )
01896     draw_pt( coordinates, i, CUBIT_BLUE_INDEX );
01897 /*  
01898     // draw new points
01899   for ( i = 0; i < (int)new_pts.size(); i += 3 )
01900     draw_pt( new_pts, i, CUBIT_RED_INDEX ); 
01901 */  
01902     // concatenate point lists
01903   for ( i = 0; i < (int)new_pts.size(); i++ )
01904     coordinates.push_back(new_pts[i]);
01905   
01906     // make reverse mapping for dudded facets
01907   std::vector<bool> dead(connections.size()/3);
01908   for ( i = 0; i < (int)dead.size(); i++ )
01909     dead[i] = false;
01910   for ( i = 0; i < (int)dudded.size(); i++ )
01911     dead[dudded[i]] = true;
01912     
01913     // make polyline sequence list (this won't work quite
01914     // right for self-intersecting polylines -- *shrug*)
01915   std::vector<int> sequence(coordinates.size() / 3);
01916   for ( i = 0; i < (int)sequence.size(); i++ )
01917     sequence[i] = 0;
01918   j = 0;
01919   for ( i = 0; i < (int)polyline.size(); i++ )
01920     sequence[polyline[i]] = ++j;
01921   
01922   
01923     // draw orginal facets
01924   for ( i = 0; i < (int)dead.size(); i++ )
01925     if( ! dead[i] )
01926       draw_tri( coordinates, connections, 3*i, sequence, CUBIT_BLUE_INDEX );
01927   
01928     // draw new facets
01929   for ( i = 0; i < (int)new_facets.size(); i += 3 )
01930     draw_tri( coordinates, new_facets, i, sequence, CUBIT_RED_INDEX );
01931   
01932     // draw polyline
01933   for ( i = 0; i < (int)polyline.size()-1; i++)
01934   {
01935     int count = polyline[i];
01936     if ( count < 2 )
01937       PRINT_ERROR("Invalid polyline length %d at polyline[%d]\n", count, i );
01938     int first = ++i;
01939     for ( j = 0; j < count-1; j++, i++ )
01940     {
01941       draw_edge( coordinates, polyline[i], polyline[i+1], CUBIT_WHITE_INDEX );
01942     }
01943     if( polyline[first] == polyline[i] )
01944       PRINT_INFO("Polyline is closed.\n");
01945   }
01946 
01947   GfxDebug::flush();
01948 #endif
01949   cleanup();  
01950 }
01951 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines