cgma
GeomMeasureTool.cpp
Go to the documentation of this file.
00001 //----------------------------------------------------------
00002 //Class: GeomMeasureTool
00003 //Description: Does Geometric measureing and statistical
00004 //             gathering on various topologies.
00005 //Creator: David R. White
00006 //Date: 7/9/2002
00007 //----------------------------------------------------------
00008 
00009 
00010 #include "GeomMeasureTool.hpp"
00011 #include "RefEntity.hpp"
00012 #include "RefVertex.hpp"
00013 #include "RefEdge.hpp"
00014 #include "RefFace.hpp"
00015 #include "RefVolume.hpp"
00016 #include "Body.hpp"
00017 #include "RefGroup.hpp"
00018 #include "RefEntityFactory.hpp"
00019 #include "CoEdge.hpp"
00020 #include "Loop.hpp"
00021 #include "Shell.hpp"
00022 #include "Chain.hpp"
00023 #include "CastTo.hpp"
00024 #include "GMem.hpp"
00025 #include "Curve.hpp"
00026 #include "Surface.hpp"
00027 #include "GeometryQueryEngine.hpp"
00028 #include "GeomSeg.hpp"
00029 #include "GeomPoint.hpp"
00030 #include "RTree.hpp"
00031 #include "AbstractTree.hpp"
00032 #include "IntersectionTool.hpp"
00033 #include "CpuTimer.hpp"
00034 #include "GeometryQueryTool.hpp"
00035 #include "GeometryModifyTool.hpp"
00036 #include "MergeTool.hpp"
00037 #include "GeomDataObserver.hpp"
00038 #include "Lump.hpp"
00039 #include "ModelQueryEngine.hpp"
00040 #include "ProgressTool.hpp"
00041 #include "AppUtil.hpp"
00042 
00043 #include <set>
00044 #include <map>
00045 #include <stack>
00046 #include <algorithm>
00047 
00048 void GeomMeasureTool::get_edges_from_list(DLIList <RefEntity*> &entity_list,
00049                                           DLIList <RefEdge*> &ref_edges )
00050 {
00051     //assume the mark is zero.
00052   int ii,jj;
00053   RefEntity *entity_ptr;
00054   RefEdge *ref_edge;
00055   
00056   TopologyEntity *topo_ent;
00057   DLIList <RefEdge*> tmp_edges;
00058   for ( ii = entity_list.size(); ii > 0; ii-- )
00059   {
00060     tmp_edges.clean_out();
00061     entity_ptr = entity_list.get_and_step();
00062     topo_ent = CAST_TO(entity_ptr, TopologyEntity);
00063     if ( topo_ent == NULL )
00064       continue;
00065     topo_ent->ref_edges(tmp_edges);
00066     for ( jj = tmp_edges.size(); jj > 0; jj-- )
00067     {
00068       ref_edge = tmp_edges.pop();
00069       if ( !ref_edge->marked() )
00070       {
00071         ref_edge->marked(1);
00072         ref_edges.append(ref_edge);
00073       }
00074     }
00075   }
00076   for ( ii = ref_edges.size(); ii > 0; ii-- )
00077     ref_edges.get_and_step()->marked(0);
00078 }
00079 
00080 void GeomMeasureTool::get_bodies_from_list(DLIList <RefVolume*> &entity_list,
00081                                            DLIList <Body*> &ref_bodies )
00082 {
00083     //assume the mark is zero.
00084   int ii, jj;
00085   RefVolume *entity_ptr;
00086   Body *ref_body;
00087   
00088   TopologyEntity *topo_ent;
00089   DLIList <Body*> tmp_bodies;
00090   for ( ii = entity_list.size(); ii > 0; ii-- )
00091   {
00092     tmp_bodies.clean_out();
00093     entity_ptr = entity_list.get_and_step();
00094     topo_ent = CAST_TO(entity_ptr, TopologyEntity);
00095     if ( topo_ent == NULL )
00096       continue;
00097     topo_ent->bodies(tmp_bodies);
00098     for ( jj = tmp_bodies.size(); jj > 0; jj-- )
00099     {
00100       ref_body = tmp_bodies.pop();
00101       if ( !ref_body->marked() )
00102       {
00103         ref_body->marked(1);
00104         ref_bodies.append(ref_body);
00105       }
00106     }
00107   }
00108   for ( ii = ref_bodies.size(); ii > 0; ii-- )
00109     ref_bodies.get_and_step()->marked(0);
00110 }
00111 
00112 void GeomMeasureTool::get_faces_from_list(DLIList <RefEntity*> &entity_list,
00113                                           DLIList <RefFace*> &ref_faces )
00114 {
00115     //assume the mark is zero.
00116   int ii,jj;
00117   RefEntity *entity_ptr;
00118   RefFace *ref_face;
00119   
00120   TopologyEntity *topo_ent;
00121   DLIList <RefFace*> tmp_faces;
00122   for ( ii = entity_list.size(); ii > 0; ii-- )
00123   {
00124     tmp_faces.clean_out();
00125     entity_ptr = entity_list.get_and_step();
00126     topo_ent = CAST_TO(entity_ptr, TopologyEntity);
00127     topo_ent->ref_faces(tmp_faces);
00128     for ( jj = tmp_faces.size(); jj > 0; jj-- )
00129     {
00130       ref_face = tmp_faces.pop();
00131       if ( !ref_face->marked() )
00132       {
00133         ref_face->marked(1);
00134         ref_faces.append(ref_face);
00135       }
00136     }
00137   }
00138   for ( ii = ref_faces.size(); ii > 0; ii-- )
00139     ref_faces.get_and_step()->marked(0);
00140 }
00141 void GeomMeasureTool::get_volumes_from_list(DLIList <RefEntity*> &entity_list,
00142                                             DLIList <RefVolume*> &ref_volumes )
00143 {
00144     //assume the mark is zero.
00145   int ii,jj;
00146   RefEntity *entity_ptr;
00147   RefVolume *ref_volume;
00148   
00149   TopologyEntity *topo_ent;
00150   DLIList <RefVolume*> tmp_volumes;
00151   for ( ii = entity_list.size(); ii > 0; ii-- )
00152   {
00153     tmp_volumes.clean_out();
00154     entity_ptr = entity_list.get_and_step();
00155     topo_ent = CAST_TO(entity_ptr, TopologyEntity);
00156     topo_ent->ref_volumes(tmp_volumes);
00157     for ( jj = tmp_volumes.size(); jj > 0; jj-- )
00158     {
00159       ref_volume = tmp_volumes.pop();
00160       if ( !ref_volume->marked() )
00161       {
00162         ref_volume->marked(1);
00163         ref_volumes.append(ref_volume);
00164       }
00165     }
00166   }
00167   for ( ii = ref_volumes.size(); ii > 0; ii-- )
00168     ref_volumes.get_and_step()->marked(0);
00169 }
00170 void GeomMeasureTool::measure_curve_length(DLIList <RefEntity*> &entity_list,
00171                                            double &smallest,
00172                                            RefEdge *&smallest_edge,
00173                                            double &largest,
00174                                            RefEdge *&largest_edge,
00175                                            double &average,
00176                                            double &sum)
00177 {
00178   DLIList<RefEdge*> ref_edges;
00179   get_edges_from_list(entity_list, ref_edges);
00180   measure_curve_length(ref_edges, smallest, smallest_edge,
00181                        largest, largest_edge, average, sum);
00182   
00183 }
00184 
00185 void GeomMeasureTool::measure_curve_length(DLIList <RefEdge*> &ref_edges,
00186                                            double &smallest,
00187                                            RefEdge *&smallest_edge,
00188                                            double &largest,
00189                                            RefEdge *&largest_edge,
00190                                            double &average,
00191                                            double &sum)
00192 {
00193   int ii;
00194   double curr_length, total_length = 0.0;
00195   smallest = CUBIT_DBL_MAX;
00196   largest = 0.0;
00197   RefEdge *curr_edge;
00198   for ( ii = ref_edges.size(); ii > 0; ii-- )
00199   {
00200     curr_edge = ref_edges.get_and_step();
00201     if( curr_edge->geometry_type() == POINT_CURVE_TYPE )
00202       continue;
00203     curr_length = curr_edge->measure();
00204     if ( curr_length < smallest )
00205     {
00206       smallest = curr_length;
00207       smallest_edge = curr_edge;
00208     }
00209     if ( curr_length > largest )
00210     {
00211       largest = curr_length;
00212       largest_edge = curr_edge;
00213     }
00214     total_length += curr_length;
00215   }
00216   if ( ref_edges.size() != 0 )
00217     average = total_length/ref_edges.size();
00218   sum = total_length;
00219   
00220 }
00221 void GeomMeasureTool::measure_face_curves( RefFace *ref_face,
00222                                            double &min_curve_length,
00223                                            double &max_curve_ratio,
00224                                            RefEdge *&min_ref_edge)
00225 {
00226   DLIList<DLIList<CoEdge*> > co_edge_loops;
00227   ref_face->co_edge_loops(co_edge_loops);
00228   CoEdge *curr_co_edge, *next_co_edge;
00229   min_curve_length = CUBIT_DBL_MAX;
00230   double curve_length, next_curve_length, ratio;
00231   max_curve_ratio = -CUBIT_DBL_MAX;
00232   CubitBoolean ratio_set = CUBIT_FALSE;
00233   int ii, jj;
00234   for ( ii = co_edge_loops.size(); ii > 0; ii-- )
00235   {
00236     DLIList <CoEdge*> &co_edge_list = co_edge_loops.get_and_step();
00237     curve_length = co_edge_list.get()->get_ref_edge_ptr()->measure();
00238     for ( jj = co_edge_list.size(); jj > 0; jj-- )
00239     {
00240       curr_co_edge = co_edge_list.get_and_step();
00241       next_co_edge = co_edge_list.get();
00242       next_curve_length = next_co_edge->get_ref_edge_ptr()->measure();
00243         //find the smallest curve
00244       if ( curve_length < min_curve_length )
00245       {
00246         min_curve_length = curve_length;
00247         min_ref_edge = curr_co_edge->get_ref_edge_ptr();
00248       }
00249       if ( co_edge_list.size() > 1 )
00250       {
00251           //find the biggest change in curve lengths.
00252         if ( curve_length > next_curve_length )
00253         {
00254           if ( next_curve_length < CUBIT_RESABS )
00255             ratio = CUBIT_DBL_MAX;
00256           else
00257             ratio = curve_length / next_curve_length;
00258         }
00259         else
00260         {
00261           if ( curve_length < CUBIT_RESABS )
00262             ratio = CUBIT_DBL_MAX;
00263           else
00264             ratio = next_curve_length /curve_length;
00265         }
00266         curve_length = next_curve_length;
00267         if ( ratio > max_curve_ratio )
00268         {
00269           ratio_set = CUBIT_TRUE;
00270           max_curve_ratio = ratio;
00271         }
00272       }
00273     }
00274   }
00275   if ( !ratio_set )
00276     max_curve_ratio = -1.0;
00277 }
00278 
00279         
00280   
00281   
00282 CubitStatus GeomMeasureTool::measure_face_loops( RefFace *face,
00283                                                  double &min_distance_between_loops,
00284                                                  double &min_distance_in_one_loop,
00285                                                  double &min_angle,
00286                                                  double &max_angle,
00287                                                  double tolerance)
00288 {
00289     //This is the most tricky of these functions.
00290     //To do this I'm going to facet the edges, then use an AbstractTree to
00291     //find the minimum distance between the loops and between a single
00292     //loop.
00293   PointLoopList boundary_point_loops;
00294   CubitStatus stat;
00295   stat = get_boundary_points(face, boundary_point_loops, tolerance);
00296   if ( stat != CUBIT_SUCCESS )
00297     return CUBIT_FAILURE;
00298   
00299   SegLoopList boundary_seg_loops;
00300   stat = convert_to_lines( boundary_point_loops,
00301                            boundary_seg_loops,
00302                            face);
00303   if ( stat != CUBIT_SUCCESS )
00304     return stat;
00305 
00306     //Now add the points to the AbstractTree.
00307   DLIList<AbstractTree<GeomSeg*>*> atree_list;
00308   AbstractTree<GeomSeg*> *curr_tree;
00309   SegList *seg_list;
00310   GeomSeg  *curr_seg, *next_seg;
00311   int ii,jj, kk;
00312   double angle;
00313   min_angle = CUBIT_DBL_MAX;
00314   max_angle = 0.0;
00315   double creation_time=0.0;
00316   CpuTimer creation_timer;
00317   const double angle_convert = 180.0/CUBIT_PI;
00318   boundary_seg_loops.reset();
00319   int num_segs = 0;
00320   for (ii = 0; ii < boundary_seg_loops.size(); ii++ )
00321   {
00322     seg_list = boundary_seg_loops.get_and_step();
00323     creation_timer.cpu_secs();
00324     curr_tree = new RTree<GeomSeg*>(GEOMETRY_RESABS);
00325     creation_time += creation_timer.cpu_secs();
00326     for ( jj = seg_list->size(); jj > 0; jj-- )
00327     {
00328         //build the r-tree.
00329       num_segs++;
00330       creation_timer.cpu_secs();
00331       curr_seg = seg_list->get_and_step();
00332       curr_tree->add(curr_seg);
00333       creation_time += creation_timer.cpu_secs();
00334         //calculate the interior angles.
00335       next_seg = seg_list->get();
00336 
00337       RefEntity* temp_ent = curr_seg->get_end()->owner();
00338         //only measure the angle at vertices.  This
00339         //should make things slightly faster.
00340       if ( CAST_TO(temp_ent, RefVertex) )
00341       {
00342         stat = interior_angle(curr_seg, next_seg,
00343                               face, angle);
00344         if ( stat != CUBIT_SUCCESS )
00345         {
00346           min_angle = 0.0;
00347           max_angle = 360.0;
00348         }
00349         else if ( angle < min_angle )
00350         {
00351           min_angle = angle;
00352         }
00353         else if ( angle > max_angle )
00354         {
00355           max_angle = angle;
00356         }
00357       }
00358     }
00359     atree_list.append(curr_tree);
00360   }
00361   min_angle *= angle_convert;
00362   max_angle *= angle_convert;
00363 
00364   if ( atree_list.size() <= 1 )
00365     min_distance_between_loops = -1.0;
00366   else
00367   {
00368       //determine the minimum distance between the loops.
00369     CpuTimer atree_time;
00370     atree_time.cpu_secs();
00371     PointList *curr_points;
00372     GeomPoint *curr_point;
00373     DLIList <GeomSeg*> nearest_neighbors;
00374     
00375     CubitVector curr_loc;
00376     double closest_dist;
00377     min_distance_between_loops = CUBIT_DBL_MAX;
00378     atree_list.reset();
00379     for ( ii = 0; ii < atree_list.size(); ii++ )
00380     {
00381       curr_points = boundary_point_loops.get_and_step();
00382       for ( jj = ii+1; jj < atree_list.size(); jj++ )
00383       {
00384         curr_tree = atree_list.next(jj);
00385           //now for every point in curr_points, find it's closest_point
00386           //in the curr_tree.
00387         for ( kk = 0; kk < curr_points->size(); kk++ )
00388         {
00389           curr_point = curr_points->get_and_step();
00390           curr_loc.set(curr_point->coordinates());
00391           nearest_neighbors.clean_out();
00392           curr_tree->k_nearest_neighbor(curr_loc,
00393                                         1, closest_dist, nearest_neighbors,
00394                                         dist_sq_point_data);
00395           if ( nearest_neighbors.size() == 0)
00396           {
00397               //hmm, not sure what is wrong here.
00398             PRINT_ERROR("Can't find closest point between loops.\n");
00399             return CUBIT_FAILURE;
00400           }
00401           if (closest_dist < min_distance_between_loops)
00402             min_distance_between_loops = closest_dist;
00403         }
00404       }
00405     }
00406     double time =  atree_time.cpu_secs();
00407     PRINT_DEBUG_139("Time to create atree: %f\n", creation_time);
00408     PRINT_DEBUG_139("Time for atree nn search: %f\n", time);
00409   }
00410   if ( min_distance_between_loops != -1.0 )
00411   min_distance_between_loops = sqrt(min_distance_between_loops);
00412     //Now find the min distance inside a loop.  This is tricky
00413     //in that we want to find distances that are across some area
00414     //and not just between adjacent curves.
00415   atree_list.reset();
00416   GeomSeg *tmp_seg;
00417   double dist;
00418   min_distance_in_one_loop = CUBIT_DBL_MAX;
00419   int k=5;
00420   CubitVector curr_loc;
00421   DLIList <GeomSeg*> nearest_neighbors;
00422   RefEntity *ent_1, *ent_2;
00423   TopologyEntity *top_1, *top_2;
00424   
00425   
00426   for ( ii = atree_list.size(); ii > 0; ii-- )
00427   {
00428     curr_tree = atree_list.get_and_step();
00429     seg_list = boundary_seg_loops.get_and_step();
00430     for ( jj = seg_list->size(); jj > 0; jj-- )
00431     {
00432       nearest_neighbors.clean_out();
00433       curr_seg = seg_list->get_and_step();
00434         //get the k closest segments.
00435       curr_loc.set(curr_seg->get_start()->coordinates());
00436       curr_tree->k_nearest_neighbor(curr_loc,
00437                                     k, dist, nearest_neighbors,
00438                                     dist_sq_point_data);
00439         //go through these nearest_neighbors and through out
00440         //the segments that are attached to curr_seg.
00441       for ( kk = 0; kk < nearest_neighbors.size(); kk++ )
00442       {
00443         tmp_seg = nearest_neighbors.get();
00444         if ( tmp_seg == curr_seg ||
00445              tmp_seg == curr_seg->get_prev() ||
00446              tmp_seg == curr_seg->get_next() )
00447           nearest_neighbors.remove();
00448       }
00449       if ( nearest_neighbors.size() == 0 )
00450         continue;
00451         //The min distance is the distance between the curr_loc
00452         //and the top nearest_neighbor.
00453       nearest_neighbors.reset();
00454       tmp_seg = nearest_neighbors.get();
00455       ent_1 = tmp_seg->owner();
00456       ent_2 = curr_seg->get_start()->owner();
00457       top_1 = CAST_TO(ent_1, TopologyEntity);
00458       top_2 = CAST_TO(ent_2, TopologyEntity);
00459       
00460       if ( top_1->is_directly_related(top_2) )
00461         continue;
00462       
00463         //we'll need to recompute it...
00464       dist = dist_sq_point_data(curr_loc, tmp_seg);
00465       if ( dist < min_distance_in_one_loop )
00466         min_distance_in_one_loop = dist;
00467     }
00468   }
00469   if ( min_distance_in_one_loop == CUBIT_DBL_MAX )
00470   {
00471     ent_1 = (RefEntity*)face;
00472     DLIList<RefEntity*> entities;
00473     entities.append(ent_1);
00474     double smallest, largest, ave, sum;
00475     RefEdge *small_e, *large_e;
00476     measure_curve_length(entities, smallest, small_e,
00477                          largest, large_e, ave, sum);
00478     min_distance_in_one_loop = smallest;
00479   }
00480   else
00481     min_distance_in_one_loop = sqrt(min_distance_in_one_loop);
00482     //clean up the memory.
00483   int ll,mm;
00484   atree_list.reset();
00485   boundary_seg_loops.reset();
00486   for ( ll = boundary_seg_loops.size(); ll > 0; ll-- )
00487   {
00488     delete atree_list.pop();
00489     seg_list = boundary_seg_loops.pop();
00490     for ( mm = seg_list->size(); mm > 0; mm-- )
00491       delete seg_list->pop();
00492     delete seg_list;
00493   }
00494         
00495   PointList *point_list;
00496   for ( ll = boundary_point_loops.size(); ll > 0; ll-- )
00497   {
00498     point_list = boundary_point_loops.pop();
00499     for ( mm = point_list->size(); mm > 0; mm-- )
00500     {
00501       GeomPoint *tmp_point = point_list->pop();
00502         //delete the CubitPointData
00503       delete tmp_point;
00504     }
00505     delete point_list;
00506   }
00507   return CUBIT_SUCCESS;
00508 }
00509 //-------------------------------------------------------------------
00510 //Function: dist_point_data
00511 //Computes the distance between the line segment and the vector position.
00512 //This function is used to determine the nearest neigbhor search and
00513 //gets passed to the AbstractTree.
00514 //-------------------------------------------------------------------
00515 double GeomMeasureTool::dist_sq_point_data( CubitVector &curr_point,
00516                                             GeomSeg *&curr_seg )
00517 {
00518   double node[3];
00519   double point_0[3], point_1[3];
00520   node[0] = curr_point.x();
00521   node[1] = curr_point.y();
00522   node[2] = curr_point.z();
00523   CubitVector start_v = curr_seg->get_start()->coordinates();
00524   CubitVector end_v = curr_seg->get_end()->coordinates();
00525   point_0[0] = start_v.x();
00526   point_0[1] = start_v.y();
00527   point_0[2] = start_v.z();
00528   point_1[0] = end_v.x();
00529   point_1[1] = end_v.y();
00530   point_1[2] = end_v.z();
00531 
00532   IntersectionTool new_tool(GEOMETRY_RESABS);
00533   double t;
00534   double dist = new_tool.distance_point_line(node, point_0,
00535                                              point_1, t);
00536   if ( dist < 0.0 )
00537   {
00538       //This happens if the point is beyond the line segment, in which
00539       //case we just want the closer end point.
00540     if ( t < 0.0 )
00541     {
00542       dist = (curr_point -start_v).length_squared();
00543     }
00544     else
00545     {
00546       dist = (curr_point - end_v).length_squared();
00547     }
00548   }
00549   else
00550   {
00551       //I hate to do this but everything else is in
00552       //terms of distance squared.
00553     dist = dist*dist;
00554   }
00555   return dist;
00556 }
00557 
00558 CubitStatus GeomMeasureTool::convert_to_lines( PointLoopList &boundary_point_loops,
00559                                                SegLoopList &boundary_line_loops,
00560                                                RefFace *ref_face )
00561 {
00562     //This does more than just converts the boundary points to line segments.  It also
00563     //Sets up a doubily linked list through the ImprintLineSegment data structure.
00564     //Calling codes need to make sure all this data is deleted.
00565     //This function also sets the PointType for each of the boundary points.  This
00566     //will depend also on the boundary_1 flag.  If we are doing boundary_1 then,
00567     //the flag will be true, otherwise it will be false.
00568   int ii, jj;
00569   PointList *point_list;
00570   SegList *segment_list;
00571   GeomSeg *new_line, *prev = NULL;
00572   SegList allocated_segs;
00573   RefEntity *entity_start, *entity_end, *seg_owner;
00574   boundary_point_loops.reset();
00575   for ( ii = boundary_point_loops.size(); ii > 0; ii-- )
00576   {
00577     point_list = boundary_point_loops.get_and_step();
00578     segment_list = new SegList;
00579     prev = NULL;
00580     new_line = NULL;
00581     seg_owner = NULL;
00582     for ( jj = point_list->size(); jj > 0; jj-- )
00583     {
00584       entity_start = point_list->get()->owner();
00585       entity_end = point_list->next()->owner();
00586       seg_owner =  CAST_TO(entity_start, RefEdge);
00587       if ( seg_owner == NULL )
00588       {
00589         seg_owner = CAST_TO(entity_end, RefEdge);
00590         if ( seg_owner == NULL )
00591         {
00592           RefVertex *ref_vert_end = CAST_TO(entity_end, RefVertex);
00593           RefVertex *ref_vert_start = CAST_TO(entity_start, RefVertex);
00594           DLIList <RefEdge*> common_edges;
00595           ref_vert_start->common_ref_edges(ref_vert_end,
00596                                            common_edges,
00597                                            ref_face);
00598           if ( common_edges.size() == 0  )
00599           {
00600             PRINT_ERROR("GeomMeasureTool::convert_to_lines having problem finding owner \n" 
00601                         "       of boundary.  Could be indicative of corrupt geometry.\n");
00602             int ll;
00603             for ( ll = boundary_line_loops.size(); ll > 0; ll-- )
00604               delete boundary_line_loops.pop();
00605             for ( ll = allocated_segs.size(); ll > 0; ll-- )
00606               delete allocated_segs.pop();
00607             delete segment_list;
00608             return CUBIT_FAILURE;
00609           }
00610           if ( common_edges.size() == 1 )
00611             seg_owner = common_edges.get();
00612           else
00613           {
00614               //Now we have to decide which of these edges is the one.
00615               //Lets take a mid point on this segment and check the distances.
00616               //First find the mid_point of the line segment.
00617             CubitVector start = point_list->get()->coordinates();
00618             CubitVector end = point_list->next()->coordinates();
00619             CubitVector mid = (start+end)/2.0;
00620             CubitVector closest_point;
00621             
00622               //Now test the edges to see which is closest.
00623             RefEdge *closest_edge = NULL, *curr_ref_edge;
00624             double min_dist = CUBIT_DBL_MAX, dist;
00625             int kk;
00626             for ( kk = common_edges.size(); kk > 0; kk-- )
00627             {
00628               curr_ref_edge = common_edges.get_and_step();
00629               curr_ref_edge->closest_point(mid, closest_point);
00630               dist = (mid - closest_point).length_squared();
00631               if ( dist < min_dist )
00632               {
00633                 min_dist = dist;
00634                 closest_edge = curr_ref_edge;
00635               }
00636             }
00637             if ( closest_edge == NULL )
00638             {
00639               PRINT_ERROR("Problems determining segment ownwership.\n"
00640                           "Internal problem with virtual imprinting.\n");
00641               int ll;
00642               for ( ll = boundary_line_loops.size(); ll > 0; ll-- )
00643                 delete boundary_line_loops.pop();
00644               for ( ll = allocated_segs.size(); ll > 0; ll-- )
00645                 delete allocated_segs.pop();
00646               delete segment_list;
00647               return CUBIT_FAILURE;
00648             }
00649             seg_owner = closest_edge;
00650           }
00651         }
00652       }
00653       new_line = new GeomSeg( point_list->get(), point_list->next(), seg_owner);
00654       allocated_segs.append(new_line);
00655       point_list->step();
00656       segment_list->append(new_line);
00657       if ( prev != NULL )
00658       {
00659         prev->set_next(new_line);
00660         new_line->set_prev(prev);
00661       }
00662       prev = new_line;
00663     }
00664     if ( prev != NULL )
00665     {
00666       segment_list->reset();
00667       assert(prev == segment_list->prev());
00668       segment_list->prev()->set_next(segment_list->get());
00669       segment_list->get()->set_prev(segment_list->prev());
00670       assert(segment_list->get()->get_next() != NULL );
00671     }
00672     if ( segment_list->size() )
00673     {
00674       boundary_line_loops.append(segment_list);
00675     }
00676     else
00677       delete segment_list;
00678   }
00679   return CUBIT_SUCCESS;
00680 }
00681 
00682 CubitStatus GeomMeasureTool::get_boundary_points( RefFace *ref_face,
00683                                                   PointLoopList &boundary_point_loops,
00684                                                   double seg_length_tol)
00685 {
00686   DLIList<DLIList<CoEdge*> > co_edge_loops;
00687   ref_face->co_edge_loops(co_edge_loops);
00688   int ii, jj;
00689   PointList *new_point_loop_ptr, tmp_point_list;
00690   RefEdge *ref_edge_ptr;
00691   CoEdge *co_edge_ptr;
00692   CubitStatus stat;
00693   CubitSense sense;
00694   
00695   for ( ii = co_edge_loops.size(); ii > 0; ii-- )
00696   {
00697     DLIList <CoEdge*> &co_edge_list_ptr = co_edge_loops.get_and_step();
00698     new_point_loop_ptr = new PointList;
00699     for ( jj = co_edge_list_ptr.size(); jj > 0; jj-- )
00700     {
00701       co_edge_ptr = co_edge_list_ptr.get_and_step();
00702       ref_edge_ptr = co_edge_ptr->get_ref_edge_ptr();
00703       tmp_point_list.clean_out();
00704       stat = get_curve_facets( ref_edge_ptr, tmp_point_list,
00705                                seg_length_tol );
00706       if ( stat != CUBIT_SUCCESS )
00707         return CUBIT_FAILURE;
00708       if ( tmp_point_list.size() == 0 )
00709         continue;
00710       tmp_point_list.reset();
00711         //the points are in order from start vertex to end vertex.
00712         //append them now according to the loop.
00713       sense = co_edge_ptr->get_sense();
00714       if ( CUBIT_FORWARD != sense )
00715         tmp_point_list.reverse();
00716         //Now take off the last point as it is a duplicate with the
00717         //other list...
00718       tmp_point_list.reset();
00719       delete tmp_point_list.pop();
00720       (*new_point_loop_ptr) += tmp_point_list;
00721     }
00722     boundary_point_loops.append(new_point_loop_ptr);
00723   }
00724       
00725   return CUBIT_SUCCESS;
00726 }
00727 
00728 CubitStatus GeomMeasureTool::get_curve_facets( RefEdge* curve,
00729                                                PointList &segments,
00730                                                double seg_length_tol ) 
00731 {
00732 //  const double COS_ANGLE_TOL =  0.965925826289068312213715; // cos(15)
00733 //  const double COS_ANGLE_TOL =  0.984807753012208020315654; // cos(10)
00734   //const double COS_ANGLE_TOL =  0.996194698091745545198705; // cos(5)
00735   GMem curve_graphics;
00736     //make this tol bigger than seg_length_tol, to
00737     //make sure the segments are larger than the tolerance.
00738   const double dist_tol = seg_length_tol;// + .05*seg_length_tol;
00739   //const double dist_tol_sqr = dist_tol*dist_tol;
00740   if ( curve->get_curve_ptr() == NULL )
00741   {
00742     PRINT_ERROR("Curve %d had no geometry, this is a bug!\n", curve->id());
00743     return CUBIT_FAILURE;
00744   }
00745   else if ( curve->get_curve_ptr()->get_geometry_query_engine() == NULL )
00746   {
00747     PRINT_ERROR("Curve %d had no geometry engine, this is a bug!\n",
00748                 curve->id());
00749     return CUBIT_FAILURE;
00750   }
00751   curve->get_graphics( curve_graphics );
00752   
00753   GPoint* gp = curve_graphics.point_list();
00754   if ( gp == NULL )
00755   {
00756     if ( curve->measure() < GEOMETRY_RESABS ) 
00757     {
00758       if( curve->geometry_type() != POINT_CURVE_TYPE )
00759       {
00760         PRINT_INFO("Curve %d has a zero length, and won't be used for loop measurements.\n",
00761                     curve->id());
00762       }
00763     }
00764     else
00765       PRINT_ERROR("Curve %d had not faceting.  Ignoring this curve.\n",
00766                   curve->id());
00767     return CUBIT_SUCCESS;
00768   }
00769   GeomPoint *last = new GeomPoint( gp[0].x, gp[0].y, gp[0].z, curve );
00770   CubitVector lastv = last->coordinates();
00771   int num_points = curve_graphics.pointListCount;
00772   segments.append( last );
00773   int ii;
00774   CubitBoolean remove_second_to_end = CUBIT_FALSE;
00775   for ( ii = 1; ii < num_points; ii++ )
00776   {
00777     CubitVector pos(  gp[ii].x, gp[ii].y, gp[ii].z );
00778     CubitVector step1 = (pos - lastv);
00779     double len1 = step1.length();
00780     if( len1 < dist_tol && ii != num_points - 1) 
00781       continue;
00782     else if ( len1 < dist_tol && ii == num_points-1 )
00783     {
00784       remove_second_to_end = CUBIT_TRUE;
00785     }
00786     last = new GeomPoint(pos,curve);
00787     segments.append( last );
00788     lastv = last->coordinates();
00789   }
00790     // Now check if the segment list is reversed wrt the curve direction.
00791   segments.reset();
00792   if ( remove_second_to_end )
00793   {
00794     if ( segments.size() == 2 )
00795     {
00796         //just leave it in, let it be small...
00797     }
00798     else
00799     {
00800         //Remove the second to last one.  To do
00801         //this efficiently (don't do remove), pop
00802         //the last one, then save that and
00803         //re-add it after poping the second one.
00804       GeomPoint *temp = segments.pop();
00805       segments.pop();
00806       segments.append(temp);
00807     }
00808   }
00809   segments.reset();
00810   if( curve->start_vertex() != curve->end_vertex() )
00811   {
00812     CubitVector start_vec, end_vec;
00813     start_vec = curve->start_vertex()->coordinates();
00814     end_vec = curve->end_vertex()->coordinates();
00815     CubitVector start_seg = segments.get()->coordinates();
00816     double dist_1 = (start_seg - start_vec).length_squared();
00817     double dist_2 = (start_seg - end_vec).length_squared();
00818     if ( dist_1 > dist_2 )
00819       segments.reverse();
00820   }
00821   else
00822   {
00823     double u1, u2;
00824     u1 = curve->u_from_position( (segments.next(1)->coordinates()) );
00825     u2 = curve->u_from_position( (segments.next(2)->coordinates()) );    
00826     if( (u2 < u1) && (curve->start_param() <= curve->end_param()) )
00827       segments.reverse();
00828   }
00829     //clean up the periodic curve case (last seg may be too small.)
00830   if ( curve->start_vertex() == curve->end_vertex() )
00831   {
00832     segments.reset();
00833     CubitVector start_v = segments.get()->coordinates();
00834     CubitVector last_v = segments.prev()->coordinates();
00835     double dist = (start_v - last_v).length();
00836     if ( dist < dist_tol )
00837     {
00838         //remove the last one.
00839       segments.pop();
00840     }
00841   }
00842   if ( segments.size() < 2 )
00843   {
00844     PRINT_ERROR("Problem faceting boundary.\n");
00845     return CUBIT_FAILURE;
00846   }
00847   segments.get()->owner(curve->start_vertex());
00848   segments.prev()->owner(curve->end_vertex());
00849   return CUBIT_SUCCESS;
00850 }
00851 CubitStatus GeomMeasureTool::interior_angle(GeomSeg *first_seg,
00852                                             GeomSeg *next_seg,
00853                                             RefFace *face,
00854                                             double &angle )
00855 {
00856   if ( first_seg->get_end() != next_seg->get_start() )
00857   {
00858     PRINT_ERROR("GeomMeasureTool::interior_angle expects segments\n"
00859                 "to be connected (Users: this is a bug, please submit.\n");
00860     return CUBIT_FAILURE;
00861   }
00862   CubitVector point = first_seg->get_end()->coordinates();
00863   CubitVector to_prev = first_seg->get_start()->coordinates();
00864   to_prev -= point;
00865   CubitVector to_next = next_seg->get_end()->coordinates();
00866   to_next -= point;
00867   CubitVector surf_norm = face->normal_at(point);
00868   if ( surf_norm.length_squared() < CUBIT_RESABS )
00869   {
00870       //Try getting it at one of the other nodes...
00871     surf_norm = face->normal_at(first_seg->get_start()->coordinates() );
00872     if (surf_norm.length_squared() < CUBIT_RESABS )
00873     {
00874         //Try getting it at one of the other nodes...
00875       surf_norm = face->normal_at(next_seg->get_end()->coordinates() );
00876       if (surf_norm.length_squared() < CUBIT_RESABS )
00877       {
00878         PRINT_ERROR("Trying to get normal at surf %d.\n"
00879                     "       Normal length being returned equals zero.\n",
00880                     face->id() );
00881         angle = CUBIT_DBL_MAX;
00882         return CUBIT_FAILURE;
00883       }
00884     }
00885   }
00886   angle = surf_norm.vector_angle ( to_next, to_prev );
00887   return CUBIT_SUCCESS;
00888 }
00889 
00890 void GeomMeasureTool::measure_face_area (DLIList <RefEntity*> &entity_list,
00891                                          double &smallest,
00892                                          RefFace *&smallest_face,
00893                                          double &largest,
00894                                          RefFace *&largest_face,
00895                                          double &average,
00896                                          double &sum)
00897 {
00898   DLIList<RefFace*> ref_faces;
00899   get_faces_from_list(entity_list, ref_faces);
00900   measure_face_area(ref_faces, smallest, smallest_face,
00901                     largest, largest_face, average, sum);
00902   
00903 }
00904 
00905 void GeomMeasureTool::measure_face_area (DLIList <RefFace*> &ref_faces,
00906                                          double &smallest,
00907                                          RefFace *&smallest_face,
00908                                          double &largest,
00909                                          RefFace *&largest_face,
00910                                          double &average,
00911                                          double &sum)
00912 {
00913   int ii;
00914   double curr_area, total_area = 0.0;
00915   smallest = CUBIT_DBL_MAX;
00916   largest = 0.0;
00917   RefFace *curr_face;
00918   for( ii = ref_faces.size(); ii > 0; ii-- )
00919   {
00920     curr_face = ref_faces.get_and_step();
00921     curr_area = measure_area(curr_face);
00922     if ( curr_area < smallest )
00923     {
00924       smallest = curr_area;
00925       smallest_face = curr_face;
00926     }
00927     if( curr_area > largest )
00928     {
00929       largest = curr_area;
00930       largest_face = curr_face;
00931     }
00932     total_area += curr_area;
00933   }
00934   if (ref_faces.size() != 0 )
00935     average = total_area/ref_faces.size();
00936   sum = total_area;
00937   
00938 }
00939 
00940 void GeomMeasureTool::measure_volume_volume (DLIList <RefEntity*> &entity_list,
00941                                              double &smallest,
00942                                              RefVolume *&smallest_volume,
00943                                              double &largest,
00944                                              RefVolume *&largest_volume,
00945                                              double &average,
00946                                              double &sum)
00947 {
00948   DLIList<RefVolume*> ref_volumes;
00949   get_volumes_from_list(entity_list, ref_volumes);
00950   measure_volume_volume(ref_volumes, smallest, smallest_volume,
00951                         largest, largest_volume, average, sum);
00952   
00953 }
00954 
00955 void GeomMeasureTool::measure_volume_volume (DLIList <RefVolume*> &ref_volumes,
00956                                              double &smallest,
00957                                              RefVolume *&smallest_volume,
00958                                              double &largest,
00959                                              RefVolume *&largest_volume,
00960                                              double &average,
00961                                              double &sum)
00962 {
00963   int ii;
00964   double curr_volume, total_volume = 0.0;
00965   smallest = CUBIT_DBL_MAX;
00966   largest = 0.0;
00967   RefVolume *curr_solid;
00968   for ( ii = ref_volumes.size(); ii > 0; ii-- )
00969   {
00970     curr_solid = ref_volumes.get_and_step();
00971     curr_volume = curr_solid->measure();
00972 
00973     if ( curr_volume < smallest )
00974     {
00975       smallest = curr_volume;
00976       smallest_volume = curr_solid;
00977     }
00978     if ( curr_volume > largest )
00979     {
00980       largest = curr_volume;
00981       largest_volume = curr_solid;
00982     }
00983     total_volume += curr_volume;
00984   }
00985   if (ref_volumes.size() != 0 )
00986     average = total_volume/ref_volumes.size();
00987   sum = total_volume;
00988   
00989 }
00990 
00991 void GeomMeasureTool::measure_face_area_and_hydraulic_radius(RefFace *curr_face,
00992                                                              double &face_area,
00993                                                              double &face_hydraulic_radius)
00994 {
00995   int ii;
00996   double total_edge_length = 0.0, curr_edge_length;
00997   RefEdge *curr_edge;
00998   DLIList <CoEdge*> my_co_edges;
00999 
01000   face_area = measure_area(curr_face);
01001   curr_face->co_edges(my_co_edges);
01002   for( ii = my_co_edges.size(); ii > 0; ii-- )
01003   {
01004     curr_edge = my_co_edges.get_and_step()->get_ref_edge_ptr();
01005     curr_edge_length = curr_edge->measure();
01006     total_edge_length += curr_edge_length;
01007   }
01008   if (total_edge_length != 0)
01009     face_hydraulic_radius = 4.0*(face_area / total_edge_length);
01010 }
01011 
01012 void GeomMeasureTool::measure_volume_volume_and_hydraulic_radius(RefVolume *curr_volume,
01013                                                                  double &volume_volume,
01014                                                                  double &volume_hydraulic_radius)
01015 {
01016   int ii;
01017   double curr_face_area, total_surface_area = 0.0;
01018   DLIList <RefFace*> ref_face_list;
01019   RefFace *curr_face;
01020 
01021   volume_volume = curr_volume->measure();
01022   curr_volume->ref_faces(ref_face_list);
01023   
01024   for ( ii = ref_face_list.size(); ii > 0; ii-- )
01025   {
01026     curr_face = ref_face_list.get_and_step();
01027     curr_face_area = measure_area(curr_face);
01028     total_surface_area += curr_face_area;
01029   }
01030   if (total_surface_area != 0)
01031     volume_hydraulic_radius = 6.0*(volume_volume / total_surface_area);
01032  
01033 }
01034 
01035 void GeomMeasureTool::angles_between_volume_surfaces(RefVolume *curr_volume,
01036                                                      double &min_angle,
01037                                                      double &max_angle,
01038                                                      RefFace *&face_min_1,
01039                                                      RefFace *&face_min_2,
01040                                                      RefFace *&face_max_1,
01041                                                      RefFace *&face_max_2)
01042 {
01043   int ii, jj, shared_angles = 0;
01044   double common_angle, sum_of_angles = 0.0;
01045   double smallest = CUBIT_DBL_MAX, largest = 0.0;
01046   DLIList <RefFace*> face_list_1, face_list_2;
01047   RefFace *curr_face_1, *curr_face_2;
01048   RefEdge *common_edge;
01049 
01050   curr_volume->ref_faces(face_list_1);
01051   face_list_2 = face_list_1;
01052 
01053   for (ii = face_list_1.size(); ii > 0; ii--)
01054   {
01055     curr_face_1 = face_list_1.get_and_step();
01056     for(jj = face_list_2.size(); jj > 0; jj--)
01057     {
01058       curr_face_2 = face_list_2.get_and_step();
01059 
01060       if(curr_face_1 == curr_face_2)
01061         continue;
01062       else
01063       {
01064         common_edge = curr_face_1->common_ref_edge(curr_face_2);
01065         if( common_edge == NULL )
01066           continue;
01067         else
01068         {
01069           shared_angles++;
01070           
01071           if ( common_edge->get_curve_ptr()->geometry_type() != STRAIGHT_CURVE_TYPE )
01072           {
01073             int kk;
01074             double frac_pos = 0.00;
01075             for ( kk = 0; kk < 4; kk++ )
01076             {
01077               common_angle =  GeometryQueryTool::instance()->surface_angle(curr_face_1, curr_face_2,
01078                                                                            common_edge, curr_volume,
01079                                                                            frac_pos);
01080               frac_pos += .25;
01081               if( common_angle < smallest )
01082               {
01083                 min_angle = common_angle * (180 / CUBIT_PI);
01084                 smallest = common_angle;
01085                 face_min_1 = curr_face_1;
01086                 face_min_2 = curr_face_2;
01087               }
01088               if( common_angle > largest )
01089               {
01090                 max_angle = common_angle * (180 / CUBIT_PI);
01091                 largest = common_angle;
01092                 face_max_1 = curr_face_1;
01093                 face_max_2 = curr_face_2;
01094               }
01095             }
01096           }
01097           else
01098           {
01099             
01100             common_angle =  GeometryQueryTool::instance()->surface_angle(curr_face_1, curr_face_2,
01101                                                                          common_edge, curr_volume );
01102             if( common_angle < smallest )
01103             {
01104               min_angle = common_angle * (180 / CUBIT_PI);
01105               smallest = common_angle;
01106               face_min_1 = curr_face_1;
01107               face_min_2 = curr_face_2;
01108             }
01109             if( common_angle > largest )
01110             {
01111               max_angle = common_angle * (180 / CUBIT_PI);
01112               largest = common_angle;
01113               face_max_1 = curr_face_1;
01114               face_max_2 = curr_face_2;
01115             }
01116           }
01117           
01118           sum_of_angles += common_angle;
01119         }
01120       }
01121     }
01122   }
01123 }
01124 
01125 void GeomMeasureTool::merged_unmerged_surface_ratio(DLIList <RefVolume*> &ref_volumes,
01126                                                     int &merged, int &unmerged,
01127                                                     double &ratio)
01128   
01129 {
01130   int ii, jj;
01131   RefVolume *curr_volume_1;
01132   DLIList <RefFace*> surface_list_1, surface_list_2;
01133   RefFace *curr_face_1;
01134   int merged_surfaces = 0, unmerged_surfaces = 0;
01135   
01136   for( ii = ref_volumes.size(); ii > 0; ii-- )
01137   {
01138     curr_volume_1 = ref_volumes.get_and_step();
01139     surface_list_1.clean_out();
01140     
01141     curr_volume_1->ref_faces(surface_list_1);
01142     for( jj = surface_list_1.size(); jj > 0; jj-- )
01143     {
01144       curr_face_1 = surface_list_1.get_and_step();
01145       if ( curr_face_1->marked() )
01146         continue;
01147       else
01148       {
01149         curr_face_1->marked(1);
01150         surface_list_2.append(curr_face_1);
01151       }
01152       if ( curr_face_1->num_ref_volumes() == 2 )
01153         merged_surfaces++;
01154       else
01155         unmerged_surfaces++;
01156     }
01157   }
01158   for ( ii = surface_list_2.size(); ii > 0; ii-- )
01159   {
01160     surface_list_2.get_and_step()->marked(0);
01161   }
01162   merged = merged_surfaces;
01163   unmerged = unmerged_surfaces;
01164   if ( unmerged_surfaces == 0 )
01165   {
01166     ratio = 100.0;
01167   }
01168   else
01169     ratio = (double) ((double)merged_surfaces / (double)unmerged_surfaces);
01170 }
01171 
01172 void GeomMeasureTool::report_intersected_volumes(DLIList <RefVolume*> &ref_vols,
01173                                                  DLIList <RefVolume*> &intersection_list)
01174 {
01175   DLIList <RefVolume*> results;
01176   RefVolume *curr_vol;
01177   int i, j;
01178   ProgressTool *progress_ptr = NULL;
01179   int total_volumes = ref_vols.size();
01180   if (total_volumes > 5)
01181   {
01182     progress_ptr = AppUtil::instance()->progress_tool();
01183     assert(progress_ptr != NULL);
01184     progress_ptr->start(0, 100, "Overlapping Volumes Progress", 
01185       NULL, CUBIT_TRUE, CUBIT_TRUE);
01186   }
01187   double curr_percent = 0.0;
01188 
01189   for( i = 0; i < ref_vols.size(); i++ )
01190   {
01191     curr_vol = ref_vols.next(i);
01192 
01193     //is this body a multi-volume-body?
01194     Body *single_volume_body = curr_vol->body();
01195     DLIList<RefVolume*> body_vols;
01196     single_volume_body->ref_volumes( body_vols );
01197     if( body_vols.size() > 1 )
01198       single_volume_body = NULL;
01199       
01200     for( j = (i + 1); j < ref_vols.size(); j++ )
01201     {
01202       RefVolume *curr_vol2 = ref_vols.next(j);
01203       if ( AppUtil::instance()->interrupt() )
01204       {
01205           //interrpt.  We need to exit.
01206         if ( progress_ptr != NULL )
01207         {
01208           progress_ptr->end();
01209         }
01210           //just leave what has been calculated...
01211         return;
01212       }
01213 
01214       Body *single_volume_body2 = curr_vol2->body();
01215       DLIList<RefVolume*> body_vols2;
01216       single_volume_body2->ref_volumes( body_vols2 );
01217       if( body_vols2.size() > 1 )
01218         single_volume_body2 = NULL;
01219 
01220       //update the progress..
01221       if ( progress_ptr != NULL )
01222       {
01223         curr_percent = ((double)(i+1))/((double)(total_volumes));
01224         progress_ptr->percent(curr_percent);
01225       }
01226     
01227       //if both are single-volume-bodies
01228       if( single_volume_body && single_volume_body2 )
01229       {
01230         if( GeometryQueryTool::instance()->bodies_overlap( single_volume_body,
01231                                                        single_volume_body2 ) )
01232         {
01233           intersection_list.append( curr_vol );
01234           intersection_list.append( curr_vol2 );
01235         }
01236       }
01237       else if( GeometryQueryTool::instance()->volumes_overlap( curr_vol, curr_vol2 ) )
01238       {
01239        intersection_list.append( curr_vol );
01240        intersection_list.append( curr_vol2 );
01241       }
01242     }
01243   }
01244   
01245   if ( progress_ptr != NULL )
01246   {
01247     progress_ptr->end();
01248   }
01249 }
01250 
01251 void GeomMeasureTool::report_intersected_bodies(DLIList <Body*> &ref_bodies,
01252                                                 DLIList <Body*> &intersection_list)
01253 {
01254   Body *curr_body, *curr_body_2;
01255   int ii, jj;
01256   ProgressTool *progress_ptr = NULL;
01257   int total_bodies = ref_bodies.size();
01258   if (total_bodies > 5)
01259   {
01260     progress_ptr = AppUtil::instance()->progress_tool();
01261     assert(progress_ptr != NULL);
01262     progress_ptr->start(0, 100, "Overlapping Volumes Progress", 
01263       NULL, CUBIT_TRUE, CUBIT_TRUE);
01264   }
01265   double curr_percent = 0.0;
01266  
01267   
01268   for( ii = 0; ii < ref_bodies.size(); ii++ )
01269   {
01270     curr_body = ref_bodies.next(ii);
01271 
01272     for( jj = (ii + 1); jj < ref_bodies.size(); jj++ )
01273     {
01274       curr_body_2 = ref_bodies.next(jj);
01275       if ( AppUtil::instance()->interrupt() )
01276       {
01277           //interrpt.  We need to exit.
01278         if ( progress_ptr != NULL )
01279         {
01280           progress_ptr->end();
01281         }
01282           //just leave what has been calculated...
01283         return;
01284       }
01285 
01286       if (GeometryQueryTool::instance()->bodies_overlap(curr_body,
01287                                                         curr_body_2) )
01288       {
01289         intersection_list.append(curr_body);
01290         intersection_list.append(curr_body_2);
01291       }
01292     }
01293        //update the progress..
01294     if ( progress_ptr != NULL )
01295     {
01296       curr_percent = ((double)(ii+1))/((double)(total_bodies));
01297       progress_ptr->percent(curr_percent);
01298     }
01299   }
01300   if ( progress_ptr != NULL )
01301   {
01302     progress_ptr->end();
01303   }
01304 }
01305 
01306 void GeomMeasureTool::face_list_from_volume_list( DLIList <RefVolume*> &input_vols,
01307                                                   DLIList <RefFace*> &all_faces )
01308 {
01309   DLIList <RefFace*> curr_faces;
01310   RefVolume *curr_vol;
01311   RefFace *curr_face;
01312   int ii, jj;
01313 
01314     //get the all RefFaces from input_vols.
01315     //put them into all_faces list.
01316   
01317   for ( ii = 0; ii < input_vols.size(); ii++ )
01318   {
01319     curr_vol = input_vols.get_and_step();
01320     curr_faces.clean_out();
01321     curr_vol->ref_faces(curr_faces);
01322     
01323     for( jj = 0; jj < curr_faces.size(); jj++ )
01324     {
01325       curr_face = curr_faces.get_and_step();
01326       if( curr_face->marked() )
01327         continue;
01328       else
01329       {
01330         curr_face->marked(1);
01331         all_faces.append(curr_face);
01332       }
01333     }
01334   }
01335   for ( ii = all_faces.size();ii >0; ii-- )
01336     all_faces.get_and_step()->marked(0);
01337 }
01338 
01339 RefFace* GeomMeasureTool::valid_start( DLIList <RefFace*> &all_faces )
01340 {
01341   int ii;
01342   RefFace *start_face;
01343 
01344   for( ii = 0; ii < all_faces.size(); ii++ )
01345   {
01346     start_face = all_faces.get_and_step();
01347 
01348     if( !start_face->marked() && start_face->num_ref_volumes()== 1 )
01349       return start_face;
01350   }
01351 
01352   return NULL;
01353   
01354 }
01355 
01356 void GeomMeasureTool::find_shells( DLIList <RefVolume*> &input_vols,
01357                                    RefGroup *&owner_groups,
01358                                    int &number_of_shells)
01359 {
01360   DLIList <DLIList<RefEntity*>*> list_of_lists;
01361   DLIList <RefEntity*> *new_list;
01362   DLIList <RefFace*> all_faces, stack, related_faces;
01363   DLIList <RefEdge*> curr_edges;
01364   RefFace *start_face, *curr_face, *related_face;
01365   RefEdge *curr_edge;
01366   int ii, jj;
01367   RefEntity *tmp_ent;
01368   number_of_shells = 0;
01369   DLIList <RefFace*> marked_faces;
01370   
01371   face_list_from_volume_list( input_vols, all_faces );
01372   
01373     // all_faces all marked - need unmarked
01374   for (ii = 0; ii < all_faces.size(); ii++ )
01375   {
01376     curr_face = all_faces.get_and_step();
01377     curr_face->marked(0);
01378   }
01379   while ( (start_face = valid_start(all_faces) ) != NULL )
01380   {
01381       //create a list.
01382     stack.clean_out();
01383     new_list = new DLIList<RefEntity*>;
01384     stack.append(start_face);
01385     start_face->marked(1);
01386     marked_faces.append(start_face);
01387 
01388     while (stack.size() > 0 )
01389     {
01390       curr_edges.clean_out();
01391 
01392       curr_face = stack.pop(); // get last face added
01393       tmp_ent = CAST_TO(curr_face, RefEntity);
01394       new_list->append(tmp_ent); // append that face to the newest list
01395       curr_face->ref_edges(curr_edges); // gets list of face's edges
01396 
01397       for( ii = 0; ii < curr_edges.size(); ii++ )  // loops edges
01398       {
01399         related_faces.clean_out();
01400 
01401         curr_edge = curr_edges.get_and_step(); // chooses next edge
01402         curr_edge->ref_faces(related_faces); // finds edge's related faces
01403         for( jj = 0; jj < related_faces.size(); jj++ ) // loops faces
01404         {
01405           related_face = related_faces.get_and_step(); // gets face
01406           if(!related_face->marked() && related_face->num_ref_volumes() == 1)
01407           {
01408             stack.append(related_face);
01409             related_face->marked(1);
01410             marked_faces.append(related_face);
01411           }
01412         }
01413       }
01414     }
01415     list_of_lists.append(new_list);
01416     number_of_shells++;
01417   }
01418 
01419   for (ii = 0; ii < marked_faces.size(); ii++ )
01420   {
01421     curr_face = marked_faces.get_and_step();
01422     curr_face->marked(0);
01423   }
01424   
01425   RefGroup *curr_group;
01426   owner_groups = RefEntityFactory::instance()->construct_RefGroup("volume_shells");
01427   for( ii = list_of_lists.size(); ii > 0; ii-- )
01428   {
01429     new_list = list_of_lists.pop();
01430     curr_group = RefEntityFactory::instance()->construct_RefGroup(*new_list);
01431     tmp_ent = CAST_TO(curr_group, RefEntity);
01432     owner_groups->add_ref_entity(tmp_ent);
01433     delete new_list;
01434   } 
01435 }
01436 
01437 void GeomMeasureTool::print_surface_measure_summary( DLIList <RefFace*> &ref_faces )
01438                                                   
01439 {
01440   
01441   double min_dist_loops=CUBIT_DBL_MAX, min_dist_loop=CUBIT_DBL_MAX;
01442   double curr_min_dist_loops=CUBIT_DBL_MAX, curr_min_dist_loop=CUBIT_DBL_MAX;
01443   double min_angle=CUBIT_DBL_MAX, max_angle=-CUBIT_DBL_MAX;
01444   double curr_min_angle=CUBIT_DBL_MAX, curr_max_angle=-CUBIT_DBL_MAX;
01445   double curr_face_area, curr_face_hydraulic_radius;
01446   double min_face_area=CUBIT_DBL_MAX, min_face_hydraulic_radius=CUBIT_DBL_MAX;
01447   double max_face_area=-CUBIT_DBL_MAX;
01448   int counter = 0;
01449   double total_area= 0.0;
01450   double min_curve_length = CUBIT_DBL_MAX;
01451   double max_curve_ratio = -CUBIT_DBL_MAX;
01452   double curve_length, curve_ratio;
01453   int ii;
01454        
01455   RefFace *curr_face;
01456   RefFace *face_min_dist_loops=NULL, *face_min_dist_loop=NULL;
01457   RefFace *face_min_angle=NULL, *face_max_angle=NULL;
01458   RefFace *face_min_area=NULL, *face_min_hydraulic_radius=NULL;
01459   RefFace *face_max_area=NULL;
01460   RefFace *face_min_curve_length=NULL, *face_max_curve_ratio=NULL;
01461   RefEdge *smallest_edge=NULL, *curr_small_edge=NULL;
01462        
01463        
01464   for( ii = 0; ii < ref_faces.size(); ii++ )
01465   {
01466     curr_face = ref_faces.get_and_step();
01467     if ( curr_face->num_ref_edges() > 0 )
01468     {
01469       GeomMeasureTool::measure_face_loops(curr_face,
01470                                           curr_min_dist_loops,
01471                                           curr_min_dist_loop,
01472                                           curr_min_angle,
01473                                           curr_max_angle,
01474                                           GEOMETRY_RESABS*500);
01475       GeomMeasureTool::measure_face_area_and_hydraulic_radius(curr_face,
01476                                                               curr_face_area,
01477                                                               curr_face_hydraulic_radius);
01478       GeomMeasureTool::measure_face_curves(curr_face, curve_length,
01479                                            curve_ratio, curr_small_edge);
01480     }
01481     else
01482       curr_face_area = curr_face->measure();
01483          
01484     total_area += curr_face_area;
01485     counter++;
01486     if ( curr_face->num_ref_edges() > 0 )
01487     {
01488       if ( curve_length < min_curve_length )
01489       {
01490         min_curve_length = curve_length;
01491         smallest_edge = curr_small_edge;
01492         face_min_curve_length = curr_face;
01493       }
01494       if ( curve_ratio != -1.0 && curve_ratio > max_curve_ratio )
01495       {
01496         max_curve_ratio = curve_ratio;
01497         face_max_curve_ratio = curr_face;
01498       }
01499       if ( curr_min_dist_loops != -1.0 && curr_min_dist_loops < min_dist_loops )
01500       {
01501         face_min_dist_loops = curr_face;
01502         min_dist_loops = curr_min_dist_loops;
01503       }
01504       if ( curr_min_dist_loop < min_dist_loop )
01505       {
01506         face_min_dist_loop = curr_face;
01507         min_dist_loop = curr_min_dist_loop;
01508       }
01509       if ( curr_min_angle < min_angle )
01510       {
01511         face_min_angle = curr_face;
01512         min_angle = curr_min_angle;
01513       }
01514       if ( curr_max_angle > max_angle )
01515       {
01516         face_max_angle = curr_face;
01517         max_angle = curr_max_angle;
01518       }
01519 
01520       if ( curr_face_hydraulic_radius < min_face_hydraulic_radius )
01521       {
01522         min_face_hydraulic_radius = curr_face_hydraulic_radius;
01523         face_min_hydraulic_radius = curr_face;
01524       }
01525 
01526     }
01527     if ( curr_face_area < min_face_area )
01528     {
01529       min_face_area = curr_face_area;
01530       face_min_area = curr_face;
01531     } 
01532     if ( curr_face_area > max_face_area )
01533     {
01534       max_face_area = curr_face_area;
01535       face_max_area = curr_face;
01536     }
01537   }
01538   if ( counter > 0 )
01539   {
01540     PRINT_INFO("\n******Summary of Surface Measure Info*******\n\n");
01541     PRINT_INFO("Measured %d surfaces: Average Area: %f \n\n", counter, total_area/counter );
01542     PRINT_INFO("Minimum Surface Area is: %f on surface %d\n\n",
01543                min_face_area, face_min_area->id());
01544     PRINT_INFO("Maximum Surface Area is: %f on surface %d\n\n",
01545                max_face_area, face_max_area->id());
01546     if ( face_min_hydraulic_radius )
01547       PRINT_INFO("Minimum Hydraulic Radius (Area/Perimeter) is: %f on surface %d\n\n",
01548                  min_face_hydraulic_radius, face_min_hydraulic_radius->id());
01549     if ( smallest_edge )
01550       PRINT_INFO("Minimum Curve Length is: %f for curve %d on surface %d\n\n",
01551                  min_curve_length, smallest_edge->id(), face_min_curve_length->id());
01552     if ( face_max_curve_ratio )
01553       PRINT_INFO("Maximum Ratio of Adjacent Curve Lengths is: %f on surface %d\n\n",
01554                  max_curve_ratio, face_max_curve_ratio->id());
01555     if ( face_min_dist_loops && face_min_dist_loops->num_loops() > 1 )
01556     {
01557       PRINT_INFO("Minimum Distance Between Loops is: %f on surface %d\n\n",
01558                  min_dist_loops, face_min_dist_loops->id());
01559     }
01560     if ( face_min_dist_loop )
01561       PRINT_INFO("Minimum Distance Between a Single Loop is: %f on surface %d\n\n",
01562                  min_dist_loop, face_min_dist_loop->id());
01563     if ( face_min_angle )
01564       PRINT_INFO("Minimum Interior Angle Between Curves is: %f Degrees on surface %d\n\n",
01565                  min_angle, face_min_angle->id());
01566     if ( face_max_angle )
01567       PRINT_INFO("Maximum Interior Angle Between Curves is: %f Degrees on surface %d\n\n",
01568                  max_angle, face_max_angle->id());
01569     center(ref_faces);
01570     PRINT_INFO("************End of Summary***********\n\n");
01571   }
01572 
01573 }
01574 void GeomMeasureTool::find_adjacent_face_ratios(RefVolume *curr_volume, double &max_face_ratio,
01575                                                 RefFace *&big_face,
01576                                                 RefFace *&small_face)
01577 {
01578     //Loop over all the faces in the model.  Find the maximum face area ratio.
01579     //(larger face area / smaller face area)
01580   
01581   DLIList <Shell*> shells;
01582   curr_volume->shells(shells);
01583   DLIList <RefFace*> ref_face_stack, adj_faces;
01584   Shell* curr_shell;
01585   RefFace *curr_ref_face, *adj_face;
01586     //AreaHashTuple is defined in the .hpp file.
01587   DLIList <AreaHashTuple*> *hash_table;
01588   max_face_ratio = -CUBIT_DBL_MAX;
01589   
01590   
01591   int ii, jj, table_size;
01592   int hash_val;
01593   double tmp_area;
01594   double tmp_ratio;
01595   
01596   AreaHashTuple* curr_tuple, *adj_tuple;
01597   
01598   for ( ii = shells.size(); ii > 0; ii-- )
01599   {
01600     curr_shell = shells.get_and_step();
01601     ref_face_stack.clean_out();
01602     curr_shell->ref_faces(ref_face_stack);
01603     table_size = ref_face_stack.size();
01604       //build a hash table with the most simple form
01605       //of hashing based on the ref face id.
01606       //To get this more reliable, make table size a prime
01607       //number...
01608       //Do this so that we don't have to calculate the
01609       //areas more than once...
01610     hash_table = new DLIList<AreaHashTuple*> [table_size];
01611     for ( jj = 0; jj < ref_face_stack.size(); jj++ )
01612     {
01613       curr_ref_face = ref_face_stack.get_and_step();
01614       tmp_area = curr_ref_face->measure();
01615       curr_tuple = new AreaHashTuple();
01616       curr_tuple->myFace = curr_ref_face;
01617       curr_tuple->myArea = tmp_area;
01618       hash_val = curr_ref_face->id() % table_size;
01619       hash_table[hash_val].append(curr_tuple);
01620     }
01621     
01622     while ( ref_face_stack.size() )
01623     {
01624       curr_ref_face = ref_face_stack.pop();
01625         //Now get the adjacent faces.
01626       get_adjacent_faces(curr_ref_face, adj_faces);
01627       find_index(hash_table, table_size, curr_ref_face, curr_tuple);
01628       assert(curr_tuple != NULL);
01629       if ( curr_tuple == NULL )
01630         return;
01631       CubitBoolean curr_is_big = CUBIT_TRUE;
01632       for ( jj = 0; jj < adj_faces.size(); jj++ )
01633       {
01634         adj_face = adj_faces.get_and_step();
01635         find_index(hash_table, table_size, adj_face, adj_tuple);
01636         if ( adj_tuple == NULL )
01637           continue;
01638         curr_is_big = CUBIT_TRUE;
01639         if ( curr_tuple->myArea < adj_tuple->myArea )
01640         {
01641           tmp_ratio = adj_tuple->myArea / curr_tuple->myArea;
01642           curr_is_big = CUBIT_FALSE;
01643         }
01644         else
01645           tmp_ratio = curr_tuple->myArea / adj_tuple->myArea;
01646         if ( tmp_ratio > max_face_ratio )
01647         {
01648           max_face_ratio = tmp_ratio;
01649           if ( curr_is_big )
01650           {
01651             big_face = curr_ref_face;
01652             small_face = adj_face;
01653           }
01654           else
01655           {
01656             big_face = adj_face;
01657             small_face = curr_ref_face;
01658           }
01659         }
01660       }
01661     }
01662     for ( jj = 0; jj < table_size; jj++ )
01663     {
01664         //delete the area tuples.
01665       while ( hash_table[jj].size() > 0 )
01666         delete hash_table[jj].pop();
01667     }
01668       //delete the hash_table.
01669     delete [] hash_table;
01670   }
01671 }
01672 void GeomMeasureTool::find_index(DLIList<AreaHashTuple*> *hash_table,
01673                                  int table_size,
01674                                  RefFace *ref_face,
01675                                  AreaHashTuple *&curr_tuple )
01676 {
01677   int ii;
01678   int hash_val = ref_face->id() % table_size;
01679   if ( hash_table[hash_val].size() == 0 )
01680   {
01681       //this ref face isn't store here.
01682     curr_tuple = NULL;
01683     return;
01684   }
01685   int curr_size = hash_table[hash_val].size();
01686     //resolve collisions by separate chaining method...
01687   for ( ii = 0;ii < curr_size; ii++ )
01688   {
01689     if ( hash_table[hash_val].next(ii)->myFace == ref_face )
01690     {
01691       curr_tuple = hash_table[hash_val].next(ii);
01692       return;
01693     }
01694   }
01695   curr_tuple = NULL;
01696   return;
01697 }
01698   
01699 void GeomMeasureTool::get_adjacent_faces(RefFace* curr_face,
01700                                          DLIList<RefFace*> &adjacent_faces)
01701 {
01702   DLIList <RefEdge*> ref_edges;
01703   DLIList <RefFace*> tmp_ref_faces;
01704   curr_face->ref_edges(ref_edges);
01705   RefEdge *ref_edge;
01706   RefFace *tmp_ref_face;
01707   int ii, jj;
01708   for ( ii = 0; ii < ref_edges.size(); ii++ )
01709   {
01710     ref_edge = ref_edges.get_and_step();
01711     tmp_ref_faces.clean_out();
01712     ref_edge->ref_faces(tmp_ref_faces);
01713     for ( jj = 0; jj < tmp_ref_faces.size(); jj++ )
01714     {
01715       tmp_ref_face = tmp_ref_faces.get_and_step();
01716       if ( !tmp_ref_face->marked() && tmp_ref_face != curr_face )
01717       {
01718         tmp_ref_face->marked(1);
01719         adjacent_faces.append(tmp_ref_face);
01720       }
01721       else
01722         continue;
01723     }
01724   }
01725   for ( ii = adjacent_faces.size(); ii > 0; ii-- )
01726     adjacent_faces.get_and_step()->marked(0);
01727 }
01728 void GeomMeasureTool::print_volume_measure_summary(DLIList <RefVolume*> &ref_volumes)
01729 {
01730   int ii;
01731   RefVolume *curr_volume;
01732   RefFace *curr_face;
01733   DLIList <RefFace*> ref_faces, tmp_faces;
01734     //First get a list of all the ref_faces.
01735     //Get them uniquely.
01736   for ( ii = 0; ii< ref_volumes.size(); ii++ )
01737   {
01738     curr_volume = ref_volumes.get_and_step();
01739     tmp_faces.clean_out();
01740     curr_volume->ref_faces(tmp_faces);
01741     int jj;
01742     for ( jj = 0; jj < tmp_faces.size(); jj++ )
01743     {
01744       curr_face = tmp_faces.get_and_step();
01745       if ( curr_face->marked() )
01746         continue;
01747       else
01748       {
01749         curr_face->marked(1);
01750         ref_faces.append(curr_face);
01751       }
01752     }
01753   }
01754     //reset the mark.
01755   for ( ii = 0; ii< ref_faces.size(); ii++ )
01756     ref_faces.get_and_step()->marked(0);
01757 
01758   double volume_volume, volume_hydraulic_radius;
01759   double min_angle=CUBIT_DBL_MAX, max_angle=-CUBIT_DBL_MAX;
01760   double total_volume = 0.0, min_volume_hydraulic_radius = CUBIT_DBL_MAX;
01761   double min_volume=CUBIT_DBL_MAX; //, max_volume=-CUBIT_DBL_MAX;
01762   double curr_min_angle = 361.0, curr_max_angle =361.;
01763   double max_face_ratio = -CUBIT_DBL_MAX, curr_face_ratio=CUBIT_DBL_MAX;
01764          
01765 
01766   RefVolume *volume_min_volume = NULL, *volume_min_hydraulic_radius=NULL;
01767   RefVolume *volume_min_angle=NULL, *volume_max_angle=NULL;
01768   RefVolume *volume_face_ratio= NULL;
01769   RefFace *max_small_face=NULL, *max_big_face=NULL;
01770   RefFace *small_face=NULL, *big_face=NULL;
01771   RefFace *face_min_1=NULL, *face_min_2=NULL;
01772   RefFace *face_max_1=NULL, *face_max_2=NULL;
01773   RefFace *face_min_angle_1=NULL, *face_min_angle_2=NULL;
01774   RefFace *face_max_angle_1=NULL, *face_max_angle_2=NULL;
01775   
01776   
01777   int counter = 0;
01778        
01779   for( ii = 0; ii< ref_volumes.size(); ii++ )
01780   {
01781     curr_volume = ref_volumes.get_and_step();
01782     GeomMeasureTool::measure_volume_volume_and_hydraulic_radius(curr_volume,
01783                                                                 volume_volume,
01784                                                                 volume_hydraulic_radius);
01785     if ( curr_volume->num_ref_faces() > 1 )
01786     {
01787       GeomMeasureTool::angles_between_volume_surfaces(curr_volume,
01788                                                       curr_min_angle,
01789                                                       curr_max_angle,
01790                                                       face_min_1, face_min_2,
01791                                                       face_max_1, face_max_2);
01792       GeomMeasureTool::find_adjacent_face_ratios(curr_volume,
01793                                                  curr_face_ratio,
01794                                                  big_face, small_face);
01795     }
01796     total_volume += volume_volume;
01797     counter++;
01798     if ( volume_volume < min_volume )
01799     {
01800       volume_min_volume = curr_volume;
01801       min_volume = volume_volume;
01802     }
01803     if ( volume_hydraulic_radius < min_volume_hydraulic_radius )
01804     {
01805       volume_min_hydraulic_radius = curr_volume;
01806       min_volume_hydraulic_radius = volume_hydraulic_radius;
01807     }
01808     if ( curr_volume->num_ref_faces() > 1 &&  curr_min_angle < min_angle )
01809     {
01810       volume_min_angle = curr_volume;
01811       face_min_angle_1 = face_min_1;
01812       face_min_angle_2 = face_min_2;
01813       min_angle = curr_min_angle;
01814     }
01815     if ( curr_volume->num_ref_faces() > 1 && curr_max_angle > max_angle )
01816     {
01817       volume_max_angle = curr_volume;
01818       max_angle = curr_max_angle;
01819       face_max_angle_1 = face_max_1;
01820       face_max_angle_2 = face_max_2;
01821     }
01822     if ( curr_volume->num_ref_faces() > 1 && curr_face_ratio > max_face_ratio )
01823     {
01824       volume_face_ratio = curr_volume;
01825       max_face_ratio = curr_face_ratio;
01826       max_small_face = small_face;
01827       max_big_face = big_face;
01828     }
01829   }
01830   if ( counter > 0 )
01831   {
01832     PRINT_INFO("*************Volume Summary************\n\n");
01833     PRINT_INFO("Measured %d Volumes: Average Volume was %f.\n\n", 
01834                counter, total_volume/counter);
01835     PRINT_INFO("Minimum Volume was %f in volume %d\n\n",
01836                min_volume, volume_min_volume->id());
01837     PRINT_INFO("Minimum Volume Hydraulic Radius (volume/surface area) was %f in volume %d\n\n",
01838                min_volume_hydraulic_radius, volume_min_hydraulic_radius->id());
01839     if ( volume_min_angle )
01840       PRINT_INFO("Minimum Interior Angle Between Two Surfaces was %f,\n"
01841                  "\tsurfaces %d and %d in volume %d\n\n",
01842                  min_angle, face_min_angle_1->id(), face_min_angle_2->id(), volume_min_angle->id());
01843     if ( volume_max_angle )
01844       PRINT_INFO("Maximum Interior Angle Between Two Surfaces was %f,\n"
01845                  "\tsurfaces %d and %d in volume %d\n\n",
01846                  max_angle, face_max_angle_1->id(), face_max_angle_2->id(), volume_max_angle->id());
01847     if ( volume_face_ratio )
01848       PRINT_INFO("Maximum Area Ratios Between Adjacent Faces was %f,\n"
01849                  "\tsurfaces %d and %d in volume %d.\n\n",
01850                  max_face_ratio, max_big_face->id(), max_small_face->id(), volume_face_ratio->id());
01851     int merged, unmerged;
01852     double ratio;     
01853     GeomMeasureTool::merged_unmerged_surface_ratio(ref_volumes, merged, unmerged, ratio);
01854     PRINT_INFO("Total number of merged surfaces:\t%d\n"
01855                "Total number of unmerged surfaces:\t%d\n"
01856                "Ratio of merged to unmerged surfaces:\t%f\n\n", merged, unmerged, ratio);
01857     RefGroup *group_of_shells;
01858     int number_shells;
01859     GeomMeasureTool::find_shells(ref_volumes, group_of_shells, number_shells);
01860     if ( number_shells > 0 )
01861     {
01862       if ( number_shells == 1 )
01863       {
01864         PRINT_INFO("Found %d shell of connected surfaces.\n", number_shells);
01865         PRINT_INFO("These shell is contained in the group %s\n",
01866                    group_of_shells->entity_name().c_str());
01867       }
01868       else
01869       {
01870         PRINT_INFO("Found %d shells of connected surfaces.\n", number_shells);
01871         PRINT_INFO("These shells are contained in the group %s\n",
01872                    group_of_shells->entity_name().c_str());
01873       }
01874     }
01875     PRINT_INFO("\n\n");
01876     PRINT_INFO("Now Printing Summaries of Surfaces contained by these volumes.\n");
01877     GeomMeasureTool::print_surface_measure_summary(ref_faces);
01878   }
01879   return;
01880 
01881 }
01882 
01883 // Find all of the surfaces in the given volumes that have narrow regions.
01884 void GeomMeasureTool::find_surfs_with_narrow_regions(DLIList <RefVolume*> &ref_vols,
01885                                           double tol,
01886                                           DLIList <RefFace*> &surfs_with_narrow_regions)
01887 {
01888   int j;
01889 
01890   int ii, jj;
01891   DLIList <RefFace*> ref_faces, temp_faces;
01892   RefVolume *ref_vol;
01893   RefFace *curr_face;
01894   for ( ii = 0; ii < ref_vols.size(); ii++ )
01895   {
01896     DLIList<RefFace*> faces;
01897     ref_vol = ref_vols.get_and_step();
01898     ref_vol->ref_faces(faces);
01899     for ( jj = faces.size(); jj > 0; jj-- )
01900     {
01901       curr_face = faces.get_and_step();
01902       curr_face->marked(0);
01903       temp_faces.append(curr_face);
01904     }
01905   }
01906 
01907   //uniquely add the faces.
01908   for ( jj = temp_faces.size(); jj > 0; jj-- )
01909   {
01910     curr_face = temp_faces.get_and_step();
01911     if ( curr_face->marked()== 0 )
01912     {
01913       curr_face->marked(1);
01914       ref_faces.append(curr_face);
01915     }
01916   }
01917 
01918   int num_faces = ref_faces.size();
01919 
01920   ProgressTool *progress_ptr = NULL;
01921   if (num_faces > 20)
01922   {
01923     progress_ptr = AppUtil::instance()->progress_tool();
01924     assert(progress_ptr != NULL);
01925     progress_ptr->start(0, 100, "Finding Surfaces with Narrow Regions", 
01926       NULL, CUBIT_TRUE, CUBIT_TRUE);
01927   }
01928 
01929   int total_faces = 0;
01930   double curr_percent = 0.0;
01931 
01932   for(j=0; j<num_faces; j++)
01933   {
01934     DLIList<CubitVector> split_pos1_list;
01935     DLIList<CubitVector> split_pos2_list;
01936     RefFace *cur_face = ref_faces.get_and_step();
01937     total_faces++;
01938     if ( progress_ptr != NULL )
01939     {
01940       curr_percent = ((double)(total_faces))/((double)(num_faces));
01941       progress_ptr->percent(curr_percent);
01942     }
01943 
01944     if ( AppUtil::instance()->interrupt() )
01945     {
01946         //interrpt.  We need to exit.
01947       if ( progress_ptr != NULL )
01948         progress_ptr->end();
01949         //just leave what has been calculated...
01950       return;
01951     }
01952     find_split_points_for_narrow_regions(cur_face,
01953       tol, split_pos1_list, split_pos2_list); 
01954     if(split_pos1_list.size() > 0)
01955       surfs_with_narrow_regions.append_unique(cur_face);
01956   }
01957 
01958   if ( progress_ptr != NULL )
01959     progress_ptr->end();
01960 }
01961 
01962 void GeomMeasureTool::get_narrow_regions(DLIList <RefVolume*> &ref_vols,
01963                                           double tol,
01964                                           DLIList <RefFace*> &surfs_with_narrow_regions)
01965 {
01966   int ii, jj;
01967   DLIList <RefFace*> ref_faces, temp_faces;
01968   RefVolume *ref_vol;
01969   RefFace *curr_face;
01970   for ( ii = 0; ii < ref_vols.size(); ii++ )
01971   {
01972     DLIList<RefFace*> faces;
01973     ref_vol = ref_vols.get_and_step();
01974     ref_vol->ref_faces(faces);
01975     for ( jj = faces.size(); jj > 0; jj-- )
01976     {
01977       curr_face = faces.get_and_step();
01978       curr_face->marked(0);
01979       temp_faces.append(curr_face);
01980     }
01981   }
01982 
01983   //uniquely add the faces.
01984   for ( jj = temp_faces.size(); jj > 0; jj-- )
01985   {
01986     curr_face = temp_faces.get_and_step();
01987     if ( curr_face->marked()== 0 )
01988     {
01989       curr_face->marked(1);
01990       ref_faces.append(curr_face);
01991     }
01992   }
01993 
01994   int num_faces = ref_faces.size();
01995 
01996   ProgressTool *progress_ptr = NULL;
01997   if (num_faces > 20)
01998   {
01999     progress_ptr = AppUtil::instance()->progress_tool();
02000     assert(progress_ptr != NULL);
02001     progress_ptr->start(0, 100, "Finding Surfaces with Narrow Regions", 
02002       NULL, CUBIT_TRUE, CUBIT_TRUE);
02003   }
02004 
02005   int total_faces = 0;
02006   double curr_percent = 0.0;
02007 
02008   for(jj=0; jj<num_faces; jj++)
02009   {
02010     RefFace *cur_face = ref_faces.get_and_step();
02011     total_faces++;
02012     if ( progress_ptr != NULL )
02013     {
02014       curr_percent = ((double)(total_faces))/((double)(num_faces));
02015       progress_ptr->percent(curr_percent);
02016     }
02017 
02018     if ( AppUtil::instance()->interrupt() )
02019     {
02020         //interrpt.  We need to exit.
02021       if ( progress_ptr != NULL )
02022         progress_ptr->end();
02023         //just leave what has been calculated...
02024       return;
02025     }
02026     if(cur_face->num_loops() == 1)
02027     {
02028       DLIList<CubitVector> split_pos1_list;
02029       DLIList<CubitVector> split_pos2_list;
02030       find_split_points_for_narrow_regions(cur_face,
02031         tol, split_pos1_list, split_pos2_list); 
02032       if(split_pos1_list.size() > 0)
02033         surfs_with_narrow_regions.append(cur_face);
02034     }
02035     else if(cur_face->num_loops() > 1)
02036     {
02037       DLIList <RefEdge*> tmp_close_edges;
02038       DLIList <double> tmp_small_lengths;
02039       find_close_loops( cur_face, tmp_close_edges,
02040                         tmp_small_lengths, tol);
02041       if ( tmp_close_edges.size() > 0 )
02042       {
02043         surfs_with_narrow_regions.append(cur_face);
02044       }
02045     }
02046   }
02047 
02048   if ( progress_ptr != NULL )
02049     progress_ptr->end();
02050 }
02051 
02052 bool GeomMeasureTool::is_surface_narrow(RefFace *face, double small_curve_size)
02053 {
02054   bool ret = true;
02055   DLIList<RefEdge*> edges;
02056   face->ref_edges(edges);
02057   RefVolume *vol = face->ref_volume();
02058   int i, j;
02059   double dist_sq = small_curve_size*small_curve_size;
02060   double proj_dist = 1.2 * small_curve_size;
02061   for(i=edges.size(); i>0 && ret == true; i--)
02062   {
02063     RefEdge *cur_edge = edges.get_and_step();
02064     double edge_length = cur_edge->measure();
02065     if(edge_length > small_curve_size)
02066     {
02067       int num_incs = (int)(edge_length/small_curve_size) + 1;
02068       double start, end;
02069       cur_edge->get_param_range(start, end);
02070       double t = start;
02071       bool one_bad = false;
02072       for(j=0; j<num_incs && ret == true; j++)
02073       {
02074         CubitVector pos1, tangent;
02075         cur_edge->position_from_u(t, pos1);
02076         cur_edge->tangent(pos1, tangent, face);
02077         CubitVector norm = face->normal_at(pos1, vol);
02078         CubitVector indir = norm * tangent;
02079         indir.normalize();
02080         CubitVector new_pos = pos1 + proj_dist * indir;
02081         CubitVector pt_on_surf;
02082         face->get_surface_ptr()->closest_point_trimmed(new_pos, pt_on_surf);
02083         if((pt_on_surf-pos1).length_squared() < dist_sq)
02084         {
02085           one_bad = false;
02086         }
02087         else // we found one out of small_curve range
02088         {
02089           if(one_bad)  // if we had already found one out of range this makes two in a row
02090             ret = false;
02091           else  // this is the first one out of range
02092             one_bad = true;
02093         }
02094       }
02095     }
02096   }
02097 
02098   return ret;
02099 }
02100 
02101 void GeomMeasureTool::find_split_points_for_narrow_regions(RefFace *face,
02102                                                           double size, 
02103                                                           DLIList<CubitVector> &split_pos1_list,
02104                                                           DLIList<CubitVector> &split_pos2_list)
02105 {
02106   int k, i, j;
02107   double size_sq = size*size;
02108   DLIList<RefEdge*> edges;
02109   face->ref_edges(edges);
02110   while(edges.size() > 1)
02111   {
02112     RefEdge *cur_edge = edges.extract();
02113     for(k=edges.size(); k--;)
02114     {
02115       RefEdge *other_edge = edges.get_and_step();
02116       DLIList<CubitVector*> e1_pos_list, e2_pos_list;
02117       DLIList<RefVertex*> e1_vert_list, e2_vert_list;
02118       if(narrow_region_exists(cur_edge, other_edge, face, size,
02119         e1_pos_list, e2_pos_list, e1_vert_list, e2_vert_list))
02120       {
02121         e1_pos_list.reset();
02122         e2_pos_list.reset();
02123         e1_vert_list.reset();
02124         e2_vert_list.reset();
02125 
02126         // Loop through each pair of positions defining a split.
02127         for(i=e1_pos_list.size(); i--;)
02128         {
02129           RefVertex *e1_vert = e1_vert_list.get_and_step();
02130           RefVertex *e2_vert = e2_vert_list.get_and_step();
02131           CubitVector *e1_pos = e1_pos_list.get_and_step();
02132           CubitVector *e2_pos = e2_pos_list.get_and_step();
02133 
02134           // Snap to existing vertices if we are not already at at vertex.
02135           if(!e1_vert)
02136           {
02137             if((cur_edge->start_vertex()->coordinates() - *e1_pos).length_squared() < size_sq)
02138               e1_vert = cur_edge->start_vertex();
02139             else if((cur_edge->end_vertex()->coordinates() - *e1_pos).length_squared() < size_sq)
02140               e1_vert = cur_edge->end_vertex();
02141           }
02142           if(!e2_vert)
02143           {
02144             if((other_edge->start_vertex()->coordinates() - *e2_pos).length_squared() < size_sq)
02145               e2_vert = other_edge->start_vertex();
02146             else if((other_edge->end_vertex()->coordinates() - *e2_pos).length_squared() < size_sq)
02147               e2_vert = other_edge->end_vertex();
02148           }
02149 
02150           // We may have multiple edges separating these two vertices but the accumulated
02151           // length of them may still be within our narrow size so check this.  If this
02152           // is the case we will not want to consider this as a place to split and 
02153           // will as a result discard these positions.
02154           if(e1_vert && e2_vert)
02155           {
02156             double dist = size*sqrt(2.0);
02157             RefVertex *cur_vert = e1_vert;
02158             RefEdge *edge = cur_edge;
02159             double length = 0.0;
02160             while(edge && cur_vert != e2_vert && length <= dist)
02161             {
02162               edge = edge->get_other_curve(cur_vert, face);
02163               if(edge)
02164               {
02165                 length += edge->get_arc_length();
02166                 cur_vert = edge->other_vertex(cur_vert);
02167               }
02168             }
02169             if(length > dist)
02170             {
02171               // We want to keep this split.
02172               split_pos1_list.append(e1_vert->coordinates());
02173               split_pos2_list.append(e2_vert->coordinates());
02174             }
02175           }
02176           else
02177           {
02178             // We want to keep this split.
02179             split_pos1_list.append(*e1_pos);
02180             split_pos2_list.append(*e2_pos);
02181           }
02182         }
02183       }
02184       while(e1_pos_list.size())
02185         delete e1_pos_list.pop();
02186       while(e2_pos_list.size())
02187         delete e2_pos_list.pop();
02188     }
02189   }
02190 
02191   // Make splits unique
02192   DLIList<CubitVector> unique_list1, unique_list2;
02193   split_pos1_list.reset();
02194   split_pos2_list.reset();
02195   for(i=split_pos1_list.size(); i--;)
02196   {
02197     CubitVector p1 = split_pos1_list.get_and_step();
02198     CubitVector p2 = split_pos2_list.get_and_step();
02199     int unique = 1;
02200     for(j=unique_list1.size(); j>0 && unique; j--)
02201     {
02202       CubitVector u1 = unique_list1.get_and_step();
02203       CubitVector u2 = unique_list2.get_and_step();
02204       if((p1.about_equal(u1) && p2.about_equal(u2)) ||
02205         (p1.about_equal(u2) && p2.about_equal(u1)))
02206       {
02207           unique = 0;
02208       }
02209     }
02210     if(unique)
02211     {
02212       unique_list1.append(p1);
02213       unique_list2.append(p2);
02214     }
02215   }
02216   split_pos1_list = unique_list1;
02217   split_pos2_list = unique_list2;
02218 }
02219 
02220 // Checks to see if at the given position the two edges are close together and 
02221 // have the same tangent.
02222 int GeomMeasureTool::is_narrow_region_at_point(RefEdge *e1,
02223                                                RefFace *face,
02224                                                const CubitVector &pt_on_e1,
02225                                                RefEdge *e2,
02226                                                const double &tol_sq,
02227                                                CubitVector &closest)
02228 {
02229   int ret = 0;
02230 
02231   CubitVector tan_1, tan_2;
02232   e2->closest_point_trimmed(pt_on_e1, closest);
02233   
02234   double dist_sq = CUBIT_DBL_MAX;
02235 
02236   // If the surface is not a plane lets project
02237   // the midpoint to the surface and fit a circular
02238   // arc through the 3 points and get the arc length
02239   // to approximate the acutal distance on the surface
02240   // between the 2 original points.
02241   if(PLANE_SURFACE_TYPE != face->geometry_type())
02242   {
02243     CubitVector closest_mid;
02244     CubitVector vec1 = closest-pt_on_e1;
02245     double straight_length_sq = vec1.length_squared();
02246     CubitVector mid = pt_on_e1 + 0.5 * vec1;
02247     face->get_surface_ptr()->closest_point(mid, &closest_mid);
02248     CubitVector mid_vec = mid-closest_mid;
02249     double mid_length_sq = mid_vec.length_squared();
02250     // If we are dealing with very small distances or
02251     // it doesn't look like there is much curvature just use
02252     // the straight line distance.
02253     if(straight_length_sq < 1e-10 || mid_length_sq < 0.01*straight_length_sq)
02254     {
02255       dist_sq = straight_length_sq;
02256     }
02257     else
02258     {
02259       double mid_length = sqrt(mid_length_sq);
02260       double half_length = sqrt(straight_length_sq)/2.0;
02261       // theta is the angle between the vector from the orig point to its projection
02262       // on the other edge and the vector from the orig point to the projected mid point.
02263       double theta = atan(mid_length/half_length); 
02264       // beta is the angle between the vector going from the circle center to
02265       // the orig point and the vector going from the circle center to the
02266       // projected mid point.
02267       double beta = 2.0*theta;
02268       // Some trig to get the radius of the circle that goes through all 3 points
02269       double radius = half_length/sin(beta);
02270       dist_sq = 2.0*radius*beta;
02271       dist_sq *= dist_sq;
02272     }
02273   }
02274   else
02275   {
02276     dist_sq = (pt_on_e1-closest).length_squared(); 
02277   }
02278 
02279   if(dist_sq < tol_sq)
02280   { 
02281     DLIList<CoEdge*> coes;
02282     e1->tangent(pt_on_e1, tan_1);
02283     e2->tangent(closest, tan_2);
02284     e1->get_co_edges(coes, face);
02285     if(coes.size() == 1)
02286     {
02287       if(coes.get()->get_sense() == CUBIT_REVERSED)
02288         tan_1 *= -1.0;
02289       coes.clean_out();
02290       e2->get_co_edges(coes, face);
02291       if(coes.size() == 1)
02292       {
02293         if(coes.get()->get_sense() == CUBIT_REVERSED)
02294           tan_2 *= -1.0;
02295         tan_1.normalize();
02296         tan_2.normalize();
02297         if(tan_1 % tan_2 < -0.9)
02298           ret = 1;
02299       }
02300     }
02301   }
02302   return ret;
02303 }
02304 
02305 int GeomMeasureTool::narrow_region_exists(RefFace *face,
02306                                           const double &tol)
02307 {
02308   int k, ret = 0;
02309   DLIList<RefEdge*> edges;
02310   face->ref_edges(edges);
02311   int num_curves = edges.size();
02312 
02313   // if the surface has no bounding edges then we won't consider
02314   // it narrow.  This can happen with complete sphere surfaces.
02315   if (0 == num_curves)
02316     return 0;
02317 
02318   int num_small_curves = 0;
02319   while(edges.size() > 1 && !ret)
02320   {
02321     // Remove the current edge each time so that we aren't
02322     // doing redundant comparisons.
02323     RefEdge *cur_edge = edges.extract();
02324     if(cur_edge->get_arc_length() < tol)
02325       num_small_curves++;
02326 
02327     // Compare this edge with the remaining edges on the face.
02328     for(k=edges.size(); k && !ret; k--)
02329     {
02330       RefEdge *other_edge = edges.get_and_step();
02331 
02332       DLIList<CubitVector*> e1_pos_list, e2_pos_list;
02333       DLIList<RefVertex*> e1_vert_list, e2_vert_list;
02334       ret = narrow_region_exists(cur_edge, other_edge, face, tol,
02335         e1_pos_list, e2_pos_list, e1_vert_list, e2_vert_list);
02336       while(e1_pos_list.size())
02337         delete e1_pos_list.pop();
02338       while(e2_pos_list.size())
02339         delete e2_pos_list.pop();
02340     }
02341   }
02342   if(!ret)
02343   {
02344     if(edges.size() == 1 && edges.get()->get_arc_length() < tol)
02345       num_small_curves++;
02346   }
02347 
02348   ret = ret || (num_small_curves == num_curves);
02349 
02350   return ret;
02351 }
02352 
02353 int GeomMeasureTool::narrow_region_exists(RefFace *face,
02354                                           RefEdge *edge,
02355                                           const double &tol)
02356 {
02357   int k, ret = 0;
02358   DLIList<RefEdge*> edges;
02359   face->ref_edges(edges);
02360   if(edges.move_to(edge))
02361   {
02362     edges.extract();
02363 
02364     // Compare this edge with the remaining edges on the face.
02365     for(k=edges.size(); k && !ret; k--)
02366     {
02367       RefEdge *other_edge = edges.get_and_step();
02368 
02369       DLIList<CubitVector*> e1_pos_list, e2_pos_list;
02370       DLIList<RefVertex*> e1_vert_list, e2_vert_list;
02371       ret = narrow_region_exists(edge, other_edge, face, tol,
02372         e1_pos_list, e2_pos_list, e1_vert_list, e2_vert_list);
02373       while(e1_pos_list.size())
02374         delete e1_pos_list.pop();
02375       while(e2_pos_list.size())
02376         delete e2_pos_list.pop();
02377     }
02378   }
02379   return ret;
02380 }
02381 
02382 // Finds the start and stop locations between two "close" curves
02383 // that meet the narrow criteria.  The lists that are returned
02384 // mark the beginning and ends of the narrow regions.
02385 int GeomMeasureTool::narrow_region_exists(
02386                                             RefEdge *e1,
02387                                             RefEdge *e2,
02388                                             RefFace *face,
02389                                             const double &tol,
02390                                             DLIList<CubitVector*> &e1_pos_list,
02391                                             DLIList<CubitVector*> &e2_pos_list,
02392                                             DLIList<RefVertex*> &e1_vert_list,
02393                                             DLIList<RefVertex*> &e2_vert_list)
02394 {
02395   int ret = 0;
02396   double tol_sq = tol*tol;
02397   double max_dist_sq = 0.0;
02398   RefVertex *e1_start_vert = e1->start_vertex();
02399   RefVertex *e1_end_vert = e1->end_vertex();
02400 
02401   CubitVector closest;
02402   DLIList<RefVertex*> e1_verts, e2_verts;
02403 
02404   e1->ref_vertices(e1_verts);
02405   e2->ref_vertices(e2_verts);
02406   e1_verts.intersect_unordered(e2_verts);
02407   int num_shared_verts = e1_verts.size();
02408   RefVertex *shared_vert = NULL;
02409   RefEdge *edge1 = NULL;
02410   RefEdge *edge2 = NULL;
02411   if(num_shared_verts == 1)
02412   {
02413     shared_vert = e1_verts.get();
02414     DLIList<CoEdge*> coes;
02415     e1->get_co_edges(coes, face);
02416     if(coes.size() == 1)
02417     {
02418       RefVolume *vol = face->ref_volume();
02419       CubitSense facevol_sense = face->sense(vol);
02420       if((coes.get()->start_vertex() == shared_vert) ==
02421         (facevol_sense == CUBIT_FORWARD))
02422       {
02423         edge1 = e1;
02424         edge2 = e2;
02425       }
02426       else
02427       {
02428         edge1 = e2;
02429         edge2 = e1;
02430       }
02431     }
02432   }
02433 
02434   // Project cur endpoints onto other.
02435 
02436   // First check the endpoints of e1 agaisnt e2
02437   int do_narrow_region_check = 1;
02438   if(num_shared_verts == 1 && shared_vert == e1_start_vert)
02439   {
02440     // Edges are next to each other.  Check the angle between them
02441     // before doing anything else.
02442     double interior_angle = edge1->angle_between(edge2, face);
02443     if(interior_angle > CUBIT_PI/4.0)
02444       do_narrow_region_check = 0;
02445   }
02446   if(do_narrow_region_check &&
02447      is_narrow_region_at_point(e1, face, e1_start_vert->coordinates(), e2, tol_sq, closest))
02448   {
02449     max_dist_sq = (closest - e1_start_vert->coordinates()).length_squared();
02450     e1_pos_list.append(new CubitVector(e1_start_vert->coordinates()));
02451     e2_pos_list.append(new CubitVector(closest));
02452     e1_vert_list.append(e1_start_vert);
02453     e2_vert_list.append(NULL);
02454   }
02455   do_narrow_region_check = 1;
02456   if(num_shared_verts == 1 && shared_vert == e1_end_vert)
02457   {
02458     // Edges are next to each other.  Check the angle between them
02459     // before doing anything else.
02460     double interior_angle = edge1->angle_between(edge2, face);
02461     if(interior_angle > CUBIT_PI/4.0)
02462       do_narrow_region_check = 0;
02463   }
02464   if(do_narrow_region_check &&
02465      is_narrow_region_at_point(e1, face, e1_end_vert->coordinates(), e2, tol_sq, closest))
02466   {
02467     double cur_dist_sq = (closest - e1_end_vert->coordinates()).length_squared();
02468     max_dist_sq = max_dist_sq > cur_dist_sq ? max_dist_sq : cur_dist_sq;
02469     e1_pos_list.append(new CubitVector(e1_end_vert->coordinates()));
02470     e2_pos_list.append(new CubitVector(closest));
02471     e1_vert_list.append(e1_end_vert);
02472     e2_vert_list.append(NULL);
02473   }
02474 
02475   if(e1_pos_list.size() < 2)
02476   {
02477     RefVertex *e2_start_vert = e2->start_vertex();
02478     RefVertex *e2_end_vert = e2->end_vertex();
02479     do_narrow_region_check = 1;
02480     // Now check the end points of e2 against e1
02481     if(num_shared_verts == 1 && shared_vert == e2_start_vert)
02482     {
02483       // Edges are next to each other.  Check the angle between them
02484       // before doing anything else.
02485       double interior_angle = edge1->angle_between(edge2, face);
02486       if(interior_angle > CUBIT_PI/4.0)
02487         do_narrow_region_check = 0;
02488     }
02489     if(do_narrow_region_check &&
02490        is_narrow_region_at_point(e2, face, e2_start_vert->coordinates(), e1, tol_sq, closest))
02491     {
02492       if(e1_pos_list.size() == 0 || !closest.about_equal(*e1_pos_list.get()))
02493       {
02494         double cur_dist_sq = (closest - e2_start_vert->coordinates()).length_squared();
02495         max_dist_sq = max_dist_sq > cur_dist_sq ? max_dist_sq : cur_dist_sq;
02496         e2_pos_list.append(new CubitVector(e2_start_vert->coordinates()));
02497         e1_pos_list.append(new CubitVector(closest));
02498         e2_vert_list.append(e2_start_vert);
02499         e1_vert_list.append(NULL);
02500       }
02501     }
02502     if(e1_pos_list.size() < 2)
02503     {
02504       do_narrow_region_check = 1;
02505       if(num_shared_verts == 1 && shared_vert == e2_end_vert)
02506       {
02507         // Edges are next to each other.  Check the angle between them
02508         // before doing anything else.
02509         double interior_angle = edge1->angle_between(edge2, face);
02510         if(interior_angle > CUBIT_PI/4.0)
02511           do_narrow_region_check = 0;
02512       }
02513       if(do_narrow_region_check &&
02514          is_narrow_region_at_point(e2, face, e2_end_vert->coordinates(), e1, tol_sq, closest))
02515       {
02516         if(e1_pos_list.size() == 0 || !closest.about_equal(*e1_pos_list.get()))
02517         {
02518           double cur_dist_sq = (closest - e2_end_vert->coordinates()).length_squared();
02519           max_dist_sq = max_dist_sq > cur_dist_sq ? max_dist_sq : cur_dist_sq;
02520           e2_pos_list.append(new CubitVector(e2_end_vert->coordinates()));
02521           e1_pos_list.append(new CubitVector(closest));
02522           e2_vert_list.append(e2_end_vert);
02523           e1_vert_list.append(NULL);
02524         }
02525       }
02526     }
02527   }
02528 
02529   // At this point we have projected each curves endpoints onto the other curve 
02530   // and if the projection was "close" we have recorded these as in the lists 
02531   // as close positions.
02532 
02533   // If we have two sets of narrow points check to see if the region
02534   // between them is all narrow.
02535   if(e1_pos_list.size() == 2)
02536   {
02537     int w;
02538     int all_good = 1;
02539  //   double dist_tol = sqrt(max_dist_sq)*5.0;
02540     e1_pos_list.reset();
02541     e2_pos_list.reset();
02542     CubitVector *cur1 = e1_pos_list.get_and_step();
02543     CubitVector *cur2 = e1_pos_list.get();
02544     CubitVector *other1 = e2_pos_list.get_and_step();
02545     CubitVector *other2 = e2_pos_list.get();
02546 
02547     double len1 = e1->get_arc_length(*cur1, *cur2);
02548     //We have to check for periodic edges because the narrow points
02549     //will be the same when one or both of the edges is periodic.
02550     if(len1 > tol || (e1->is_periodic()|| e2->is_periodic()))
02551     {
02552       double len2 = e2->get_arc_length(*other1, *other2);
02553       if(len2 > tol || (e1->is_periodic()|| e2->is_periodic()))
02554       {
02555         double cur_param1 = e1->u_from_position(*cur1);
02556         double cur_param2 = e1->u_from_position(*cur2);
02557         int num_divisions = 2;
02558         CubitVector cur_pos;
02559         double param_step = (cur_param2-cur_param1)/num_divisions;
02560         double cur_param = cur_param1 + param_step;
02561         for(w=1; w<num_divisions && all_good; w++)
02562         {
02563           e1->position_from_u(cur_param, cur_pos);
02564           cur_param += param_step;
02565           if(is_narrow_region_at_point(e1, face, cur_pos, e2, tol_sq, closest))
02566           {
02567             // Sanity check to make sure we aren't splitting off negative space.
02568             CubitVector mid = (cur_pos + closest)/2.0;
02569             CubitVector tmp_pt;
02570             face->get_surface_ptr()->closest_point_trimmed(mid, tmp_pt);
02571             // If it didn't move we are fine but if it did we need to check
02572             // a little more.
02573             if(!mid.about_equal(tmp_pt))
02574             {
02575               // It is possible that the midpoint between the two curves does not lie on the
02576               // surface.  This would be the case in a blend surface.  Check to see that
02577               // the projected point is just an offset along the normal of the surface
02578               // rather than a laterl movement because the space between the two 
02579               // curves was actually negative space and not part of the surface.
02580               CubitVector norm = face->normal_at(tmp_pt);
02581               CubitVector dir(tmp_pt - mid);
02582               dir.normalize();
02583               if(fabs(norm % dir) < .9)
02584                 all_good = 0;
02585             }
02586           }
02587         }
02588       }
02589       else
02590         all_good = 0;
02591     }
02592     else 
02593       all_good = 0;
02594 
02595     if(all_good)
02596       ret = 1;
02597     else
02598     {
02599       e1_pos_list.remove(cur1);
02600       e1_pos_list.remove(cur2);
02601       e2_pos_list.remove(other1);
02602       e2_pos_list.remove(other2);
02603       delete cur1;
02604       delete cur2;
02605       delete other1;
02606       delete other2;
02607     }
02608   }
02609 
02610   // If we didn't find a continous region that was narrow we may curves
02611   // that start narrow at one or more of the ends and then move away from each other.
02612   // We need to find if there is a significant region that is narrow and identify
02613   // where the two curves start diverging.  Add additional points to the lists
02614   // at the points where things start diverging.
02615   if(!ret && e1_pos_list.size() > 0)
02616   {
02617     int i;
02618     e1_pos_list.reset();
02619     e2_pos_list.reset();
02620     e1_vert_list.reset();
02621     e2_vert_list.reset();
02622     for(i=e1_pos_list.size(); i--;)
02623     {
02624       RefVertex *e1_vert = e1_vert_list.get_and_step();
02625       RefVertex *e2_vert = e2_vert_list.get_and_step();
02626 
02627       RefVertex *cur_vert = NULL;
02628       RefEdge *cur_edge = NULL, *other_edge = NULL;
02629       if(e1_vert)
02630       {
02631         cur_edge = e1;
02632         other_edge = e2;
02633         cur_vert = e1_vert;
02634       }
02635       else if(e2_vert)
02636       {
02637         cur_edge = e2;
02638         other_edge = e1;
02639         cur_vert = e2_vert;
02640       }
02641       if(cur_vert)
02642       {
02643         CubitVector next_pos;
02644         CubitVector prev_pos = cur_vert->coordinates();
02645         int num_incs = 20;
02646         double step = cur_edge->get_arc_length()/(double)num_incs;
02647         int still_good = 1;
02648         int cntr = 0;
02649         double direction = (cur_vert == cur_edge->start_vertex() ? 1.0 : -1.0);
02650         // Do coarse traversal along curve to see where we start deviating from
02651         // narrow.
02652         while(still_good)
02653         {
02654           cur_edge->point_from_arc_length(prev_pos, step*direction, next_pos);
02655           if(is_narrow_region_at_point(cur_edge, face, next_pos, other_edge, tol_sq, closest) &&
02656                       cntr < num_incs)
02657           {
02658             prev_pos = next_pos;
02659             cntr++;
02660           }
02661           else
02662             still_good = 0;
02663         }
02664         if(cntr < num_incs)
02665         {
02666           cntr = 0;
02667           double cur_arc_length = cur_edge->get_arc_length(prev_pos, next_pos);
02668           // Do bisection on remaining interval to zero in on point where
02669           // we go from narrow to non-narrow.
02670           CubitVector mid_pos;
02671           double close_enough = tol/20.0;
02672           while(cur_arc_length > close_enough && cntr < 100)
02673           {
02674             cntr++;  // prevent infinite looping
02675             cur_edge->point_from_arc_length(prev_pos, cur_arc_length*direction/2.0, mid_pos);
02676             if(is_narrow_region_at_point(cur_edge, face, mid_pos, other_edge, tol_sq, closest))
02677               prev_pos = mid_pos;
02678             else
02679               next_pos = mid_pos;
02680             cur_arc_length = cur_edge->get_arc_length(prev_pos, next_pos);
02681           }
02682           if(cur_edge->get_arc_length(cur_vert->coordinates(), prev_pos) > tol)
02683           {
02684             // end up with the position that guarantees a new split curve that
02685             // is smaller than the small curve size.
02686             other_edge->closest_point_trimmed(prev_pos, closest);
02687             if(cur_edge == e1)
02688             {
02689               e1_pos_list.append(new CubitVector(prev_pos));
02690               e2_pos_list.append(new CubitVector(closest));
02691             }
02692             else
02693             {
02694               e2_pos_list.append(new CubitVector(prev_pos));
02695               e1_pos_list.append(new CubitVector(closest));
02696             }
02697             e1_vert_list.append(NULL);
02698             e2_vert_list.append(NULL);
02699             ret = 1;
02700           }
02701         }
02702       }
02703     }
02704   }
02705 
02706   // Now remove and regions of close curves between which 
02707   // there is negative space (not actually part of the surface).
02708   if(ret)
02709   {
02710     int i, j;
02711     e1_pos_list.reset();
02712     e2_pos_list.reset();
02713     e1_vert_list.reset();
02714     e2_vert_list.reset();
02715     for(i=e1_pos_list.size(); i--;)
02716     {
02717       CubitVector *e1_pos = e1_pos_list.get();
02718       CubitVector *e2_pos = e2_pos_list.get();
02719       int num_divisions = 6;
02720       CubitVector step = (*e2_pos - *e1_pos)/num_divisions;
02721       CubitVector cur_pos = *e1_pos + step;
02722       int removed = 0;
02723       for(j=1; j<num_divisions; j++)
02724       {
02725         CubitVector tmp_pt;
02726         face->get_surface_ptr()->closest_point_trimmed(cur_pos, tmp_pt);
02727         if(!cur_pos.about_equal(tmp_pt))
02728         {
02729           CubitVector norm = face->normal_at(tmp_pt);
02730           CubitVector dir(tmp_pt - cur_pos);
02731           dir.normalize();
02732           if(fabs(norm % dir) < .9)
02733           {
02734             removed = 1;
02735             j=num_divisions;
02736             delete e1_pos_list.remove();
02737             delete e2_pos_list.remove();
02738             e1_vert_list.remove();
02739             e2_vert_list.remove();
02740           }
02741         }
02742         else
02743           cur_pos += step;
02744       }
02745       if(!removed)
02746       {
02747         e1_pos_list.step();
02748         e2_pos_list.step();
02749         e1_vert_list.step();
02750         e2_vert_list.step();
02751       }
02752     }
02753 
02754     if(e1_pos_list.size() == 0)
02755       ret = 0;
02756   }
02757 
02758   return ret;
02759 }
02760 
02761 void GeomMeasureTool::find_small_curves( DLIList <RefVolume*> &ref_vols,
02762                                          double tol,
02763                                          DLIList <RefEdge*> &small_curves,
02764                                          DLIList <double> &small_lengths)
02765 {
02766   int ii, jj;
02767   DLIList <RefEdge*> ref_edges, temp_edges;
02768   RefVolume *ref_vol;
02769   RefEdge *curr_edge;
02770   for ( ii = 0; ii < ref_vols.size(); ii++ )
02771   {
02772     ref_vol = ref_vols.get_and_step();
02773     ref_vol->ref_edges(temp_edges);
02774       //uniquely add the edges.
02775     for ( jj = temp_edges.size(); jj > 0; jj-- )
02776     {
02777       curr_edge = temp_edges.pop();
02778       if ( curr_edge->marked()== 0 )
02779       {
02780         curr_edge->marked(1);
02781         ref_edges.append(curr_edge);
02782       }
02783     }
02784   }
02785 
02786   int num_curves = ref_edges.size();
02787   ProgressTool *progress_ptr = NULL;
02788   if (num_curves> 20)
02789   {
02790     progress_ptr = AppUtil::instance()->progress_tool();
02791     assert(progress_ptr != NULL);
02792     progress_ptr->start(0, 100, "Finding Small Curves", 
02793       NULL, CUBIT_TRUE, CUBIT_TRUE);
02794   }
02795 
02796   int total_curves = 0;
02797   double curr_percent = 0.0;
02798   double length;
02799     //find the small curves and reset the marked flag.
02800   for ( ii = ref_edges.size(); ii > 0; ii-- )
02801   {
02802     total_curves++;
02803     if ( progress_ptr != NULL )
02804     {
02805       curr_percent = ((double)(total_curves))/((double)(num_curves));
02806       progress_ptr->percent(curr_percent);
02807     }
02808     if ( AppUtil::instance()->interrupt() )
02809     {
02810         //interrpt.  We need to exit.
02811       if ( progress_ptr != NULL )
02812         progress_ptr->end();
02813         //just leave what has been calculated...
02814       return;
02815     }
02816 
02817     curr_edge = ref_edges.get_and_step();
02818     curr_edge->marked(0);
02819 
02820     //skip point curves
02821     if( curr_edge->geometry_type() == POINT_CURVE_TYPE )
02822       continue;
02823 
02824     length = curr_edge->measure();
02825     if ( length <= tol )
02826     {
02827       small_curves.append(curr_edge);
02828       small_lengths.append(length);
02829     }
02830   }
02831 
02832   if ( progress_ptr != NULL )
02833     progress_ptr->end();
02834 
02835   return;
02836 }
02837 
02838 RefEdge* GeomMeasureTool::find_first_small_curve(RefVolume* vol,
02839                                          double tol)
02840 {
02841   RefEdge *ret = NULL;
02842   int j;
02843   DLIList <RefEdge*> ref_edges;
02844   vol->ref_edges(ref_edges);
02845   for(j=ref_edges.size(); j > 0 && !ret; j--)
02846   {
02847     RefEdge *curr_edge = ref_edges.get_and_step();
02848     if(curr_edge->measure() <= tol)
02849       ret = curr_edge;
02850   }
02851   return ret;
02852 }
02853 
02854 void GeomMeasureTool::find_narrow_faces(DLIList<RefVolume*> &ref_vols,
02855                                         double small_curve_size,
02856                                         DLIList<RefFace*> &narrow_faces,
02857                                         DLIList<RefFace*> &surfs_to_ignore)
02858 {
02859   int ii, jj;
02860   DLIList <RefFace*> ref_faces, temp_faces;
02861   RefVolume *ref_vol;
02862   RefFace *curr_face;
02863 
02864   for ( ii = 0; ii < ref_vols.size(); ii++ )
02865   {
02866     DLIList<RefFace*> faces;
02867     ref_vol = ref_vols.get_and_step();
02868     ref_vol->ref_faces(faces);
02869     for ( jj = faces.size(); jj > 0; jj-- )
02870     {
02871       curr_face = faces.get_and_step();
02872       curr_face->marked(0);
02873       temp_faces.append(curr_face);
02874     }
02875   }
02876 
02877   //uniquely add the faces.
02878   for ( jj = temp_faces.size(); jj > 0; jj-- )
02879   {
02880     curr_face = temp_faces.get_and_step();
02881     if ( curr_face->marked()== 0 )
02882     {
02883       curr_face->marked(1);
02884       ref_faces.append(curr_face);
02885     }
02886   }
02887 
02888   int num_faces = ref_faces.size();
02889   ProgressTool *progress_ptr = NULL;
02890   if (num_faces > 20)
02891   {
02892     progress_ptr = AppUtil::instance()->progress_tool();
02893     assert(progress_ptr != NULL);
02894     progress_ptr->start(0, 100, "Finding Narrow Surfaces", 
02895       NULL, CUBIT_TRUE, CUBIT_TRUE);
02896   }
02897 
02898   int total_faces = 0;
02899   double curr_percent = 0.0;
02900 
02901   for ( ii = ref_faces.size(); ii > 0; ii-- )
02902   {
02903     total_faces++;
02904     if ( progress_ptr != NULL )
02905     {
02906       curr_percent = ((double)(total_faces))/((double)(num_faces));
02907       progress_ptr->percent(curr_percent);
02908     }
02909 
02910     if ( AppUtil::instance()->interrupt() )
02911     {
02912         //interrpt.  We need to exit.
02913       if ( progress_ptr != NULL )
02914         progress_ptr->end();
02915         //just leave what has been calculated...
02916       return;
02917     }
02918 
02919     curr_face = ref_faces.get_and_step();
02920     if(!surfs_to_ignore.is_in_list(curr_face))
02921     {
02922       if(narrow_region_exists(curr_face, small_curve_size))
02923       {
02924         DLIList<CubitVector> split_pos1_list;
02925         DLIList<CubitVector> split_pos2_list;
02926         find_split_points_for_narrow_regions(curr_face,
02927           small_curve_size, split_pos1_list, split_pos2_list); 
02928         if(split_pos1_list.size() == 0)
02929           narrow_faces.append_unique(curr_face);
02930       }
02931     }
02932   }
02933 
02934   if ( progress_ptr != NULL )
02935     progress_ptr->end();
02936 
02937   return;
02938 }
02939 
02940 void GeomMeasureTool::find_small_faces( DLIList <RefVolume*> &ref_vols,
02941                                         double tol,
02942                                         DLIList <RefFace*> &small_faces)
02943 {
02944   int ii, jj;
02945   DLIList <RefFace*> ref_faces, temp_faces;
02946   RefVolume *ref_vol;
02947   RefFace *curr_face;
02948   for ( ii = 0; ii < ref_vols.size(); ii++ )
02949   {
02950     DLIList<RefFace*> faces;
02951     ref_vol = ref_vols.get_and_step();
02952     ref_vol->ref_faces(faces);
02953     for ( jj = faces.size(); jj > 0; jj-- )
02954     {
02955       curr_face = faces.get_and_step();
02956       curr_face->marked(0);
02957       temp_faces.append(curr_face);
02958     }
02959   }
02960 
02961   //uniquely add the faces.
02962   for ( jj = temp_faces.size(); jj > 0; jj-- )
02963   {
02964     curr_face = temp_faces.get_and_step();
02965     if ( curr_face->marked()== 0 )
02966     {
02967       curr_face->marked(1);
02968       ref_faces.append(curr_face);
02969     }
02970   }
02971 
02972   int num_faces = ref_faces.size();
02973   ProgressTool *progress_ptr = NULL;
02974   if (num_faces > 20)
02975   {
02976     progress_ptr = AppUtil::instance()->progress_tool();
02977     assert(progress_ptr != NULL);
02978     progress_ptr->start(0, 100, "Finding Small Surfaces", 
02979       NULL, CUBIT_TRUE, CUBIT_TRUE);
02980   }
02981 
02982   int total_faces = 0;
02983   double curr_percent = 0.0;
02984   double area;
02985     //find the small curves and reset the marked flag.
02986   for ( ii = ref_faces.size(); ii > 0; ii-- )
02987   {
02988     total_faces++;
02989     if ( progress_ptr != NULL )
02990     {
02991       curr_percent = ((double)(total_faces))/((double)(num_faces));
02992       progress_ptr->percent(curr_percent);
02993     }
02994 
02995     if ( AppUtil::instance()->interrupt() )
02996     {
02997         //interrpt.  We need to exit.
02998       if ( progress_ptr != NULL )
02999         progress_ptr->end();
03000         //just leave what has been calculated...
03001       return;
03002     }
03003 
03004     curr_face = ref_faces.get_and_step();
03005     area = measure_area(curr_face);
03006     if ( area <= tol )
03007       small_faces.append(curr_face);
03008   }
03009 
03010   if ( progress_ptr != NULL )
03011     progress_ptr->end();
03012 
03013   return;
03014 }
03015 
03016 void GeomMeasureTool::find_closed_narrow_faces( DLIList <RefVolume*> &ref_vols,
03017                                         double tol,
03018                                         DLIList <RefFace*> &face_list)
03019 {
03020   int ii, jj;
03021   DLIList <RefFace*> ref_faces, temp_faces;
03022   RefVolume *ref_vol;
03023   RefFace *curr_face;
03024   for ( ii = 0; ii < ref_vols.size(); ii++ )
03025   {
03026     DLIList<RefFace*> faces;
03027     ref_vol = ref_vols.get_and_step();
03028     ref_vol->ref_faces(faces);
03029     for ( jj = faces.size(); jj > 0; jj-- )
03030     {
03031       curr_face = faces.get_and_step();
03032       curr_face->marked(0);
03033       temp_faces.append(curr_face);
03034     }
03035   }
03036 
03037   //uniquely add the faces.
03038   for ( jj = temp_faces.size(); jj > 0; jj-- )
03039   {
03040     curr_face = temp_faces.get_and_step();
03041     if ( curr_face->marked()== 0 )
03042     {
03043       curr_face->marked(1);
03044       ref_faces.append(curr_face);
03045     }
03046   }
03047 
03048   int num_faces = ref_faces.size();
03049   ProgressTool *progress_ptr = NULL;
03050   if (num_faces > 20)
03051   {
03052     progress_ptr = AppUtil::instance()->progress_tool();
03053     assert(progress_ptr != NULL);
03054     progress_ptr->start(0, 100, "Finding Small Surfaces", 
03055       NULL, CUBIT_TRUE, CUBIT_TRUE);
03056   }
03057 
03058   int total_faces = 0;
03059   double curr_percent = 0.0;
03060   for ( ii = ref_faces.size(); ii > 0; ii-- )
03061   {
03062     total_faces++;
03063     if ( progress_ptr != NULL )
03064     {
03065       curr_percent = ((double)(total_faces))/((double)(num_faces));
03066       progress_ptr->percent(curr_percent);
03067     }
03068 
03069     if ( AppUtil::instance()->interrupt() )
03070     {
03071         //interrpt.  We need to exit.
03072       if ( progress_ptr != NULL )
03073         progress_ptr->end();
03074         //just leave what has been calculated...
03075       return;
03076     }
03077 
03078     curr_face = ref_faces.get_and_step();
03079     if(curr_face->get_surface_ptr()->is_closed_in_U() || curr_face->get_surface_ptr()->is_closed_in_V())
03080     {
03081       if(curr_face->num_loops() == 2)
03082       {
03083         DLIList<RefEdge*> ref_edge_list;
03084         curr_face->ref_edges(ref_edge_list);
03085         if(ref_edge_list.size() == 2)
03086         {
03087           int num_pts = 5;
03088           double start, end;
03089           RefEdge *e1 = ref_edge_list.get_and_step();
03090           RefEdge *e2 = ref_edge_list.get();
03091           e1->get_param_range(start, end);
03092           double dt = (end-start)/(double)(num_pts-1);
03093           double t = start;
03094           bool one_bad = false;
03095           double dist_sq = tol*tol;
03096           for(jj=0; jj<num_pts && one_bad == false; jj++)
03097           {
03098             CubitVector pos1, pos2;
03099             e1->position_from_u(t, pos1);
03100             e2->closest_point_trimmed(pos1, pos2);
03101             if((pos1-pos2).length_squared() > dist_sq)
03102               one_bad = true;
03103             t += dt;
03104           }
03105           if(one_bad == false)
03106           {
03107             face_list.append(curr_face);
03108           }
03109         }
03110       }
03111     }
03112   }
03113 
03114   if ( progress_ptr != NULL )
03115     progress_ptr->end();
03116 
03117   return;
03118 }
03119 
03120 void GeomMeasureTool::find_small_faces_hydraulic_radius( DLIList <RefVolume*> &ref_vols,
03121                                                          double tol,
03122                                                          DLIList <RefFace*> &small_faces,
03123                                                          DLIList <double> &small_hyd_rad,
03124                                                          double &percent_planar,
03125                                                          double &percent_pl_co)
03126 {
03127   int ii, jj;
03128   DLIList <RefFace*> ref_faces, temp_faces;
03129   RefVolume *ref_vol;
03130   RefFace *curr_face;
03131   for ( ii = 0; ii < ref_vols.size(); ii++ )
03132   {
03133     ref_vol = ref_vols.get_and_step();
03134     ref_vol->ref_faces(temp_faces);
03135       //uniquely add the faces.
03136     for ( jj = temp_faces.size(); jj > 0; jj-- )
03137     {
03138       curr_face = temp_faces.pop();
03139       if ( curr_face->marked()== 0 )
03140       {
03141         curr_face->marked(1);
03142         ref_faces.append(curr_face);
03143       }
03144     }
03145   }
03146   double area;
03147   double length;
03148   int total_faces = 0;
03149   int total_plane = 0;
03150   int total_cone = 0;
03151   double area_plane = 0.0;
03152   double total_area = 0.0;
03153     //find the small curves and reset the marked flag.
03154   DLIList <CoEdge*> co_edges;
03155   int num_faces = ref_faces.size();
03156   ProgressTool *progress_ptr = NULL;
03157   if (num_faces > 20)
03158   {
03159     progress_ptr = AppUtil::instance()->progress_tool();
03160     assert(progress_ptr != NULL);
03161     progress_ptr->start(0, 100, "Small Surface Progress", 
03162       NULL, CUBIT_TRUE, CUBIT_TRUE);
03163   }
03164   double curr_percent = 0.0;
03165   for ( ii = ref_faces.size(); ii > 0; ii-- )
03166   {
03167     curr_face = ref_faces.get_and_step();
03168     curr_face->marked(0);
03169     area = measure_area(curr_face);
03170     total_area += area;
03171     total_faces++;
03172       //update the progress..
03173     if ( progress_ptr != NULL )
03174     {
03175       curr_percent = ((double)(total_faces))/((double)(num_faces));
03176       progress_ptr->percent(curr_percent);
03177     }
03178     if ( AppUtil::instance()->interrupt() )
03179     {
03180         //interrpt.  We need to exit.
03181       if ( progress_ptr != NULL )
03182         progress_ptr->end();
03183         //just leave what has been calculated...
03184       return;
03185     }
03186     if ( curr_face->geometry_type() == PLANE_SURFACE_TYPE )
03187     {
03188       total_plane++;
03189       area_plane += area;
03190     }
03191     else if ( curr_face->geometry_type() == CONE_SURFACE_TYPE )
03192       total_cone++;
03193     co_edges.clean_out();
03194     curr_face->co_edges(co_edges);
03195     length = 0.0;
03196     for ( jj = co_edges.size(); jj > 0; jj-- )
03197     {
03198       length += co_edges.get_and_step()->get_ref_edge_ptr()->measure();
03199     }
03200       //reset the mark.
03201     curr_face->marked(0);
03202       //compute the hydraulic radius 4*(A/P).
03203     if ( length <= CUBIT_RESABS )
03204     {
03205       PRINT_INFO("Total Perimeter Length of Surface %d is less than tolerance.\n",
03206                  curr_face->id());
03207       
03208       continue;
03209     }
03210     area = 4.0*(area/length);
03211     if ( area <= tol )
03212     {
03213       small_faces.append(curr_face);
03214       small_hyd_rad.append(area);
03215     }
03216   }
03217   if ( total_faces > 0 )
03218   {
03219     percent_planar = 100.0*area_plane/total_area;
03220     percent_pl_co = 100.0*((double)total_plane+(double)total_cone)/(double)total_faces;
03221   }
03222   else
03223   {
03224     percent_planar = 0.;
03225     percent_pl_co = 0.;
03226   }
03227   if ( progress_ptr != NULL )
03228   {
03229     progress_ptr->end();
03230   }
03231   
03232   return;
03233 }
03234 void GeomMeasureTool::find_small_volumes( DLIList <RefVolume*> &ref_vols,
03235                                           double tol,
03236                                           DLIList <RefVolume*> &small_volumes)
03237 {
03238   int ii;
03239   RefVolume *ref_vol;
03240   double volume;
03241   
03242   for ( ii = 0; ii < ref_vols.size(); ii++ )
03243   {
03244     ref_vol = ref_vols.get_and_step();
03245     volume = ref_vol->measure();
03246     if ( volume <= tol )
03247       small_volumes.append(ref_vol);
03248   }
03249   return;
03250 }
03251 void GeomMeasureTool::find_small_volumes_hydraulic_radius( DLIList <RefVolume*> &ref_vols,
03252                                                            double tol,
03253                                                            DLIList <RefVolume*> &small_volumes,
03254                                                            DLIList <double> &small_hyd_rad)
03255 {
03256   int ii, jj;
03257   DLIList <RefFace*> ref_faces, temp_faces;
03258   RefVolume *ref_vol;
03259   RefFace *curr_face;
03260   double area;
03261   double hyd_volume;
03262   
03263   for ( ii = 0; ii < ref_vols.size(); ii++ )
03264   {
03265     ref_vol = ref_vols.get_and_step();
03266     ref_vol->ref_faces(temp_faces);
03267     area = 0.0;
03268     for ( jj = temp_faces.size(); jj > 0; jj-- )
03269     {
03270       curr_face = temp_faces.pop();
03271       area += measure_area(curr_face);
03272     }
03273     if ( area <= CUBIT_RESABS )
03274       continue;
03275     hyd_volume = 6.0*(ref_vol->measure() / area);
03276     if ( hyd_volume <= tol )
03277     {
03278       small_hyd_rad.append(hyd_volume);
03279       small_volumes.append(ref_vol);
03280     }
03281   }
03282   return;
03283 }
03298 void GeomMeasureTool::find_sharp_tangential_meets( RefVolume *ref_volume,
03299                                                    DLIList <RefEdge*> &small_angle_edge_pairs,
03300                                                    DLIList <RefFace*> &tangential_surface_pairs )
03301 {
03302     //Okay, given the small angle edge pairs.  See if there are surfaces that meet tangentially.
03303     //I'm really looking for something spefic.
03304   RefEdge *ref_edge_1, *ref_edge_2;
03305   int ii, jj;
03306   DLIList <Loop*> loops_1, loops_2;
03307   Loop *common_loop;
03308   const double rad_to_deg = 180.0/CUBIT_PI;
03309   RefFace *ref_face;
03310   
03311     //loop over for each edge pair.
03312   for ( ii = 0; ii < small_angle_edge_pairs.size(); ii += 2 )
03313   {
03314     ref_edge_1 = small_angle_edge_pairs.get_and_step();
03315     ref_edge_2 = small_angle_edge_pairs.get_and_step();
03316       //find the loop that these edges are both on.
03317     loops_1.clean_out();
03318     loops_2.clean_out();
03319     ref_edge_1->loops(loops_1);
03320     ref_edge_2->loops(loops_2);
03321     common_loop = NULL;
03322     for ( jj = loops_1.size(); jj > 0; jj-- )
03323     {
03324       if ( !loops_2.move_to(loops_1.get()) )
03325         loops_1.remove();
03326       else
03327         loops_1.step();
03328     }
03329     if ( loops_1.size() == 1 )
03330       common_loop = loops_1.get();
03331     else
03332     {
03333         //find the loop that has a co_edge from ref_edge_1
03334         //then a co_edge from ref_edge_2.
03335       DLIList <CoEdge*> co_edges;
03336       Loop *tmp_loop;
03337       int kk;
03338       CoEdge *co_edge;
03339       for ( jj = loops_1.size(); jj > 0; jj-- )
03340       {
03341         tmp_loop = loops_1.get_and_step();
03342         co_edges.clean_out();
03343         tmp_loop->ordered_co_edges(co_edges);
03344         for ( kk = co_edges.size(); kk > 0; kk-- )
03345         {
03346           co_edge = co_edges.get_and_step();
03347           if ( co_edge->get_ref_edge_ptr() == ref_edge_1 )
03348           {
03349             if ( co_edges.get()->get_ref_edge_ptr() == ref_edge_2 )
03350             {
03351               common_loop = tmp_loop;
03352               break;
03353             }
03354             else
03355                 //not this loop..
03356               break;
03357           }
03358         }
03359         if ( common_loop != NULL )
03360           break;
03361       }
03362     }
03363     assert ( common_loop != NULL );
03364       //now find surface.
03365       //Then see if there are surfaces that meet at 180 degrees that
03366       //are 90* rotated from here.
03367     ref_face = common_loop->get_ref_face_ptr();
03368       //Okay, refface is the surface on which we have the bad angle.
03369       //One of these curves will have a concave angle and the other a
03370       //convex angle.  The surface on the other side of the convex angle
03371       //we need to check to see if it comes tangentially into another
03372       //surface (dihedral angle of 180.)
03373     RefFace *ref_face_1 = ref_edge_1->other_face(ref_face, ref_volume);
03374     if( ref_face_1 == NULL )
03375       continue;
03376 
03377     double surf_angle_1, surf_angle_2;
03378     RefFace *tangential_face = NULL;
03379     RefEdge *edge_next_to = NULL;
03380     surf_angle_1 =  GeometryQueryTool::instance()->surface_angle(ref_face,
03381                                                                  ref_face_1,
03382                                                                  ref_edge_1,
03383                                                                  ref_volume);
03384     surf_angle_1 *= rad_to_deg;
03385     RefFace *ref_face_2;
03386     ref_face_2 = ref_edge_2->other_face(ref_face, ref_volume);
03387     if( ref_face_2 == NULL )
03388       continue;
03389     surf_angle_2 = GeometryQueryTool::instance()->surface_angle(ref_face,
03390                                                                 ref_face_2,
03391                                                                 ref_edge_2,
03392                                                                 ref_volume);
03393     surf_angle_2 *= rad_to_deg;
03394     if ( surf_angle_1 < 180.0 && surf_angle_2 > 180.0 )
03395     {
03396       tangential_face = ref_face_1;
03397       edge_next_to = ref_edge_1;
03398     }
03399     else if ( surf_angle_2 < 180.0 && surf_angle_1 > 180.0 )
03400     {
03401       tangential_face = ref_face_2;
03402       edge_next_to = ref_edge_2;
03403     }
03404     if ( tangential_face == NULL ) // This isn't the config we seek.
03405       continue;
03406 
03407       //Now, on this face, there is an edge, next to edge_next_to,
03408       //that will have an interior angle of 180.0.  If that is the
03409       //case then this is what we seek.
03410     RefVertex* ref_vert = ref_edge_1->common_ref_vertex(ref_edge_2,
03411                                                         ref_face);
03412     DLIList <RefEdge*> tmp_ref_edges;
03413     ref_vert->ref_edges(tmp_ref_edges);
03414       //Find the edge that is on tangential face and is not edge_next_to.
03415     RefEdge *tangential_edge=NULL, *tmp_edge;
03416     for ( jj = tmp_ref_edges.size(); jj > 0; jj-- )
03417     {
03418       tmp_edge = tmp_ref_edges.get_and_step();
03419       if ( tmp_edge != edge_next_to &&
03420            tmp_edge->is_directly_related(tangential_face) )
03421       {
03422         tangential_edge = tmp_edge;
03423         break;
03424       }
03425     }
03426     if ( tangential_edge == NULL )
03427         //not what we are looking for.
03428       continue;
03429       //Now get the other face and measure the dihedral angle.
03430     RefFace *other_face = tangential_edge->other_face(tangential_face,
03431                                                       ref_volume);
03432     double tan_angle;
03433     tan_angle = GeometryQueryTool::instance()->surface_angle(tangential_face,
03434                                                              other_face,
03435                                                              tangential_edge,
03436                                                              ref_volume);
03437     tan_angle *= rad_to_deg;
03438       //If they are within tolerance, then we have it...
03439     if ( tan_angle > 165.0 && tan_angle < 195.0 )
03440     {
03441       tangential_surface_pairs.append(tangential_face);
03442       tangential_surface_pairs.append(other_face);
03443     }
03444   }
03445 }
03446     
03447       
03448     
03449 void GeomMeasureTool::find_interior_curve_angles( RefVolume *ref_volume,
03450                                                   double upper_bound,
03451                                                   double lower_bound,
03452                                                   DLIList <RefEdge*> &large_edge_angles,
03453                                                   DLIList <RefEdge*> &small_edge_angles,
03454                                                   DLIList <double> &large_angles,
03455                                                   DLIList <double> &small_angles,
03456                                                   int &total_interior,
03457                                                   int &total_fuzzy)
03458 {
03459   int ii, jj, kk;
03460   DLIList <RefFace*> ref_faces;
03461   RefFace *curr_face;
03462   total_interior = 0;
03463   total_fuzzy = 0;
03464   ref_volume->ref_faces(ref_faces);
03465   
03466     //Okay now loop over the edge loops of the reffaces
03467     //and find the interior angles.
03468   DLIList<Loop*> loops;
03469   Loop *curr_loop;
03470   DLIList <CoEdge*> co_edges;
03471   CoEdge *curr_edge, *next_edge;
03472   double angle = 0.0;
03473   const double RAD_TO_DEG = 180.0/CUBIT_PI;
03474   for ( ii = ref_faces.size(); ii > 0; ii-- )
03475   {
03476     curr_face = ref_faces.get_and_step(); 
03477     curr_face->loops(loops);
03478       //Now loop over these loops.
03479     for ( jj = loops.size(); jj > 0; jj-- )
03480     {
03481       curr_loop = loops.pop();
03482         //now loop over the edges in this
03483       co_edges.clean_out();
03484       curr_loop->ordered_co_edges(co_edges);
03485       for ( kk = co_edges.size(); kk > 0; kk-- )
03486       {
03487         curr_edge = co_edges.get_and_step();
03488         next_edge = co_edges.get();
03489        
03490         if ( curr_edge == next_edge )
03491         {
03492             //This is a hardpoint
03493           if ( curr_edge->get_ref_edge_ptr()->geometry_type() == POINT_CURVE_TYPE )
03494           {
03495             angle = 360.0;
03496               //for now ignore this..
03497             continue;
03498           }
03499           else
03500           {
03501             RefEdge *the_edge = curr_edge->get_ref_edge_ptr();
03502               //Get points 1% from the ends in both directions.
03503             CubitVector point_start, point_end;
03504             CubitStatus stat = the_edge->position_from_fraction(0.01,
03505                                                                 point_start);
03506             if ( stat != CUBIT_SUCCESS )
03507               angle = 180.0;
03508             else
03509             {
03510               stat = the_edge->position_from_fraction(0.99,
03511                                                       point_end);
03512               if ( stat != CUBIT_SUCCESS )
03513                 angle = 180.0;
03514               else
03515               {
03516                 CubitVector normal = curr_face->normal_at(point_start, NULL);
03517                 CubitVector tangent_1, tangent_2;
03518                 stat = the_edge->tangent( point_start, tangent_1, curr_face);
03519                 CubitStatus stat2 = the_edge->tangent( point_end, tangent_2, curr_face);
03520                 if ( stat != CUBIT_SUCCESS || stat2 != CUBIT_SUCCESS )
03521                   angle = 180.0;
03522                 else
03523                 {
03524                     //one of the tangents will be pointing the wrong
03525                     //direction, so reverse it based on the sense.
03526                   if ( curr_edge->get_sense() == CUBIT_REVERSED )
03527                     tangent_1 = -tangent_1;
03528                   else
03529                     tangent_2 = -tangent_2;
03530                   angle = normal.vector_angle(tangent_2, tangent_1);
03531                   angle *= RAD_TO_DEG;
03532                 }
03533               }
03534             }
03535           }
03536         }
03537         else if ( curr_edge->get_ref_edge_ptr() == next_edge->get_ref_edge_ptr() )
03538         {
03539             //This is the end of a hard line. Set this angle because
03540             //it is tricky to measure it.
03541           angle = 360.0;
03542         }
03543         else if ( curr_edge != next_edge )
03544         {
03545           angle =
03546             GeometryQueryTool::instance()->geometric_angle(curr_edge,
03547                                                            next_edge);
03548           angle *= RAD_TO_DEG;
03549         }
03550         total_interior++;
03551           //count the number of angles that are fuzzy.
03552         if ( angle <= GEOM_END_LOWER )
03553           total_fuzzy++;
03554         else if ( angle >= GEOM_END_UPPER &&
03555                   angle <= GEOM_SIDE_LOWER )
03556           total_fuzzy++;
03557         else if ( angle >= GEOM_SIDE_UPPER &&
03558                   angle <= GEOM_CORNER_LOWER )
03559           total_fuzzy++;
03560         else if ( angle >= GEOM_CORNER_UPPER )
03561           total_fuzzy++;
03562         if ( angle <= lower_bound )
03563         {
03564           small_edge_angles.append(curr_edge->get_ref_edge_ptr());
03565           small_edge_angles.append(next_edge->get_ref_edge_ptr());
03566           small_angles.append(angle);
03567         }
03568         if ( angle >= upper_bound )
03569         {
03570           large_edge_angles.append(curr_edge->get_ref_edge_ptr());
03571           large_edge_angles.append(next_edge->get_ref_edge_ptr());
03572           large_angles.append(angle);
03573         }
03574       }
03575     }
03576   }
03577   return;
03578 }
03579 
03580 void GeomMeasureTool::find_dihedral_angles( DLIList <RefVolume*> &ref_vols,
03581                                             double upper_bound,
03582                                             double lower_bound,
03583                                             DLIList <RefFace*> &large_face_angles,
03584                                             DLIList <RefFace*> &small_face_angles,
03585                                             DLIList <double> &large_angles,
03586                                             DLIList <double> &small_angles,
03587                                             int &total_interior,
03588                                             int &total_fuzzy,
03589                                             int &total_not_flat)
03590 {
03591   int ii;
03592   RefVolume *ref_vol;
03593   total_interior = 0;
03594   total_fuzzy = 0;
03595   total_not_flat = 0;
03596   for ( ii = 0; ii < ref_vols.size(); ii++ )
03597   {
03598     ref_vol = ref_vols.get_and_step();
03599     find_dihedral_angles(ref_vol, upper_bound,
03600                          lower_bound,
03601                          large_face_angles,
03602                          small_face_angles,
03603                          large_angles, small_angles,
03604                          total_interior, total_fuzzy, total_not_flat);
03605   }
03606   return;
03607 }
03608 void GeomMeasureTool::find_dihedral_angles(RefVolume *curr_volume,
03609                                            double upper_bound,
03610                                            double lower_bound,
03611                                            DLIList <RefFace*> &large_face_angles,
03612                                            DLIList <RefFace*> &small_face_angles,
03613                                            DLIList <double> &large_angles,
03614                                            DLIList <double> &small_angles,
03615                                            int &total_interior,
03616                                            int &total_fuzzy, int &total_not_flat)
03617 {
03618   int ii, jj, shared_angles = 0;
03619   double common_angle;
03620   DLIList <RefFace*> face_list;
03621   RefFace *curr_face_1, *curr_face_2;
03622   RefFace *tmp_face;
03623   RefEdge *common_edge;
03624   DLIList <RefVolume*> ref_vols;
03625   DLIList<RefEdge*> vol_edges;
03626   curr_volume->ref_edges(vol_edges);
03627 
03628   //maybe the body is a sheet-body...don't print errors if so
03629   bool is_sheet_body = false;
03630   if( curr_volume->measure() == 0 )
03631     is_sheet_body = true;
03632 
03633   for ( ii = vol_edges.size(); ii > 0; ii-- )
03634   {
03635     common_edge = vol_edges.get_and_step();
03636 
03637     //skip point curves
03638     if( common_edge->geometry_type() == POINT_CURVE_TYPE )
03639       continue;
03640 
03641     curr_face_1 = NULL;
03642     curr_face_2 = NULL;
03643     face_list.clean_out();
03644     common_edge->ref_faces(face_list);
03645       //find the two faces that are part of this volume.
03646     for ( jj = face_list.size(); jj > 0; jj-- )
03647     {
03648       tmp_face = face_list.get_and_step();
03649       ref_vols.clean_out();
03650       tmp_face->ref_volumes(ref_vols);
03651       if ( ref_vols.move_to(curr_volume) )
03652       {
03653         if ( curr_face_1 == NULL )
03654           curr_face_1 = tmp_face;
03655         else if ( curr_face_2 == NULL )
03656         {
03657           curr_face_2 = tmp_face;
03658           break;
03659         }
03660       }
03661     }
03662     if ( curr_face_1 == NULL ||
03663          curr_face_2 == NULL )
03664     {
03665       if( is_sheet_body == false )
03666         PRINT_ERROR("Problems finding connected surfaces to a curve\n"
03667                    "\tfor measuring dihedral angles.\n");
03668       continue;
03669     }
03670     shared_angles++;
03671     common_angle =  GeometryQueryTool::instance()->surface_angle(curr_face_1,
03672                                                                  curr_face_2,
03673                                                                  common_edge,
03674                                                                  curr_volume);
03675     common_angle *= 180.0/CUBIT_PI;
03676     total_interior++;
03677       //use hard coded angles because we want these a little looser.
03678     if ( common_angle >= 200.0 ||
03679          common_angle <= 160.0 )
03680       total_not_flat++;
03681     
03682     if ( curr_face_1->geometry_type() != PLANE_SURFACE_TYPE ||
03683          curr_face_2->geometry_type() != PLANE_SURFACE_TYPE )
03684       total_fuzzy++;
03685     else if ( common_angle <= GEOM_END_LOWER )
03686       total_fuzzy++;
03687     else if ( common_angle >= GEOM_END_UPPER &&
03688               common_angle <= GEOM_SIDE_LOWER )
03689       total_fuzzy++;
03690     else if ( common_angle >= GEOM_SIDE_UPPER &&
03691               common_angle <= GEOM_CORNER_LOWER )
03692       total_fuzzy++;
03693     else if ( common_angle >= GEOM_CORNER_UPPER )
03694       total_fuzzy++;
03695     
03696     if( common_angle <= lower_bound )
03697     {
03698       small_face_angles.append(curr_face_1);
03699       small_face_angles.append(curr_face_2);
03700       small_angles.append(common_angle);
03701     }
03702     if( common_angle >= upper_bound )
03703     {
03704       large_face_angles.append(curr_face_1);
03705       large_face_angles.append(curr_face_2);
03706       large_angles.append(common_angle);
03707     }
03708   }
03709     //clean up the edge marks.
03710   for ( ii = vol_edges.size(); ii > 0; ii-- )
03711     vol_edges.get_and_step()->marked(0);
03712 }
03713 void GeomMeasureTool::find_close_loops(DLIList <RefVolume*> &ref_vols,
03714                                        DLIList <RefEdge*> &close_edges,
03715                                        DLIList <RefFace*> &close_loop_faces,
03716                                        DLIList <double> &small_lengths,
03717                                        double tol)
03718 {
03719   int ii, jj;
03720   DLIList <RefFace*> ref_faces, temp_faces;
03721   RefVolume *ref_vol;
03722   RefFace *curr_face;
03723   for ( ii = 0; ii < ref_vols.size(); ii++ )
03724   {
03725     ref_vol = ref_vols.get_and_step();
03726     ref_vol->ref_faces(temp_faces);
03727       //uniquely add the faces.
03728     for ( jj = temp_faces.size(); jj > 0; jj-- )
03729     {
03730       curr_face = temp_faces.pop();
03731       ref_faces.append(curr_face);
03732     }
03733   }
03734   DLIList <RefEdge*> tmp_close_edges;
03735   DLIList <double> tmp_small_lengths;
03736   
03737   for ( ii = ref_faces.size(); ii > 0; ii-- )
03738   {
03739     curr_face = ref_faces.get_and_step();
03740     curr_face->marked(0);
03741     if ( curr_face->num_loops() < 2 )
03742       continue;
03743     tmp_close_edges.clean_out();
03744     tmp_small_lengths.clean_out();
03745     find_close_loops( curr_face, tmp_close_edges,
03746                       tmp_small_lengths, tol);
03747     if ( tmp_close_edges.size() > 0 )
03748     {
03749       assert(tmp_close_edges.size()%2 == 0 );
03750       tmp_close_edges.reset();
03751       tmp_small_lengths.reset();
03752       for ( jj = tmp_close_edges.size()/2; jj > 0; jj-- )
03753       {
03754         close_edges.append(tmp_close_edges.get_and_step());
03755         close_edges.append(tmp_close_edges.get_and_step());
03756         small_lengths.append(tmp_small_lengths.get_and_step());
03757         close_loop_faces.append(curr_face);
03758       }
03759     }
03760   }
03761   return;
03762 }
03763 
03764   
03765 void GeomMeasureTool::find_close_loops(RefFace *face,
03766                                        DLIList <RefEdge*> &close_edges,
03767                                        DLIList <double> &small_lengths,
03768                                        double tol)
03769 {
03770     //only do faces with multiple loops since we are measuring the
03771     //distances between these loops.
03772   if ( face->num_loops() < 2 )
03773     return;
03774     //This is the most tricky of these functions.
03775     //To do this I'm going to facet the edges, then use an AbstractTree to
03776     //find the minimum distance between the loops and between a single
03777     //loop.
03778   PointLoopList boundary_point_loops;
03779   CubitStatus stat;
03780   stat = get_boundary_points(face, boundary_point_loops, GEOMETRY_RESABS*100);
03781   if ( stat != CUBIT_SUCCESS )
03782     return;
03783   
03784   SegLoopList boundary_seg_loops;
03785   stat = convert_to_lines( boundary_point_loops,
03786                            boundary_seg_loops,
03787                            face);
03788   if ( stat != CUBIT_SUCCESS )
03789     return;
03790 
03791     //Now add the points to the AbstractTree.
03792   DLIList<AbstractTree<GeomSeg*>*> atree_list;
03793   AbstractTree<GeomSeg*> *curr_tree;
03794   SegList *seg_list;
03795   GeomSeg  *curr_seg;
03796   int ii,jj, kk;
03797   //const double angle_convert = 180.0/CUBIT_PI;
03798   boundary_seg_loops.reset();
03799   int num_segs = 0;
03800   for (ii = 0; ii < boundary_seg_loops.size(); ii++ )
03801   {
03802     seg_list = boundary_seg_loops.get_and_step();
03803     curr_tree = new RTree<GeomSeg*>(GEOMETRY_RESABS);
03804     for ( jj = seg_list->size(); jj > 0; jj-- )
03805     {
03806         //build the r-tree.
03807       num_segs++;
03808       curr_seg = seg_list->get_and_step();
03809       curr_tree->add(curr_seg);
03810         //calculate the interior angles.
03811     }
03812     atree_list.append(curr_tree);
03813   }
03814     //determine the minimum distance between the loops.
03815   PointList *curr_points;
03816   GeomPoint *curr_point;
03817   DLIList <GeomSeg*> nearest_neighbors;
03818   RefEdge *closest_edge_1, *closest_edge_2;
03819   CubitVector curr_loc;
03820   double closest_dist;
03821   double min_for_loop_squared;
03822   atree_list.reset();
03823   // This searching was orignally coded to not be an n-squared
03824   // search but I changed it to be n-squared because the
03825   // way the loops are compared comparing loop 1 against loop 2 may
03826   // give a different result than comparing loop 2 against loop 1.
03827   for ( ii = 0; ii < atree_list.size(); ii++ )
03828   {
03829     curr_points = boundary_point_loops.get_and_step();
03830     for ( jj = 1; jj < atree_list.size(); jj++ )
03831     {
03832       min_for_loop_squared = CUBIT_DBL_MAX;
03833       closest_edge_1 = NULL;
03834       closest_edge_2 = NULL;
03835       curr_tree = atree_list.next(ii+jj);
03836         //now for every point in curr_points, find it's closest_point
03837         //in the curr_tree.
03838       for ( kk = 0; kk < curr_points->size(); kk++ )
03839       {
03840         curr_point = curr_points->get_and_step();
03841         curr_loc.set(curr_point->coordinates());
03842         nearest_neighbors.clean_out();
03843         curr_tree->k_nearest_neighbor(curr_loc,
03844                                       1, closest_dist, nearest_neighbors,
03845                                       dist_sq_point_data);
03846         if ( nearest_neighbors.size() == 0)
03847         {
03848             //hmm, not sure what is wrong here.
03849           PRINT_ERROR("Can't find closest point between loops.\n");
03850         }
03851         if (closest_dist <= tol*tol && closest_dist < min_for_loop_squared )
03852         {
03853             //just store the smaller ones.  Don't store
03854           RefEntity *ref_ent = curr_point->owner();
03855           GeomSeg *near_seg = nearest_neighbors.get();
03856           RefEdge *ref_edge_1 = CAST_TO(ref_ent, RefEdge);
03857           if ( ref_edge_1 == NULL )
03858           {
03859               //assume this is a vertex.
03860             RefVertex *ref_vert = CAST_TO(ref_ent, RefVertex);
03861             if ( ref_vert == NULL )
03862             {
03863               PRINT_ERROR("problem with point ownership in closest loops.\n");
03864               continue;
03865             }
03866             DLIList <RefEdge*> ref_edges;
03867             ref_vert->ref_edges(ref_edges);
03868             int ll;
03869             RefEdge *the_edge = NULL;
03870             for ( ll = ref_edges.size(); ll > 0; ll-- )
03871             {
03872               if ( ref_edges.get()->is_directly_related(face) )
03873               {
03874                 the_edge = ref_edges.get();
03875                 break;
03876               }
03877               ref_edges.step();
03878             }
03879             if ( the_edge == NULL )
03880             {
03881               PRINT_ERROR("problem with point ownership in closest loops.\n");
03882               continue;
03883             }
03884             ref_edge_1 = the_edge;
03885           }
03886           ref_ent = near_seg->owner();
03887           RefEdge *ref_edge_2 = CAST_TO(ref_ent, RefEdge);
03888           if ( ref_edge_2 == NULL )
03889           {
03890             PRINT_ERROR("problem with point ownership in closest loops.\n");
03891             continue;
03892           }
03893           
03894           min_for_loop_squared = closest_dist;
03895           closest_edge_1 = ref_edge_1;
03896           closest_edge_2 = ref_edge_2;
03897         }
03898       }
03899       if ( closest_edge_1 != NULL )
03900       {
03901         close_edges.append(closest_edge_1);
03902         close_edges.append(closest_edge_2);
03903         small_lengths.append(sqrt(min_for_loop_squared));
03904       }
03905     }
03906   }
03907     //clean up the memory.
03908   int ll,mm;
03909   atree_list.reset();
03910   boundary_seg_loops.reset();
03911   for ( ll = boundary_seg_loops.size(); ll > 0; ll-- )
03912   {
03913     delete atree_list.pop();
03914     seg_list = boundary_seg_loops.pop();
03915     for ( mm = seg_list->size(); mm > 0; mm-- )
03916       delete seg_list->pop();
03917     delete seg_list;
03918   }
03919         
03920   PointList *point_list;
03921   for ( ll = boundary_point_loops.size(); ll > 0; ll-- )
03922   {
03923     point_list = boundary_point_loops.pop();
03924     for ( mm = point_list->size(); mm > 0; mm-- )
03925     {
03926       GeomPoint *tmp_point = point_list->pop();
03927         //delete the CubitPointData
03928       delete tmp_point;
03929     }
03930     delete point_list;
03931   }
03932   return;
03933 }
03934 double GeomMeasureTool::measure_area(RefFace* curr_face)
03935 {
03936   double area;
03937   GeomDataObserver *g_d_o = GeomDataObserver::get(curr_face);
03938   if ( g_d_o == NULL )
03939     g_d_o = GeomDataObserver::create(curr_face);
03940 
03941   if ( g_d_o == NULL )
03942   {
03943     assert(g_d_o != NULL );
03944     return curr_face->measure();
03945   }
03946   if ( !g_d_o->is_measure_set() )
03947   {
03948     area  = (curr_face->get_surface_ptr())->measure();
03949     g_d_o->set_measure(area);
03950   }
03951   area = g_d_o->get_measure();
03952   return area;
03953 }
03954 
03960 void GeomMeasureTool::find_bad_geometry(RefVolume *volume,
03961                                         DLIList <RefEntity*> &bad_ents)
03962 {
03963   Lump *lump_ptr = volume->get_lump_ptr();
03964   DLIList <TopologyEntity*> bad_top_ents;
03965   lump_ptr->validate(volume->entity_name(), bad_top_ents);
03966     //Now go through and get a RefEnttiy for each entity.
03967   RefEntity *ref_ent;
03968   TopologyEntity *topo_ent;
03969   GroupingEntity *gr_ent;
03970   SenseEntity *se_ent;
03971   
03972   int ii;
03973   for (ii = bad_top_ents.size(); ii > 0; ii-- )
03974   {
03975     topo_ent = bad_top_ents.pop();
03976     ref_ent = CAST_TO(topo_ent, RefEntity);
03977     if ( ref_ent != NULL )
03978     {
03979       bad_ents.append(ref_ent);
03980       continue;
03981     }
03982     gr_ent = CAST_TO(topo_ent, GroupingEntity);
03983     if ( gr_ent != NULL )
03984     {
03985       Chain *ch = CAST_TO(gr_ent, Chain);
03986       if ( ch != NULL )
03987       {
03988         RefVertex *start_v = ch->start_vertex();
03989         RefVertex *end_v = ch->end_vertex();
03990         RefEdge *ref_edge = start_v->common_ref_edge(end_v);
03991         bad_ents.append(ref_edge);
03992         continue;
03993       }
03994       Loop *lo = CAST_TO(gr_ent, Loop);
03995       if ( lo != NULL )
03996       {
03997         RefFace *rf = lo->get_ref_face_ptr();
03998         bad_ents.append(rf);
03999         continue;
04000       }
04001       Shell *sh = CAST_TO(gr_ent, Shell);
04002       if (sh != NULL )
04003       {
04004         RefVolume *rv = sh->get_ref_volume_ptr();
04005         bad_ents.append(rv);
04006         continue;
04007       }
04008       else
04009       {
04010           //this is incase there is some other entity we don't
04011           //know about, just get the body...
04012         Body *body = topo_ent->body();
04013         bad_ents.append(body);
04014         continue;
04015       }
04016     }
04017     se_ent = CAST_TO(topo_ent, SenseEntity);
04018     if ( se_ent != NULL )
04019     {
04020       BasicTopologyEntity *be = se_ent->get_basic_topology_entity_ptr();
04021       bad_ents.append(be);
04022       continue;
04023     }
04024       //If we are here we didn't get the right entity.  Just get
04025       //the body associated with this topology entity.
04026     Body *b = topo_ent->body();
04027     bad_ents.append(b);
04028   }
04029   return;
04030 }
04031 
04037 void GeomMeasureTool::find_irregular_valence( DLIList <RefVolume*> &ref_volumes,
04038                                               DLIList <RefVertex*> &irregular_vertices)
04039 {
04040   RefVolume *ref_vol;
04041   int ii;
04042   for ( ii = 0; ii < ref_volumes.size(); ii++ )
04043   {
04044     ref_vol = ref_volumes.get_and_step();
04045     find_irregular_valence(ref_vol, irregular_vertices);
04046   }
04047 }
04053 void GeomMeasureTool::find_irregular_valence( RefVolume *ref_volume,
04054                                               DLIList <RefVertex*> &irregular_vertices)
04055 {
04056   RefVertex *ref_vert;
04057   DLIList <RefVertex*> ref_vertices;
04058   ref_volume->ref_vertices(ref_vertices);
04059   int ii;
04060   for ( ii = ref_vertices.size(); ii > 0; ii-- )
04061   {
04062     ref_vert = ref_vertices.get_and_step();
04063     if ( ref_vert->num_ref_edges() > 4 )
04064       irregular_vertices.append(ref_vert);
04065   }
04066 }
04067 
04071 #define FACE_BLEND      1
04072 #define VERTEX_BLEND    2
04073 #define BLEND_PROCESSED 3
04074 
04075 void GeomMeasureTool::find_blends( RefVolume *ref_volume,
04076                                    DLIList <RefFace*> &blend_faces,
04077                                    DLIList <DLIList<RefFace*>*> &blend_groups )
04078 {
04079     //Assume for now that blends have cartesian angles between
04080     //all the attached faces and itself.
04081     //Also a blend cannot be planar.  And directly accross from
04082     //one blend edge, there is usually another blended edge but
04083     //the angle between the two normals must be orthogonal. In other
04084     //words their must be some sort of transition.
04085   
04086     //mark all the faces as 0.  The mark is used here to find the
04087     // blend groups.  The function is_vertex_blend() marks faces
04088     // as being a vertex blend.  The group of faces is then traversed
04089     // until the adjacent surface (across the cross curve) is either
04090     // a vertex blend or not a blend.
04091   DLIList <RefFace*> ref_faces;
04092   DLIList <RefFace*> vertex_blends, face_blends;
04093   ref_volume->ref_faces(ref_faces);
04094   int ii, jj;
04095   for ( ii =0; ii < ref_faces.size(); ii++ )
04096     ref_faces.get_and_step()->marked(0);
04097 
04098     //First go through each face and each edge on each face.
04099   DLIList <RefEdge*> ref_edges;
04100   RefFace *ref_face;
04101   RefEdge *ref_edge, *other_edge;
04102   int num_vertex_blends = 0;
04103   int num_face_blends = 0;
04104   for ( ii = ref_faces.size(); ii > 0; ii-- )
04105   {
04106     ref_face = ref_faces.get_and_step();
04107       //Test the face to see if it is a face or vertex blend.
04108     other_edge = NULL;
04109     ref_edge = NULL;
04110     if ( is_face_blend(ref_face, ref_volume,
04111                        ref_edge, other_edge ) )
04112     {
04113       face_blends.append(ref_face);
04114       ref_face->marked(1);
04115       num_face_blends++;
04116     }
04117     else if ( is_vertex_blend(ref_face, ref_volume) )
04118     {
04119       vertex_blends.append(ref_face);
04120       ref_face->marked(2);
04121       num_vertex_blends++;  // keep track of the number of vertex blends
04122     }
04123   }
04124   if ( face_blends.size() + vertex_blends.size() == 0 )
04125   {
04126     return;
04127   }
04128     //Find out how many different groups of surfaces there
04129     //are that share curves.
04130   DLIList <RefFace*> *blend_group = NULL;
04131   DLIList <RefFace*> stack;
04132   RefFace *start_face = NULL, *other_face;
04133   
04134   // continue while there are blends to process
04135   while ( vertex_blends.size() + face_blends.size() > 0 || 
04136           start_face != NULL)
04137   {
04138     // if we don't have a start_face, get a new one.
04139     if (!start_face) 
04140     {
04141       // this is the start of a new group
04142       blend_group = new DLIList <RefFace*>;
04143       blend_groups.append(blend_group);
04144 
04145       // always prefer to start at a vertex blend if one exists
04146       if (vertex_blends.size() > 0 )
04147       {
04148         start_face = vertex_blends.pop();
04149       }
04150       else 
04151       {
04152         start_face = face_blends.pop();
04153       }
04154     }
04155 
04156     blend_group->append(start_face);  // add ref_face to group
04157     blend_faces.append(start_face);   // add the ref_face to the returned list
04158     start_face->marked(BLEND_PROCESSED);
04159 
04160     ref_edges.clean_out();
04161     start_face->ref_edges(ref_edges);
04162     for ( jj = 0; jj < ref_edges.size(); jj++ )
04163     {
04164       ref_edge = ref_edges.get_and_step();
04165       other_face = ref_edge->other_face(start_face, ref_volume);
04166       if (other_face == NULL)
04167       {
04168         start_face = NULL;
04169         break;
04170       }
04171 
04172       // if this blend has been processed and isn't in the current group
04173       if (other_face->marked() == BLEND_PROCESSED && 
04174           !blend_group->is_in_list(other_face) )
04175       {
04176         start_face = NULL;
04177         break;  // reached an end of this loop
04178       }
04179       else if (other_face->marked() == VERTEX_BLEND)
04180       {
04181         blend_group->append(other_face);  // add the ref_face to the group
04182         blend_faces.append(other_face);   // add the ref_face to the returned list
04183         vertex_blends.remove(other_face);
04184         other_face->marked(BLEND_PROCESSED);
04185         start_face = NULL;
04186         break;  // a vertex blend is also end of the line
04187       }
04188       else if (other_face->marked() == FACE_BLEND)
04189       {
04190         start_face = other_face;
04191         face_blends.remove(other_face);
04192         break;  // continue using this face as the new start face
04193       }
04194     }
04195     // if we traversed through all the edges of this blend without
04196     // finding another blend attached, this is the end of the chain 
04197     // and we need to find another starting place.
04198     if ( jj >= ref_edges.size() ) 
04199     {
04200       start_face = NULL;
04201     }
04202     
04203   }
04204 
04205     //cleanup marks.
04206   for ( ii =0; ii < ref_faces.size(); ii++ )
04207     ref_faces.get_and_step()->marked(0);
04208   for ( ii =0; ii < blend_faces.size(); ii++ )
04209     blend_faces.get_and_step()->marked(0);
04210     //Okay we should have it now.
04211 
04212     //REMEMBER TO DELETE the dllists in blend_faces.
04213   return;
04214 }
04215 
04216 // struct type object private to the next function
04217 class NextBlend
04218 {
04219   public:
04220   RefEdge* ref_edge;
04221   RefFace* ref_face;
04222 };
04223 
04224 void GeomMeasureTool::find_blends_from_edge( RefVolume* ref_volume, 
04225                                                     RefFace *start_face, 
04226                                                     RefEdge* start_edge,
04227                                                     std::vector <RefFace*> &blend_faces)
04228 {
04229   // initialize the mark on all surfaces in the volume
04230   DLIList <RefFace*> ref_faces;
04231   RefEdge *next_edge, *other_edge;
04232   ref_volume->ref_faces(ref_faces);
04233   int ii;
04234   for ( ii =0; ii < ref_faces.size(); ii++ )
04235     ref_faces.get_and_step()->marked(0);
04236 
04237   std::stack<NextBlend> blend_stack;
04238 
04239   NextBlend blend;
04240   blend.ref_edge = start_edge;
04241   blend.ref_face = start_face;
04242   blend_faces.push_back(start_face);
04243 
04244   blend_stack.push(blend);
04245 
04246   while (blend_stack.size() > 0)
04247   {
04248     // this is an oddity with std::stack.  Get the top of the stack
04249     // and then you have to remove it from the stack with pop
04250     NextBlend next_blend = blend_stack.top();
04251     blend_stack.pop();
04252     RefEdge* ref_edge = next_blend.ref_edge;
04253     RefFace* ref_face = next_blend.ref_face;
04254 
04255     RefFace* next_face;
04256     next_face = ref_edge->other_face(ref_face, ref_volume);
04257 
04258     // the the next face exists and it's a face blend, save the current blend,
04259     // create a new blend object and add it to the stack
04260     if ( next_face && is_face_blend(next_face, ref_volume, next_edge, other_edge ) )
04261     {
04262       // if the next face is already processed, the chain loops back on itself.
04263       if (next_face->marked() == BLEND_PROCESSED)
04264       {
04265         for ( ii =0; ii < ref_faces.size(); ii++ )
04266           ref_faces.get_and_step()->marked(0);
04267         return;
04268       }
04269 
04270       blend_faces.push_back(next_face);
04271       next_face->marked(BLEND_PROCESSED);
04272 
04273       // the is_face_blend function assumes (poorly) rectangular surfaces
04274       // without any holes.  It returns the spring curves and we want
04275       // the cross curves. So get all four edges and remove the two
04276       // spring curves returned by is_face_blend
04277       DLIList <RefEdge*> ref_edges;
04278       next_face->ref_edges(ref_edges);
04279       ref_edges.remove(next_edge);
04280       ref_edges.remove(other_edge);
04281 
04282       if (ref_edges.get() != ref_edge)
04283       {
04284         next_blend.ref_edge = ref_edges.get();
04285       }
04286       else
04287       {
04288         next_blend.ref_edge = ref_edges.next(1);
04289       }
04290       next_blend.ref_face = next_face;
04291       blend_stack.push(next_blend);
04292     }
04293     else if ( next_face && is_vertex_blend(next_face, ref_volume) )
04294     {
04295       // we will stop the chain at a vertex blend
04296       blend_faces.push_back(next_face);
04297       next_face->marked(BLEND_PROCESSED);
04298     }
04299   }
04300 
04301   // clean up the marks
04302   std::vector<RefFace*>::iterator iter;
04303   for ( iter = blend_faces.begin(); iter != blend_faces.end(); iter++)
04304     (*iter)->marked(0);
04305 }
04306 
04307 // given a starting blend surface find a chain of blends from
04308 // that surface.  
04309 //
04310 // Note that this function intentionally does _not_
04311 // clear the blend_face list so that additional chains can be added.
04312 CubitStatus GeomMeasureTool::find_blend_chains( RefFace *start_face,
04313                                 std::vector<std::vector< RefFace*> > &blend_chains)
04314 {
04315 
04316   if (start_face == NULL)
04317   {
04318     return CUBIT_FAILURE;
04319   }
04320 
04321   std::vector <RefFace*> blend_faces;
04322 
04323   // get the owning volume of this blend
04324   DLIList<RefEntity*> entity_list;
04325   RefVolume* ref_volume;
04326   int ii;
04327 
04328   start_face->get_parent_ref_entities(entity_list);
04329 
04330   // this indicates merged enitites and potential problems
04331   if (entity_list.size() > 1)
04332   {
04333     return CUBIT_FAILURE;
04334   }
04335 
04336   // make sure we're at the beginning of the list and get the first
04337   // and only item on the list and cast it to a RefVolume
04338   entity_list.reset();
04339   ref_volume = CAST_TO(entity_list.get(), RefVolume);
04340  
04341   if (!ref_volume)
04342   {
04343     return CUBIT_FAILURE;
04344   }
04345 
04346   // initialize the mark on all surfaces in the volume
04347   DLIList <RefFace*> ref_faces;
04348   ref_volume->ref_faces(ref_faces);
04349 
04350   RefEdge *spring_curve1, *spring_curve2;
04351   if ( is_face_blend(start_face, ref_volume, spring_curve1, spring_curve2 ) )
04352   {
04353     // the is_face_blend function assumes (poorly) rectangular surfaces
04354     // without any holes.  It returns the spring curves and we want
04355     // the cross curves. So get all four edges and remove the two
04356     // spring curves returned by is_face_blend
04357     DLIList <RefEdge*> ref_edges;
04358     start_face->ref_edges(ref_edges);
04359     ref_edges.remove(spring_curve1);
04360     ref_edges.remove(spring_curve2);
04361 
04362     // there is a special case where the blend is a periodic surface
04363     // meaning that there _are_ no cross curves
04364     if ( ref_edges.size() == 0 )
04365     {
04366       blend_faces.push_back(start_face);
04367     }
04368     else
04369     {
04370       // now find additional blends extending from either side of the blend
04371       blend_faces.clear();
04372       find_blends_from_edge( ref_volume, start_face, ref_edges.get(), blend_faces);
04373       find_blends_from_edge( ref_volume, start_face, ref_edges.next(1), blend_faces);
04374 
04375       // make sure that we have a unique list (the start surface is probably here twice)
04376       std::vector<RefFace*>::iterator new_end;
04377       std::sort( blend_faces.begin(), blend_faces.end() );
04378       new_end = std::unique( blend_faces.begin(), blend_faces.end() );
04379       blend_faces.erase(new_end, blend_faces.end());
04380     }
04381 
04382     blend_chains.push_back(blend_faces);
04383   }
04384   else if ( is_vertex_blend(start_face, ref_volume) )
04385   {
04386     DLIList<RefEdge*> ref_edges;
04387     start_face->ref_edges(ref_edges);
04388 
04389     for (ii = 0; ii < ref_edges.size(); ii++)
04390     {
04391       RefEdge* start_edge = ref_edges.get_and_step();
04392       blend_faces.clear();
04393       find_blends_from_edge( ref_volume, start_face, start_edge, blend_faces);
04394 
04395       blend_chains.push_back(blend_faces);
04396     }
04397   }
04398   else
04399   {
04400     // the given face is not a blend
04401     return CUBIT_FAILURE;
04402   }
04403 
04404   return CUBIT_SUCCESS;
04405 }
04406 
04407 //--------------------------------------------------------------------
04408 //Function: Public, is_face_blend
04409 //Description: Determines if a face is a blend surface, returns the
04410 //   ref_edge on one side of the blend and other_edge on the opposite side.
04411 //   For this type of blend, only ref_edge must be tangentially meeting
04412 //   with another surface.  Other_edge must be oriented orthogonally to
04413 //   it and may or may not blend with another surface.  This assumes
04414 //   a rectangular blend surface, without holes.
04415 //---------------------------------------------------------------------
04416 CubitBoolean GeomMeasureTool::is_face_blend(RefFace *ref_face,
04417                                             RefVolume *ref_volume,
04418                                             RefEdge *&ref_edge,
04419                                             RefEdge *&other_edge)
04420 {
04421     //first we know that blend surfaces are not planar.
04422     //so remove these first.
04423     //Also, don't look at faces with more than 2 loops.
04424   if ( ref_face->geometry_type() == PLANE_SURFACE_TYPE ||
04425        ref_face->num_loops() > 2 )
04426     return CUBIT_FALSE;
04427   
04428   CubitBoolean is_cartesian;
04429   DLIList<RefEdge*> ref_edges;
04430   ref_edge = NULL;
04431   other_edge = NULL;
04432   RefFace *other_face;
04433   int jj;
04434   double angle;
04435   ref_face->ref_edges(ref_edges);
04436   for ( jj = ref_edges.size(); jj > 0; jj-- )
04437   {
04438     ref_edge = ref_edges.get_and_step();
04439 
04440     //Weed-out case where one edge is shared between more 
04441     //than 2 surfaces of the same volume
04442     DLIList<RefFace*> tmp_faces;
04443     ref_edge->ref_faces( tmp_faces );
04444 
04445     if( tmp_faces.size() > 2 )
04446     {
04447       int kk;
04448       for(kk=tmp_faces.size(); kk--;)
04449       {
04450         if( !tmp_faces.get()->is_child( ref_volume ) )
04451           tmp_faces.change_to(NULL);
04452         tmp_faces.step();
04453       }
04454       tmp_faces.remove_all_with_value( NULL );
04455       if( tmp_faces.size() > 2 )
04456         //this isn't the type of surface we are looking for...
04457         continue;
04458     }
04459 
04460     other_face = ref_edge->other_face(ref_face, ref_volume);
04461     if ( other_face == NULL )
04462     {
04463         //this isn't the type of surface we are looking for...
04464       break;
04465     }
04466     angle = GeometryQueryTool::instance()->surface_angle(ref_face,
04467                                                          other_face,
04468                                                          ref_edge,
04469                                                          ref_volume);
04470     angle *= 180.0/CUBIT_PI;
04471     is_cartesian = CUBIT_TRUE;
04472     if ( angle <= GEOM_SIDE_LOWER ||
04473          angle >= GEOM_SIDE_UPPER )
04474       is_cartesian = CUBIT_FALSE;
04475       //Okay, we have one major criteria achieved, this edge is a cartesian meet.
04476 
04477     if ( !is_cartesian )
04478       continue;
04479       //Now we need to check the radius of curvatures between these
04480       //two surfaces. I'm not totally sure about this but I think we
04481       //don't want the same radius of curvature.
04482     double k1_s1, k2_s1, k1_s2, k2_s2;
04483     CubitVector mid_point;
04484     ref_edge->mid_point(mid_point);
04485     ref_face->get_principal_curvatures( mid_point, k1_s1, k2_s1, ref_volume);
04486     other_face->get_principal_curvatures( mid_point, k1_s2, k2_s2, ref_volume);
04487     if (( is_equal(k1_s1, k1_s2) || is_equal(k1_s1, k2_s2) ) &&
04488         ( is_equal(k2_s1, k1_s2) || is_equal(k2_s1, k2_s2) ) )
04489         //try a different edge.
04490       continue;
04491 
04492       //Okay, this looks pretty good.
04493       //Now try and find an edge that is on the "opposite" side of
04494       // this surface. Sort of assume this is like a mapable surface.
04495       //check the oposite side.
04496     other_edge = NULL;
04497     CubitVector opp_point;
04498     if ( !find_opposite_edge( ref_edge, ref_face, other_edge, opp_point) ||
04499          other_edge == NULL )
04500         //if we can't find an opposite edge, we need to not
04501         //continue checking this surface.
04502       return CUBIT_FALSE;
04503     
04504       //Okay, we have the point opposite.  If this is a blend surface
04505       //the normal of this point should be pointing 90 degrees from
04506       //us.
04507     CubitVector normal_1 = ref_face->normal_at(mid_point, ref_volume);
04508     CubitVector normal_2 = ref_face->normal_at(opp_point, ref_volume);
04509     normal_1.normalize();
04510     normal_2.normalize();
04511     double dot_product = normal_1 % normal_2;
04512     if ( dot_product >= .5 || dot_product <= -.5  )
04513         //try a different edge.
04514       continue;
04515     CubitVector tangent_1, tangent_2;
04516     ref_edge->tangent(mid_point, tangent_1, ref_face);
04517     other_edge->tangent(opp_point, tangent_2, ref_face);
04518     tangent_1.normalize();      
04519     tangent_2.normalize();      
04520     dot_product = tangent_1 % tangent_2;
04521       //basically these edges should be somewhat parallel at
04522       //this point.
04523     if ( dot_product < .5 && dot_product > -.5 )
04524       continue;
04525 //This was called but never checked, not sure why.
04526 //      //check the curvature opposite.
04527 //    double k1_s3, k2_s3;
04528 //    ref_face->get_principal_curvatures( opp_point, k1_s3, k2_s3, ref_volume);
04529       //we are done with this face. It is a blend face.
04530     return CUBIT_TRUE;
04531   }
04532   return CUBIT_FALSE;
04533 }
04534 
04535 //--------------------------------------------------------------------
04536 //Function: Public, is_vertex_blend
04537 //Description: Determines if a face is a vertex blend surface.
04538 //   For this type of blend, all ref_edges must be meet tangentially
04539 //   with another surface.  This assumes blend surface with no holes. 
04540 //---------------------------------------------------------------------
04541 CubitBoolean GeomMeasureTool::is_vertex_blend(RefFace *ref_face,
04542                                               RefVolume* ref_volume)
04543 {
04544     //first we know that blend surfaces are not planar.
04545     //so remove these first.
04546     //Also, don't look at faces with more than 2 loops.
04547   if ( ref_face->geometry_type() == PLANE_SURFACE_TYPE ||
04548        ref_face->num_loops() > 2 )
04549     return CUBIT_FALSE;
04550   
04551   CubitBoolean is_cartesian;
04552   DLIList<RefEdge*> ref_edges;
04553   RefFace *other_face;
04554   RefEdge *ref_edge = NULL;
04555   int jj;
04556   double angle;
04557   ref_face->ref_edges(ref_edges);
04558   for ( jj = ref_edges.size(); jj > 0; jj-- )
04559   {
04560     ref_edge = ref_edges.get_and_step();
04561 
04562     //Weed-out case where one edge is shared between more 
04563     //than 2 surfaces of the same volume
04564     DLIList<RefFace*> tmp_faces;
04565     ref_edge->ref_faces( tmp_faces );
04566 
04567     if( tmp_faces.size() > 2 )
04568     {
04569       int kk;
04570       for(kk=tmp_faces.size(); kk--;)
04571       {
04572         if( !tmp_faces.get()->is_child( ref_volume ) )
04573           tmp_faces.change_to(NULL);
04574         tmp_faces.step();
04575       }
04576       tmp_faces.remove_all_with_value( NULL );
04577       if( tmp_faces.size() > 2 )
04578         //this isn't the type of surface we are looking for...
04579         continue;
04580     }
04581 
04582     other_face = ref_edge->other_face(ref_face, ref_volume);
04583     if ( other_face == NULL )
04584     {
04585         //this isn't the type of surface we are looking for...
04586       break;
04587     }
04588     angle = GeometryQueryTool::instance()->surface_angle(ref_face,
04589                                                          other_face,
04590                                                          ref_edge,
04591                                                          ref_volume);
04592     angle *= 180.0/CUBIT_PI;
04593     is_cartesian = CUBIT_TRUE;
04594     if ( angle <= GEOM_SIDE_LOWER ||
04595          angle >= GEOM_SIDE_UPPER )
04596       is_cartesian = CUBIT_FALSE;
04597       //Okay, we have one major criteria achieved, this edge is a cartesian meet.
04598 
04599     if ( !is_cartesian )
04600       return CUBIT_FALSE;
04601       //Now we need to check the radius of curvatures between these
04602       //two surfaces. I'm not totally sure about this but I think we
04603       // want the same radius of curvature.
04604     double k1_s1, k2_s1, k1_s2, k2_s2;
04605     CubitVector mid_point;
04606     ref_edge->mid_point(mid_point);
04607     ref_face->get_principal_curvatures( mid_point, k1_s1, k2_s1, ref_volume);
04608     other_face->get_principal_curvatures( mid_point, k1_s2, k2_s2, ref_volume);
04609     if (( is_equal(k1_s1, k1_s2) || is_equal(k1_s1, k2_s2) ) &&
04610         ( is_equal(k2_s1, k1_s2) || is_equal(k2_s1, k2_s2) ) )
04611         //try a different edge.
04612         continue;
04613     else
04614       return CUBIT_FALSE;
04615   }
04616 
04617   // if all edges are tangent and share curvatures then it must be a 
04618   // vertex blend.
04619   return CUBIT_TRUE;
04620 }
04621 
04622 CubitBoolean GeomMeasureTool::find_opposite_edge( RefEdge* ref_edge,
04623                                                   RefFace* ref_face,
04624                                                   RefEdge *&other_edge,
04625                                                   CubitVector &closest_point)
04626 {
04627     //Here we want to find the edge that is "opposite" to ref_edge.
04628     //Opposite is defined by assuming the surface is a square or cylinder and
04629     //finding the edge on the other side, so that there is a point closest
04630     //to the mid point on the ref_edge.
04631   //double need_to_turn = CUBIT_PI;
04632   int ii;
04633   DLIList <Loop*> loops;
04634   DLIList <CoEdge*> co_edges;
04635   Loop *curr_loop;
04636   CoEdge *co_edge, *this_co_edge;
04637   ref_edge->get_co_edges(co_edges, ref_face);
04638   if ( co_edges.size() != 1 )
04639     return CUBIT_FALSE;
04640   this_co_edge = co_edges.get();
04641   co_edges.clean_out();
04642   ref_face->loops(loops);
04643   if ( loops.size() > 2 )
04644     return CUBIT_FALSE;
04645   CubitBoolean found_loop = CUBIT_FALSE;
04646     //find the loop where this_co_edge is on.
04647   for(ii = loops.size(); ii > 0; ii-- )
04648   {
04649     curr_loop = loops.get_and_step();
04650     co_edges.clean_out();
04651     curr_loop->ordered_co_edges(co_edges);
04652     if ( co_edges.move_to(this_co_edge) )
04653     {
04654       found_loop = CUBIT_TRUE;
04655       break;
04656     }
04657   }
04658   if ( !found_loop )
04659     return CUBIT_FALSE;
04660   CoEdge *opp_co_edge = NULL;
04661   if ( loops.size() == 1 )
04662   {
04663       //This is the normal case, I hope.
04664       //this co_edge, go around the interior until
04665       //we turn twice.  If there are angles greater than CUBIT_PI, this isn't
04666       //a traditional blend.  Everything should be convex.
04667     if ( co_edges.size() == 4 )
04668     {
04669         //Just get the one opposite...
04670       opp_co_edge = co_edges.next(2);
04671     }
04672     else
04673     {
04674       //double total_angle = 0.0;
04675       int turn_counter = 0;
04676       CoEdge *next_co_edge = NULL;
04677       for ( ii = co_edges.size(); ii > 0; ii-- )
04678       {
04679         co_edge = co_edges.get_and_step();
04680         next_co_edge = co_edges.get();
04681         double angle = GeometryQueryTool::instance()->geometric_angle(co_edge, next_co_edge);
04682         angle *= 180.0/CUBIT_PI;
04683           //we don't deal with small or concave angles.
04684         if ( angle < 10.0 )
04685           return CUBIT_FALSE;
04686         else if ( angle > 185.0 )
04687           return CUBIT_FALSE;
04688         else if ( angle < 140.0 )
04689           turn_counter++;
04690         if ( turn_counter == 2 )
04691         {
04692           opp_co_edge = next_co_edge;
04693           break;
04694         }
04695       }
04696     }
04697   }
04698   else
04699   {
04700       //Just get the other loop and get one coedge in it.
04701     co_edges.clean_out();
04702     curr_loop = loops.get();
04703     curr_loop->ordered_co_edges(co_edges);
04704     opp_co_edge = co_edges.get();
04705   }
04706   if ( opp_co_edge == NULL )
04707     return CUBIT_FALSE;
04708     //This is just here to test correctness to this point, it can be taken out
04709     //if need be.
04710   if ( !co_edges.move_to(opp_co_edge) )
04711   {
04712     PRINT_ERROR("logic messed-up in find_opposite edge\n");
04713     assert(0);
04714     return CUBIT_FALSE;
04715   }
04716     //okay, we now have the co_edges in the correct position and
04717     //we have a canidate co_edge.
04718     //Search from the co_edge with the point that has closest point
04719     //to the mid point of this_co_edge.
04720   other_edge = NULL;
04721   CubitVector mid_point;
04722   ref_edge->mid_point(mid_point);
04723   CubitVector tmp_closest_point;
04724   double min_dist_sq = CUBIT_DBL_MAX;
04725   double dist_sq;
04726   RefEdge *test_ref_edge;
04727   CoEdge *prev_co_edge;
04728   CubitBoolean first = CUBIT_TRUE;
04729   for ( ii = co_edges.size(); ii > 0; ii-- )
04730   {
04731     co_edge = co_edges.get_and_step();
04732     prev_co_edge = co_edges.prev(2);
04733     test_ref_edge = co_edge->get_ref_edge_ptr();
04734       //We don't want edges that are connected to ref edge.
04735       //If we have them, we have gone to far around the surface.
04736     if ( test_ref_edge->common_ref_vertex(ref_edge) )
04737       break;
04738     if ( !first )
04739     {
04740         //break out of the for loop if we turn again. We just want the closest
04741         //point opposite.
04742       double angle = GeometryQueryTool::instance()->
04743         geometric_angle(prev_co_edge, co_edge);
04744       angle *= 180.0/CUBIT_PI;
04745       if ( angle > 195.0 || angle < 145.0 )
04746         break;
04747     }
04748     else
04749       first = CUBIT_FALSE;
04750       //now find the closest point to the mid point of ref_edge.
04751     test_ref_edge->closest_point_trimmed( mid_point,
04752                                           tmp_closest_point);
04753     dist_sq = (mid_point - closest_point).length_squared();
04754       //find and keep the closest one.
04755     if ( dist_sq < min_dist_sq )
04756     {
04757       other_edge = test_ref_edge;
04758       min_dist_sq = dist_sq;
04759       closest_point = tmp_closest_point;
04760     }
04761   }
04762   if ( other_edge == NULL )
04763     return CUBIT_FALSE;
04764 
04765   return CUBIT_TRUE;
04766 }
04767 CubitBoolean GeomMeasureTool::is_equal(double v1, double v2)
04768 {
04769   double smallv = v1-v2;
04770   if ( smallv < 0.0 )
04771   {
04772     if ( smallv >= -CUBIT_RESABS )
04773       return CUBIT_TRUE;
04774     else
04775       return CUBIT_FALSE;
04776   }
04777   else
04778   {
04779     if ( smallv <= CUBIT_RESABS )
04780       return CUBIT_TRUE;
04781     else
04782       return CUBIT_FALSE;
04783   }
04784 }
04785 CubitStatus GeomMeasureTool::get_centroid( RefFace *ref_face, CubitVector &centroid, double &tot_area )
04786 {
04787   GMem g_mem;
04788   unsigned short norm_tol = 5;
04789   double dist_tol = -1.0;
04790 
04791   ref_face->get_geometry_query_engine()->
04792       get_graphics(ref_face->get_surface_ptr(), &g_mem, norm_tol, dist_tol );
04793 
04794   if(g_mem.fListCount < 1)
04795   {
04796       // Decrease tolerance and try again (we can get this for small features)
04797       norm_tol /= 2;
04798       ref_face->get_geometry_query_engine()->
04799           get_graphics(ref_face->get_surface_ptr(), &g_mem, norm_tol, dist_tol );
04800   }
04801 
04802   if(g_mem.fListCount < 1)
04803   {
04804       // Lets give up 
04805       PRINT_ERROR( "Unable to find the center of a surface\n" );
04806       return CUBIT_FAILURE;
04807   }
04808 
04809   // Initialize
04810   tot_area = 0.0;
04811   centroid.set( 0.0, 0.0, 0.0 );
04812 
04813   // Loop through the triangles
04814   double tri_area;
04815   GPoint p[3];
04816   GPoint* plist = g_mem.point_list();
04817   int* facet_list = g_mem.facet_list();
04818   int c = 0;
04819   for( ; c < g_mem.fListCount; )
04820   {
04821       p[0] = plist[facet_list[++c]];
04822       p[2] = plist[facet_list[++c]];
04823       p[1] = plist[facet_list[++c]]; 
04824       c++;
04825 
04826       // Get centroid
04827       CubitVector p1( p[0].x, p[0].y, p[0].z );
04828       CubitVector p2( p[2].x, p[2].y, p[2].z );
04829       CubitVector p3( p[1].x, p[1].y, p[1].z );
04830 
04831       CubitVector center = (p1 + p2 + p3)/3.0;
04832 
04833       //Calculate area
04834       CubitVector vec1 = (p2 - p1);
04835       CubitVector vec2 = (p3 - p1);
04836 
04837       CubitVector cross = vec1 * vec2;
04838         
04839       tri_area = 0.5* sqrt(cross.x() * cross.x() + cross.y() * cross.y() + cross.z() * cross.z());
04840         
04841         
04842       centroid += (tri_area * center);
04843 
04844       tot_area += tri_area;
04845         
04846     }
04847   if( tot_area == 0 )
04848     return CUBIT_FAILURE;
04849 
04850   centroid /= tot_area;
04851   return CUBIT_SUCCESS;
04852 }    
04853 
04854 CubitStatus
04855 GeomMeasureTool::center( DLIList<RefFace*> ref_faces )
04856 {
04857   int ii,id;
04858   double surf_area;
04859   double tot_area = 0.0;
04860   CubitVector surf_centroid;
04861   CubitVector centroid;
04862   centroid.set( 0.0, 0.0, 0.0 );
04863 
04864   for ( ii = ref_faces.size(); ii > 0; ii-- )
04865   {
04866     RefFace *ref_face = ref_faces.get_and_step();
04867     if( GeomMeasureTool::get_centroid( ref_face, surf_centroid, surf_area) != CUBIT_SUCCESS )
04868       return CUBIT_FAILURE;
04869     centroid += (surf_area * surf_centroid);
04870     tot_area += surf_area;
04871     id = ref_face->id();
04872   }
04873   
04874   centroid /= tot_area;
04875 
04876   if( ref_faces.size() > 1 )
04877     PRINT_INFO("Centroid of Composite surface located at: ( %f , %f , %f )\n\n",
04878                  centroid.x(), centroid.y(), centroid.z());
04879   else
04880     PRINT_INFO("Centroid of Surface %i located at: ( %f , %f , %f )\n\n",
04881                   id, centroid.x(), centroid.y(), centroid.z());
04882   return CUBIT_SUCCESS;
04883 }
04884 
04885 CubitStatus GeomMeasureTool::find_near_coincident_vertices( 
04886                             DLIList<RefVolume*> &ref_volumes,
04887                             DLIList<RefVertex*> &ref_vertices_out,
04888                             DLIList<double> &distances,
04889                             double low_tol,
04890                             double high_tol,
04891                             bool filter_same_volume_cases)
04892 {
04893   DLIList<RefVertex*> tmp_vert_list;
04894   DLIList<RefVertex*> ref_verts; 
04895   int i,j;
04896   for( i=ref_volumes.size(); i--; )
04897   {
04898     RefVolume *tmp_vol = ref_volumes.get_and_step();
04899     tmp_vert_list.clean_out();
04900     tmp_vol->ref_vertices( tmp_vert_list );
04901     ref_verts += tmp_vert_list;
04902   }
04903 
04904   //put all the vertices in a tree 
04905   AbstractTree <RefVertex*> *a_tree = new RTree<RefVertex*>( high_tol ); 
04906   for (i=ref_verts.size(); i--;)
04907     a_tree->add(ref_verts.get_and_step());
04908 
04909   std::multimap<double, dist_vert_struct> distance_vertex_map; 
04910 
04911   //for each vertex
04912   for (i=ref_verts.size(); i--;)
04913   {
04914     RefVertex *tmp_vert = ref_verts.get_and_step();
04915     RefVolume *v1 = tmp_vert->ref_volume();
04916     CubitVector vert_xyz = tmp_vert->coordinates();
04917 
04918     //get all close vertices
04919     DLIList<RefVertex*> close_verts;
04920     a_tree->find(tmp_vert->bounding_box(), close_verts);
04921 
04922     //if any vertex is between low_tol and high_tol
04923     //add it to the list
04924     DLIList<RefVertex*> near_coincident_verts;
04925     for( j=close_verts.size(); j--; )
04926     {
04927       RefVertex *close_vert = close_verts.get_and_step();
04928       if( close_vert == tmp_vert ) 
04929         continue;
04930 
04931       RefVolume *v2 = close_vert->ref_volume();
04932       bool check_distance = true;
04933       if(filter_same_volume_cases && v1 && v2 && v1 == v2)
04934         check_distance = false;
04935       if(check_distance)
04936       {
04937         double distance = vert_xyz.distance_between( close_vert->coordinates() );
04938         if( distance >= low_tol && distance <= high_tol )
04939         {
04940           dist_vert_struct tmp_struct;
04941           tmp_struct.dist = distance;
04942           tmp_struct.v1 = tmp_vert;
04943           tmp_struct.v2 = close_vert;
04944           distance_vertex_map.insert( std::multimap<double, dist_vert_struct>::
04945                                       value_type( distance, tmp_struct ));
04946         }
04947       }
04948     }
04949 
04950     a_tree->remove( tmp_vert );
04951   }
04952 
04953   std::multimap<double, dist_vert_struct>::reverse_iterator iter;
04954   
04955   iter = distance_vertex_map.rbegin();
04956   for(; iter!=distance_vertex_map.rend(); iter++ )
04957   {
04958     distances.append( (*iter).second.dist );
04959     ref_vertices_out.append( (*iter).second.v1 );
04960     ref_vertices_out.append( (*iter).second.v2 );
04961   }
04962 
04963   delete a_tree;
04964 
04965   return CUBIT_SUCCESS;
04966 }
04967 
04968 // This function is similar to find_near_coincident_vertices except for the
04969 // fact that it will only find the closest vertex in a given volume for
04970 // a vertex in another volume to be close to.  This tries to exclude the case where
04971 // you would attempt to merge one vertex from one volume to two different
04972 // vertices in another volume.
04973 CubitStatus GeomMeasureTool::find_near_coincident_vertices_unique( 
04974                             DLIList<RefVolume*> &ref_volumes,
04975                             double high_tol,
04976                             std::map <RefVertex*, DLIList<dist_vert_struct*>*> &vert_dist_map)
04977 {
04978   DLIList<RefVertex*> tmp_vert_list;
04979   DLIList<RefVertex*> ref_verts; 
04980   int i,j;
04981   for( i=ref_volumes.size(); i--; )
04982   {
04983     RefVolume *tmp_vol = ref_volumes.get_and_step();
04984     tmp_vert_list.clean_out();
04985     tmp_vol->ref_vertices( tmp_vert_list );
04986     ref_verts += tmp_vert_list;
04987   }
04988 
04989   //put all the vertices in a tree 
04990   AbstractTree <RefVertex*> *a_tree = new RTree<RefVertex*>( high_tol ); 
04991   for (i=ref_verts.size(); i--;)
04992     a_tree->add(ref_verts.get_and_step());
04993 
04994   //for each vertex
04995   for (i=ref_verts.size(); i--;)
04996   {
04997     RefVertex *tmp_vert = ref_verts.get_and_step();
04998     RefVolume *vol1 = tmp_vert->ref_volume();
04999     CubitVector vert_xyz = tmp_vert->coordinates();
05000 
05001     //get all close vertices
05002     DLIList<RefVertex*> close_verts;
05003     a_tree->find(tmp_vert->bounding_box(), close_verts);
05004 
05005     //if any vertex is between low_tol and high_tol
05006     //add it to the list
05007     DLIList<dist_vert_struct*> *near_coincident_verts = NULL;
05008     for( j=close_verts.size(); j--; )
05009     {
05010       RefVertex *close_vert = close_verts.get_and_step();
05011       if( close_vert == tmp_vert ) 
05012         continue;
05013 
05014       RefVolume *vol2 = close_vert->ref_volume();
05015       if(vol1 && vol2 && vol1 != vol2)
05016       {
05017         if(!near_coincident_verts)
05018         {
05019           near_coincident_verts = new DLIList<dist_vert_struct*>;
05020           vert_dist_map[tmp_vert] = near_coincident_verts;
05021         }
05022         double distance = vert_xyz.distance_between( close_vert->coordinates() );
05023         int h;
05024         bool found_entry_with_same_vol = false;
05025         for(h=near_coincident_verts->size(); h>0 && !found_entry_with_same_vol; h--)
05026         {
05027           dist_vert_struct* vds = near_coincident_verts->get_and_step();
05028           if(vds->vol2 == vol2)
05029           {
05030             found_entry_with_same_vol = true;
05031             if(distance < vds->dist)
05032             {
05033               vds->dist = distance;
05034               vds->vol2 = vol2;
05035               vds->v2 = close_vert;
05036             }
05037           }
05038         }
05039         if(!found_entry_with_same_vol)
05040         {
05041           dist_vert_struct *new_vds = new dist_vert_struct;
05042           new_vds->dist = distance;
05043           new_vds->v2 = close_vert;
05044           new_vds->vol2 = vol2;
05045           near_coincident_verts->append(new_vds);
05046         }
05047       }
05048     }
05049     a_tree->remove( tmp_vert );
05050   }
05051 
05052   delete a_tree;
05053 
05054   return CUBIT_SUCCESS;
05055 }
05056 
05057 struct dist_vert_vert_struct
05058 {
05059   double dist;
05060   RefVertex *vert1;
05061   RefVertex *vert2;
05062 };
05063 
05064 struct dist_vert_curve_struct
05065 {
05066   double dist;
05067   RefVertex *vert;
05068   RefEdge *edge;
05069 //  bool operator<( const dist_vert_curve_struct& b ) const
05070 //  {
05071 //    return this->dist < b.dist;
05072  // }
05073 };
05074 
05075 struct vert_curve_dist_sort
05076 {
05077   bool operator()( const dist_vert_curve_struct& a, const dist_vert_curve_struct& b ) const
05078   {
05079     return a.dist < b.dist;
05080   }
05081 };
05082 
05083 struct vert_curve_dist_sort_ptr
05084 {
05085   bool operator()( dist_vert_curve_struct *a, dist_vert_curve_struct *b ) const
05086   {
05087     return a->dist < b->dist;
05088   }
05089 };
05090 
05091 struct vert_vert_dist_sort_ptr
05092 {
05093   bool operator()( dist_vert_vert_struct *a, dist_vert_vert_struct *b ) const
05094   {
05095     return a->dist < b->dist;
05096   }
05097 };
05098 
05099 CubitStatus GeomMeasureTool::find_closest_vertex_curve_pairs(
05100                                   DLIList<RefVolume*> &vol_list,
05101                                   int &num_to_return,
05102                                   DLIList<RefVertex*> &vert_list,
05103                                   DLIList<RefEdge*> &curve_list,
05104                                   DLIList<double> &distances)
05105 {
05106   DLIList<RefFace*> surfs;
05107 
05108   int i, total_num_entries = 0;
05109   for( i=vol_list.size(); i--; )
05110   {
05111     RefVolume *tmp_vol = vol_list.get_and_step();
05112     tmp_vol->ref_faces( surfs );
05113   }
05114 
05115   std::set<dist_vert_curve_struct*,vert_curve_dist_sort_ptr> distance_vertex_curve_set; 
05116 
05117   for(i=surfs.size(); i>0; i--)
05118   {
05119     RefFace *surf = surfs.get_and_step();
05120     DLIList<RefVertex*> surf_verts;
05121     surf->ref_vertices(surf_verts);
05122     DLIList<RefEdge*> surf_curves;
05123     surf->ref_edges(surf_curves);
05124 
05125     int j;
05126     for(j=surf_verts.size(); j>0; j--)
05127     {
05128       RefVertex *tmp_vert = surf_verts.get_and_step();
05129       CubitVector vert_xyz = tmp_vert->coordinates();
05130       CubitVector closest_pt;
05131       int k;
05132       for(k=surf_curves.size(); k>0; k--)
05133       {
05134         RefEdge *cur_edge = surf_curves.get_and_step();
05135         if(cur_edge->start_vertex() != tmp_vert &&
05136           cur_edge->end_vertex() != tmp_vert)
05137         {
05138           cur_edge->closest_point_trimmed(vert_xyz, closest_pt);
05139           if(!closest_pt.about_equal(cur_edge->start_coordinates()) &&
05140              !closest_pt.about_equal(cur_edge->end_coordinates()))
05141           {
05142             double dist_sq = vert_xyz.distance_between_squared(closest_pt);
05143             dist_vert_curve_struct *tmp_struct = new dist_vert_curve_struct;
05144             tmp_struct->dist = dist_sq;
05145             tmp_struct->vert = tmp_vert;
05146             tmp_struct->edge = cur_edge; 
05147             distance_vertex_curve_set.insert( tmp_struct );
05148             total_num_entries++;
05149           }
05150         }
05151       }
05152     }
05153   }
05154 
05155   std::set<dist_vert_curve_struct*, vert_curve_dist_sort_ptr>::iterator iter, upper_iter; 
05156   
05157   int local_num_to_return = num_to_return;
05158   if(num_to_return == -1)
05159   {
05160     local_num_to_return = total_num_entries;
05161   }
05162   int cntr = 0;
05163   iter = distance_vertex_curve_set.begin();
05164   for(; iter!=distance_vertex_curve_set.end() && cntr < local_num_to_return; iter++ )
05165   {
05166     distances.append( sqrt((*iter)->dist) );
05167     vert_list.append( (*iter)->vert );
05168     curve_list.append( (*iter)->edge );
05169     cntr++;
05170   }
05171 
05172   iter = distance_vertex_curve_set.begin();
05173   for(; iter!=distance_vertex_curve_set.end(); iter++ )
05174   {
05175     delete *iter;
05176   }
05177 
05178   return CUBIT_SUCCESS;
05179 }
05180 
05181 CubitStatus GeomMeasureTool::find_closest_vertex_vertex_pairs(
05182                                   DLIList<RefVolume*> &vol_list,
05183                                   int &num_to_return,
05184                                   DLIList<RefVertex*> &vert_list1,
05185                                   DLIList<RefVertex*> &vert_list2,
05186                                   DLIList<double> &distances)
05187 {
05188   std::set<dist_vert_vert_struct*,vert_vert_dist_sort_ptr> distance_vertex_vertex_set; 
05189 
05190   int i, total_num_entries = 0;
05191   for( i=vol_list.size(); i--; )
05192   {
05193     RefVolume *tmp_vol = vol_list.get_and_step();
05194     DLIList<RefVertex*> vol_verts;
05195     tmp_vol->ref_vertices(vol_verts);
05196     while(vol_verts.size() > 1)
05197     {
05198       RefVertex *vert1 = vol_verts.pop();
05199       CubitVector vert1_xyz = vert1->coordinates();
05200       int j;
05201       for(j=vol_verts.size(); j>0; j--)
05202       {
05203         RefVertex *vert2 = vol_verts.get_and_step();
05204         double dist_sq = vert2->coordinates().distance_between_squared(vert1_xyz);
05205         dist_vert_vert_struct *tmp_struct = new dist_vert_vert_struct;
05206         tmp_struct->dist = dist_sq;
05207         tmp_struct->vert1 = vert1;
05208         tmp_struct->vert2 = vert2; 
05209         distance_vertex_vertex_set.insert( tmp_struct );
05210         total_num_entries++;
05211       }
05212     }
05213   }
05214 
05215   std::set<dist_vert_vert_struct*, vert_vert_dist_sort_ptr>::iterator iter, upper_iter; 
05216   
05217   int local_num_to_return = num_to_return;
05218   if(num_to_return == -1)
05219     local_num_to_return = total_num_entries;
05220   int cntr = 0;
05221   iter = distance_vertex_vertex_set.begin();
05222   for(; iter!=distance_vertex_vertex_set.end() && cntr < local_num_to_return; iter++ )
05223   {
05224     distances.append( sqrt((*iter)->dist) );
05225     vert_list1.append( (*iter)->vert1 );
05226     vert_list2.append( (*iter)->vert2 );
05227     cntr++;
05228   }
05229 
05230   iter = distance_vertex_vertex_set.begin();
05231   for(; iter!=distance_vertex_vertex_set.end(); iter++ )
05232   {
05233     delete *iter;
05234   }
05235 
05236   return CUBIT_SUCCESS;
05237 }
05238 
05239 CubitStatus GeomMeasureTool::find_near_coincident_vertex_curve_pairs( 
05240                                 DLIList<RefVolume*> &ref_vols,
05241                                 DLIList<RefEdge*> &ref_edges,
05242                                 DLIList<RefVertex*> &ref_verts,
05243                                 DLIList<double> &distances,
05244                                 double low_tol,
05245                                 double high_tol,
05246                                 bool filter_same_volume_cases)
05247 {
05248   //get all the curves and vertices of volumes in list
05249   DLIList<RefVertex*> verts;
05250   DLIList<RefEdge*> curves;
05251 
05252   RTree<RefEdge*> a_tree(high_tol);
05253 
05254   int i,j;
05255   for( i=ref_vols.size(); i--; )
05256   {
05257     RefVolume *tmp_vol = ref_vols.get_and_step();
05258     tmp_vol->ref_vertices( verts );
05259     
05260     curves.clean_out();
05261     tmp_vol->ref_edges( curves );
05262     for( j=curves.size(); j--; )
05263     {
05264       RefEdge *tmp_edge = curves.get_and_step();
05265       a_tree.add( tmp_edge );
05266     }
05267   }
05268 
05269   ProgressTool *progress_ptr = NULL;
05270   int total_verts = verts.size();
05271   if (total_verts > 5)
05272   {
05273     progress_ptr = AppUtil::instance()->progress_tool();
05274     assert(progress_ptr != NULL);
05275     progress_ptr->start(0, 100, "Finding Near Coincident Vertex-Curve Pairs", 
05276       NULL, CUBIT_TRUE, CUBIT_TRUE);
05277   }
05278     
05279   double curr_percent = 0.0;
05280   int processed_verts = 0;
05281   std::multimap<double, dist_vert_curve_struct> distance_vertex_curve_map; 
05282 
05283   //for each vertex
05284   for( i=verts.size(); i--; )
05285   {
05286     processed_verts++;
05287     if ( progress_ptr != NULL )
05288     {
05289       curr_percent = ((double)(processed_verts))/((double)(total_verts));
05290       progress_ptr->percent(curr_percent);
05291     }
05292 
05293     RefVertex *tmp_vert = verts.get_and_step();
05294     RefVolume *v1 = tmp_vert->ref_volume();
05295     CubitBox vertex_box ( tmp_vert->coordinates(), tmp_vert->coordinates() );
05296     DLIList<RefEdge*> close_curves;
05297     a_tree.find( vertex_box, close_curves );
05298 
05299     CubitVector vertex_xyz = tmp_vert->coordinates(); 
05300 
05301     for( j=close_curves.size(); j--; )
05302     {
05303       RefEdge *tmp_edge = close_curves.get_and_step();
05304       RefVolume *v2 = tmp_edge->ref_volume();
05305 
05306       bool check_distance = true;
05307       if(filter_same_volume_cases && v1 && v2 && v1 == v2)
05308         check_distance = false;
05309       if(check_distance)
05310       {
05311         CubitVector closest_location;
05312         tmp_edge->closest_point_trimmed( vertex_xyz, closest_location );
05313         double distance = closest_location.distance_between( vertex_xyz );
05314 
05315         if( distance >= low_tol && distance <= high_tol )
05316         {
05317           dist_vert_curve_struct tmp_struct;
05318           tmp_struct.dist = distance;
05319           tmp_struct.vert = tmp_vert;
05320           tmp_struct.edge = tmp_edge; 
05321 
05322           distance_vertex_curve_map.insert( std::multimap<double, dist_vert_curve_struct>::
05323                                             value_type( distance, tmp_struct ));
05324         }
05325       }
05326     }
05327   }
05328 
05329   if ( progress_ptr != NULL )
05330     progress_ptr->end();
05331 
05332   //std::set<dist_vert_curve_struct, vert_curve_dist_sort>::iterator iter, upper_iter; 
05333   std::multimap<double, dist_vert_curve_struct>::reverse_iterator iter;
05334   
05335   iter = distance_vertex_curve_map.rbegin();
05336   for(; iter!=distance_vertex_curve_map.rend(); iter++ )
05337   {
05338     distances.append( (*iter).second.dist );
05339     ref_verts.append( (*iter).second.vert );
05340     ref_edges.append( (*iter).second.edge );
05341   }
05342 
05343   return CUBIT_SUCCESS;
05344 }
05345 
05346 
05347 struct dist_vert_surf_struct
05348 {
05349   double dist;
05350   RefVertex *vert;
05351   RefFace *face;
05352 };
05353 
05354 struct vert_surf_dist_sort
05355 {
05356   bool operator()( dist_vert_surf_struct a, dist_vert_surf_struct b ) const
05357   {
05358     if( a.dist < b.dist )
05359       return true;
05360     else if( a.dist > b.dist )
05361       return false;
05362     else 
05363     {
05364       if( a.vert < b.vert )
05365         return true;
05366       else if( a.vert > b.vert )
05367         return false;
05368       else if( a.face < b.face )
05369         return true;
05370       else if( a.face > b.face )
05371         return false;
05372     }
05373     return false;
05374   }
05375 };
05376 
05377 
05378 CubitStatus GeomMeasureTool::find_near_coincident_vertex_surface_pairs( 
05379                                 DLIList<RefVolume*> &ref_vols,
05380                                 DLIList<RefFace*> &ref_faces,
05381                                 DLIList<RefVertex*> &ref_verts,
05382                                 DLIList<double> &distances,
05383                                 double low_tol,
05384                                 double high_tol,
05385                                 bool filter_same_volume_cases)
05386 {
05387   //get all the curves and vertices of volumes in list
05388   DLIList<RefVertex*> verts;
05389   DLIList<RefFace*> faces;
05390 
05391   AbstractTree<RefFace*> *a_tree = new RTree<RefFace*>( high_tol );
05392 
05393   int i,j;
05394   for( i=ref_vols.size(); i--; )
05395   {
05396     RefVolume *tmp_vol = ref_vols.get_and_step();
05397     tmp_vol->ref_vertices( verts );
05398     
05399     faces.clean_out();
05400     tmp_vol->ref_faces( faces );
05401     // Populate the Surface AbstractTree
05402     for( j=faces.size(); j--; )
05403     {
05404       RefFace *tmp_face = faces.get_and_step();
05405       a_tree->add( tmp_face );
05406     }
05407   }
05408 
05409   ProgressTool *progress_ptr = NULL;
05410   int total_verts = verts.size();
05411   if (total_verts > 50)
05412   {
05413     progress_ptr = AppUtil::instance()->progress_tool();
05414     assert(progress_ptr != NULL);
05415     progress_ptr->start(0, 100, "Finding Near Coincident Vertex-Surface Pairs", 
05416       NULL, CUBIT_TRUE, CUBIT_TRUE);
05417   }
05418   double curr_percent = 0.0;
05419   int processed_verts = 0;
05420 
05421   std::multimap<double, dist_vert_surf_struct> distance_vertex_surface_map; 
05422 
05423   //for each vertex
05424   for( i=verts.size(); i--; )
05425   {
05426     processed_verts++;
05427     if ( progress_ptr != NULL )
05428     {
05429       curr_percent = ((double)(processed_verts))/((double)(total_verts));
05430       progress_ptr->percent(curr_percent);
05431     }
05432 
05433     RefVertex *tmp_vert = verts.get_and_step();
05434     RefVolume *v1 = tmp_vert->ref_volume();
05435     CubitBox vertex_box ( tmp_vert->coordinates(), tmp_vert->coordinates() );
05436     DLIList<RefFace*> close_faces;
05437     a_tree->find( vertex_box, close_faces);
05438 
05439     CubitVector vertex_xyz = tmp_vert->coordinates(); 
05440 
05441     for( j=close_faces.size(); j--; )
05442     {
05443       RefFace *tmp_face = close_faces.get_and_step();
05444       RefVolume *v2 = tmp_face->ref_volume();
05445 
05446       bool check = true;
05447       if(filter_same_volume_cases && v1 && v2 && v1 == v2)
05448         check = false;
05449 
05450       if(check)
05451       {
05452         DLIList<RefVertex*> tmp_verts;
05453         tmp_face->ref_vertices( tmp_verts );
05454         if( tmp_verts.is_in_list( tmp_vert ) ) 
05455           continue;
05456 
05457         CubitVector closest_location;
05458         tmp_face->find_closest_point_trimmed( vertex_xyz, closest_location );
05459         double distance = closest_location.distance_between( vertex_xyz );
05460 
05461         if( distance > low_tol && distance < high_tol )
05462         {
05463           dist_vert_surf_struct tmp_struct;
05464           tmp_struct.dist = distance;
05465           tmp_struct.vert = tmp_vert;
05466           tmp_struct.face = tmp_face; 
05467           distance_vertex_surface_map.insert( std::multimap<double, dist_vert_surf_struct>::
05468                                       value_type( distance, tmp_struct ));
05469         }
05470       }
05471     }
05472   }
05473 
05474   if ( progress_ptr != NULL )
05475     progress_ptr->end();
05476 
05477   std::multimap<double, dist_vert_surf_struct>::reverse_iterator iter;
05478 
05479   iter = distance_vertex_surface_map.rbegin();
05480   for(; iter!=distance_vertex_surface_map.rend(); iter++ )
05481   {
05482     distances.append( (*iter).second.dist );
05483     ref_verts.append( (*iter).second.vert );
05484     ref_faces.append( (*iter).second.face);
05485   }
05486 
05487   delete a_tree;
05488 
05489   return CUBIT_SUCCESS;
05490 }
05491 
05492 
05493 
05494 
05495 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines