MeshKit  1.0
arc.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <iomanip> // for setprecision
00003 #include <limits>  // for double min/max
00004 #include <assert.h>
00005 #include <vector>
00006 #include "moab/Range.hpp"
00007 #include "moab/AdaptiveKDTree.hpp"
00008 #include "meshkit/arc.hpp"
00009 #include "meshkit/zip.hpp"
00010 #include "meshkit/gen.hpp"
00011 
00012 #include "moab/GeomTopoTool.hpp"
00013 #include "moab/FileOptions.hpp"
00014 
00015 using namespace moab;
00016 namespace arc {
00017 
00018   ErrorCode orient_edge_with_tri( const EntityHandle edge, const EntityHandle tri ) {
00019     ErrorCode result;
00020     // get the connected vertices, properly ordered
00021     const EntityHandle *tri_conn;
00022     int n_verts;
00023     result = MBI()->get_connectivity( tri, tri_conn, n_verts );
00024     assert(MB_SUCCESS == result);
00025     assert( 3 == n_verts );
00026 
00027     // get the endpoints of the edge
00028     const EntityHandle *edge_conn;
00029     result = MBI()->get_connectivity( edge, edge_conn, n_verts );
00030     assert(MB_SUCCESS == result);
00031     assert( 2 == n_verts );
00032 
00033     // if the edge is backwards, reverse it
00034     if (( edge_conn[0]==tri_conn[0] && edge_conn[1]==tri_conn[2] ) ||
00035         ( edge_conn[0]==tri_conn[1] && edge_conn[1]==tri_conn[0] ) ||
00036         ( edge_conn[0]==tri_conn[2] && edge_conn[1]==tri_conn[1] ) ) {
00037       EntityHandle new_conn[2];
00038       new_conn[0] = edge_conn[1];
00039       new_conn[1] = edge_conn[0];
00040       result = MBI()->set_connectivity( edge, new_conn, 2 );
00041       assert(MB_SUCCESS == result);
00042     }
00043     return result;
00044   } 
00045 
00046   // Degenerate edges (same topological endpts) are caused by a prior step in which
00047   // coincident verts are merged.
00048   ErrorCode remove_degenerate_edges( Range &edges, const bool debug ) {
00049     Range::iterator i = edges.begin();
00050     while (i!=edges.end()) {
00051       // get the endpoints of the edge
00052       ErrorCode rval;
00053       const EntityHandle *endpts;
00054       int n_verts;
00055       rval = MBI()->get_connectivity( *i, endpts, n_verts );
00056       if(gen::error(MB_SUCCESS!=rval,"could not get connectivity"))
00057         return rval;
00058 
00059       // remove the edge if degenerate
00060       if(2==n_verts && endpts[0]!=endpts[1]) {
00061         ++i;
00062       } else if( (2==n_verts && endpts[0]==endpts[1]) ||
00063                  (1==n_verts                        ) ) {
00064         if(debug) {
00065           std::cout << "remove_degenerate_edges: deleting degenerate edge and tris " 
00066                     << std::endl;
00067         }
00068         rval = zip::delete_adj_degenerate_tris( endpts[0] );
00069         if(gen::error(MB_SUCCESS!=rval,"could not delete degenerate tris")) return rval;
00070         rval = MBI()->delete_entities( &(*i), 1 );
00071         if(gen::error(MB_SUCCESS!=rval,"could not delete degenerate edge")) return rval;
00072         i = edges.erase(i);
00073       } else {
00074         std::cout << "remove_degenerate_edge: wrong edge connectivity size" << std::endl;
00075         return MB_FAILURE;
00076       }
00077     
00078     }
00079     return MB_SUCCESS;
00080   }
00081 
00082 
00083   // Given a range of edges, remove pairs that have vertices (a,b) (b,a)
00084   ErrorCode remove_opposite_pairs_of_edges( Range &edges, const bool debug ) {
00085 
00086     // do this in O(n) by using adjacencies instead of O(n^2)
00087     ErrorCode result;
00088     //for(Range::iterator i=edges.begin(); i!=edges.end(); i++ ) {
00089     for(unsigned int i=0; i<edges.size(); i++) {
00090       EntityHandle the_edge = edges[i];
00091 
00092       // get endpoint verts
00093       Range two_verts;
00094       result = MBI()->get_adjacencies( &the_edge, 1, 0, false, two_verts);
00095       if(MB_SUCCESS != result) {
00096         std::cout << "result=" << result << " could not get adjacencies of edge" << std::endl;
00097         return result;
00098       }
00099 
00100       // get adjacent edges, but only keep the edges adjacent to both verts
00101       Range adj_edges;
00102       result = MBI()->get_adjacencies( two_verts, 1, false, adj_edges, Interface::INTERSECT);
00103       assert(MB_SUCCESS == result);
00104       // remove the original edge
00105       //adj_edges.erase( *i );
00106 
00107       // if any other edges exist, they are opposite the original edge and should be
00108       // removed from the skin
00109       if ( 1<adj_edges.size() ) {
00110         if(debug) {
00111           std::cout << adj_edges.size() 
00112                     << " opposite edges will be removed from the surface skin " 
00113                     << adj_edges[0] << " " << adj_edges[1] << std::endl;
00114         }
00115         //gen::print_range_of_edges( adj_edges );
00116         //gen::print_range_of_edges( edges );
00117         //edges = edges.subtract( adj_edges );
00118         //edges.erase( *i );
00119         edges = subtract( edges, adj_edges );
00120         result = MBI()->delete_entities( adj_edges );
00121         assert(MB_SUCCESS == result);
00122         i--;
00123       }
00124     }
00125     return MB_SUCCESS;
00126   }
00127 
00128   ErrorCode remove_opposite_pairs_of_edges_fast( Range &edges, const bool debug) {
00129     // special case
00130     ErrorCode rval;
00131     if(1==edges.size()) {
00132       std::cout << "cannot remove pairs: only one input edge" << std::endl;
00133       return MB_FAILURE;
00134     }
00135 
00136     // populate edge array, used only for searching
00137     unsigned n_orig_edges = edges.size();
00138     gen::edge *my_edges = new gen::edge[n_orig_edges];
00139     unsigned j = 0;
00140     for(Range::const_iterator i=edges.begin(); i!=edges.end(); ++i) {
00141       // get the endpoints of the edge
00142       const EntityHandle *endpts;
00143       int n_verts;
00144       rval = MBI()->get_connectivity( *i, endpts, n_verts );
00145       if(gen::error(MB_SUCCESS!=rval || 2!=n_verts,"could not get connectivity"))
00146         return rval;
00147     
00148       // store the edges
00149       my_edges[j].e    = *i;
00150       my_edges[j].v0   = endpts[0];
00151       my_edges[j].v1   = endpts[1];
00152 
00153       // sort edge by handle
00154       if(my_edges[j].v1 < my_edges[j].v0) {
00155         EntityHandle temp = my_edges[j].v0;
00156         my_edges[j].v0 = my_edges[j].v1;
00157         my_edges[j].v1 = temp;
00158       }
00159       ++j;
00160     }
00161 
00162     // sort edge array
00163     qsort(my_edges, n_orig_edges, sizeof(struct gen::edge), gen::compare_edge);
00164 
00165     // find duplicate edges
00166     j=0;
00167     Range duplicate_edges;
00168     for(unsigned i=1; i<n_orig_edges; ++i) {
00169       // delete edge if a match exists
00170       if(my_edges[j].v0==my_edges[i].v0 && my_edges[j].v1==my_edges[i].v1) {
00171         duplicate_edges.insert( my_edges[j].e );
00172         // find any remaining matches
00173         while( my_edges[j].v0==my_edges[i].v0 && 
00174                my_edges[j].v1==my_edges[i].v1 &&
00175                i<n_orig_edges) {
00176           duplicate_edges.insert( my_edges[i].e );
00177           ++i;
00178         }
00179         // delete the matches
00180         edges = subtract( edges, duplicate_edges );
00181         rval = MBI()->delete_entities( duplicate_edges );
00182         if(gen::error(MB_SUCCESS!=rval,"cannot delete edge")) {
00183           delete[] my_edges;
00184           return rval;
00185         }
00186         if(debug) {
00187           std::cout << "remove_opposite_edges: deleting " << duplicate_edges.size() 
00188                     << " edges" << std::endl;
00189         }
00190         duplicate_edges.clear();
00191       }
00192       j = i;
00193     }
00194     
00195     delete[] my_edges;
00196     return MB_SUCCESS;
00197   }
00198 
00199   ErrorCode get_next_oriented_edge( const Range edges, const EntityHandle edge,
00200                                       EntityHandle &next_edge ) {
00201 
00202     // get the back vertex
00203     ErrorCode result;
00204     const EntityHandle *end_verts;
00205     int n_verts;
00206     result = MBI()->get_connectivity( edge, end_verts, n_verts );
00207     assert(MB_SUCCESS==result);
00208     assert( 2 == n_verts );
00209 
00210     // get the edges adjacent to the back vertex
00211     Range adj_edges;
00212     result = MBI()->get_adjacencies( &(end_verts[1]), 1, 1, false, adj_edges );
00213     assert(MB_SUCCESS==result);
00214 
00215     // keep the edges that are part of the input range
00216     adj_edges = intersect( adj_edges, edges );
00217     // don't want the input edge
00218     adj_edges.erase( edge );
00219 
00220     // make sure the edge is oriented correctly
00221     for(Range::iterator i=adj_edges.begin(); i!=adj_edges.end(); i++) {
00222       const EntityHandle *adj_end_verts;
00223       result = MBI()->get_connectivity( *i, adj_end_verts, n_verts );
00224       if(MB_SUCCESS != result) {
00225         MBI()->list_entity(*i);
00226         std::cout << "result=" << result 
00227                   << " could not get connectivity of edge" << std::endl;
00228         return result;
00229         //gen::print_edge( *i );
00230       }
00231       assert(MB_SUCCESS==result);
00232       assert( 2 == n_verts );
00233       if ( end_verts[1]!=adj_end_verts[0] ) i = adj_edges.erase(i) - 1;
00234     }
00235 
00236     /* The next edge could be ambiguous if more than one exists.This happens in 
00237        surfaces that are ~1D, and in surfaces that have pinch points 
00238        (mod13surf881). Although I didn't handle this case yet, if it occurs:
00239        -Remember that a pinch point could have not just 2, but multiple ears.
00240        -Select the next edge as the edge that shares the same triangle.
00241        -This approach will only work if the pinch point also exists in the 
00242         geometric curves.
00243        -If the pinch point does not exist in the geometric curves (~1D surfs),
00244         there is no robust way to handle it. There should be an input assumption
00245         that this never happens. "The faceting skin cannot have pinch points
00246         unless they also occur in the surface's geometric curves."
00247      */
00248     if ( 0==adj_edges.size() ) {
00249       next_edge = 0;
00250     } else if ( 1==adj_edges.size() ) {
00251       next_edge = adj_edges.front();
00252     } else {
00253       std::cout << "get_next_oriented_edge: " << adj_edges.size() <<
00254         " possible edges indicates a pinch point." << std::endl;
00255       result = MBI()->list_entity( end_verts[1] );
00256       //assert(MB_SUCCESS == result);
00257       //return MB_MULTIPLE_ENTITIES_FOUND;
00258       next_edge = adj_edges.front();
00259     }
00260     return MB_SUCCESS;
00261   }
00262 
00263   ErrorCode create_loops_from_oriented_edges_fast( Range edges,
00264                                                      std::vector< std::vector<EntityHandle> > &loops_of_edges,
00265                                                      const bool debug ) {
00266     // place all edges in map
00267     std::multimap<EntityHandle,gen::edge> my_edges;
00268     ErrorCode rval;
00269     for(Range::const_iterator i=edges.begin(); i!=edges.end(); ++i) {
00270       // get the endpoints of the edge
00271       const EntityHandle *endpts;
00272       int n_verts;
00273       rval = MBI()->get_connectivity( *i, endpts, n_verts );
00274       if(gen::error(MB_SUCCESS!=rval || 2!=n_verts,"could not get connectivity"))
00275         return rval;
00276     
00277       // store the edges
00278       gen::edge temp;
00279       temp.e    = *i;
00280       temp.v0   = endpts[0];
00281       temp.v1   = endpts[1];
00282       my_edges.insert( std::pair<EntityHandle,gen::edge>(temp.v0,temp) );
00283     }
00284     std::cout << "error: function not complete" << std::endl;
00285     return MB_FAILURE;
00286 
00287     return MB_SUCCESS;
00288   }
00289 
00290   // This function should be rewritten using multimaps or something to avoid
00291   // upward adjacency searching. Vertices are searched for their adjacent edges.
00292   ErrorCode create_loops_from_oriented_edges( Range edges,
00293                                                 std::vector< std::vector<EntityHandle> > &loops_of_edges,
00294                                                 const bool debug ) {
00295 
00296     // conserve edges
00297     ErrorCode result;
00298     unsigned int n_edges_in  = edges.size();
00299     unsigned int n_edges_out = 0;
00300     //gen::print_range_of_edges( edges );
00301     // there could be several arcs for each surface
00302     while ( 0!= edges.size() ) {
00303       std::vector<EntityHandle> loop_of_edges;
00304       // pick initial edge and point
00305       EntityHandle edge = edges.front();
00306 
00307       // 20091201 Update: Pinch points may not be important. If not, there is no
00308       // purpose detecting them. Instead assume that pinch points coincide with 
00309       // the endpoints of geometric curves. Also assume that the loop creation at
00310       // pinch points does not matter. Pinch points can result in one or more 
00311       // loops, depending upon the path of traversal through the point.
00312 
00313       // Check to make sure the beginning endpt of the first edge is not a pinch 
00314       // point. If it is a pinch point the loop is ambiguous. Maybe--see watertightness notes for 20091201
00315       {
00316         const EntityHandle *end_verts;
00317         int n_verts;
00318         result = MBI()->get_connectivity( edge, end_verts, n_verts );
00319         assert(MB_SUCCESS==result);
00320         assert( 2 == n_verts );
00321         // get the edges adjacent to the back vertex
00322         Range adj_edges;
00323         result = MBI()->get_adjacencies( &(end_verts[0]), 1, 1, false, adj_edges );
00324         assert(MB_SUCCESS==result);
00325         // keep the edges that are part of the input range
00326         adj_edges = intersect( adj_edges, edges );
00327         if(2!=adj_edges.size() && debug) {
00328           std::cout << "  create_loops: adj_edges.size()=" << adj_edges.size() << std::endl;
00329           std::cout << "  create_loops: pinch point exists" << std::endl;
00330           result = MBI()->list_entity( end_verts[0] );
00331           assert(MB_SUCCESS == result);
00332           //return MB_MULTIPLE_ENTITIES_FOUND;
00333         }
00334       }
00335 
00336       // add it to the loop
00337       loop_of_edges.push_back( edge );
00338       if(debug) std::cout << "push_back: " << edge << std::endl;
00339       n_edges_out++;
00340       edges.erase( edge );
00341 
00342       // find connected edges and add to the loop
00343       EntityHandle next_edge = 0;
00344       while (true) {
00345 
00346         // get the next vertex and next edge
00347         result = get_next_oriented_edge( edges, edge, next_edge );
00348         if(MB_ENTITY_NOT_FOUND == result) {
00349           return result;
00350         } else if(MB_SUCCESS != result) {
00351           gen::print_arc_of_edges( loop_of_edges );
00352           return result;
00353         }
00354         //assert( MB_SUCCESS == result );
00355 
00356         // if the next edge was found
00357         if ( 0!=next_edge ) {
00358           // add it to the loop
00359           loop_of_edges.push_back( next_edge );
00360           if(debug) std::cout << "push_back: " << next_edge << std::endl;
00361           n_edges_out++;
00362 
00363           // remove the edge from the possible edges
00364           edges.erase( next_edge );
00365 
00366           // set the new reference vertex
00367           //vert = next_vert;
00368           edge = next_edge;
00369 
00370           // if another edge was not found
00371         } else {
00372           break;
00373 
00374         }
00375       }
00376 
00377       // check to ensure the arc is closed
00378       Range first_edge;
00379       first_edge.insert( loop_of_edges.front() );
00380       result = get_next_oriented_edge( first_edge, loop_of_edges.back(), next_edge );
00381       assert(MB_SUCCESS == result);
00382       if(next_edge != first_edge.front()) {
00383         std::cout << "create_loops: loop is not closed" << std::endl;
00384         gen::print_arc_of_edges(loop_of_edges);
00385         return MB_FAILURE;
00386       }
00387 
00388       // add the current arc to the vector of arcs
00389       loops_of_edges.push_back(loop_of_edges);
00390     }
00391 
00392     // check to make sure that we have the same number of verts as we started with
00393     if(gen::error(n_edges_in!=n_edges_out,"edges not conserved")) return MB_FAILURE;
00394     assert( n_edges_in == n_edges_out );
00395 
00396     return MB_SUCCESS;
00397   }
00398 
00399   // return a set of ordered_verts and remaining unordered_edges
00400   ErrorCode order_verts_by_edge( Range unordered_edges,
00401                                    std::vector<EntityHandle> &ordered_verts ) {
00402     if(unordered_edges.empty()) return MB_SUCCESS;
00403  
00404     // get the endpoints of the curve. It should have 2 endpoints, unless is it a circle.
00405     Range end_verts;
00406     Skinner tool(MBI());
00407     ErrorCode result;
00408     result = tool.find_skin( 0 , unordered_edges, 0, end_verts, false );
00409     if(MB_SUCCESS != result) gen::print_range_of_edges( unordered_edges );
00410     assert(MB_SUCCESS == result);
00411 
00412     // start with one endpoint
00413     EntityHandle vert, edge;
00414     if(2 == end_verts.size()) {
00415       vert = end_verts.front();
00416     } else if (0 == end_verts.size()) {
00417       result = MBI()->get_adjacencies( &unordered_edges.front(), 1, 0, false, end_verts );
00418       assert(MB_SUCCESS == result);
00419       assert(2 == end_verts.size());
00420       vert = end_verts.front();
00421     } else return MB_FAILURE;
00422     
00423     // build the ordered set of verts. It will be as large as the number
00424     // of edges, plus one extra endpoint.
00425     ordered_verts.clear();
00426     ordered_verts.push_back( vert );
00427    
00428     // this cannot be used if multiple loops exist
00429     while(!unordered_edges.empty()) {
00430       // get an edge of the vert
00431       Range adj_edges;
00432       result = MBI()->get_adjacencies( &vert, 1, 1, false, adj_edges );
00433       assert(MB_SUCCESS == result);
00434       adj_edges = intersect( adj_edges, unordered_edges );
00435       //assert(!adj_edges.empty());
00436       if(adj_edges.empty()) {
00437         std::cout << "    order_verts_by_edgs: adj_edges is empty" << std::endl;
00438         return MB_FAILURE;
00439       }
00440       edge = adj_edges.front(); 
00441       unordered_edges.erase( edge );
00442 
00443       // get the next vert
00444       end_verts.clear();
00445       result = MBI()->get_adjacencies( &edge, 1, 0, false, end_verts );
00446       assert(MB_SUCCESS == result);
00447       if(2 != end_verts.size()) {
00448         std::cout << "end_verts.size()=" << end_verts.size() << std::endl;
00449         gen::print_edge( edge );
00450       }
00451       assert(2 == end_verts.size());
00452       vert = end_verts.front()==vert ? end_verts.back() : end_verts.front();
00453       ordered_verts.push_back( vert );
00454     }
00455     return MB_SUCCESS;
00456   }
00457      
00458   ErrorCode get_meshset( const EntityHandle set, std::vector<EntityHandle> &vec) {
00459     ErrorCode result;
00460     vec.clear();
00461     result = MBI()->get_entities_by_handle( set, vec );
00462     assert(MB_SUCCESS == result);
00463     return result;
00464   }
00465 
00466   ErrorCode set_meshset( const EntityHandle set, const std::vector<EntityHandle> vec) {
00467     ErrorCode result;
00468     result = MBI()->clear_meshset( &set, 1 );
00469     assert(MB_SUCCESS == result);
00470     result = MBI()->add_entities( set, &vec[0], vec.size() );       
00471     assert(MB_SUCCESS == result);
00472     return result;
00473   }
00474 
00475   ErrorCode merge_curves( Range curve_sets, const double facet_tol,
00476                             Tag id_tag, Tag merge_tag, const bool debug ) {
00477     // find curve endpoints to add to kd tree
00478     ErrorCode result;
00479     const double SQR_TOL = facet_tol*facet_tol;
00480     Range endpoints;
00481     for(Range::iterator i=curve_sets.begin(); i!=curve_sets.end(); i++) {
00482       std::vector<EntityHandle> curve;
00483       result = get_meshset( *i, curve );
00484       assert(MB_SUCCESS == result);
00485       //if(2 > curve.size()) continue;
00486       assert(1 < curve.size());
00487       EntityHandle front_endpt = curve[0];
00488       EntityHandle back_endpt  = curve[curve.size()-1];
00489       // ADD CODE TO HANDLE SPECIAL CASES!!
00490       if(front_endpt == back_endpt) { // special case
00491         if(0 == gen::length(curve)) { // point curve
00492         } else {                      // circle
00493         }
00494       } else {                        // normal curve
00495         endpoints.insert( front_endpt );
00496         endpoints.insert( back_endpt );
00497       }
00498     }
00499 
00500     // add endpoints to kd-tree. Tree must track ownership to know when verts are
00501     // merged away (deleted).  
00502 
00503     AdaptiveKDTree kdtree(MBI()); //, 0, MESHSET_TRACK_OWNER);
00504     EntityHandle root;
00505 
00506     //set tree options
00507     const char settings[]="MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;MESHSET_FLAGS=0x1;TAG_NAME=0";
00508     FileOptions fileopts(settings);
00509 
00510     /* Old settings for the KD Tree                                            
00511     // initialize settings of the KD Tree
00512     AdaptiveKDTree::Settings settings;
00513     // sets the tree to split any leaves with more than 1 entity                
00514     settings.maxEntPerLeaf = 1;                                 
00515     // tells the tree how many candidate planes to consider for each dim of the tree                    
00516     settings.candidateSplitsPerDir = 1;                           
00517     // planes are set to be at evenly spaced intervals
00518     settings.candidatePlaneSet = AdaptiveKDTree::SUBDIVISION;
00519     */
00520     
00521 
00522 
00523 
00524     // initialize the tree and pass the root entity handle back into root
00525     result = kdtree.build_tree( endpoints, &root, &fileopts);            
00526     assert(MB_SUCCESS == result);
00527     // create tree iterator                                         
00528     AdaptiveKDTreeIter tree_iter;
00529     kdtree.get_tree_iterator( root, tree_iter );     
00530 
00531     // search for other endpoints that match each curve's endpoints
00532     for(Range::iterator i=curve_sets.begin(); i!=curve_sets.end(); i++) {
00533       std::vector<EntityHandle> curve_i_verts;
00534       result = get_meshset( *i, curve_i_verts );
00535       assert(MB_SUCCESS == result);
00536       double curve_length = gen::length( curve_i_verts );
00537       //if(2 > curve.size()) continue; // HANDLE SPECIAL CASES (add logic)     
00538       if(curve_i_verts.empty()) continue;
00539       EntityHandle endpts[2] = { curve_i_verts.front(), curve_i_verts.back() };
00540       CartVect endpt_coords;
00541       //if( endpts[0] == endpts[1]) continue; // special case of point curve or circle
00542       std::vector<EntityHandle> leaves;
00543 
00544       // initialize an array which will contain matched of front points in [0] and 
00545       // matches for back points in [1]
00546       Range adj_curves[2];
00547       // match the front then back endpts                        
00548       for(unsigned int j=0; j<2; j++) {                                  
00549         result = MBI()->get_coords( &endpts[j], 1, endpt_coords.array());
00550         assert(MB_SUCCESS == result);
00551         // takes all leaves of the tree within a distance (facet_tol) of the coordinates
00552         // passed in by endpt_coords.array() and returns them in leaves
00553         result = kdtree.distance_search( endpt_coords.array(), 
00554                                                 facet_tol, leaves, root );
00555         assert(MB_SUCCESS == result);
00556         for(unsigned int k=0; k<leaves.size(); k++) {           
00557           // retrieve all information about vertices in leaves
00558           std::vector<EntityHandle> leaf_verts;
00559           result = MBI()->get_entities_by_type( leaves[k], MBVERTEX, leaf_verts);
00560           assert(MB_SUCCESS == result);
00561           for(unsigned int l=0; l<leaf_verts.size(); l++) {               
00562             double sqr_dist;                                             
00563             result = gen::squared_dist_between_verts( endpts[j], leaf_verts[l], sqr_dist);        
00564             assert(MB_SUCCESS == result);
00565             /* Find parent curves. There will be no parent curves if the curve has
00566                already been merged and no longer exists. */     
00567             if(SQR_TOL >= sqr_dist) {
00568               Range temp_curves;
00569               // get the curves for all vertices that are within the squared distance of each other
00570               result = MBI()->get_adjacencies( &leaf_verts[l], 1, 4, false, temp_curves);
00571               assert(MB_SUCCESS == result);
00572               // make sure the sets are curve sets before adding them to the list
00573               // of candidates.
00574               temp_curves = intersect(temp_curves, curve_sets);
00575               adj_curves[j].merge( temp_curves );
00576             }
00577           }
00578         }
00579       }
00580 
00581       // now find the curves that have matching endpts
00582       Range candidate_curves;
00583       // select curves that do not have coincident front AND back points
00584       // place them into candidate curves
00585       candidate_curves =intersect( adj_curves[0], adj_curves[1] );
00586       if(candidate_curves.empty()) continue;
00587 
00588       // subtract the current curve
00589       //int n_before = candidate_curves.size();
00590       candidate_curves.erase( *i );
00591       //int n_after = candidate_curves.size();
00592       
00593 
00594       // now find curves who's interior vertices are also coincident and merge them
00595       for(Range::iterator j=candidate_curves.begin(); j!=candidate_curves.end(); j++) {
00596         std::vector<EntityHandle> curve_j_verts;
00597         result = get_meshset( *j, curve_j_verts );
00598         assert(MB_SUCCESS == result);
00599         double j_curve_length = gen::length( curve_j_verts );
00600 
00601         int i_id, j_id;
00602         result = MBI()->tag_get_data( id_tag, &(*i), 1, &i_id );                         
00603         assert(MB_SUCCESS == result);
00604         result = MBI()->tag_get_data( id_tag, &(*j), 1, &j_id );                         
00605         assert(MB_SUCCESS == result);
00606         if(debug) {
00607           std::cout << "curve i_id=" << i_id << " j_id=" << j_id 
00608                     << " leng0=" << curve_length << " leng1=" << j_curve_length << std::endl;
00609         }
00610 
00611         // reject curves with significantly different length (for efficiency)
00612         if( facet_tol < abs(curve_length - j_curve_length)) continue;
00613 
00614         // align j_curve to be the same as i_curve
00615         bool reversed;
00616         if(gen::dist_between_verts(curve_i_verts.front(), curve_j_verts.front()) >
00617            gen::dist_between_verts(curve_i_verts.front(), curve_j_verts.back())) {
00618           reverse( curve_j_verts.begin(), curve_j_verts.end() );
00619           reversed = true;
00620         } else {
00621           reversed = false;
00622         }
00623 
00624         // Reject curves if the average distance between them is greater than 
00625         // facet_tol.
00626         double dist;
00627         result = gen::dist_between_arcs( debug, curve_i_verts, curve_j_verts, dist );
00628         assert(MB_SUCCESS == result);
00629         if( facet_tol < dist ) continue;
00630 
00631         // THE CURVE WILL BE MERGED
00632         if (debug)
00633           {
00634             std::cout << "  merging curve " << j_id << " to curve " << i_id 
00635                       << ", dist_between_curves=" << dist << " cm" << std::endl;
00636           }
00637 
00638         // Merge the endpts of the curve to preserve topology. Merging (and deleting)
00639         // the endpoints will also remove them from the KDtree so that the merged
00640         // curve cannot be selected again. Update the curves when merging to avoid
00641         // stale info.
00642         if(curve_i_verts.front() != curve_j_verts.front()) { 
00643           result = zip::merge_verts( curve_i_verts.front(), curve_j_verts.front(), 
00644                                      curve_i_verts, curve_j_verts );
00645           if(MB_SUCCESS != result) std::cout << result << std::endl;
00646           assert(MB_SUCCESS == result);
00647         }
00648         if(curve_i_verts.back() != curve_j_verts.back()) {
00649           result = zip::merge_verts( curve_i_verts.back(), curve_j_verts.back(),
00650                                      curve_i_verts, curve_j_verts );
00651           if(MB_SUCCESS != result) std::cout << result << std::endl;
00652           assert(MB_SUCCESS == result);
00653         }
00654 
00655         // Tag the curve that is merged away. We do not delete it so that its 
00656         // parent-child links are preserved. Later, a surface's facets will be
00657         // deleted if all of its curves are deleted or merged with themselves.
00658         // Only after this can the merged away curves be deleted.
00659         result = MBI()->tag_set_data( merge_tag, &(*j), 1, &(*i) );
00660         assert(MB_SUCCESS == result);
00661 
00662         // clear the sets contents
00663         curve_j_verts.clear();
00664         result = set_meshset( *j, curve_j_verts ); 
00665         assert(MB_SUCCESS == result);
00666 
00667         // reverse the curve-surf senses if the merged curve was found to be opposite the 
00668         // curve we keep
00669         if(reversed) {
00670           std::vector<EntityHandle> surfs;
00671           std::vector<int> senses;
00672           GeomTopoTool gt(MBI(), false);
00673           result = gt.get_senses( *j, surfs, senses );
00674           if(gen::error(MB_SUCCESS!=result,"failed to get senses")) return result;
00675           for(unsigned k=0; k<surfs.size(); ++k) {
00676             //forward to reverse
00677             if(SENSE_FORWARD==senses[k])
00678               senses[k] = SENSE_REVERSE;
00679             //reverse to forward
00680             else if(SENSE_REVERSE==senses[k])
00681               senses[k] = SENSE_FORWARD;
00682             //unknown to unknown 
00683             else if(SENSE_UNKNOWN==senses[k])
00684               senses[k] = SENSE_UNKNOWN;
00685             else
00686               if(gen::error(true,"unrecognized sense")) return MB_FAILURE;
00687           }   
00688           result = gt.set_senses( *j, surfs, senses );
00689           if(gen::error(MB_SUCCESS!=result,"failed to set senses")) return result;
00690         }
00691 
00692       }
00693     }
00694     return MB_SUCCESS;
00695   }
00696 
00697 
00698 
00699 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines