cgma
Faceter.cpp
Go to the documentation of this file.
00001 //-------------------------------------------------------------------------
00002 // Filename      : Faceter.cpp
00003 //
00004 // Purpose       : Implementation of Faceter
00005 //
00006 // Special Notes : 
00007 //
00008 // Creator       : David White
00009 //
00010 // Creation Date : 03/01/02
00011 //-------------------------------------------------------------------------
00012 
00013 #include "Faceter.hpp"
00014 #include "CubitPoint.hpp"
00015 #include "CubitFacet.hpp"
00016 #include "FaceterPointData.hpp"
00017 #include "FaceterFacetData.hpp"
00018 #include "CubitVector.hpp"
00019 
00020 #include "DLIList.hpp"
00021 
00022 #include "GeometryDefines.h"
00023 #include "CubitDefines.h"
00024 #include "GfxDebug.hpp"
00025 #include "GeometryQueryEngine.hpp"
00026 #include "GMem.hpp"
00027 
00028 #include "RefFace.hpp"
00029 #include "RefEdge.hpp"
00030 #include "CoEdge.hpp"
00031 #include "RefVertex.hpp"
00032 #include "CubitPlane.hpp"
00033 #include "Curve.hpp"
00034 #include "PointGridSearch.hpp"
00035 
00036 //-------------------------------------------------------------------------
00037 // Purpose       : Constructor
00038 //
00039 // Special Notes : 
00040 //
00041 // Creator       : David White
00042 //
00043 // Creation Date : 03/01/02
00044 //-------------------------------------------------------------------------
00045 Faceter::Faceter( RefFace *face ) 
00046     : thisRefFacePtr(face), gridCellScale(2.5)
00047 {
00048   globalPointList = new DLIList <CubitPoint*>;
00049   gridSearchPtr = NULL;
00050   avoidedOverlap = CUBIT_FALSE;
00051 }
00052 //-------------------------------------------------------------------------
00053 // Purpose       : Destructor
00054 //
00055 // Special Notes : 
00056 //
00057 // Creator       : David White
00058 //
00059 // Creation Date : 03/01/02
00060 //-------------------------------------------------------------------------
00061 Faceter::~Faceter()
00062 {
00063   if (gridSearchPtr)
00064     delete gridSearchPtr;
00065   if ( globalPointList )
00066     delete globalPointList;
00067 }
00068 //-------------------------------------------------------------------------
00069 // Purpose       : facet_surface: facets the surface.
00070 //
00071 // Special Notes : 
00072 //
00073 // Creator       : David White
00074 //
00075 // Creation Date : 03/01/02
00076 //-------------------------------------------------------------------------
00077 CubitStatus Faceter::facet_surface(DLIList <CubitFacet*> &results,
00078                                    DLIList <CubitPoint*> &point_list)
00079 {
00080   if ( DEBUG_FLAG(129) )
00081   {
00082     GfxDebug::clear();
00083     GfxDebug::draw_ref_face_edges(thisRefFacePtr);
00084     GfxDebug::flush();
00085     int debug = 0;
00086     if ( debug )
00087     {
00088       GfxDebug::mouse_xforms();
00089       GfxDebug::flush();
00090     }
00091   }
00092   if ( thisRefFacePtr->number_of_Loops() > 1 )
00093     return CUBIT_FAILURE;
00094     //Get the ordered boundary loops.
00095   int ii, jj;
00096   DLIList <DLIList<CubitPoint*>*>  boundary_point_loops;
00097   DLIList <CubitPoint*> *tmp_list_ptr;
00098   CubitStatus stat = get_boundary_points( boundary_point_loops );
00099   if ( stat != CUBIT_SUCCESS )
00100   {
00101       //clean up the data...
00102     for ( ii = boundary_point_loops.size(); ii > 0; ii-- )
00103     {
00104       tmp_list_ptr = boundary_point_loops.pop();
00105       for ( jj = tmp_list_ptr->size(); jj > 0; jj-- )
00106         delete tmp_list_ptr->pop();
00107       delete tmp_list_ptr;
00108     }
00109     return stat;
00110   }
00111     //Set up the gridsearch.
00112   double ratio = gridCellScale, cell_size = 0.0;
00113   max_min_edge_ratio(boundary_point_loops, ratio, cell_size);
00114   if (ratio <= gridCellScale) {
00115     ratio = gridCellScale;
00116   }
00117     //Get all of the points into a single list.
00118   for ( ii = boundary_point_loops.size(); ii > 0; ii-- )
00119   {
00120     tmp_list_ptr = boundary_point_loops.get_and_step();
00121     for ( jj = tmp_list_ptr->size(); jj > 0; jj-- )
00122     {
00123       globalPointList->append(tmp_list_ptr->get_and_step());
00124     }
00125   }
00126   gridSearchPtr = new PointGridSearch(*globalPointList,
00127                                       cell_size,
00128                                       ratio);
00129       //fill in the grid...
00130   for ( ii = globalPointList->size(); ii > 0; ii-- )
00131     gridSearchPtr->add_point(globalPointList->get_and_step());
00132 
00133     //Now start faceting.
00134   stat = facet_loop( boundary_point_loops.get(), results );
00135     //clean up the data...
00136   for ( ii = boundary_point_loops.size(); ii > 0; ii-- )
00137     delete boundary_point_loops.pop();
00138   if ( stat != CUBIT_SUCCESS )
00139   {
00140       //clean the data and return..
00141     for ( ii = results.size(); ii > 0; ii-- )
00142       delete results.pop();
00143     for ( ii = globalPointList->size(); ii > 0; ii-- )
00144       delete globalPointList->pop();
00145     return stat;
00146   }
00147     //Didn't add any points...
00148   point_list += *globalPointList;
00149   return CUBIT_SUCCESS;
00150 }
00151 
00152 
00153 CubitStatus Faceter::facet_loop(DLIList <CubitPoint*> *loop_ptr,
00154                                 DLIList <CubitFacet*> &results)
00155 {
00156     //Assume there is just this loop and facet it as if there was nothing
00157     //in the middle.  This is about the most stupid, simplistic method
00158     //of faceting I could think of.  I'm sure anyone reading this could
00159     //think of something more sophisticated and better.
00160 
00161     //First get the inteior angle for all of the points.
00162   DLIList <FaceterPointData*> ordered_point_list;
00163   int ii;
00164   double angle;
00165   CubitStatus stat;
00166   CubitPoint *curr_point;
00167   FaceterPointData *curr_faceter_data, *prev, *next;
00168   loop_ptr->reset();
00169   int debug = 0;
00170   for ( ii = 0; ii < loop_ptr->size(); ii++ )
00171   {
00172     stat = interior_angle(loop_ptr, angle);
00173     if ( stat == CUBIT_FAILURE )
00174       return CUBIT_FAILURE;
00175     curr_point = loop_ptr->get_and_step();
00176     if ( debug )
00177     {
00178       PRINT_INFO("%d\n",curr_point->id());
00179       GfxDebug::draw_point( curr_point->coordinates(), CUBIT_ORANGE_INDEX);
00180     }
00181     curr_faceter_data = (FaceterPointData*) curr_point;
00182     prev = (FaceterPointData*) loop_ptr->prev(2);
00183     next = (FaceterPointData*) loop_ptr->get();
00184 
00185     if ( angle < CUBIT_RESABS || angle > DEGREES_TO_RADIANS( 359.95 ))
00186     {
00187       
00188       PRINT_ERROR("Zero degree angle in loop\n");
00189       RefEdge *owner_e = (RefEdge*) curr_faceter_data->owner();
00190       RefEdge *prev_e = (RefEdge*) prev->owner();
00191       RefEdge *next_e = (RefEdge*) next->owner();
00192       PRINT_ERROR("On edge %d, with next %d and prev %d\n", owner_e->id(),
00193                   next_e->id(), prev_e->id());
00194       return CUBIT_FAILURE;
00195     }
00196     curr_faceter_data->set_interior_angle(angle);
00197     curr_faceter_data->set_prev(prev);
00198     curr_faceter_data->set_next(next);
00199     ordered_point_list.append(curr_faceter_data);
00200   }
00201   ordered_point_list.sort(FaceterPointData::sort_by_angle);
00202 
00203     //Now go through and just cut off the smallest triangles...
00204     //The list should be ordered largest angle to smallest.
00205     //so go from the end to the front.
00206   CubitFacet *new_facet = NULL;
00207   int safe_size = 4* ordered_point_list.size();
00208   int safe_count = 0;
00209   int redo_count = ordered_point_list.size();
00210   while (ordered_point_list.size() >= 3 &&
00211          safe_count < safe_size )
00212   {
00213     curr_faceter_data = ordered_point_list.pop();
00214     new_facet = NULL;
00215     stat = facet_factory(curr_faceter_data, new_facet, ordered_point_list);
00216     if ( new_facet != NULL )
00217       results.append(new_facet);
00218     if ( stat != CUBIT_SUCCESS )
00219       return stat;
00220     if ( safe_count == redo_count )
00221     {
00222       if ( avoidedOverlap )
00223       {
00224         CubitBoolean reorder = CUBIT_FALSE;
00225         for ( ii = ordered_point_list.size(); ii > 0; ii-- )
00226         {
00227           curr_faceter_data = ordered_point_list.get_and_step();
00228           angle = curr_faceter_data->get_interior_angle();
00229           if ( angle > 2*CUBIT_PI )
00230           {
00231             angle = angle - 2*CUBIT_PI;
00232             curr_faceter_data->set_interior_angle(angle);
00233             reorder = CUBIT_TRUE;
00234           }
00235         }
00236         if ( reorder )
00237           ordered_point_list.sort(FaceterPointData::sort_by_angle);
00238           //set the next redo to be later.
00239         avoidedOverlap = CUBIT_FALSE;
00240       }
00241       redo_count += redo_count;
00242     }
00243     safe_count++;
00244   }
00245   if ( ordered_point_list.size() < 3 )
00246     return CUBIT_SUCCESS;
00247   else
00248     return CUBIT_FAILURE;
00249 }
00250 CubitStatus Faceter::facet_factory(FaceterPointData *curr_faceter_data,
00251                                    CubitFacet *&new_facet,
00252                                    DLIList <FaceterPointData*> &order_list)
00253 {
00254   double angle;
00255   DLIList<CubitPoint*> loop_list;
00256   new_facet = NULL;
00257   if ( avoid_facet(curr_faceter_data, order_list) )
00258     return CUBIT_SUCCESS;
00259   FaceterPointData *prev = curr_faceter_data->get_prev();
00260   FaceterPointData *next = curr_faceter_data->get_next();
00261     //create the new facet.
00262     //create a triangle from these three points.
00263   new_facet = (CubitFacet*) new FaceterFacetData((CubitPoint*)prev,
00264                                                  (CubitPoint*)curr_faceter_data,
00265                                                  (CubitPoint*)next);
00266   CubitVector facet_normal = new_facet->normal();
00267   CubitVector surface_normal = thisRefFacePtr->normal_at(curr_faceter_data->coordinates());
00268   double dot = facet_normal%surface_normal;
00269   if ( dot < 0 )
00270   {
00271       //new need to try and do this in a different way...
00272     delete new_facet;
00273     new_facet = NULL;
00274     if ( order_list.size() < 3 )
00275     {
00276       PRINT_DEBUG_129("Last triangle is inverted. (surface %d)\n",
00277                       thisRefFacePtr->id());
00278       return CUBIT_FAILURE;
00279     }
00280       //Try and create facets on either side of this facet.
00281     return CUBIT_FAILURE;
00282   }
00283   
00284   int debug = 1;
00285   if ( debug )
00286     GfxDebug::draw_facet(new_facet, CUBIT_RED_INDEX);
00287     //Now fix up the ordered_list.
00288     //First update the angles.
00289   prev->set_next(next);
00290   next->set_prev(prev);
00291   loop_list.clean_out();
00292   loop_list.append(prev->get_prev());
00293   loop_list.append(prev);
00294   loop_list.append(next);
00295   loop_list.step();
00296   interior_angle(&loop_list, angle);
00297   prev->set_interior_angle(angle);
00298   loop_list.clean_out();
00299   loop_list.append(prev);
00300   loop_list.append(next);
00301   loop_list.append(next->get_next());
00302   loop_list.step();
00303   interior_angle(&loop_list, angle);
00304   next->set_interior_angle(angle);
00305     //Now do basically a sort again.  There is probably
00306     //a faster way to do this but there are also much slower
00307     //ways in worst case too. At least this has guaranteed O(nlogn).
00308   order_list.sort(FaceterPointData::sort_by_angle);
00309   CubitPoint *tmp_point = (CubitPoint*) curr_faceter_data;
00310   gridSearchPtr->remove_point(tmp_point);
00311   return CUBIT_SUCCESS;
00312 }
00313 CubitBoolean Faceter::avoid_facet(FaceterPointData *curr_faceter_data,
00314                                   DLIList<FaceterPointData*> &order_list)
00315 {
00316   
00317   int ii;
00318   CubitPoint *tmp_point;
00319   FaceterPointData *prev = curr_faceter_data->get_prev();
00320   FaceterPointData *next = curr_faceter_data->get_next();
00321     //get the closest edges to this point.
00322   CubitVector curr_v = curr_faceter_data->coordinates();
00323   CubitVector prev_v = prev->coordinates();
00324   CubitVector next_v = next->coordinates();
00325   double r1 = (curr_v - prev_v).length();
00326   double r2 = (curr_v - next_v).length();
00327   double radius;
00328   double angle;
00329   if ( r1 > r2 )
00330     radius = r1;
00331   else
00332     radius = r2;
00333   gridSearchPtr->set_neighborhood_bounds(curr_v, radius);
00334   DLIList<CubitPoint*> neighborhood_points;
00335   DLIList<CubitPoint*> close_points, loop_list;
00336   gridSearchPtr->get_neighborhood_points(neighborhood_points);
00337     //Loop through and see if there are nodes in here closer
00338     //than the r1 or r2.
00339   FaceterPointData *tmp_faceter;
00340   CubitVector tmp_v;
00341   radius *= radius;
00342   double angle1, angle2;
00343   for (ii = neighborhood_points.size(); ii > 0; ii-- )
00344   {
00345     tmp_point = neighborhood_points.get_and_step();
00346     tmp_faceter = (FaceterPointData*) tmp_point;
00347     if ( tmp_faceter == curr_faceter_data ||
00348          tmp_faceter == prev ||
00349          tmp_faceter == next ||
00350          tmp_faceter == next->get_next() ||
00351          tmp_faceter == prev->get_prev() )
00352       continue;
00353     tmp_v = tmp_faceter->coordinates();
00354     if ( (curr_v-tmp_v).length_squared() < radius )
00355       close_points.append(tmp_point);
00356   }
00357     //Now of these close points if any of them make
00358     //good trinagles for the other nodes then don't do
00359     //this triangle...
00360   CubitBoolean abort_facet = CUBIT_FALSE;
00361   for ( ii = close_points.size(); ii > 0; ii-- )
00362   {
00363     tmp_point = close_points.get_and_step();
00364     loop_list.clean_out();
00365     loop_list.append(tmp_point);
00366     loop_list.append((CubitPoint*)prev);
00367     loop_list.append((CubitPoint*)curr_faceter_data);
00368     loop_list.step();
00369     interior_angle(&loop_list, angle1);
00370     loop_list.clean_out();
00371     loop_list.append((CubitPoint*)curr_faceter_data);
00372     loop_list.append((CubitPoint*)next);
00373     loop_list.append(tmp_point);
00374     loop_list.step();
00375     interior_angle(&loop_list, angle2);
00376     if ( angle1 < DEGREES_TO_RADIANS(160) &&
00377          angle2 < DEGREES_TO_RADIANS(160))
00378     {
00379       int tmp_debug = 1;
00380       if ( tmp_debug )
00381       {
00382         GfxDebug::draw_point(prev->coordinates(), CUBIT_GREEN_INDEX);
00383         GfxDebug::draw_point(curr_faceter_data->coordinates(), CUBIT_RED_INDEX);
00384         GfxDebug::draw_point(tmp_point->coordinates(), CUBIT_BLUE_INDEX);
00385         GfxDebug::draw_point(next->coordinates(), CUBIT_MAGENTA_INDEX);
00386         loop_list.clean_out();
00387         loop_list.append(tmp_point);
00388         loop_list.append((CubitPoint*)prev);
00389         loop_list.append((CubitPoint*)curr_faceter_data);
00390         loop_list.step();
00391         interior_angle(&loop_list, angle1);
00392       }
00393       abort_facet = CUBIT_TRUE;
00394       break;
00395     }
00396   }
00397   if ( abort_facet )
00398   {
00399       //We need to add curr_facerter_data back to the list.
00400     order_list.insert(curr_faceter_data);
00401       //Set the angle above 360...
00402     angle = curr_faceter_data->get_interior_angle();
00403     curr_faceter_data->set_interior_angle(angle + 2*CUBIT_PI);
00404     avoidedOverlap = CUBIT_TRUE;
00405     return CUBIT_TRUE;
00406   }
00407   return CUBIT_FALSE;
00408 }
00409 
00410 CubitStatus Faceter::interior_angle(DLIList <CubitPoint*> *loop_ptr,
00411                                     double &my_angle)
00412 {
00413   CubitVector point = loop_ptr->get()->coordinates();
00414   CubitVector to_prev = loop_ptr->prev()->coordinates();
00415   to_prev  -= point;
00416   CubitVector to_next = loop_ptr->next()->coordinates();
00417   to_next -= point;
00418   CubitVector surf_norm = thisRefFacePtr->normal_at ( point );
00419   if ( surf_norm.length_squared() < CUBIT_RESABS )
00420   {
00421       //Try getting it at one of the other nodes...
00422     surf_norm = thisRefFacePtr->normal_at(loop_ptr->next()->coordinates() );
00423     if (surf_norm.length_squared() < CUBIT_RESABS )
00424     {
00425         //Try getting it at one of the other nodes...
00426       surf_norm = thisRefFacePtr->normal_at(loop_ptr->prev()->coordinates() );
00427       if (surf_norm.length_squared() < CUBIT_RESABS )
00428       {
00429         PRINT_ERROR("Trying to get normal at point %d on surf %d.\n"
00430                     "       Normal length being returned equals zero.\n",
00431                     loop_ptr->get()->id(), thisRefFacePtr->id() );
00432         return CUBIT_FAILURE;
00433       }
00434     }
00435   }
00436   double angle = surf_norm.vector_angle ( to_next, to_prev );
00437 
00438     //Now if the angle is very small (zero) we don't know if this
00439     //should be 0 or 2PI, so do some tests.
00440   const double NEAR_ZERO = DEGREES_TO_RADIANS( 0.15 );
00441   const double NEAR_360  = DEGREES_TO_RADIANS( 359.95 );
00442   if( NEAR_ZERO < angle && angle < NEAR_360 )
00443   {
00444     my_angle = angle;
00445     return CUBIT_SUCCESS;
00446   }
00447     //If we are close to zero or 2PI then take some other loop nodes
00448     //into consideration to see if we can get a better approximation
00449     //as to we are zero or 2 PI.
00450   CubitPoint *prev_node = loop_ptr->prev();
00451   CubitPoint *next_node = loop_ptr->next();
00452   int stop_it = loop_ptr->size();
00453   int count = 0;
00454   while (( angle < NEAR_ZERO || angle > NEAR_360 ) && ( count < stop_it ) )
00455   {
00456     prev_node = loop_ptr->prev( count + 2 );
00457     next_node = loop_ptr->next( count + 2 );
00458     to_prev = prev_node->coordinates() - point;
00459     to_next = next_node->coordinates() - point;
00460     angle = surf_norm.vector_angle ( to_next, to_prev );
00461     count++;
00462   }
00463   const double TWO_PI = 2.0 * CUBIT_PI;
00464   if ( angle < NEAR_ZERO )
00465     my_angle = TWO_PI;
00466   else if( angle < CUBIT_PI )
00467     my_angle = 0.0;
00468   else
00469     my_angle = TWO_PI;
00470   return CUBIT_SUCCESS;
00471 }
00472 void Faceter::max_min_edge_ratio( DLIList <DLIList <CubitPoint*>*> &boundary_point_loops,
00473                                   double &ratio,
00474                                   double &cell_size)
00475 {
00476   double max_edge_length = CUBIT_DBL_MIN;
00477   double min_edge_length = CUBIT_DBL_MAX;
00478   double sum = 0.0;
00479   int total_count = 0;
00480   DLIList<CubitPoint*> *loopPtr;
00481   CubitPoint *node1, *node2;
00482   CubitVector edge;
00483   for (int i = boundary_point_loops.size(); i > 0; i--) {
00484     loopPtr = boundary_point_loops.get_and_step();
00485     node2 = loopPtr->prev();
00486     for (int j = loopPtr->size(); j > 0; j--) {
00487       node1 = node2;
00488       node2 = loopPtr->get_and_step();
00489       CubitVector p1 = node1->coordinates();
00490       CubitVector p2 = node2->coordinates();
00491       edge = p2-p1;
00492       double edge_length = edge.length_squared();
00493       if (edge_length > max_edge_length) max_edge_length = edge_length;
00494       if (edge_length < min_edge_length) min_edge_length = edge_length;
00495       total_count++;
00496       sum += edge_length;
00497     }
00498   }
00499 
00500   if (min_edge_length > CUBIT_RESABS) {
00501     ratio = sqrt(max_edge_length/min_edge_length);
00502   }
00503   else {
00504     ratio = sqrt(max_edge_length);
00505   }
00506   if ( total_count > 0 && sum > 0 )
00507     cell_size = sqrt(sum/total_count);
00508   else
00509     cell_size = 0.0;
00510 }
00511 
00512 
00513 CubitStatus Faceter::get_boundary_points( DLIList <DLIList<CubitPoint*>*> &boundary_point_loops ) const
00514 {
00515   DLIList<DLIList<CoEdge*> > co_edge_loops;
00516   thisRefFacePtr->co_edge_loops(co_edge_loops);
00517   int ii, jj;
00518   DLIList <CubitPoint*> *new_point_loop_ptr, tmp_point_list;
00519   RefEdge *ref_edge_ptr;
00520   CoEdge *co_edge_ptr;
00521   CubitStatus stat;
00522   CubitSense sense;
00523   
00524   for ( ii = co_edge_loops.size(); ii > 0; ii-- )
00525   {
00526     DLIList <CoEdge*> &co_edge_list_ptr = co_edge_loops.get_and_step();
00527     new_point_loop_ptr = new DLIList <CubitPoint*>;
00528     for ( jj = co_edge_list_ptr.size(); jj > 0; jj-- )
00529     {
00530       co_edge_ptr = co_edge_list_ptr.get_and_step();
00531       ref_edge_ptr = co_edge_ptr->get_ref_edge_ptr();
00532       tmp_point_list.clean_out();
00533       stat = get_curve_facets( ref_edge_ptr, tmp_point_list );
00534       PRINT_DEBUG_129("curve %d has %d points\n",
00535                       ref_edge_ptr->id(),
00536                       tmp_point_list.size());
00537       if ( !stat )
00538         return CUBIT_FAILURE;
00539       tmp_point_list.reset();
00540         //the points are in order from start vertex to end vertex.
00541         //append them now according to the loop.
00542       sense = co_edge_ptr->get_sense();
00543       if ( CUBIT_FORWARD != sense )
00544         tmp_point_list.reverse();
00545         //Now take off the last point as it is a duplicate with the
00546         //other list...
00547       tmp_point_list.reset();
00548       delete tmp_point_list.pop();
00549       (*new_point_loop_ptr) += tmp_point_list;
00550     }
00551     CubitVector curr, prev;
00552     for ( jj = new_point_loop_ptr->size(); jj > 0; jj-- )
00553     {
00554       prev = new_point_loop_ptr->prev()->coordinates();
00555       curr = new_point_loop_ptr->get_and_step()->coordinates();
00556       if ( prev.about_equal(curr) )
00557       {
00558         PRINT_DEBUG_129("Points within tolerance in boundaryloop.\n");
00559         new_point_loop_ptr->back();
00560         delete new_point_loop_ptr->remove();
00561       }
00562     }
00563     boundary_point_loops.append(new_point_loop_ptr);
00564   }
00565   
00566   return CUBIT_SUCCESS;
00567 }
00568 
00569 CubitStatus Faceter::get_curve_facets( RefEdge* curve, DLIList<CubitPoint*>& segments ) const
00570 {
00571 //const double COS_ANGLE_TOL =  0.965925826289068312213715; // cos(15)
00572   const double COS_ANGLE_TOL =  0.984807753012208020315654; // cos(10)
00573 //const double COS_ANGLE_TOL =  0.996194698091745545198705; // cos(5)
00574   GMem curve_graphics;
00575   const double dist_tol = GEOMETRY_RESABS;
00576   const double dist_tol_sqr = dist_tol*dist_tol;
00577   Curve* curve_ptr = curve->get_curve_ptr();
00578   curve_ptr->get_geometry_query_engine()->get_graphics( curve_ptr, &curve_graphics );
00579   
00580   GPoint* gp = curve_graphics.point_list();
00581   CubitPoint* last = (CubitPoint*) new FaceterPointData( gp->x, gp->y, gp->z );
00582   ((FaceterPointData*)last)->owner(dynamic_cast<RefEntity*>(curve));
00583   CubitVector lastv = last->coordinates();
00584   segments.append( last );
00585   GPoint* end = gp + curve_graphics.pointListCount - 1;
00586   
00587   for( gp++; gp < end; gp++ )
00588   {
00589     CubitVector pos(  gp->x, gp->y, gp->z );
00590     CubitVector step1 = (pos - lastv);
00591     double len1 = step1.length();
00592     if( len1 < dist_tol ) continue;
00593     
00594     GPoint* np = gp + 1;
00595     CubitVector next( np->x, np->y, np->z );
00596     CubitVector step2 = next - pos;
00597     double len2 = step2.length();
00598     if( len2 < dist_tol ) continue;
00599     
00600     double cosine = (step1 % step2) / (len1 * len2);
00601     if( cosine > COS_ANGLE_TOL ) continue;
00602     
00603     last = new FaceterPointData( pos );
00604     ((FaceterPointData*)last)->owner(dynamic_cast<RefEntity*>(curve));
00605     segments.append( last );
00606     lastv = last->coordinates();
00607   }
00608   
00609   CubitVector last_pos( gp->x, gp->y, gp->z );
00610   segments.last();
00611   while( (last_pos - (segments.get()->coordinates())).length_squared() < dist_tol_sqr )
00612   {
00613     delete segments.pop();
00614     segments.last();
00615   }
00616   CubitPoint *tmp_point = (CubitPoint*) new FaceterPointData( last_pos );
00617   segments.append( tmp_point );
00618   ((FaceterPointData*)tmp_point)->owner( dynamic_cast<RefEntity*>(curve) );
00619     
00620   // Now check if the segment list is reversed wrt the curve direction.
00621   segments.reset();
00622   double u1, u2;
00623   if( segments.size() > 2 )
00624   {
00625     u1 = curve->u_from_position( (segments.next(1)->coordinates()) );
00626     u2 = curve->u_from_position( (segments.next(2)->coordinates()) );    
00627   }
00628   else
00629   {
00630     u1 = curve->u_from_position( (segments.get()->coordinates() ) );
00631     u2 = curve->u_from_position( (segments.next()->coordinates()) );
00632   }
00633   if( (u2 < u1) && (curve->start_param() <= curve->end_param()) )
00634     segments.reverse();
00635 
00636     //Make sure we don't have duplicate points.
00637   int jj;
00638   CubitVector curr, prev;
00639   for ( jj = segments.size(); jj > 0; jj-- )
00640   {
00641     prev = segments.prev()->coordinates();
00642     curr = segments.get_and_step()->coordinates();
00643     if ( prev.about_equal(curr) )
00644     {
00645       PRINT_DEBUG_129("Points on curve %d within tolerance...\n", curve->id());
00646       segments.back();
00647       delete segments.remove();
00648     }
00649   }
00650   return CUBIT_SUCCESS;
00651 }
00652 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines