cgma
FacetProjectTool Class Reference

#include <FacetProjectTool.hpp>

List of all members.

Public Member Functions

 FacetProjectTool ()
 ~FacetProjectTool ()
CubitStatus project (DLIList< CubitVector * > &segments, const std::vector< double > &coordinates, const std::vector< int > &connections, std::vector< int > &duddedFacets, std::vector< int > &newFacets, std::vector< int > &newFacetsIndex, std::vector< double > &newPoints, std::vector< int > &edges, std::vector< int > &segmentPoints, const double *tolerance_length=NULL)
void debug (DLIList< CubitVector * > &segments, std::vector< double > &coordinates, std::vector< int > &connections)

Protected Member Functions

void cleanup ()
PST_Pointproject (const CubitVector &point, const double *tolerance_length=NULL)
CubitStatus get_result_set (DLIList< PST_CoEdge * > &edges)
CubitStatus do_projection (const DLIList< CubitVector * > &segments, CubitBoolean &point_changed, const double *tolerance_length=NULL)
CubitStatus project (PST_Point *&start, PST_Point *&end)
PST_Faceclosest_face (const CubitVector &position)
PST_Edgeclosest_edge (const CubitVector &position)
PST_Pointclosest_point (const CubitVector &position)
PST_Pointinsert_in_face (PST_Face *face, const CubitVector &position)
PST_Pointsplit_edge (PST_Edge *edge, double t)
CubitStatus populate_data (const DLIList< PST_Face * > &facets)
CubitStatus next_around_point (PST_Point *&point, const CubitVector &segment_end, PST_Face *&closest_face, PST_Edge *&closest_edge, CubitBoolean &is_boundary_edge, PST_Edge *last_closest_edge)

Static Protected Member Functions

static bool project_to_face (PST_Face *face, const CubitVector &position, CubitVector &result)
static bool inside_face (PST_Face *face, const CubitVector &position)
static void closest_on_lines (const CubitVector &base1, const CubitVector &direction1, const CubitVector &base2, const CubitVector &direction2, double &t1, double &t2)
static double closest_on_line (const CubitVector &base, const CubitVector &direction, const CubitVector &from_position)
static bool within_tolerance (PST_Edge *edge, const CubitVector &point, double tolerance)

Private Member Functions

CubitStatus find_closest_int_point_with_boundary (CubitVector &start_pt, CubitVector &end_pt, CubitVector &closest_pt)
int point_on_facets (CubitVector *pos, double tolerance_sq)
void finalize_results ()

Private Attributes

DLIList< PST_Face * > facetList
DLIList< PST_Edge * > edgeList
DLIList< PST_Point * > pointList
DLIList< PST_Point * > segPoints
DLIList< PST_Edge * > boundaryList
DLIList< PST_CoEdge * > imprintResult

Detailed Description

Project a polyline onto a manifold patch of (triangular) facets.

See also:
FacetProjectTool::project

Definition at line 34 of file FacetProjectTool.hpp.


Constructor & Destructor Documentation

Definition at line 26 of file FacetProjectTool.cpp.

{ }

Definition at line 471 of file FacetProjectTool.cpp.

{
  cleanup();
}

Member Function Documentation

PST_Edge * FacetProjectTool::closest_edge ( const CubitVector position) [protected]

Definition at line 384 of file FacetProjectTool.cpp.

{
  PST_Edge* result = 0;
  double dist_sqr = 0;
  
  for( int i = edgeList.size(); i--; )
  {
    PST_Edge* edge = edgeList.step_and_get();
    
    double t = edge->closest_on_line( position );
    if( (t > -CUBIT_RESABS) && (t < (1.0 + CUBIT_RESABS)) )
    {
      double d = (edge->position(t) - position).length_squared();
      if( !result || d < dist_sqr )
      {
        result = edge;
        dist_sqr = d;
      }
    }
  }

  return result;
}
PST_Face * FacetProjectTool::closest_face ( const CubitVector position) [protected]

Definition at line 349 of file FacetProjectTool.cpp.

{
  PST_Face* result = 0;
  double dist_sqr = 0;
  CubitVector projected;
  
  for( int i = facetList.size(); i--; )
  {
    PST_Face* facet = facetList.step_and_get();
    if( project_to_face( facet, position, projected ) )
    {
      double d = (projected - position).length_squared();
      if( !result || d < dist_sqr )
      {
        result = facet;
        dist_sqr = d;
      }
    }
  }
  
  return result;
}
double FacetProjectTool::closest_on_line ( const CubitVector base,
const CubitVector direction,
const CubitVector from_position 
) [static, protected]

Definition at line 1735 of file FacetProjectTool.cpp.

{ 
  const double dist_tol_sqr = GEOMETRY_RESABS*GEOMETRY_RESABS;
  double D = M.length_squared();
  return (D < dist_tol_sqr ) ? 0.0 : (M % (P - B)) / D;
}
void FacetProjectTool::closest_on_lines ( const CubitVector base1,
const CubitVector direction1,
const CubitVector base2,
const CubitVector direction2,
double &  t1,
double &  t2 
) [static, protected]

Definition at line 1745 of file FacetProjectTool.cpp.

{
  CubitVector N = M2 * (M2 * M1);
  double d = N % M1;
  if( fabs(d) < CUBIT_RESABS ) //parallel lines
    t1 = 0.0;
  else
    t1 = ((N % B2) - (N % B1)) / d;

  t2 = closest_on_line( B2, M2, B1 + t1 * M1 );
}
PST_Point * FacetProjectTool::closest_point ( const CubitVector position) [protected]

Definition at line 417 of file FacetProjectTool.cpp.

{
  PST_Point* result = 0;
  double dist_sqr = 0;
  
  for( int i = pointList.size(); i--; )
  {
    PST_Point* pt = pointList.step_and_get();
    double d = (position - *pt).length_squared();
    if( !result || d < dist_sqr )
    {
      dist_sqr = d;
      result = pt;
    }
  }
  
  return result;
}
void FacetProjectTool::debug ( DLIList< CubitVector * > &  segments,
std::vector< double > &  coordinates,
std::vector< int > &  connections 
)

Definition at line 1809 of file FacetProjectTool.cpp.

{
#if 1 /* test with old interface */  
  DLIList<PST_Edge*> edges;
  PST_Edge::make_facets(coordinates, connections, GEOMETRY_RESABS, edges);
  
  DLIList<PST_Face*> faces;
  PST_Edge::faces( edges, faces );
  populate_data(faces);
  
  CubitBoolean point_changed;
  do_projection( segments, point_changed );
  
  PST_Edge::debug_draw_edges( edgeList, CUBIT_BLUE_INDEX );
  
  DLIList<PST_CoEdge*> coedges;
  finalize_results();
  while( get_result_set( coedges ) ) {
    edges.clean_out();
    coedges.reverse();
    while (coedges.size()) edges.append(coedges.pop()->edge());
    PST_Edge::debug_draw_edges( edges, CUBIT_WHITE_INDEX );
    PST_Edge::debug_draw_points( edges, CUBIT_WHITE_INDEX );
  }

#else /* test with new interface */

  int i, j;

    // do projection
  std::vector<double> new_pts;
  std::vector<int> dudded, new_facets, facet_index, polyline;
  if( !project( segments, coordinates, connections, dudded, 
                new_facets, facet_index, new_pts, polyline ) )
  {
    PRINT_ERROR("FacetProjectTool::project returned CUBIT_FAILURE\n");
    return;
  } 
  
  assert( connections.size() % 3 == 0 );
  assert( coordinates.size() % 3 == 0 );
  assert( new_facets.size() % 3 == 0 );
  assert( new_pts.size() % 3 == 0 );
  assert( dudded.size() == facet_index.size() );
  int num_old_tri = connections.size() / 3;
  int num_new_tri = new_facets.size() / 3;
  int num_old_pts = coordinates.size() / 3;
  int num_new_pts = new_pts.size() / 3;
  int num_tot_tri = num_new_tri + num_old_tri;
  int num_tot_pts = num_new_pts + num_old_pts;
  
  PRINT_INFO("%d input triangles using %d points.\n", num_old_tri, num_old_pts );
  PRINT_INFO("%d new triangles and %d new poitns.\n", num_new_tri, num_new_pts );
  PRINT_INFO("%d total triangles using %d total points.\n", num_tot_tri, num_tot_pts );
  PRINT_INFO("connections.size() = %d\n", (int)connections.size());
  PRINT_INFO("coordinates.size() = %d\n", (int)coordinates.size());
  PRINT_INFO("dudded.size() = %d\n", (int)dudded.size());
  PRINT_INFO("new_facets.size() = %d\n", (int)new_facets.size());
  PRINT_INFO("facet_index.size() = %d\n", (int)facet_index.size());
  PRINT_INFO("new_pts.size() = %d\n", (int)new_pts.size());
  PRINT_INFO("polyline.size() = %d\n", (int)polyline.size());
  
  for ( i = 0; i < (int)new_facets.size(); i++ )
    if ( new_facets[i] >= num_tot_pts || new_facets[i] < 0 )
      PRINT_ERROR("Invalid value %d at new_facets[%d]\n", new_facets[i], i);
  for ( i = 0; i < (int)dudded.size(); i++ )
    if ( dudded[i] >= num_old_tri || dudded[i] < 0 )
      PRINT_ERROR("Invalid value %d at dudded[%d]\n", dudded[i], i);
  for ( i = 0; i < (int)polyline.size()-1; )
  {
    int count = polyline[i];
    if ( count < 2 )
      PRINT_ERROR("Invalid polyline length %d at polyline[%d]\n", count, i );
    i++;
    for ( j = 0; j < count; j++, i++ )
    {
      if( i >= (int)polyline.size() )
        PRINT_ERROR("Ran out of indices in polyline.\n");
      if ( polyline[i] >= num_tot_pts || polyline[i] < 0 )
        PRINT_ERROR("Invalid value %d at dudded[%d]\n", dudded[i], i);
    }
  }
  
    // draw original points
  for ( i = 0; i < (int)coordinates.size(); i += 3 )
    draw_pt( coordinates, i, CUBIT_BLUE_INDEX );
/*  
    // draw new points
  for ( i = 0; i < (int)new_pts.size(); i += 3 )
    draw_pt( new_pts, i, CUBIT_RED_INDEX ); 
*/  
    // concatenate point lists
  for ( i = 0; i < (int)new_pts.size(); i++ )
    coordinates.push_back(new_pts[i]);
  
    // make reverse mapping for dudded facets
  std::vector<bool> dead(connections.size()/3);
  for ( i = 0; i < (int)dead.size(); i++ )
    dead[i] = false;
  for ( i = 0; i < (int)dudded.size(); i++ )
    dead[dudded[i]] = true;
    
    // make polyline sequence list (this won't work quite
    // right for self-intersecting polylines -- *shrug*)
  std::vector<int> sequence(coordinates.size() / 3);
  for ( i = 0; i < (int)sequence.size(); i++ )
    sequence[i] = 0;
  j = 0;
  for ( i = 0; i < (int)polyline.size(); i++ )
    sequence[polyline[i]] = ++j;
  
  
    // draw orginal facets
  for ( i = 0; i < (int)dead.size(); i++ )
    if( ! dead[i] )
      draw_tri( coordinates, connections, 3*i, sequence, CUBIT_BLUE_INDEX );
  
    // draw new facets
  for ( i = 0; i < (int)new_facets.size(); i += 3 )
    draw_tri( coordinates, new_facets, i, sequence, CUBIT_RED_INDEX );
  
    // draw polyline
  for ( i = 0; i < (int)polyline.size()-1; i++)
  {
    int count = polyline[i];
    if ( count < 2 )
      PRINT_ERROR("Invalid polyline length %d at polyline[%d]\n", count, i );
    int first = ++i;
    for ( j = 0; j < count-1; j++, i++ )
    {
      draw_edge( coordinates, polyline[i], polyline[i+1], CUBIT_WHITE_INDEX );
    }
    if( polyline[first] == polyline[i] )
      PRINT_INFO("Polyline is closed.\n");
  }

  GfxDebug::flush();
#endif
  cleanup();  
}
CubitStatus FacetProjectTool::do_projection ( const DLIList< CubitVector * > &  segments,
CubitBoolean point_changed,
const double *  tolerance_length = NULL 
) [protected]

Definition at line 784 of file FacetProjectTool.cpp.

{
  const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;
  point_changed = CUBIT_FALSE;
  DLIList<CubitVector*> segments( passed_list );  
  
  bool first_and_last_points_coincident = false;  
  if( segments[0]->distance_between_squared( *(segments[ segments.size()-1]) ) < TOL_SQR )
    first_and_last_points_coincident = true; 

  segments.reset();
  CubitVector end = *(segments.get_and_step());
  PST_Point* last_point = 0;

  bool did_any_projection = false;
  
  int where_is_last_pt;
  int where_is_end_pt;
  
  bool debug = false;  
  
  for( int i = segments.size(); i > 1; i-- )
  {
    CubitVector start(end);
    end = *(segments.get_and_step());
    
    if( debug )
    {
      GfxDebug::clear();       
      for( int j=facetList.size(); j--; )
        facetList.get_and_step()->debug_draw( CUBIT_CYAN_INDEX+j, false );        
      GfxDebug::draw_point( start, CUBIT_GREEN_INDEX );
      GfxDebug::draw_point( end, CUBIT_RED_INDEX );
      GfxDebug::flush();
      GfxDebug::mouse_xforms();      
    }

    if( DEBUG_FLAG( VG_FACET_DEBUG ) )
    {
      GfxDebug::draw_line( 
        (float)start.x(), (float)start.y(), (float)start.z(),
        (float)end.x()  , (float)end.y()  , (float)end.z(),
        VG_FACET_SEGMENT_COLOR );     
      GfxDebug::flush();
      GfxDebug::mouse_xforms();
    }

    double facet_dist_tol = start.distance_between_squared( end );
    facet_dist_tol *= 0.0001;

    //determine if start is on/off/on-boundary of facets
    if( !last_point )
    {
      where_is_last_pt = point_on_facets( &start, facet_dist_tol );

      if( 1 == where_is_last_pt || 
          2 == where_is_last_pt  )
      {
        last_point = project( start, tolerance_length );
        
        segPoints.append(last_point); 
        start = *last_point;

        if( first_and_last_points_coincident )
          *(segments[segments.size()-1]) = start;
      }
      else
        segPoints.append(NULL);  //it's not on the facets
    } 

    where_is_end_pt = point_on_facets( &end, facet_dist_tol );

    //both points off the facets
    //or one-on, one-off
    PST_Point *end_pt = NULL;

    //both points off
    //one-point-on && one-point-off
    if( (0 == where_is_last_pt && 0 == where_is_end_pt ) || 
        (0 == where_is_last_pt && 1 == where_is_end_pt ) ||
        (1 == where_is_last_pt && 0 == where_is_end_pt ) )
    {
      //see if segment intersects boundary
      CubitVector int_pt;
      CubitStatus my_status = find_closest_int_point_with_boundary( start, end, int_pt );
      if( CUBIT_FAILURE == my_status )      
      {
        segPoints.append( NULL );
        continue;
      }      

      if( 0 == where_is_last_pt )
      {
        last_point = project( int_pt, tolerance_length );

        if( 0 == where_is_end_pt )
        {
          my_status = find_closest_int_point_with_boundary( end, *last_point, int_pt );
          if( CUBIT_SUCCESS == my_status )
          {
            end_pt = project( int_pt, tolerance_length );
            i++;
            segments.back();
          }          
        } 
        else
        {
          end_pt = project( end, tolerance_length );
          segPoints.append( end_pt );
        }
      }
      else if( 0 == where_is_end_pt )
      {
        end_pt = project( int_pt, tolerance_length );
        i++;
        segments.back();
      }      
    }
    else if( (1 == where_is_last_pt || 2 == where_is_last_pt) && 
             (1 == where_is_end_pt || 2 == where_is_end_pt) )
    {
      end_pt = project( end, tolerance_length );
      segPoints.append( end_pt );
    }         
    else if( (0 == where_is_last_pt && 2 == where_is_end_pt ) || 
             (2 == where_is_last_pt && 0 == where_is_end_pt ) )
    {
      if( 0 == where_is_last_pt && 2 == where_is_end_pt )
      {
        //see if segment intersects boundary
        CubitVector int_pt;
        CubitStatus my_status = find_closest_int_point_with_boundary( start, end, int_pt );
        if( CUBIT_SUCCESS == my_status )
        {
          end_pt = project( end, tolerance_length );
          last_point = project( int_pt, tolerance_length );
          segPoints.append( end_pt );
        }
        else
        {
          segPoints.append(NULL);
          continue;
        }
      }
      else
      {
        //see if segment intersects boundary
        CubitVector int_pt;
        CubitStatus my_status = find_closest_int_point_with_boundary( end, start, int_pt );      

        if( CUBIT_SUCCESS == my_status )
        {
          end_pt = project( int_pt, tolerance_length );                        
          segPoints.append( NULL );
        }
        else
        {
          segPoints.append( NULL );
          continue;
        }
      }
    }

        
    PST_Point* start_point = last_point;
    PST_Point* end_point = end_pt;

    if( debug )
    {
      GfxDebug::draw_line( *last_point, *end_pt, CUBIT_MAGENTA_INDEX );
      GfxDebug::flush();
      GfxDebug::mouse_xforms();
    }

    if(project( last_point, end_pt ) == CUBIT_SUCCESS)
    {
      did_any_projection = true;
      if ((*start_point - *last_point).length_squared() > TOL_SQR)
      {
        segPoints.move_to(start_point);
        segPoints.change_to(last_point);
        point_changed = CUBIT_TRUE;
      }
      if ((*end_pt - *end_point).length_squared() > TOL_SQR)
      {
        segPoints.move_to(end_point);
        segPoints.change_to(end_pt);
        point_changed = CUBIT_TRUE;
      }    
    }

    last_point = end_pt;
    where_is_last_pt = where_is_end_pt;    
  
  } // end for( i = segments )
  
  if( !did_any_projection )
    return CUBIT_FAILURE;

  return CUBIT_SUCCESS;
}

Definition at line 1669 of file FacetProjectTool.cpp.

{
  imprintResult.reset();
  for ( int i = imprintResult.size(); i--; )
  {
    PST_Edge* edge = imprintResult.step_and_get()->edge();
    if (edge->mark)
      imprintResult.change_to(0);
    else
      edge->mark = 1;
  }
  imprintResult.remove_all_with_value(0);
  imprintResult.reset();
}

Definition at line 1120 of file FacetProjectTool.cpp.

{
  CubitVector seg_direction = end_pt - start_pt;   
  double seg_length_sq = seg_direction.length_squared();
  double shortest_length = CUBIT_DBL_MAX;    
  
  const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;

  bool debug = false;
  
  //find the farthest boundary edge that is on the line segment    
  CubitVector split_position;
  for( int k=boundaryList.size(); k--; )
  {
    PST_Edge *tmp_edge = boundaryList.get_and_step();

    //if lines are parallel
    CubitVector tmp_vec1 = end_pt - start_pt;
    CubitVector tmp_vec2 = tmp_edge->end_coord() - tmp_edge->start_coord();
    tmp_vec1 *= tmp_vec2;
    if( tmp_vec1.length_squared() < TOL_SQR )
      continue;

    double t1 = tmp_edge->closest_on_line( start_pt, end_pt - start_pt );
    
    //ignore points that are coincident with start_pt or end_pt
    CubitVector position_on_edge = tmp_edge->position( t1 );
    double dist_sq_start = position_on_edge.distance_between_squared( start_pt );
    double dist_sq_end = position_on_edge.distance_between_squared( end_pt );
    if( dist_sq_start < TOL_SQR || dist_sq_end < TOL_SQR )      
      continue;

    if( t1 >= 0 && t1 <= 1 )    
    {       
      if( debug )
      {
        GfxDebug::draw_line( start_pt, end_pt, CUBIT_MAGENTA_INDEX );      
        PRINT_INFO("t1 = %f\n", t1 );
        tmp_edge->debug_draw( CUBIT_YELLOW_INDEX ); 
      }
  
      CubitVector tmp_pos = tmp_edge->position( t1 );

      //This checks if the intersection between the edge and segment
      //is not within the segment.
      if( tmp_pos.distance_between_squared( start_pt ) > seg_length_sq || 
          tmp_pos.distance_between_squared( end_pt ) > seg_length_sq )      
        continue;              
      
      if( debug )
      {
        GfxDebug::draw_point( tmp_pos, CUBIT_RED_INDEX );        
        GfxDebug::flush();
        GfxDebug::mouse_xforms(); 
      }

      double tmp_dist = tmp_pos.distance_between( start_pt );
      
      if( tmp_dist < shortest_length )
      {
        shortest_length = tmp_dist;        
        closest_pt = tmp_pos;
      }    
    }
  }

  if( shortest_length == CUBIT_DBL_MAX )
    return CUBIT_FAILURE;

  return CUBIT_SUCCESS;
}

Definition at line 1696 of file FacetProjectTool.cpp.

{
  if ( !imprintResult.size() )
    return CUBIT_FAILURE;
  
  imprintResult.reset();
  PST_CoEdge* prev = imprintResult.get();
  imprintResult.change_to(0);
  coedges.append(prev);
  
  for ( int i = imprintResult.size() - 1; i--; )
  {
    PST_CoEdge* coedge = imprintResult.step_and_get();
    PST_Point* pt = prev->end_point();
    if (pt != coedge->start_point())
      break;
    
    int count = 1;
    for (PST_Edge* pt_edge = pt->next(coedge->edge());
         pt_edge != coedge->edge();
         pt_edge = pt->next(pt_edge))
    {
      if (pt_edge->mark)
        count++;
    }
     
    if (count != 2)
      break;
    
    coedges.append(coedge);
    prev = coedge;
    imprintResult.change_to(0);
  }
  
  imprintResult.remove_all_with_value(0);
  return CUBIT_SUCCESS;
}
PST_Point * FacetProjectTool::insert_in_face ( PST_Face face,
const CubitVector position 
) [protected]

Definition at line 508 of file FacetProjectTool.cpp.

{
  PST_Point* new_point = new PST_Point( position );
  PST_CoEdge* ce1 = face->first_coedge();
  PST_CoEdge* ce2 = ce1->next();
  PST_CoEdge* ce3 = ce2->next();

    // insert non-mainifold edge in face after ce1
  PST_Edge* edge1 = PST_Edge::insert_in_face( new_point, ce1 );
  if( ! edge1 )
  {
    delete new_point;
    return 0;
  }
  new_point->sequence = pointList.size();
  pointList.append(new_point);
  edgeList.append( edge1 );
  if( DEBUG_FLAG( VG_FACET_DEBUG ) )
    edge1->debug_draw( VG_FACET_FACET_COLOR );
  
    // connect non-manifold edge to opposite side of the face
    // to make manifold topology.
  PST_Edge* edge2 = PST_Edge::split_face( ce2, edge1->start_point() == new_point ?
                                               edge1->reverse() : edge1->forward() );   
  if( ! edge2 ) return new_point;
  edgeList.append( edge2 );
  PST_Face* new_face = edge2->other(face);
  new_face->sequence = facetList.size();
  new_face->parent = face->parent;
  facetList.append( new_face );
  if( face->owner() )
    face->owner()->notify_split( face, new_face );
    
  if( DEBUG_FLAG( VG_FACET_DEBUG ) )
    edge2->debug_draw( VG_FACET_FACET_COLOR );
  
    // Split face so that all faces are triangles.
  PST_Face* split_face = ce3->face();
  PST_Edge* edge3 = PST_Edge::split_face( new_point, ce3->end_point(), split_face );
  if( edge3 )
  {
    edgeList.append( edge3 );
    new_face = edge3->other(split_face);
    new_face->sequence = facetList.size();
    new_face->parent = face->parent;
    facetList.append( new_face );
    if( split_face->owner() )
      split_face->owner()->notify_split( split_face, new_face );
  }
  
  return new_point;
}
bool FacetProjectTool::inside_face ( PST_Face face,
const CubitVector position 
) [static, protected]

Definition at line 1762 of file FacetProjectTool.cpp.

{      
    // This code checks if the position is in the interior of 
    // the prism formed by sweeping the facet infinetly in the
    // direction of, and opposite direction of the facet normal.
  PST_CoEdge* coe = face->first_coedge();
  do {
    CubitVector vect( pos );
    vect -= *coe->start_point();
    vect *= coe->direction();
    if( vect % face->normal() > -CUBIT_RESABS )
      return false;
    coe = coe->next();
  } while( coe != face->first_coedge() );
  return true;
}
CubitStatus FacetProjectTool::next_around_point ( PST_Point *&  point,
const CubitVector segment_end,
PST_Face *&  closest_face,
PST_Edge *&  closest_edge,
CubitBoolean is_boundary_edge,
PST_Edge last_closest_edge 
) [protected]

Definition at line 1435 of file FacetProjectTool.cpp.

{
  DLIList<PST_Face*> face_list;
  DLIList<PST_Edge*> boundary_edges;
  is_boundary_edge = CUBIT_FALSE; 
  double shortest_dist_sqr = CUBIT_DBL_MAX;
  closest_face = 0;
  closest_edge = 0;
  
  //const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;
  
  const CubitVector tangent = seg_end - *start;
   
    // To traverse the edges in clockwise order about the point,
    // it is necessary to start with the appropriate boundary edge
    // if there is one.  If this point is not on the boundary of 
    // the facet patch, any edge is acceptable.

    // First mark all faces adjacent to the vertex.  This is to make
    // sure disjoint sets of faces (sets of faces with no edges in
    // common) are not missed when traversing edge-face adjacencies.
  start->faces( &face_list );
  for (int i = face_list.size(); i--; )  
    face_list.step_and_get()->mark = true;
  
    // Loop until all faces in face_list have been un-marked
  boundary_edges.clean_out();
  closest_edge = 0;
  closest_face = 0; 
  while (face_list.size())
  {
    PST_Face* face = face_list.pop();
    if (!face->mark)
      continue;

    PST_CoEdge* coedge = face->first_coedge();
    if (coedge->edge()->other(start) == 0)
      coedge = coedge->next();

      // search for the 'first' coedge in the clockwise list.
    PST_CoEdge* first = coedge;
    while (coedge->face())
    {
      if (coedge->end_point() == start)
        coedge = coedge->other();
      else
        coedge = coedge->previous();

      if (coedge == first) // no start, point is not on boundary
        break;
    }

      // keep boundary edges encountered for later
    if (!coedge->face())
      boundary_edges.append(coedge->edge());

      // Of the two edges on the face that are adjacent to the
      // point, begin with the one that is first in clockwise
      // order around the point.
    if (coedge->start_point() == start)
      coedge = coedge->other();
    first = coedge;

      // Traverse adjacent facets in clockwise order
      // around point.
    bool clockwise_of_prev = false;
    while ( (face = coedge->face()) )
    {
      face->mark = false;

        // get vectors for edges of face pointing outwards
        // from vertex.
      assert(coedge->end_point() == start && coedge->next()->start_point() == start);
      CubitVector prev = coedge->other()->direction();
      CubitVector next = coedge->next()->direction();

        // calculate if segment is to the 'inside' of each
        // adjacent edge.
      CubitVector normal = next * prev;  // ccw 
      bool inside_prev = ((tangent * prev) % normal) >= 0.0;
      bool inside_next = ((next * tangent) % normal) >= 0.0;

        // If edge is inside face, check distance from face
      if( (inside_prev && inside_next) )
      {
        // calculate distance from end point to plane of face
        const double plane_coeff = normal % *start;
        const double end_coeff = normal % seg_end;
        const double numerator = end_coeff - plane_coeff;
        const double dist_sqr = (numerator * numerator) / normal.length_squared();         
        
        if (dist_sqr < shortest_dist_sqr)
        {
          closest_face = face;
          closest_edge = 0;
          shortest_dist_sqr = dist_sqr;
        }

        clockwise_of_prev = false;
      }
        // If the edge is above a peak formed by the shared edge of two faces
      else if (inside_next && clockwise_of_prev)
      {
          // calculate distance from end of segment to line of facet edge
        const double dot_prod = tangent % prev;
        const double dist_sqr = 
          tangent.length_squared() - dot_prod * dot_prod / prev.length_squared();

        if (dist_sqr <= shortest_dist_sqr && coedge->edge() != last_closest_edge)
        {
          closest_face = 0;
          closest_edge = coedge->edge();
          shortest_dist_sqr = dist_sqr;
        }
        clockwise_of_prev = false;
      }
        // If clockwise of the face, note for next iteration
      else if (inside_prev)
      {
        clockwise_of_prev = true;
      }
      else
      {
        clockwise_of_prev = false;
      }  

      coedge = coedge->next()->other();
      if (coedge == first)
        break;
    } // while(coedge->face())
    
    if (!coedge->face())
      boundary_edges.append(coedge->edge());
  
  } // while(face_list.size())
  
    // If a closest entity was found, then done.
  if (closest_face || closest_edge)
    return CUBIT_SUCCESS;
    

  //Here we try to intersect the segment with the boundary 
  //facet edge closest to 'start' (farthest from 'end')
  bool debug = false;
  if( !closest_edge )
  {
    CubitVector seg_direction = seg_end - *start;   
    double seg_length = seg_direction.length();
    double longest_length = 0;            
    
    if( debug )
    {
      GfxDebug::clear(); 
      start->debug_draw(CUBIT_GREEN_INDEX);
      GfxDebug::draw_point( seg_end, CUBIT_RED_INDEX );
      for( int j=facetList.size(); j--; )
        facetList.get_and_step()->debug_draw( CUBIT_CYAN_INDEX );
      GfxDebug::mouse_xforms(); 
    }
    
    //find the farthest boundary edge that is on the line segment    
    CubitVector split_position;
    for( int k=boundaryList.size(); k--; )
    {
      PST_Edge *tmp_edge = boundaryList.get_and_step();

      if( tmp_edge->start_point() == start ||
          tmp_edge->end_point() == start )
          continue;      

      double t1 = tmp_edge->closest_on_line( *start, seg_end - *start );
      
      if( t1 > 1e-6 && t1 < 1 )
      { 
        CubitVector tmp_pos = tmp_edge->position( t1 );

        if( debug )
        {
          GfxDebug::draw_line( *start, seg_end, CUBIT_MAGENTA_INDEX );
          PRINT_INFO("t1 = %f\n", t1 );
          tmp_edge->debug_draw( CUBIT_YELLOW_INDEX );
          GfxDebug::draw_point( tmp_pos, CUBIT_RED_INDEX );        
          GfxDebug::flush();
          GfxDebug::mouse_xforms(); 
        }

        double tmp_dist = tmp_pos.distance_between( seg_end );

        if( tmp_dist > longest_length && tmp_dist < seg_length )
        {
          longest_length = tmp_dist;
          closest_edge = tmp_edge;
          split_position = tmp_pos;
        }
      }
    }
    
    if( closest_edge )
    {
      start = project( split_position );    

      CubitStatus my_stat = next_around_point( start, seg_end, closest_face, 
        closest_edge, is_boundary_edge, last_closest_edge );

      if( CUBIT_FAILURE == my_stat )
        return CUBIT_FAILURE;
     
      return closest_face ? CUBIT_SUCCESS : CUBIT_FAILURE;      
    }   
  } 
  
  return closest_edge ? CUBIT_SUCCESS : CUBIT_FAILURE;
}
int FacetProjectTool::point_on_facets ( CubitVector pos,
double  tolerance_sq 
) [private]

Definition at line 1001 of file FacetProjectTool.cpp.

{
  //Return Values:
  //0 is outside of facets
  //1 is inside of facets
  //2 is on boundary of facets

  PST_Point* point = closest_point( *pos );
  PST_Edge *  edge = closest_edge ( *pos );
  PST_Face *  face = closest_face ( *pos );

  // find closest facet
  double face_dist = CUBIT_DBL_MAX;
  CubitVector face_pos;
  if( face )
  {
    project_to_face( face, *pos, face_pos );
    face_dist = (*pos - face_pos).length_squared();
  }

  // find closest edge
  double edge_dist = CUBIT_DBL_MAX;
  double edge_pos  = CUBIT_DBL_MAX;
  if( edge )
  {
    edge_pos = edge->closest_on_edge( *pos );
    edge_dist = (*pos - edge->position(edge_pos)).length_squared();
  }

  // find closest point
  double point_dist = CUBIT_DBL_MAX;
  if( point )
  {
    point_dist = (*pos - *point).length_squared();
  }  

  // if the position is closer to a facet than a point or edge...
  if( face && (face_dist < point_dist) && (face_dist < edge_dist) )
  {
    // check if projected position is within tolerance of
    // any bounding edge of the face.
    PST_CoEdge* coe = face->first_coedge();
    bool edge_closer = false;
    do {
      PST_Edge* e = coe->edge();
      CubitVector p = e->position( e->closest_on_line( face_pos ) );
      double d = (p - face_pos).length_squared();

      if( d < tolerance_sq ) 
      {
        edge = e;
        edge_pos = edge->closest_on_edge( *pos );
        face_dist = d;
        edge_dist = d;
        edge_closer = true;
      }
      coe = coe->next();
    } while( coe != face->first_coedge() );
    
    if( edge_closer )
    {      
      if( boundaryList.is_in_list( edge ) )
        return 2;
      else
        return 1;    
    }      
    else
      return 1;    
    
    return 0;
  }

  // If the point was closer to an edge than a point
  if( edge && (edge_dist < point_dist) )
  {
    if( edge_dist < tolerance_sq )
    {
      if( boundaryList.is_in_list( edge ) )
        return 2;
      else
        return 1;
    }
    else
      return 0;
  }

  PST_Edge *first_edge = point->edge();
  PST_Edge *tmp_edge = first_edge;
  if( point_dist < tolerance_sq )
  {
    do
    {
      if( boundaryList.is_in_list( tmp_edge ) )
        return 2;
      tmp_edge = point->next( tmp_edge );
    }
    while( first_edge != tmp_edge );

    return 1;
  }
  return 0;
}
CubitStatus FacetProjectTool::populate_data ( const DLIList< PST_Face * > &  facets) [protected]

Definition at line 238 of file FacetProjectTool.cpp.

{
  facetList = facets;
  int i;
  
    // Mark all connected edges with a 3, and initialize 
    // parent number to sequence number.
  facetList.last();
  for( i = facetList.size(); i--; )
  {
    PST_Face* face = facetList.step_and_get();
    face->parent = face->sequence;
    PST_CoEdge* coe = face->first_coedge();
    do
    {
      coe->edge()->mark = 3;
      coe = coe->next();
    } while( coe != face->first_coedge() );
  }
  
    // Decrement the edge mark by one for each
    // face that contains the edge.
  for( i = facetList.size(); i--; )
  {
    PST_Face* face = facetList.step_and_get();
    PST_CoEdge* coe = face->first_coedge();
    do
    {
      coe->edge()->mark--;
      coe = coe->next();
    } while( coe != face->first_coedge() );
  }
  
    // Populate edgeList and boundaryList with edges
  for( i = facetList.size(); i--; )
  {
    PST_Face* face = facetList.step_and_get();
    PST_CoEdge* coe = face->first_coedge();
    do
    {
      PST_Edge* edge = coe->edge();
        // If mark is 2, edge was only in one face in the
        // passed list, and is thus on the boundary of the
        // passed patch of faces.
      if( edge->mark == 2 )
      {
        boundaryList.append( edge );
        edgeList.append( edge );
        edge->mark = 0;
      }
        // If mark is 1, the edge is an interior edge in the
        // patch of faces
      else if( edge->mark == 1 )
      {
        edgeList.append( edge );
        edge->mark = 0;
      }
        // If the mark on the edge is zero, we have already
        // added it to the list(s).
      else
      {
        assert( edge->mark == 0 );  
      }
      
      coe = coe->next();
    } while( coe != face->first_coedge() );
  }
  
    // Get all the points
  PST_Edge::points( edgeList, pointList );
  
    // Clear marks on all faces connected to our edge
    // list (the passed faces and any faces connected
    // at the boundary edges.)
  for( i = edgeList.size(); i--; )
  {
    PST_Edge* edge = edgeList.step_and_get();
    edge->mark = 0;
    if( edge->forward()->face() )
      edge->forward()->face()->mark = 0;
    if( edge->reverse()->face() )
      edge->reverse()->face()->mark = 0;
  }

    // clear marks on faces
  for( i = facetList.size(); i--; )
    facetList.step_and_get()->mark = 0;
    // Set mark to 1 on boundary edges
  for( i = boundaryList.size(); i--; )
    boundaryList.step_and_get()->mark = 1;
    
  if( DEBUG_FLAG( VG_FACET_DEBUG ) )
  {
//    PST_Face::draw_faces( facetList, VG_FACET_FACET_COLOR, false );
    PST_Edge::debug_draw_edges( boundaryList, VG_FACET_BOUNDARY_COLOR );
  }

  return CUBIT_SUCCESS;
} 
CubitStatus FacetProjectTool::project ( DLIList< CubitVector * > &  segments,
const std::vector< double > &  coordinates,
const std::vector< int > &  connections,
std::vector< int > &  duddedFacets,
std::vector< int > &  newFacets,
std::vector< int > &  newFacetsIndex,
std::vector< double > &  newPoints,
std::vector< int > &  edges,
std::vector< int > &  segmentPoints,
const double *  tolerance_length = NULL 
)

Project a polyline onto a manifold patch of (triangular) facets.

Parameters:
segmentsInput list of ordered polyline vertices.
coordinatesInput list of facet vertex coordinate triples ordered as { {x,y,z}, {x,y,z}, ...}
connectionsInput list of facets specified as { {p1,p2,p3}, ... } where each value p in the list corresponds to the vertex at coordinates[3*p]
duddedFacetsOutput list of facets that were modified during the operation where each value i in the list corresponds to the facet at connections[3*i]
newFacetsOutput list of new facets replacing those specified in duddedFacets where each value p corresponds to the vertex at coordinates[3*p] or the vertex at newPoints[3*p - coordinates.size()].
newFacetsIndexThis list is populated with indices specifying which facets in newFacets replace each facet in duddedFacets. The facets in newFacets are grouped by the dudded facet they replace. For each dudded facet specified at duddedFacets[i], the new facets replacing that facet begin at newFacets[a] and end at newFacets[b-1] where a = newFacetsIndex[i-1] and b = newFacetsIndex[i]. Thus this list will always be the same size as duddedFacets and the last value will always be newFacets.size().
newPointsOutput list of new facet vertices created during the operation ordered as { {x,y,z}, {x,y,z}, ...}
edgesOutput list of facet vertices specifying the projection of the passed polyline onto the facets. The list may contain multiple polylines if the projection was self-intersecting or intersected the boundary of the facet patch in more than two places. The format is: {n1, {p1.1, p1.2, .., p1.n1}, n2, {p2.1, p2.2, .., p2.n2), ... 0 } The list contains a count 'n' followed by n vertex indices. The last value is always zero. Each vertex index p corresponds the vertex at coordinates[3*p] or the vertex at newPoints[3*p - coordinates.size()].
segmentPointsThe list of facet points corresponding to input segment points. Point indices are the specified in the same fashion as those in 'edges'. A value of -1 indicates that the segment point did not project into the interior of the bounded surface.
Returns:
A return value of CUBIT_FAILURE indicates an internal error in FacetProjectTool. If the passed polyline does not project onto any part of the facet patch, the return value will be CUBIT_SUCCESS and the output lists will be empty.
See also:
FSInt

Definition at line 51 of file FacetProjectTool.cpp.

{
  int i;
//   const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS; 
                     
  assert(connections.size()%3 == 0);  //  For triangles these must be triples
  assert(coordinates.size()%3 == 0);  //  Coordinates must come in triples   
  DLIList<PST_Edge*> facet_edges;  

    // make the PST edges and faces from the coordinates and connections.
  PST_Edge::make_facets( coordinates, connections, GEOMETRY_RESABS, facet_edges );
  DLIList<PST_Face*> faces;

  CubitStatus success = CUBIT_SUCCESS;
  if(facet_edges.size() > 0)
  {

    PST_Edge::faces( facet_edges, faces );
    populate_data(faces);

      // Order the orignal points by sequence number.  New points
      // will be appended in order.
    pointList.sort(&sequence_compare_pts);
        
      // Do the work
    CubitBoolean point_changed;
    success = do_projection(segments, point_changed, tolerance_length);
    
    bool debug = false;
    if( debug )
    {
      GfxDebug::clear();
      for( int j=0; j<facetList.size(); j++ )
        facetList.get_and_step()->debug_draw( CUBIT_CYAN_INDEX );
      GfxDebug::mouse_xforms();
    }

      // fill in segmentPoints
//    assert (segPoints.size() == segments.size());
    segmentPoints.resize( segPoints.size() );
    segPoints.reset();
    for ( i = 0; i < segPoints.size(); i++ )
    {
      PST_Point* pt = segPoints.get_and_step();     
      segmentPoints[i] = pt ? pt->sequence : -1;
      if (point_changed)
        success = CUBIT_FAILURE;
    }
  }
  else
    success = CUBIT_FAILURE;
  
  if (!success)
  {
     cleanup();
     return success;
  }
                                                                                
    //  Now put the points with sequence > coordinates.size()/3 in newPoints.
  int orig_num_points = coordinates.size()/3;
  int num_new_points = pointList.size() - orig_num_points;
  pointList.reset();
  pointList.step(orig_num_points);
  newPoints.resize(num_new_points*3);
  std::vector<double>::iterator ditor = newPoints.begin();
  while( num_new_points-- ) {
    PST_Point* pt = pointList.get_and_step();
    *ditor++ = pt->x();
    *ditor++ = pt->y();
    *ditor++ = pt->z();
  }

  if ( facetList.size() > 1 ) {  //  If only one facet, skip finding dudded facets
    // Sort (group) facets by parent.
    facetList.sort(&parent_compare_facets);

    // Fill in the duddedFacets, newFacets, and newFacetsIndex vectors.  

#ifndef NDEBUG
    int orig_num_facets = connections.size()/3;
#endif
    int new_facet_index = 0;
    DLIList<PST_Point*> facet_pts(3);

    facetList.reset();
    duddedFacets.clear();
    newFacetsIndex.clear();
    newFacets.clear();

    for ( i = facetList.size(); i > 0; ) {
      PST_Face* dudded_facet = 0;
      if( facetList.get()->parent != facetList.next()->parent ) {

          // unmodified original facet - skip it
#ifndef NDEBUG
        PST_Face* facet = facetList.get();
        assert(facet->sequence < orig_num_facets &&
               facet->sequence == facet->parent);
#endif
        facetList.step();
        i--;
      } 
      else {
          // new facets
        
          // put all new facets split from the
          // same original facet (having same parent),
          // including the orignal facet which was 
          // recycled, into newFacets
        PST_Face* facet = 0;
        do {
          facet = facetList.get_and_step();
          i--;
        
          if ( facet->parent == facet->sequence ) {
            assert(!dudded_facet);
            dudded_facet = facet;
          }
      
          facet_pts.clean_out();
          facet->append_points(facet_pts);
          facet_pts.reset();
          newFacets.push_back(facet_pts.get_and_step()->sequence);
          newFacets.push_back(facet_pts.get_and_step()->sequence);
          newFacets.push_back(facet_pts.get_and_step()->sequence);
          new_facet_index += 3;
        } while ( (facet->parent == facetList.get()->parent) && (i > 0) );
      
          // add replaced facet to duddedFacets
        assert(dudded_facet && dudded_facet->sequence < orig_num_facets);
        duddedFacets.push_back(dudded_facet->sequence);
          // add end position in new_facets for the
          // set of replacement facets to newFacetsIndex
        newFacetsIndex.push_back(new_facet_index);
      }
    }  
  }
  
  DLIList<PST_CoEdge*> coedge_list;
  edges.clear();
  finalize_results();
  while( get_result_set(coedge_list) ) {

      // add count and first point
    coedge_list.reset();
    edges.push_back(coedge_list.size() + 1);
    PST_Point* pt = coedge_list.get()->start_point();
    edges.push_back(pt->sequence);
    
      // add all other points
    for( i = coedge_list.size(); i--; ) {
      pt = coedge_list.get_and_step()->end_point();
      edges.push_back(pt->sequence);
    }
  
      // clean up for next iteration
    coedge_list.clean_out();
  }
  edges.push_back(0);  //  terminating zero for next segment length      
  
  cleanup();
  return CUBIT_SUCCESS;
}
PST_Point * FacetProjectTool::project ( const CubitVector point,
const double *  tolerance_length = NULL 
) [protected]

Definition at line 574 of file FacetProjectTool.cpp.

{
  PST_Point* point = closest_point( pos );
  PST_Edge *  edge = closest_edge ( pos );
  PST_Face *  face = closest_face ( pos );
  
    // find closest facet
  double face_dist = CUBIT_DBL_MAX;
  CubitVector face_pos;
  if( face )
  {
    project_to_face( face, pos, face_pos );
    face_dist = (pos - face_pos).length_squared();
  }
  
    // find closest edge
  double edge_dist = CUBIT_DBL_MAX;
  double edge_pos  = CUBIT_DBL_MAX;
  if( edge )
  {
    edge_pos = edge->closest_on_edge( pos );
    edge_dist = (pos - edge->position(edge_pos)).length_squared();
  }
  
    // find closest point
  double point_dist = CUBIT_DBL_MAX;
  if( point )
  {
    point_dist = (pos - *point).length_squared();
  }
  
    // if the position is closer to a facet than a point or edge...
  if( face && (face_dist < point_dist) && (face_dist < edge_dist) )
  {
      // get a tolerance based on a sample face size. I attempted
      // to do this with a representative facet size but still 
      // ended up with problems.  While this is slower is more
      // robust.  KGM
    double facet_length;

    if(tolerance_length)
    {
      facet_length = *tolerance_length;
    }
    else
    {
      facet_length = face->bounding_length();
    }

    // use the square of the facet length times an arbitrary factor to
    // determine the tolerance
    double tol_sqr = facet_length*facet_length*.00001;

      // check if projected position is within tolerance of
      // any bounding edge of the face.
    PST_CoEdge* coe = face->first_coedge();
    bool insert = true;
    do {
      PST_Edge* e = coe->edge();
      CubitVector p = e->position( e->closest_on_line( face_pos ) );
      double d = (p - face_pos).length_squared();

      // use a relative tolerance here - KGM
      if( d < tol_sqr ) 
      {
        edge = e;
        edge_pos = edge->closest_on_edge( pos );
        face_dist = d;
        edge_dist = d;
        insert = false;
      }
      coe = coe->next();
    } while( coe != face->first_coedge() );
    
      // If projected position was sufficiently far from the edges
      // of the facet, insert an interior point in the facet.  Otherwise
      // fall through and split edge.
    if( insert )
    {
      return insert_in_face( face, face_pos );
    }
  }
  
    // If the point was closer to an edge than a point
  if( edge && (edge_dist <= point_dist) )
  {
      // If within tolerance of an end point of the edge, return
      // the end point.
    CubitVector pt = edge->position( edge_pos );
    if( (pt - *edge->start_point()).length_squared() < GEOMETRY_RESABS*GEOMETRY_RESABS )
      point = edge->start_point();
    else if( (pt - *edge->end_point()).length_squared() < GEOMETRY_RESABS*GEOMETRY_RESABS )
      point = edge->end_point();
      // Otherwise split the edge at the projected position
    else point = split_edge( edge, edge_pos );
  }
  
  if( DEBUG_FLAG( VG_FACET_DEBUG ) && point )
    point->debug_draw( VG_FACET_IMPRINT_COLOR );
  
    // If this was not set in the above if-block, then it is either
    // the closest point in the facet patch or NULL if the position
    // projected outside the facet patch.
  return point;
}    
CubitStatus FacetProjectTool::project ( PST_Point *&  start,
PST_Point *&  end 
) [protected]

Definition at line 1209 of file FacetProjectTool.cpp.

{
  PST_Edge* edge,* closest_edge, *last_closest_edge=NULL;
  PST_Face* closest_face;
  PST_Point* point;
  PST_Point* start_point = start;  
  CubitStatus stat;
  
  const double TOL_SQR = GEOMETRY_RESABS*GEOMETRY_RESABS;

  while( start_point != end )
  {
      // If the points share an edge, that edge is the result.
      // Return it.
    CubitBoolean is_boundary_edge;
    if( (edge = start_point->common( end )) )
    {
      PST_CoEdge* coedge = edge->end_point() == end ? 
                           edge->forward() : edge->reverse();
      imprintResult.append( coedge );

      if( DEBUG_FLAG( VG_FACET_DEBUG ) )
      {
        edge->debug_draw( VG_FACET_IMPRINT_COLOR );        
      }
    
      return CUBIT_SUCCESS;
    }    

      // Draw the current segment
    if( DEBUG_FLAG( VG_FACET_DEBUG ) )
    {
      GfxDebug::draw_line( 
        (float)start_point->x(), (float)start_point->y(), 
        (float)start_point->z(),
        (float)end->x()  , (float)end->y()  , (float)end->z(),
        VG_FACET_SEGMENT_COLOR );
      GfxDebug::flush();
      GfxDebug::mouse_xforms();
    }
    
      // Get face or edge adjacent to start and closest to the polyline segment
    stat = next_around_point( start_point, *end, closest_face, closest_edge, 
                              is_boundary_edge, last_closest_edge );
    
    if(closest_edge)    
      last_closest_edge = closest_edge;    

    if (!stat || (closest_face && closest_edge))
    {
 //     assert(false);
      return stat;
    }

    const CubitVector seg_tan = *end - *start_point;
    double t_seg = 0.0;
    point = 0;
    
    if (closest_face)
    {
        // calculate intersection with triangle edge opposite
        // 'start_point'.
      double t_edge;
      edge = closest_face->opposite(start_point);
      PST_Point* edge_start = edge->start_point();
      PST_Point* edge_end   = edge->end_point();
      closest_on_lines( *edge_start, *edge_end - *edge_start,
                        *start_point,  seg_tan, t_edge, t_seg);
      
      const CubitVector  seg_pt = *start_point + t_seg  *  seg_tan;
      const CubitVector edge_pt = (1.0 - t_edge) * *edge_start
                                       + t_edge  * *edge_end;
      
        // if line segment intersects opposite edge of triangle...
      if ( (t_seg < 1.0) || ((seg_pt - edge_pt).length_squared() < TOL_SQR) )
      {
          // Check if intersection is near start of opposite edge
        if ((t_edge <= 0.0) || 
            ((*edge_start - edge_pt).length_squared() < TOL_SQR) ||
            ((*edge_start - seg_pt ).length_squared() < TOL_SQR) )
          point = edge_start;
          // Check if intersection is near end of opposite edge
        else if ((t_edge >= 1.0) || 
                 ((*edge_end - edge_pt).length_squared() < TOL_SQR) ||
                 ((*edge_end -  seg_pt).length_squared() < TOL_SQR))
          point = edge_end;
          // Otherwise intersection is in the interior of the edge
        else
          point = split_edge(edge, t_edge);
      }
        // otherwise segment end is inside triangle
      else
      {
        t_seg = 1.0;
        const CubitVector proj = closest_face->plane().project(*end);
        PST_Edge *edge1, *edge2;
        closest_face->two_edges (start_point, edge1, edge2);
          // Too close to start position - skip segment
        if ( (*start_point - proj).length_squared() < TOL_SQR )
          point = start_point;
          // If close to side of triangle, fall through to edge
          // handling code
        else if (within_tolerance(edge1, *end, GEOMETRY_RESABS) ||
                 within_tolerance(edge1, proj, GEOMETRY_RESABS))
          closest_edge = edge1;
          // If close to side of triangle, fall through to edge
          // handling code
        else if (within_tolerance(edge2, *end, GEOMETRY_RESABS) ||
                 within_tolerance(edge2, proj, GEOMETRY_RESABS))
          closest_edge = edge2;
          // Insert new point in triangle interior...
        else
          point = insert_in_face (closest_face, proj);
      }
    } // if(closest_face)
    
    if (closest_edge)
    {
      PST_Point *const edge_end = closest_edge->other(start_point);
        // If edge and segment ends are equal...
      if ( (*edge_end - *end).length_squared() < TOL_SQR )
      {
        point = edge_end;
        edge = closest_edge;
        if (is_boundary_edge)
          end = edge_end;
      }
        // Otherwise calculate closest point
      else
      {
          // edge and segment share a start point so just need
          // to calculate the projection of each on the other.
        const CubitVector edge_tan = *edge_end - *start_point;
        const double dot_prod = seg_tan % edge_tan;
        const double t_seg  = dot_prod /  seg_tan.length_squared();
        const double t_edge = dot_prod / edge_tan.length_squared();
        
        const CubitVector  seg_pt = *start_point + t_seg  *  seg_tan;
        const CubitVector edge_pt = *start_point + t_edge * edge_tan;
        
          // Too close to start point -- skip segment
        if (t_edge <= 0.0 || t_seg <= 0.0 ||
            (*start_point - edge_pt).length_squared() < TOL_SQR ||
            (*start_point -  seg_pt).length_squared() < TOL_SQR )
          point = start_point;
          // If segment end is near or passed edge end, 
          // then the edge end is the intersection with this triangle
        else if (t_edge > 1.0 || 
                 (*edge_end - edge_pt).length_squared() < TOL_SQR ||
                 (*edge_end -  seg_pt).length_squared() < TOL_SQR)
        {
          //make sure the facet edge is not a boundary edge.
          if (is_boundary_edge && start_point == start)
            start =  edge_end;
          point = edge_end;
        }
          // Otherwise closest point to segment end is in the interior
          // of the edge -- split the edge.
        else
        {
          if(start_point != closest_edge->start_point())
            point = split_edge( closest_edge, 1.0 - t_edge );
          else
            point = split_edge( closest_edge, t_edge );
        }
      }
    } // if(closest_edge)

    if (point == start_point) // skip perpendicular or very short segment
      continue;

    edge = start_point->common( point );
    if (!edge)
    {
      assert(false);
      continue;
    }
    
    if ( !is_boundary_edge || start_point != start )
    { 
       PST_CoEdge* coedge = edge->end_point() == point ?
                            edge->forward() : edge->reverse();
       imprintResult.append(coedge);
    }

    if( DEBUG_FLAG( VG_FACET_DEBUG ) )
    {
      edge->debug_draw( VG_FACET_IMPRINT_COLOR );
      point->debug_draw( VG_FACET_IMPRINT_COLOR );
    }    
    start_point = point;
  }
  
  return CUBIT_SUCCESS;
}
bool FacetProjectTool::project_to_face ( PST_Face face,
const CubitVector position,
CubitVector result 
) [static, protected]

Definition at line 448 of file FacetProjectTool.cpp.

{
    // First check if projection of the point will be within
    // the bounds of the triangle.  
  if( ! inside_face( triangle, position ) )
    return false;
  
    // Calculate the actual projection
  result = triangle->plane().project( position );
  return true;
}
PST_Point * FacetProjectTool::split_edge ( PST_Edge edge,
double  t 
) [protected]

Definition at line 692 of file FacetProjectTool.cpp.

{
  assert( t > 0.0 && t < 1.0 );
  
    // create point
  PST_Point* point = new PST_Point( edge->position(t) );
  
    // get connected faces 
  PST_Face* face1 = edge->forward()->face();
  PST_Face* face2 = edge->reverse()->face();
  
    // point on each triangle opposite the edge to be split
  PST_Point* face1_pt = face1 ? edge->forward()->next()->end_point() : 0;
  PST_Point* face2_pt = face2 ? edge->reverse()->next()->end_point() : 0;

  PST_Edge* face1_edge = 0, *face2_edge = 0;
  
    // split edge edge
  PST_Edge* split_edge = edge->split( point );
  if( edge->owner() )
    edge->owner()->notify_split( edge, split_edge );
    // if we split a boundary edge, put the new edge in the
    // list of boundary edges too.
  if( edge->mark )
  {
    split_edge->mark = edge->mark;
    boundaryList.append( split_edge );
  }
    // add edge to our list of edges
  edgeList.append( split_edge );
    // If the first face was not a boundary face (different
    // than the boundary of our patch of facets), then split
    // it so that it remains a triangle.
  if( face1 ) 
  {
    face1_edge = PST_Edge::split_face( point, face1_pt, face1 );
    
    if( face1_edge )
    {
      if( DEBUG_FLAG( VG_FACET_DEBUG ) )
        face1_edge->debug_draw( VG_FACET_FACET_COLOR );

      edgeList.append( face1_edge );
        PST_Face* new_face = face1_edge->other(face1);
      new_face->sequence = facetList.size();
      new_face->parent = face1->parent;
      facetList.append(new_face );
      
      if( face1->owner() )
        face1->owner()->notify_split( face1, face1_edge->other(face1) );
    }
  }
    // If the second face is not a boundary face, split it
    // so that we still have triangles.
  if( face2 )
  {
    face2_edge = PST_Edge::split_face( point, face2_pt, face2 );
    if( face2_edge )
    {
      if( DEBUG_FLAG( VG_FACET_DEBUG ) )
        face2_edge->debug_draw( VG_FACET_FACET_COLOR );

      edgeList.append( face2_edge );
        PST_Face* new_face = face2_edge->other(face2);
      new_face->sequence = facetList.size();
      new_face->parent = face2->parent;
      facetList.append( new_face );
      
      if( face2->owner() )
        face2->owner()->notify_split( face2, face2_edge->other(face2) );
    }
  }
  point->sequence = pointList.size();
  pointList.append(point);
  
  return point;
}
bool FacetProjectTool::within_tolerance ( PST_Edge edge,
const CubitVector point,
double  tolerance 
) [static, protected]

Definition at line 1415 of file FacetProjectTool.cpp.

{
  const double t_edge = edge->closest_on_edge( point );
  const CubitVector edge_pos = edge->position(t_edge);
  return (edge_pos - point).length_squared() < tolerance*tolerance;
}

Member Data Documentation

Definition at line 222 of file FacetProjectTool.hpp.

Definition at line 214 of file FacetProjectTool.hpp.

Definition at line 213 of file FacetProjectTool.hpp.

Definition at line 215 of file FacetProjectTool.hpp.

Definition at line 218 of file FacetProjectTool.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines