MeshKit  1.0
gen.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <iomanip> // for setprecision
00003 #include <limits> // for min/max values
00004 #include <assert.h>
00005 #include <math.h>
00006 #include <time.h>
00007 #include <vector>
00008 
00009 #include "moab/GeomTopoTool.hpp"
00010 #include "moab/FileOptions.hpp"
00011 #include "moab/Core.hpp"
00012 
00013 #include "meshkit/gen.hpp"
00014 #include "meshkit/zip.hpp"
00015 #include "moab/Skinner.hpp"
00016 
00017 
00018 const char GEOM_SENSE_2_TAG_NAME[] = "GEOM_SENSE_2";
00019 const char GEOM_SENSE_N_ENTS_TAG_NAME[] = "GEOM_SENSE_N_ENTS";
00020 const char GEOM_SENSE_N_SENSES_TAG_NAME[] = "GEOM_SENSE_N_SENSES"; 
00021 
00022 using namespace moab;
00023 namespace gen {
00024 
00025   bool error( const bool error_has_occured, const std::string message ) {
00026     if(error_has_occured) {
00027       if("" == message) {
00028         std::cout << "Error at " << __FILE__ << ":" << __LINE__ 
00029                   << std::endl;
00030       } else {
00031         std::cout << message << std::endl;
00032       }
00033       return true;
00034     } else {
00035       return false;
00036     }
00037   }
00038 
00039 void moab_printer(ErrorCode error_code)
00040 {
00041   if ( error_code == MB_INDEX_OUT_OF_RANGE )
00042     {
00043       std::cerr << "ERROR: MB_INDEX_OUT_OF_RANGE" << std::endl;
00044     }
00045   if ( error_code == MB_MEMORY_ALLOCATION_FAILED )
00046     {
00047       std::cerr << "ERROR: MB_MEMORY_ALLOCATION_FAILED" << std::endl;
00048     }
00049   if ( error_code == MB_ENTITY_NOT_FOUND )
00050     {
00051       std::cerr << "ERROR: MB_ENTITY_NOT_FOUND" << std::endl;
00052     }
00053   if ( error_code == MB_MULTIPLE_ENTITIES_FOUND )
00054     {
00055       std::cerr << "ERROR: MB_MULTIPLE_ENTITIES_FOUND" << std::endl;
00056     }
00057   if ( error_code == MB_TAG_NOT_FOUND )
00058     {
00059       std::cerr << "ERROR: MB_TAG_NOT_FOUND" << std::endl;
00060     }
00061   if ( error_code == MB_FILE_DOES_NOT_EXIST )
00062     {
00063       std::cerr << "ERROR: MB_FILE_DOES_NOT_EXIST" << std::endl;
00064     }    
00065   if ( error_code == MB_FILE_WRITE_ERROR )
00066     {
00067       std::cerr << "ERROR: MB_FILE_WRITE_ERROR" << std::endl;
00068     }    
00069   if ( error_code == MB_ALREADY_ALLOCATED )
00070     {
00071       std::cerr << "ERROR: MB_ALREADY_ALLOCATED" << std::endl;
00072     }    
00073   if ( error_code == MB_VARIABLE_DATA_LENGTH )
00074     {
00075       std::cerr << "ERROR: MB_VARIABLE_DATA_LENGTH" << std::endl;
00076     }  
00077   if ( error_code == MB_INVALID_SIZE )
00078     {
00079       std::cerr << "ERROR: MB_INVALID_SIZE" << std::endl;
00080     }  
00081   if ( error_code == MB_UNSUPPORTED_OPERATION )
00082     {
00083       std::cerr << "ERROR: MB_UNSUPPORTED_OPERATION" << std::endl;
00084     }  
00085   if ( error_code == MB_UNHANDLED_OPTION )
00086     {
00087       std::cerr << "ERROR: MB_UNHANDLED_OPTION" << std::endl;
00088     }  
00089   if ( error_code == MB_FAILURE )
00090     {
00091       std::cerr << "ERROR: MB_FAILURE" << std::endl;
00092     }  
00093   return;
00094 }
00095 
00096 
00097 
00098 
00099 
00100   void print_vertex_cubit( const EntityHandle vertex ) {
00101 
00102     ErrorCode result;
00103     double coords[3];
00104     int n_precision = 20;
00105     result = MBI()->get_coords( &vertex, 1, coords );
00106     if(gen::error(MB_SUCCESS!=result, "failed to get vertex coords"));
00107     assert(MB_SUCCESS == result);
00108     std::cout << "  create vertex " 
00109               << std::setprecision(n_precision)
00110               << coords[0] << " " << coords[1] << " " << coords[2]
00111               << std::endl;
00112   }
00113 
00114   void print_vertex_coords( const EntityHandle vertex ) {
00115 
00116     ErrorCode result;
00117     double coords[3];
00118     result = MBI()->get_coords( &vertex, 1, coords );
00119     if(MB_SUCCESS!=result) std::cout << "vert=" << vertex << std::endl;
00120     assert(MB_SUCCESS == result);
00121     std::cout << "    vertex " << vertex << " coords= (" 
00122               << coords[0] << "," << coords[1] << "," << coords[2] << ")" 
00123               << std::endl;
00124   }
00125 
00126   void print_triangles( const Range tris ) {
00127     for(Range::const_iterator i=tris.begin(); i!=tris.end(); i++) {
00128       print_triangle( *i, false );
00129     }
00130   }
00131   // If the edges of the tri are ambiguous, do not print edges!
00132   void print_triangle( const EntityHandle tri, bool print_edges ) {
00133     ErrorCode result;
00134     double area;
00135     result = triangle_area( tri, area );
00136     assert(MB_SUCCESS == result);
00137     std::cout << "    triangle " << tri << " area=" << area << std::endl;
00138     const EntityHandle *conn;
00139     int n_verts;
00140     result = MBI()->get_connectivity( tri, conn, n_verts );
00141     assert(MB_SUCCESS == result);
00142     assert(3 == n_verts);
00143     for(int i=0; i<3; i++) print_vertex_coords( conn[i] );
00144     
00145     if(print_edges) {
00146       Range edges;
00147       result = MBI()->get_adjacencies( &tri, 1, 1, true, edges );
00148       if(MB_SUCCESS != result) std::cout << "result=" << result << std::endl;
00149       assert(MB_SUCCESS == result);
00150       //std::cout << "      edges: ";
00151       for(Range::iterator i=edges.begin(); i!=edges.end(); i++) {
00152         //std::cout << *i << " ";
00153         print_edge( *i );
00154       }
00155       //std::cout << std::endl;
00156     }
00157   }
00158 
00159   void print_edge( const EntityHandle edge ) {
00160     const EntityHandle *conn;
00161     int n_verts;
00162     std::cout << "    edge " << edge << std::endl;
00163     ErrorCode result = MBI()->get_connectivity( edge, conn, n_verts );
00164     if(gen::error(MB_SUCCESS!=result, "failed to get edge connectivity"));
00165     assert(MB_SUCCESS == result);
00166     assert(2 == n_verts);
00167     print_vertex_coords( conn[0] );   
00168     print_vertex_coords( conn[1] ); 
00169 
00170     Range tris;
00171     result = MBI()->get_adjacencies( &edge, 1, 2, false, tris );
00172     assert(MB_SUCCESS == result);
00173     std::cout << "     tris: ";
00174     for(Range::iterator i=tris.begin(); i!=tris.end(); i++) {
00175       std::cout << *i << " ";
00176     }
00177     std::cout << std::endl; 
00178   }
00179 
00180   void print_range( const Range range ) {
00181     std::cout << "print range:" << std::endl;
00182     Range::iterator i;
00183     for(i=range.begin(); i!=range.end(); i++) {
00184       std::cout << "    " << *i << std::endl;
00185     }
00186   }
00187 
00188   void print_range_of_edges( const Range range ) {
00189     std::cout << "print range:" << std::endl;
00190     Range::const_iterator i;
00191     for(i=range.begin(); i!=range.end(); i++) {
00192       print_edge( *i );
00193     }
00194   }
00195 
00196 
00197   void print_vertex_count(const EntityHandle input_meshset) {
00198   
00199     // get the range of facets of the surface meshset
00200     ErrorCode result;
00201     Range vertices;
00202     result = MBI()->get_entities_by_type(0, MBVERTEX, vertices);
00203     if(gen::error(MB_SUCCESS!=result, "failed to get vertex entities from the mesh"));
00204     assert( MB_SUCCESS == result );
00205   
00206     std::cout<< "    " << vertices.size() << " vertices found." << std::endl;
00207   }
00208 
00209   void print_arcs( const std::vector< std::vector<EntityHandle> > arcs ) {
00210     for(unsigned int i=0; i<arcs.size(); i++) {
00211       std::cout << "arc " << i << std::endl;
00212       print_loop( arcs[i] );
00213     }
00214   }
00215 
00216   void print_arc_of_edges( const std::vector<EntityHandle> arc_of_edges ) {
00217   
00218     ErrorCode result;
00219     std::vector<EntityHandle>::const_iterator i;
00220     double dist = 0;
00221     for( i=arc_of_edges.begin(); i!=arc_of_edges.end(); i++ ) {
00222       int n_verts;
00223       const EntityHandle *conn;
00224       result = MBI()->get_connectivity( *i, conn, n_verts ); 
00225       if(gen::error(MB_SUCCESS!=result, "failed to get edge connectivity"));
00226       assert(MB_SUCCESS == result);
00227       assert( 2 == n_verts );
00228       dist += dist_between_verts( conn[0], conn[1] );
00229       print_vertex_coords( conn[0] );
00230       print_vertex_coords( conn[1] );
00231     }
00232     std::cout << "  dist= " << dist << std::endl;
00233   }
00234 
00235   void print_loop( const std::vector<EntityHandle> loop_of_verts ) {
00236    
00237     std::cout << "  size=" << loop_of_verts.size() << std::endl;
00238     double dist = 0;
00239     //std::vector<EntityHandle>::iterator i;
00240     //for( i=loop_of_verts.begin(); i!=loop_of_verts.end(); i++ ) {
00241     for(unsigned int i=0; i<loop_of_verts.size(); i++) {
00242       print_vertex_coords( loop_of_verts[i] );
00243       if(i != loop_of_verts.size()-1) {
00244         dist += dist_between_verts( loop_of_verts[i], loop_of_verts[i+1] );
00245       }
00246     }
00247     std::cout << "  dist=" << dist << std::endl;
00248   }
00249 
00250 
00254 ErrorCode find_closest_vert( const EntityHandle reference_vert,
00255                                const std::vector<EntityHandle> arc_of_verts,
00256                                unsigned &position,
00257                                const double dist_limit ) {
00258   ErrorCode rval;
00259   const bool debug = false;
00260   double min_dist_sqr = std::numeric_limits<double>::max();
00261   CartVect ref_coords;
00262   rval = MBI()->get_coords( &reference_vert, 1, ref_coords.array() );
00263   if(gen::error(MB_SUCCESS!=rval,"failed to get ref coords")) return rval;
00264   double length = 0;
00265   CartVect prev_coords;
00266 
00267   for(unsigned i=0; i<arc_of_verts.size(); ++i) {
00268     CartVect coords;
00269     rval = MBI()->get_coords( &arc_of_verts[i], 1, coords.array() );
00270     if(gen::error(MB_SUCCESS!=rval,"failed to get coords")) return rval;
00271 
00272     // use dist_limit to exit early; avoid checking the entire arc
00273     if(0!=i) {
00274       CartVect temp = prev_coords - coords;
00275       length += temp.length();
00276       if(length>dist_limit && debug) 
00277         std::cout << "length=" << length << " dist_limit=" << dist_limit << std::endl;
00278       if(length > dist_limit) return MB_SUCCESS;
00279     }
00280     prev_coords = coords;
00281 
00282     // get distance to ref_vert
00283     CartVect temp = ref_coords - coords;
00284     double dist_sqr = temp.length_squared();
00285     if(dist_sqr < min_dist_sqr) {
00286       position = i;
00287       min_dist_sqr = dist_sqr;
00288       if(debug) std::cout << "min_dist_sqr=" << min_dist_sqr << std::endl;
00289     }
00290   }
00291  
00292   return MB_SUCCESS;
00293 }
00294 
00295 
00296 
00297   // Return the closest vert and all within tol. This is needed because sometimes
00298   // the correct vert is not the closest. For example, iter_surf4010 the skin
00299   // loop has the same point in it twice, at two different locations (center of L).
00300   // This ensure that both are returned as candidates.
00301   ErrorCode find_closest_vert( const double tol,
00302                                  const EntityHandle reference_vert,
00303                                  const std::vector<EntityHandle> loop_of_verts,
00304                                  std::vector<unsigned> &positions, 
00305                                  std::vector<double> &dists) {
00306 
00307     ErrorCode rval;
00308     positions.clear();
00309     dists.clear();
00310     const double TOL_SQR = tol*tol;
00311     unsigned min_pos;
00312     double sqr_min_dist = std::numeric_limits<double>::max();
00313     for(unsigned int i=0; i<loop_of_verts.size(); i++) {
00314         double sqr_dist = std::numeric_limits<double>::max();
00315         rval = squared_dist_between_verts(reference_vert, loop_of_verts[i], sqr_dist);
00316         if(gen::error(MB_SUCCESS!=rval,"could not get dist")) return rval;
00317         if(sqr_dist < sqr_min_dist) {
00318           if(sqr_dist >= TOL_SQR) {
00319             sqr_min_dist = sqr_dist;
00320             min_pos = i;
00321           } else {
00322             sqr_min_dist = TOL_SQR;
00323             positions.push_back(i);
00324             dists.push_back( sqrt(sqr_dist) );
00325           }
00326         }
00327     }
00328  
00329     if(dists.empty()) {
00330       dists.push_back( sqrt(sqr_min_dist) );
00331       positions.push_back( min_pos );
00332     }
00333 
00334             
00335    return MB_SUCCESS;
00336   }
00337 
00338   ErrorCode merge_vertices( Range verts /* in */, const double tol /* in */ ) {
00339 
00340     ErrorCode result;
00341     const double SQR_TOL = tol*tol;
00342     // Clean up the created tree, and track verts so that if merged away they are
00343     // removed from the tree.
00344     AdaptiveKDTree kdtree(MBI()); //, true, 0, MESHSET_TRACK_OWNER);
00345     // initialize the KD Tree
00346     EntityHandle root;
00347     const char settings[]="MAX_PER_LEAF=6;MAX_DEPTH=50;SPLITS_PER_DIR=1;PLANE_SET=2;MESHSET_FLAGS=0x1;TAG_NAME=0";
00348     FileOptions fileopts(settings);
00349 
00350 
00351     /* Old KDTree settings
00352     AdaptiveKDTree::Settings settings;
00353    
00354     // tells the tree to split leaves with more than 6 entities
00355     settings.maxEntPerLeaf = 6;
00356     // tells the tree a maximum depth limit
00357     settings.maxTreeDepth  = 50;
00358     // tells the tree how many candidate split planed to consider in each dimension
00359     settings.candidateSplitsPerDir = 1;
00360     // tells the tree to use the median vertex coordinate values to set planes
00361     settings.candidatePlaneSet = AdaptiveKDTree::VERTEX_MEDIAN;
00362    
00363     */
00364 
00365     // builds the KD Tree, making the EntityHandle root the root of the tree
00366     result = kdtree.build_tree( verts, &root, &fileopts);
00367     assert(MB_SUCCESS == result);
00368     // create tree iterator to loop over all verts in the tree
00369     AdaptiveKDTreeIter tree_iter;
00370     kdtree.get_tree_iterator( root, tree_iter );
00371  
00372     //for(unsigned int i=0; i<verts.size(); i++) {
00373     for(Range::iterator i=verts.begin(); i!=verts.end(); ++i) {
00374       double from_point[3];
00375       //EntityHandle vert = *i;
00376       result = MBI()->get_coords( &(*i), 1, from_point);
00377       assert(MB_SUCCESS == result);
00378       std::vector<EntityHandle> leaves_out;
00379       result = kdtree.distance_search( from_point, tol, leaves_out, root);
00380       assert(MB_SUCCESS == result);
00381       for(unsigned int j=0; j<leaves_out.size(); j++) {
00382         std::vector<EntityHandle> leaf_verts;
00383         result = MBI()->get_entities_by_type( leaves_out[j], MBVERTEX, leaf_verts);
00384         assert(MB_SUCCESS == result);
00385         if(100 < leaf_verts.size()) std::cout << "*i=" << *i << " leaf_verts.size()=" << leaf_verts.size() << std::endl;
00386         for(unsigned int k=0; k<leaf_verts.size(); k++) {
00387           if( leaf_verts[k] == *i ) continue; 
00388           double sqr_dist;
00389           result = gen::squared_dist_between_verts( *i, leaf_verts[k], sqr_dist);
00390           assert(MB_SUCCESS == result);
00391 
00392           if(SQR_TOL >= sqr_dist) {
00393             //  std::cout << "merge_vertices: vert " << leaf_verts[k] << " merged to vert "
00394             //            << verts[i] << " dist=" << dist << " leaf=" << std::endl;
00395 
00396             // The delete_vert is automatically remove from the tree because it
00397             // uses tracking meshsets. merge_verts checks for degenerate tris.
00398             // Update the list of leaf verts to prevent stale handles.
00399             std::vector<EntityHandle> temp_arc;
00400             EntityHandle keep_vert   = *i;
00401             EntityHandle delete_vert = leaf_verts[k];
00402             result = zip::merge_verts( keep_vert, delete_vert, leaf_verts, temp_arc );
00403             assert(MB_SUCCESS == result);
00404             // Erase delete_vert from verts
00405             // Iterator should remain valid because delete_vert > keep_vert handle.
00406             verts.erase( delete_vert );
00407           }
00408         }
00409       }
00410     }   
00411     return result;
00412   }
00413 
00414   ErrorCode squared_dist_between_verts( const EntityHandle v0,
00415                                           const EntityHandle v1,
00416                                           double &d) {
00417     ErrorCode result;
00418     CartVect coords0, coords1;
00419     result = MBI()->get_coords( &v0, 1, coords0.array() );
00420     if(MB_SUCCESS != result) {
00421       std::cout << "dist_between_verts: get_coords on v0=" << v0 << " result=" 
00422                 << result << std::endl; 
00423       return result;
00424     }
00425     result = MBI()->get_coords( &v1, 1, coords1.array() );
00426     if(MB_SUCCESS != result) {
00427       std::cout << "dist_between_verts: get_coords on v1=" << v1 << " result=" 
00428                 << result << std::endl; 
00429       return result;
00430     }
00431     const CartVect diff = coords0 - coords1;
00432     d = diff.length_squared();
00433     return MB_SUCCESS;
00434   }
00435 
00436   double dist_between_verts( const CartVect v0, const CartVect v1 ) {
00437     CartVect v2 = v0 - v1;
00438     return v2.length();
00439   }
00440   ErrorCode dist_between_verts( const EntityHandle v0, const EntityHandle v1, double &d) {
00441     ErrorCode result;
00442     CartVect coords0, coords1;
00443     result = MBI()->get_coords( &v0, 1, coords0.array() );
00444     if(MB_SUCCESS != result) {
00445       std::cout << "dist_between_verts: get_coords on v0=" << v0 << " result=" 
00446                 << result << std::endl; 
00447       return result;
00448     }
00449     result = MBI()->get_coords( &v1, 1, coords1.array() );
00450     if(MB_SUCCESS != result) {
00451       std::cout << "dist_between_verts: get_coords on v1=" << v1 << " result=" 
00452                 << result << std::endl; 
00453       return result;
00454     }
00455     d = dist_between_verts( coords0, coords1 );
00456     return MB_SUCCESS;
00457   }
00458 
00459   double dist_between_verts( double coords0[], double coords1[] ) {
00460     return sqrt( (coords0[0]-coords1[0])*(coords0[0]-coords1[0]) +
00461                  (coords0[1]-coords1[1])*(coords0[1]-coords1[1]) +
00462                  (coords0[2]-coords1[2])*(coords0[2]-coords1[2]) );
00463   }
00464   double dist_between_verts( EntityHandle vert0, EntityHandle vert1 ) {
00465     double coords0[3], coords1[3];
00466     ErrorCode result;
00467     result = MBI()->get_coords( &vert0, 1, coords0 );
00468     if(MB_SUCCESS!=result) std::cout << "result=" << result << " vert="
00469                                      << vert0 << std::endl;
00470     assert(MB_SUCCESS == result);
00471     result = MBI()->get_coords( &vert1, 1, coords1 );
00472     if(MB_SUCCESS!=result) std::cout << "result=" << result << " vert="
00473                                      << vert1 << std::endl;
00474     assert(MB_SUCCESS == result);
00475     return dist_between_verts( coords0, coords1 );
00476   }
00477 
00478   // Return the length of the curve defined by MBEDGEs or ordered MBVERTEXs.
00479   double length( std::vector<EntityHandle> edges ) {
00480     if(edges.empty()) return 0;
00481 
00482     ErrorCode result;
00483     std::vector<EntityHandle>::iterator i;
00484     double dist = 0;
00485     EntityType type = MBI()->type_from_handle( edges[0] );
00486 
00487     // if vector has both edges and verts, only use edges
00488     // NOTE: The curve sets from ReadCGM do not contain duplicate endpoints for loops!
00489     EntityType end_type = MBI()->type_from_handle( edges.back() );
00490     if(type != end_type) {
00491       for(std::vector<EntityHandle>::iterator i=edges.begin(); i!=edges.end(); i++) {
00492         if(MBVERTEX == MBI()->type_from_handle( *i )) {
00493           i = edges.erase(i) - 1;
00494         }
00495       }
00496     }
00497 
00498     // determine if vector defines an arc by edges of verts
00499     type = MBI()->type_from_handle( edges[0] ); 
00500     if        (MBEDGE == type) {
00501       if(edges.empty()) return 0.0; 
00502       for( i=edges.begin(); i!=edges.end(); i++ ) {
00503         int n_verts;
00504         const EntityHandle *conn;
00505         result = MBI()->get_connectivity( *i, conn, n_verts );
00506         if( MB_SUCCESS!=result ) std::cout << "result=" << result << std::endl;
00507         assert(MB_SUCCESS == result);
00508         assert( 2 == n_verts );
00509         if(conn[0] == conn[1]) continue;
00510         dist += dist_between_verts( conn[0], conn[1] );
00511         //std::cout << "length: " << dist << std::endl;
00512       }
00513     } else if (MBVERTEX == type) {
00514       if(2 > edges.size()) return 0.0;
00515       EntityHandle front_vert = edges.front();
00516       for( i=edges.begin()+1; i!=edges.end(); i++) {
00517         dist += dist_between_verts( front_vert, *i );
00518         front_vert = *i;
00519       }
00520     } else return MB_FAILURE;
00521 
00522     return dist;
00523   }
00524 
00525   // Given a vertex and vector of edges, return the number of edges adjacent to the vertex.
00526   unsigned int n_adj_edges( EntityHandle vert, Range edges ) {
00527     ErrorCode result;
00528     Range adj_edges;
00529     result = MBI()->get_adjacencies( &vert, 1, 1, false, adj_edges );
00530     gen::error(MB_SUCCESS!=result, "could not get edges adjacent to the vertex");
00531     assert(MB_SUCCESS == result);
00532     //adj_edges = adj_edges.intersect(edges);
00533     adj_edges = intersect( adj_edges, edges );
00534     return adj_edges.size();
00535   }
00536 
00537 
00538 
00539   // Return true if the edges share a vertex. Does not check for coincident edges.
00540   bool edges_adjacent( EntityHandle edge0, EntityHandle edge1 ) {
00541     ErrorCode result;
00542     Range verts0, verts1;
00543     result = MBI()->get_adjacencies( &edge0, 1, 0, false, verts0 );
00544     gen::error(MB_SUCCESS!=result, "could not get edge0 adjacencies");
00545     assert( MB_SUCCESS == result );
00546     assert( 2 == verts0.size() );
00547     result = MBI()->get_adjacencies( &edge1, 1, 0, false, verts1 );
00548     gen::error(MB_SUCCESS!=result, "could not get edge1 adjacencies");
00549     assert( MB_SUCCESS == result );
00550     assert( 2 == verts1.size() );
00551     if      ( verts0.front() == verts1.front() ) return true;
00552     else if ( verts0.front() == verts1.back()  ) return true;
00553     else if ( verts0.back()  == verts1.back()  ) return true;
00554     else if ( verts0.back()  == verts1.front() ) return true;
00555     else                                         return false;
00556   }
00557 
00558   // get the direction unit vector from one vertex to another vertex
00559   ErrorCode get_direction( const EntityHandle from_vert, const EntityHandle to_vert,
00560                              CartVect &dir ) {
00561     // double d[3];
00562     ErrorCode result;
00563     CartVect coords0, coords1;
00564     result = MBI()->get_coords( &from_vert, 1, coords0.array() );
00565     assert(MB_SUCCESS==result);
00566     result = MBI()->get_coords( &to_vert, 1, coords1.array() );
00567     assert(MB_SUCCESS==result);
00568     dir = coords1 - coords0;
00569     if(0 == dir.length()) {
00570       CartVect zero_vector( 0.0 );
00571       dir = zero_vector;
00572       std::cout << "direction vector has 0 magnitude" << std::endl;
00573       return MB_SUCCESS;
00574     }
00575     dir.normalize();
00576     return result;
00577   }    
00578  
00579   // from http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry1
00580   double edge_point_dist( const CartVect a, const CartVect b, const CartVect c ) {
00581     CartVect ab, bc, ba, ac;
00582     ab = b - a;
00583     bc = c - b;
00584     ba = a - b;
00585     ac = c - a;
00586   
00587     // find the magnitude of the cross product and test the line
00588     CartVect cross_product = ab*ac;
00589     double dist = cross_product.length() / dist_between_verts(a,b);
00590  
00591     // test endpoint1
00592     if (ab%bc > 0) {
00593       //std::cout << "edge_point_dist=" << dist_between_verts(b,c) 
00594       //<< " at endpt1" << std::endl;
00595       return dist_between_verts(b,c);
00596     }
00597 
00598     // test endpoint0
00599     if (ba%ac > 0) {
00600       //std::cout << "edge_point_dist=" << dist_between_verts(a,c) 
00601       //<< " at endpt0" << std::endl;
00602       return dist_between_verts(a,c);
00603     }
00604  
00605     //std::cout << "edge_point_dist=" << fabs(dist) << " at middle" 
00606     //<< std::endl;      
00607     return fabs(dist);
00608   }
00609   double edge_point_dist( const EntityHandle endpt0, const EntityHandle endpt1,
00610                           const EntityHandle pt ) {
00611     ErrorCode result;
00612     CartVect a, b, c;
00613     result = MBI()->get_coords( &endpt0, 1, a.array() );
00614     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00615     assert(MB_SUCCESS==result);
00616     result = MBI()->get_coords( &endpt1, 1, b.array() );
00617     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00618     assert(MB_SUCCESS==result);
00619     result = MBI()->get_coords( &pt,     1, c.array() );
00620     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00621     assert(MB_SUCCESS==result);
00622     return edge_point_dist( a, b, c);
00623   }
00624   double edge_point_dist( const EntityHandle edge, const EntityHandle pt ) {
00625     ErrorCode result;
00626     const EntityHandle *conn;
00627     int n_verts;     
00628     result = MBI()->get_connectivity( edge, conn, n_verts );
00629     gen::error(MB_SUCCESS!=result, "could not get edge connectivity");
00630     assert(MB_SUCCESS==result);
00631     assert( 2 == n_verts );
00632     return edge_point_dist( conn[0], conn[1], pt );
00633   }
00634   /*
00635   ErrorCode  point_curve_min_dist( const std::vector<EntityHandle> curve, // of verts
00636                                      const EntityHandle pt,
00637                                      double &min_dist,
00638                                      const double max_dist_along_curve ) { 
00639   
00640     min_dist = std::numeric_limits<double>::max();
00641     double cumulative_dist = 0;
00642     bool last_edge = false;
00643   
00644     // it is a curve of verts or a curve of edges?
00645     EntityType type = MBI()->type_from_handle( curve.front() );
00646     std::vector<EntityHandle>::const_iterator i;
00647     if(MBVERTEX == type) {
00648       EntityHandle front_vert = curve.front();
00649       for( i=curve.begin()+1; i!=curve.end(); i++) {
00650         // if we are using verts, do not explicitly create the edge in MOAB
00651         cumulative_dist += gen::dist_between_verts( front_vert, *i );
00652         //std::cout << "  point_curve_min_dist: cumulative_dist="
00653         //        << cumulative_dist << std::endl;
00654         double d = edge_point_dist( front_vert, *i, pt );
00655         //std::cout << "  point_curve_min_dist: d=" << d << std::endl;        
00656         if(d < min_dist) {
00657           min_dist = d;
00658           //std::cout << "min_dist=" << min_dist << std::endl;
00659           //print_vertex_coords( front_vert );
00660           //print_vertex_coords( *i );
00661         }
00662         // check one edge past the point after max_dist_along_curve
00663         if(last_edge) return MB_SUCCESS;
00664         if(max_dist_along_curve<cumulative_dist) last_edge = true;
00665       
00666         front_vert = *i;
00667       }
00668     } //else if(MBEDGE == type) {
00669       //for( i=curve.begin(); i!=curve.end(); i++) {
00670 //      double d = edge_point_dist( *i, pt );
00671         //if(d < min_dist) min_dist = d;
00672       //}
00673      // } 
00674        else return MB_FAILURE;
00675     //std::cout << "point_curve_min_dist=" << min_dist << " curve.size()=" << curve.size() << std::endl;    
00676     return MB_SUCCESS;
00677   }
00678 
00679   ErrorCode  point_curve_min_dist( const std::vector<EntityHandle> curve, // of verts
00680                                      const EntityHandle pt,
00681                                      double &min_dist ) { 
00682     const double max_dist_along_curve = std::numeric_limits<double>::max();
00683     return point_curve_min_dist( curve, pt, min_dist, max_dist_along_curve ); 
00684   }
00685 */
00686   double triangle_area( const CartVect a, const CartVect b,
00687                         const CartVect c) {
00688     CartVect d = c - a;
00689     CartVect e = c - b;
00690     CartVect f = d*e;
00691     return 0.5*f.length();
00692   }
00693   ErrorCode triangle_area( const EntityHandle conn[], double &area ) {
00694     CartVect coords[3];
00695     ErrorCode result = MBI()->get_coords( conn, 3, coords[0].array() );
00696     gen::error(MB_SUCCESS!=result, "could not get triangle vertex coordinates");
00697     assert(MB_SUCCESS == result);
00698     area = triangle_area( coords[0], coords[1], coords[2] );
00699     return result;
00700   }
00701   ErrorCode triangle_area( const EntityHandle tri, double &area ) {
00702     ErrorCode result;
00703     const EntityHandle *conn;
00704     int n_verts;
00705     result = MBI()->get_connectivity( tri, conn, n_verts );
00706     gen::error(MB_SUCCESS!=result, "could not get trangle vertices");
00707     assert(MB_SUCCESS == result);
00708     assert(3 == n_verts);
00709     
00710     result = triangle_area( conn, area );
00711     gen::error(MB_SUCCESS!=result, "could not get triangle area");
00712     assert(MB_SUCCESS == result);
00713     return result;
00714   }
00715   double triangle_area( const Range tris ) {
00716     double a, area = 0;
00717     ErrorCode result;
00718     for(Range::iterator i=tris.begin(); i!=tris.end(); i++) {
00719       result = triangle_area( *i, a);
00720       gen::error(MB_SUCCESS!=result, "could not get triangle area");
00721       assert(MB_SUCCESS == result);
00722       area += a;
00723     }
00724     return area;
00725   }
00726 
00727   bool triangle_degenerate( const EntityHandle tri ) {
00728     ErrorCode result;
00729     const EntityHandle *conn;
00730     int n_verts;
00731     result = MBI()->get_connectivity( tri, conn, n_verts );
00732     gen::error(MB_SUCCESS!=result, "could not get triangle vertices");
00733     assert(MB_SUCCESS == result);
00734     assert(3 == n_verts);
00735     return triangle_degenerate( conn[0], conn[1], conn[2] );
00736   }
00737 
00738   bool triangle_degenerate( const EntityHandle v0, const EntityHandle v1,
00739                             const EntityHandle v2 ) {
00740     if(v0==v1 || v1==v2 || v2==v0) return true;
00741     return false;
00742   }
00743 
00744   ErrorCode triangle_normals( const Range tris, std::vector<CartVect> &normals ) {
00745     ErrorCode result;
00746     normals.clear();
00747     for(Range::const_iterator i=tris.begin(); i!=tris.end(); i++) {
00748       CartVect normal;
00749       result = triangle_normal( *i, normal );
00750       gen::error(MB_SUCCESS!=result, "could not get triangle normal vector");
00751       assert(MB_SUCCESS==result || MB_ENTITY_NOT_FOUND==result);
00752       normals.push_back( normal );
00753     }
00754     // if we've gotten here, then we have succeeded
00755     result = MB_SUCCESS;
00756     return result;
00757   }
00758 
00759   ErrorCode triangle_normal( const EntityHandle tri, CartVect &normal) {
00760     ErrorCode result;
00761     const EntityHandle *conn;
00762     int n_verts;
00763     result = MBI()->get_connectivity( tri, conn, n_verts );
00764     if(MB_ENTITY_NOT_FOUND == result) {
00765       std::cout << "triangle_normal: triangle not found" << std::endl;
00766       CartVect zero_vector( 0.0 );
00767       normal = zero_vector;
00768       return result;
00769     }else if(MB_SUCCESS != result) {
00770       return result;
00771     } else {
00772       assert(3 == n_verts);
00773       return triangle_normal( conn[0], conn[1], conn[2], normal );
00774     }
00775   }
00776 
00777   ErrorCode triangle_normal( const EntityHandle v0, const EntityHandle v1,
00778                                const EntityHandle v2, CartVect &normal ) {
00779 
00780     // if tri is degenerate return 0,0,0
00781     if( triangle_degenerate(v0, v1, v2) ) {
00782       CartVect zero_vector( 0.0 );
00783       normal = zero_vector;
00784       std::cout << "  normal=" << normal << std::endl;
00785       return MB_SUCCESS;
00786     }
00787 
00788     EntityHandle conn[3];
00789     conn[0] = v0;
00790     conn[1] = v1;
00791     conn[2] = v2;
00792     ErrorCode result;
00793     CartVect coords[3];
00794     result = MBI()->get_coords( conn, 3, coords[0].array() );
00795     gen::error(MB_SUCCESS!=result, "could not get coordinates of the triangle vertices");
00796     assert(MB_SUCCESS == result);
00797     return triangle_normal( coords[0], coords[1], coords[2], normal );
00798   }
00799 
00800   ErrorCode triangle_normal( const CartVect coords0, const CartVect coords1,
00801                                const CartVect coords2, CartVect &normal ) {
00802     CartVect edge0, edge1;
00803     edge0 = coords1-coords0;
00804     edge1 = coords2-coords0;
00805     normal = edge0*edge1;
00806 
00807     // do not normalize if magnitude is zero (avoid nans)
00808     if(0 == normal.length()) return MB_SUCCESS;
00809 
00810     normal.normalize();
00811     //if(debug) std::cout << "  normal=" << normal << std::endl;
00812     return MB_SUCCESS;
00813   }
00814     
00815 
00816 
00817   // Distance between a point and line. The line is defined by two verts.
00818   // We are using a line and not a line segment!
00819   // http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
00820   ErrorCode line_point_dist( const EntityHandle line_pt1, const EntityHandle line_pt2,
00821                                const EntityHandle pt0, double &dist ) {
00822     ErrorCode result;
00823     CartVect x0, x1, x2;
00824     result = MBI()->get_coords( &line_pt1, 1, x1.array() );
00825     assert(MB_SUCCESS == result);
00826     result = MBI()->get_coords( &line_pt2, 1, x2.array() );
00827     assert(MB_SUCCESS == result);
00828     result = MBI()->get_coords( &pt0, 1, x0.array() );
00829     assert(MB_SUCCESS == result);
00830 
00831     dist = ( ((x0-x1)*(x0-x2)).length() ) / ( (x2-x1).length() );
00832     return result;
00833   }
00834 
00835   // Project the point onto the line. Not the line segment!
00836   ErrorCode point_line_projection( const EntityHandle line_pt1,
00837                                      const EntityHandle line_pt2,
00838                                      const EntityHandle pt0 ) {
00839     CartVect projected_coords;
00840     double parameter;
00841     ErrorCode result = point_line_projection( line_pt1, line_pt2,
00842                                                 pt0, projected_coords,
00843                                                 parameter );
00844     assert(MB_SUCCESS == result);
00845     result = MBI()->set_coords( &pt0, 1, projected_coords.array() );   
00846     assert(MB_SUCCESS == result);
00847     return result;
00848   }    
00849   ErrorCode point_line_projection( const EntityHandle line_pt1,
00850                                      const EntityHandle line_pt2,
00851                                      const EntityHandle pt0,
00852                                      CartVect &projected_coords,
00853                                      double &parameter  ) {
00854 
00855     ErrorCode result;
00856     CartVect coords[3];
00857     result = MBI()->get_coords( &line_pt1, 1, coords[1].array() );
00858     assert(MB_SUCCESS == result);
00859     result = MBI()->get_coords( &line_pt2, 1, coords[2].array() );           
00860     assert(MB_SUCCESS == result);
00861     result = MBI()->get_coords( &pt0, 1, coords[0].array() );          
00862     assert(MB_SUCCESS == result);
00863 
00864     // project the t_joint between the endpts                  
00865     // http://en.wikipedia.org/wiki/Vector_projection           
00866     CartVect a = coords[0] - coords[1];
00867     CartVect b = coords[2] - coords[1];
00868     parameter    = (a%b)/(b%b);
00869     CartVect c = parameter*b;
00870     projected_coords = c     + coords[1];      
00871     return result;
00872   }
00873   ErrorCode point_line_projection( const EntityHandle line_pt1,
00874                                      const EntityHandle line_pt2,
00875                                      const EntityHandle pt0,
00876                                      double &dist_along_edge  ) {
00877 
00878     ErrorCode result;
00879     CartVect coords[3];
00880     result = MBI()->get_coords( &line_pt1, 1, coords[1].array() );
00881     assert(MB_SUCCESS == result);
00882     result = MBI()->get_coords( &line_pt2, 1, coords[2].array() );           
00883     assert(MB_SUCCESS == result);
00884     result = MBI()->get_coords( &pt0, 1, coords[0].array() );          
00885     assert(MB_SUCCESS == result);
00886 
00887     // project the t_joint between the endpts                  
00888     // http://en.wikipedia.org/wiki/Vector_projection           
00889     CartVect a = coords[0] - coords[1];
00890     CartVect b = coords[2] - coords[1];
00891     dist_along_edge = a%b / b.length();      
00892     return result;
00893   }
00894   
00895 
00896   double area2( const EntityHandle pt_a, const EntityHandle pt_b,
00897                 const EntityHandle pt_c, const CartVect plane_normal ) {
00898     //std::cout << "area2: a=" << pt_a << " b=" << pt_b << " c=" << pt_c << std::endl;
00899     ErrorCode result;
00900     CartVect a, b, c;
00901     result = MBI()->get_coords( &pt_a, 1, a.array() );
00902     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00903     assert(MB_SUCCESS == result);
00904     result = MBI()->get_coords( &pt_b, 1, b.array() );
00905     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00906     assert(MB_SUCCESS == result);
00907     result = MBI()->get_coords( &pt_c, 1, c.array() );
00908     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00909     assert(MB_SUCCESS == result);
00910     CartVect d = b - a;
00911     CartVect e = c - a;
00912     // project onto a plane defined by the plane's normal vector
00913 
00914     return (d*e)%plane_normal;
00915   }
00916 
00917   // Is point c to the left of line ab?
00918   bool left( const EntityHandle a, const EntityHandle b,
00919              const EntityHandle c, const CartVect n ) {
00920     double area_2 = area2(a,b,c,n);
00921     //std::cout << "left: a=" << a << " b=" << b << " c=" << c
00922     //          << " area2=" << area_2 << std::endl;
00923     if(area_2 > 0) return true;
00924     else return false;
00925   } 
00926 
00927   // Is point c to the left of line ab or collinear?
00928   bool left_on( const EntityHandle a, const EntityHandle b,
00929                 const EntityHandle c, const CartVect n ) {
00930     double area_2 = area2(a,b,c,n);
00931     //std::cout << "left_on: a=" << a << " b=" << b << " c=" << c
00932     //          << " area2=" << area_2 << std::endl;
00933     if(area_2 >= 0) return true;
00934     else return false;
00935   } 
00936 
00937   // Are pts a,b,c collinear?
00938   bool collinear( const EntityHandle a, const EntityHandle b,
00939                   const EntityHandle c, const CartVect n ) {
00940     double area_2 = area2(a,b,c,n);
00941     //std::cout << "collinear: a=" << a << " b=" << b << " c=" << c
00942     //          << " area2=" << area_2 << std::endl;
00943     if( area_2 ==0) return true;
00944     else return false;
00945   } 
00946 
00947   // Exclusive or: T iff exactly one argument is true
00948   bool logical_xor( const bool x, const bool y ) {
00949     return (x || y) && !(x && y);
00950   }
00951 
00952   bool intersect_prop( const EntityHandle a, const EntityHandle b,
00953                        const EntityHandle c, const EntityHandle d,
00954                        const CartVect n ) {
00955     if( collinear(a,b,c,n) ||
00956         collinear(a,b,d,n) ||
00957         collinear(c,d,a,n) ||
00958         collinear(c,d,b,n) ) {
00959       return false;
00960     } else {
00961       return logical_xor(left(a,b,c,n), left(a,b,d,n)) && 
00962         logical_xor(left(c,d,a,n), left(c,d,b,n));
00963     }
00964   }
00965       
00966   bool between( const EntityHandle pt_a, const EntityHandle pt_b,
00967                 const EntityHandle pt_c, const CartVect n) {
00968     if( !collinear(pt_a,pt_b,pt_c,n) ) return false;
00969 
00970     ErrorCode result;
00971     CartVect a, b, c;
00972     result = MBI()->get_coords( &pt_a, 1, a.array() );
00973     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00974     assert(MB_SUCCESS == result);
00975     result = MBI()->get_coords( &pt_b, 1, b.array() );
00976     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00977     assert(MB_SUCCESS == result);
00978     result = MBI()->get_coords( &pt_c, 1, c.array() );
00979     gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
00980     assert(MB_SUCCESS == result);
00981 
00982     // if ab not vertical, check betweenness on x; else on y.
00983     if(a[0] != b[0]) {
00984       return ((a[0] <= c[0]) && (c[0] <= b[0])) || ((a[0] >= c[0]) && (c[0] >= b[0]));
00985     }else if(a[1] != b[1]) {
00986       return ((a[1] <= c[1]) && (c[1] <= b[1])) || ((a[1] >= c[1]) && (c[1] >= b[1]));
00987     } else {
00988       return ((a[2] <= c[2]) && (c[2] <= b[2])) || ((a[2] >= c[2]) && (c[2] >= b[2]));
00989     }
00990   }
00991 
00992   bool intersect( const EntityHandle a, const EntityHandle b,
00993                   const EntityHandle c, const EntityHandle d,
00994                   const CartVect n ) {
00995     if(intersect_prop(a,b,c,d,n)) return true;
00996     else if( between(a,b,c,n) ||
00997              between(a,b,d,n) ||
00998              between(c,d,a,n) ||
00999              between(c,d,b,n) ) return true;
01000     else return false;
01001   }
01002 
01003   // verts is an ordered polygon of verts
01004   bool diagonalie( const EntityHandle a, const EntityHandle b,
01005                    const CartVect n,
01006                    const std::vector<EntityHandle> verts ) {
01007     for(unsigned int i=0; i<verts.size(); i++) {
01008       EntityHandle c = verts[i];
01009       EntityHandle c1;
01010       if(verts.size()-1 == i) c1 = verts[0];
01011       else                    c1 = verts[i+1];
01012 
01013       if( (c != a) && (c1 != a) &&
01014           (c != b) && (c1 != b) &&
01015           intersect( a, b, c, c1, n ) ) {
01016         //std::cout << "diagonalie a=" << a << " b=" << b << " c="
01017         //        << c << " c1=" << c1 << " result=";
01018         //std::cout << "false" << std::endl;
01019         return false;
01020       }
01021     }
01022     //std::cout << "diagonalie a=" << a << " b=" << b << " result=";
01023     //std::cout << "true" << std::endl;
01024     return true;
01025   }
01026 
01027   // verts is an ordered polygon of verts
01028   bool in_cone( const EntityHandle a, const EntityHandle b,
01029                 const CartVect n,
01030                 const std::vector<EntityHandle> verts ) {
01031     std::vector<EntityHandle>::const_iterator a_iter;
01032     a_iter = find( verts.begin(), verts.end(), a );
01033     EntityHandle a0, a1;
01034     // a0 is before a
01035     if(verts.begin() == a_iter) a0 = verts[verts.size()-1];
01036     else a0 = *(a_iter-1);
01037     // a1 is after a
01038     if(verts.end()-1 == a_iter) a1 = verts[0];
01039     else a1 = *(a_iter+1);
01040 
01041     //std::cout << "in_cone: a=" << a << " b=" << b << " a0="
01042     //        << a0 << " a1=" << a1 << std::endl;
01043     // if a is a convex vertex
01044     if(left_on(a,a1,a0,n)) return left(a,b,a0,n) && left(b,a,a1,n);
01045 
01046     // else a is reflex
01047     else return !(left_on(a,b,a1,n) && left_on(b,a,a0,n));
01048   }
01049 
01050   bool diagonal( const EntityHandle a, const EntityHandle b,
01051                  const CartVect n,
01052                  const std::vector<EntityHandle> verts ) {
01053     bool result = in_cone(a,b,n,verts) && in_cone(b,a,n,verts) && diagonalie(a,b,n,verts);
01054     //std::cout << "diagonal a=" << a << " b=" << b << " result="
01055     //          << result << std::endl;
01056     return result;
01057   }
01058 
01059 
01060   // Determine if each vertex is an ear. Input an ordered polygon of verts.
01061   ErrorCode ear_init( const std::vector<EntityHandle> verts,
01062                         const CartVect n, // plane normal vector
01063                         std::vector<bool> &is_ear ) {
01064     if(verts.size() != is_ear.size()) return MB_FAILURE;
01065     for(unsigned int i=0; i<verts.size(); i++) {
01066       EntityHandle prev, next;
01067       if(0 == i) prev = verts.back();
01068       else prev = verts[i-1];
01069       if(verts.size()-1 == i) next = verts[0];
01070       else next = verts[i+1];
01071       is_ear[i] = diagonal(prev,next,n,verts);
01072       //std::cout << "is_ear[" << i << "]=" << is_ear[i] << std::endl;
01073     }
01074     return MB_SUCCESS;
01075   }
01076 
01077 
01078 
01079   // Input an ordered polygon of verts and a normal vector of the plane
01080   // that the polygon is mostly in. The vector is required for orientation.
01081   ErrorCode ear_clip_polygon( std::vector<EntityHandle> verts,
01082                                 CartVect n,
01083                                 Range &new_tris ) {
01084 
01085     // initialize the status of ears
01086     //std::cout << "begin ear clipping----------------------" << std::endl;
01087     //for(unsigned int i=0; i<verts.size(); i++) {
01088     //  print_vertex_cubit( verts[i] );
01089     //}
01090 
01091     //print_loop( verts );
01092     ErrorCode result;
01093     std::vector<bool> is_ear( verts.size() );
01094     result = ear_init( verts, n, is_ear );
01095     assert(MB_SUCCESS == result);
01096  
01097     // if the polygon intersects itself the algorithm will not stop
01098     int counter = 0; 
01099     int n_initial_verts = verts.size();
01100 
01101     while(3 < verts.size()) {
01102       for(unsigned int i=0; i<verts.size(); i++) {
01103         if(is_ear[i]) {
01104           EntityHandle v0, v1, v2, v3, v4;
01105           if(0 == i)      v0 = verts[verts.size()-2];
01106           else if(1 == i) v0 = verts[verts.size()-1];
01107           else            v0 = verts[i-2]; 
01108           if(0 == i) v1 = verts[verts.size()-1];
01109           else       v1 = verts[i-1];
01110           v2 = verts[i];
01111           if(verts.size()-1 == i) v3 = verts[0];
01112           else                    v3 = verts[i+1];
01113           if(verts.size()-2 == i)       v4 = verts[0];
01114           else if (verts.size()-1 == i) v4 = verts[1];
01115           else                          v4 = verts[i+2];
01116 
01117           //std::cout << "ear_clip_polygon: triangle=" << std::endl;
01118           //print_vertex_coords( v1 ); 
01119           //print_vertex_coords( v2 ); 
01120           //print_vertex_coords( v3 );
01121           EntityHandle new_tri;
01122           EntityHandle conn[3] = {v1,v2,v3};
01123           result = MBI()->create_element( MBTRI, conn, 3, new_tri );
01124           assert(MB_SUCCESS == result);
01125           new_tris.insert( new_tri );
01126 
01127           // update ear status
01128           if(0 == i) is_ear[verts.size()-1] = diagonal( v0, v3, n, verts );
01129           else       is_ear[i-1]            = diagonal( v0, v3, n, verts );
01130           if(verts.size()-1 == i) is_ear[0] = diagonal( v1, v4, n, verts );
01131           else                    is_ear[i+1] = diagonal( v1, v4, n, verts );
01132 
01133           // cut off the ear at i
01134           verts.erase( verts.begin()+i );
01135           is_ear.erase( is_ear.begin()+i );
01136           //i--;
01137           break;
01138         }
01139       }
01140 
01141       // If the polygon has intersecting edges this loop will continue until it
01142       // hits this return.
01143       //std::cout << "counter=" << counter << " verts.size()=" << verts.size() << std::endl;
01144       if(counter > n_initial_verts) {
01145         //std::cout << "ear_clip_polygon: no ears found" << std::endl;
01146         result = MBI()->delete_entities( new_tris );
01147         assert(MB_SUCCESS == result);
01148         new_tris.clear();
01149         return MB_FAILURE;
01150       }
01151       counter++;
01152     } 
01153     //std::cout << "ear_clip_polygon: triangle=" << std::endl;
01154     //print_vertex_coords( verts[0] ); 
01155     //print_vertex_coords( verts[1] ); 
01156     //print_vertex_coords( verts[2] );
01157     EntityHandle new_tri;
01158     EntityHandle conn[3] = {verts[0],verts[1],verts[2]};
01159     result = MBI()->create_element( MBTRI, conn, 3, new_tri );
01160     assert(MB_SUCCESS == result);
01161     new_tris.insert( new_tri );
01162     
01163     return result; 
01164   }
01165 
01166   int geom_id_by_handle( const EntityHandle set ) {
01167     ErrorCode result;
01168     Tag id_tag;
01169     //result = MBI()->tag_create( GLOBAL_ID_TAG_NAME, sizeof(int), MB_TAG_DENSE,
01170     //                                MB_TYPE_INTEGER, id_tag, 0, true );           
01171     result = MBI()->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER,id_tag,MB_TAG_DENSE);
01172     gen::error(MB_SUCCESS!=result && MB_ALREADY_ALLOCATED != result, "could not get the tag handle");
01173     assert(MB_SUCCESS==result || MB_ALREADY_ALLOCATED==result);
01174     int id;
01175     result = MBI()->tag_get_data( id_tag, &set, 1, &id );                  
01176     //assert(MB_SUCCESS == result);
01177     return id;
01178   }
01179   
01180   ErrorCode save_normals( Range tris, Tag normal_tag ) {
01181     std::vector<CartVect> normals(tris.size());
01182     ErrorCode result;
01183     result = triangle_normals( tris, normals );
01184     gen::error(MB_SUCCESS!=result, "could not get triangle normals");
01185     assert(MB_SUCCESS == result);
01186 
01187     result = MBI()->tag_set_data(normal_tag, tris, &normals[0]);
01188     gen::error(MB_SUCCESS!=result, "could not get the tag set data");
01189     assert(MB_SUCCESS == result);
01190     return result;
01191   }
01192 
01193   ErrorCode flip(const EntityHandle tri, const EntityHandle vert0,
01194                    const EntityHandle vert2, const EntityHandle surf_set) {
01195 
01196     // get the triangles in the surface. The tri and adj_tri must be in the surface.
01197     Range surf_tris;
01198     ErrorCode result = MBI()->get_entities_by_type( surf_set, MBTRI, surf_tris);
01199     assert(MB_SUCCESS == result);
01200 
01201     // get the triangle across the edge that will be flipped
01202     Range adj_tri;
01203     EntityHandle edge[2] = {vert0, vert2};
01204     result = MBI()->get_adjacencies( edge, 2, 2, false, adj_tri );
01205     assert(MB_SUCCESS == result);
01206     adj_tri = intersect(adj_tri, surf_tris);
01207     assert(2 == adj_tri.size());
01208     adj_tri.erase( tri );
01209     print_triangle( adj_tri.front(), false );
01210     //result = MBI()->list_entity( adj_tri.front() );
01211     //assert(MB_SUCCESS == result);
01212 
01213     // get the remaining tri vert
01214     Range tri_verts;
01215     result = MBI()->get_adjacencies( &tri, 1, 0, false, tri_verts );
01216     assert(MB_SUCCESS == result);
01217     assert(3 == tri_verts.size());
01218     tri_verts.erase(vert0);
01219     tri_verts.erase(vert2);
01220     assert(1 == tri_verts.size());
01221     EntityHandle vert1 = tri_verts.front();
01222 
01223     // get the remaining adj_tri vert
01224     Range adj_tri_verts;
01225     result = MBI()->get_adjacencies( &adj_tri.front(), 1, 0, false, adj_tri_verts );
01226     assert(MB_SUCCESS == result);
01227     assert(3 == adj_tri_verts.size());
01228     adj_tri_verts.erase(vert0);
01229     adj_tri_verts.erase(vert2);
01230     assert(1 == adj_tri_verts.size());
01231     EntityHandle vert3 = adj_tri_verts.front();
01232 
01233     // original tri_conn    = {vert0, vert1, vert2}
01234     // original adj_tri_conn= {vert2, vert3, vert0}
01235  
01236     // set the new connectivity
01237     EntityHandle tri_conn[3] = {vert0, vert1, vert3};
01238     result = MBI()->set_connectivity( tri, tri_conn, 3 );
01239     assert(MB_SUCCESS == result);
01240     EntityHandle adj_tri_conn[3] = {vert1, vert2, vert3};
01241     result = MBI()->set_connectivity( adj_tri.front(), adj_tri_conn, 3 );
01242     assert(MB_SUCCESS == result);
01243     print_triangle( tri, false );
01244     print_triangle( adj_tri.front(), false );
01245     return result;
01246   }
01247 
01248   ErrorCode ordered_verts_from_ordered_edges( const std::vector<EntityHandle> ordered_edges,
01249                                                 std::vector<EntityHandle> &ordered_verts ) {
01250     ErrorCode result = MB_SUCCESS;
01251     ordered_verts.clear();
01252     ordered_verts.reserve(ordered_edges.size()+1);
01253 
01254     // Save the back of the previous edge to check for continuity.
01255     EntityHandle previous_back_vert = 0;
01256     
01257     for(std::vector<EntityHandle>::const_iterator i=ordered_edges.begin();
01258         i!=ordered_edges.end(); i++) {
01259       const EntityHandle *conn;
01260       int n_verts;
01261       result = MBI()->get_connectivity( *i, conn, n_verts);
01262       gen::error(MB_SUCCESS!=result, "could not get edge connectivity");
01263       assert(MB_SUCCESS == result);
01264       assert(2 == n_verts);
01265       if(ordered_edges.begin() == i) {
01266         ordered_verts.push_back(conn[0]);
01267       } else {
01268         gen::error(previous_back_vert!=conn[0], "edges are not adjacent");
01269         assert(previous_back_vert == conn[0]);
01270       }
01271       ordered_verts.push_back(conn[1]);
01272       previous_back_vert = conn[1];
01273     }
01274     return result;
01275   }
01276 
01277   /* Find the distance between two arcs. Assume that their endpoints are somewhat
01278      close together. */
01279   ErrorCode dist_between_arcs( bool debug,
01280                                  const std::vector<EntityHandle> arc0,
01281                                  const std::vector<EntityHandle> arc1,
01282                                  double &dist ) {
01283     dist = 0;
01284 
01285     // Special Case: arcs have no verts.
01286     if( arc0.empty() || arc1.empty() ) {
01287       std::cout << "arc has no vertices" << std::endl;
01288       return MB_FAILURE;
01289     }
01290 
01291     //print_loop(arc0);
01292     //print_loop(arc1);
01293 
01294     // for simplicity, put arcs into the same structure
01295     std::vector<EntityHandle> arcs[2] = {arc0, arc1};
01296 
01297     // Special Case: Remove duplicate vert handles
01298     for(unsigned int i=0; i<2; ++i) {
01299       if( 2>arcs[i].size() ) continue;
01300       for(std::vector<EntityHandle>::iterator j=arcs[i].begin()+1; j!=arcs[i].end(); ++j) {
01301         if(*j == *(j-1)) {
01302           if(debug) {
01303             gen::print_loop( arcs[i] );
01304             std::cout << "dist_between_arcs: found duplicate vert handle in arc" << std::endl;
01305           }
01306           j = arcs[i].erase(j) - 1;
01307         }
01308       }
01309     }
01310     
01311     // get the coords in one call per arc. For speed, do not ask MOAB again for coords.
01312     ErrorCode result;
01313     std::vector<CartVect> coords[2];
01314     for(unsigned int i=0; i<2; i++) {
01315       coords[i].resize( arcs[i].size() );
01316       result = MBI()->get_coords( &arcs[i][0], arcs[i].size(), coords[i][0].array());
01317       gen::error(MB_SUCCESS!=result, "could not get vertex coordinates");
01318       assert(MB_SUCCESS == result);
01319     }
01320 
01321     // Special case: arc has 1 vert or a length of zero
01322     bool point_arc_exists = false;
01323     unsigned int point_arc_index;
01324     for(unsigned int i=0; i<2; ++i) {
01325       if( 1==arcs[i].size() ) {
01326         point_arc_exists = true;
01327         point_arc_index  = i;
01328         break;
01329       } else if ( 0==length(arcs[i]) ) {
01330         point_arc_exists = true;
01331         point_arc_index  = i;
01332         break;
01333       }
01334     }
01335     // If the special case exists, we can still get an average distance
01336     if(point_arc_exists) {
01337       unsigned int other_arc_index = (0==point_arc_index) ? 1 : 0;
01338       // Both are point arcs
01339       if(1==arcs[other_arc_index].size()) {
01340         dist = dist_between_verts( arcs[point_arc_index].front(), arcs[other_arc_index].front());
01341         return MB_SUCCESS;
01342         // The other arc has more than one point
01343       } else {
01344         double area = 0.0;
01345         for(unsigned int i=0; i<arcs[other_arc_index].size()-1; ++i) {
01346           area += fabs(triangle_area( coords[other_arc_index][i], 
01347                                       coords[other_arc_index][i+1], 
01348                                       coords[point_arc_index].front() ));
01349         }
01350         dist = area / gen::length(arcs[other_arc_index]);
01351         return MB_SUCCESS;
01352       }
01353     }
01354 
01355     // get the arc length and parametric coordinates. The parametrize the arcs by
01356     // length from 0 to 1.
01357     double arc_len[2] = {0.0};
01358     std::vector<double> params[2];
01359     for(unsigned int i=0; i<2; i++) {
01360       params[i].resize( arcs[i].size() );
01361       params[i][0] = 0.0;
01362       for(unsigned int j=0; j<coords[i].size()-1; j++) {
01363         double d = dist_between_verts( coords[i][j], coords[i][j+1] );
01364         arc_len[i] += d;
01365         params[i][1+j] = arc_len[i];
01366       }
01367       for(unsigned int j=0; j<coords[i].size(); j++) {
01368         params[i][j] /= arc_len[i];
01369         //std::cout << "params[" << i << "][" << j << "]=" << params[i][j] << std::endl;
01370       }
01371     }
01372 
01373     // Merge the two sets of parameters into one ordered set without duplicates.
01374     std::vector<double> mgd_params;
01375     mgd_params.reserve( params[0].size() + params[1].size() );
01376     unsigned int a=0, b=0;
01377     while(a<params[0].size() && b<params[1].size()) {
01378       if       (params[0][a] < params[1][b]) {
01379         mgd_params.push_back( params[0][a] );
01380         a++;
01381       } else if(params[0][a] > params[1][b]) {
01382         mgd_params.push_back( params[1][b] );
01383         b++;
01384       } else {
01385         mgd_params.push_back( params[0][a] );
01386         a++;
01387         b++;
01388       }
01389     }
01390 
01391     for(unsigned int i=0; i<mgd_params.size(); i++) {
01392       //std::cout << "mgd_params[" << i << "]=" << mgd_params[i] << std::endl;
01393     }
01394 
01395     // Insert new points to match the other arcs points, by parameter.
01396     for(unsigned int i=0; i<2; i++) {
01397       for(unsigned int j=0; j<mgd_params.size(); j++) {
01398         //std::cout << "params[" << i << "][" << j << "]=" << params[i][j] 
01399         //          << " mgd_params[" << j << "]=" << mgd_params[j] << std::endl;
01400         if(params[i][j] > mgd_params[j]) {
01401           double ratio = (mgd_params[j]-params[i][j-1]) / (params[i][j]-params[i][j-1]);
01402           //std::cout << "j=" << j << " ratio=" << ratio << std::endl;
01403           CartVect pt = coords[i][j-1] + ratio*(coords[i][j]-coords[i][j-1]);
01404           coords[i].insert( coords[i].begin()+j, pt);
01405           params[i].insert( params[i].begin()+j, mgd_params[j]); 
01406         }
01407       }
01408     }
01409 
01410 
01411     // Each arc should now have the same number of points
01412     if(coords[0].size() != coords[1].size()) {
01413       for(unsigned int i=0; i<2; i++) {
01414         for(unsigned int j=0; j<coords[i].size(); j++) {
01415           std::cout << "coords[" << i << "][" << j << "]=" << coords[i][j] << std::endl;
01416         }
01417       }
01418     }
01419     assert( coords[0].size() == coords[1].size() );
01420 
01421     // Get the area between arcs. Use absolute value to prevent cancelling.
01422     double area = 0;
01423     for(unsigned int i=0; i<coords[0].size()-1; i++) {
01424       area += fabs(triangle_area( coords[0][i], coords[0][i+1], coords[1][i+1] ));
01425       //std::cout << "area0=" << area << std::endl;
01426       area += fabs(triangle_area( coords[0][i], coords[1][i+1], coords[1][i] ));
01427       //std::cout << "area1=" << area << std::endl;
01428     }
01429 
01430     // Divide the area by the average length to get the average distance between arcs.
01431     dist = fabs(2*area / (arc_len[0] + arc_len[1] ));
01432     //std::cout << "dist_between_arcs=" << dist << std::endl;
01433     return MB_SUCCESS;
01434   }
01435 
01436 
01437     // qsort struct comparision function        
01438     // If the first handle is the same, compare the second                                        
01439     int compare_edge(const void *a, const void *b) {                                 
01440       struct edge *ia = (struct edge *)a;           
01441       struct edge *ib = (struct edge *)b;
01442       if(ia->v0 == ib->v0) {        
01443         return (int)(ia->v1 - ib->v1);                   
01444       } else {
01445         return (int)(ia->v0 - ib->v0);                   
01446       }
01447       /* float comparison: returns negative if b > a           
01448       and positive if a > b. We multiplied result by 100.0        
01449       to preserve decimal fraction */                
01450     }  
01451 
01452   // WARNING: This skinner goes 10x faster by assuming that no edges already exist
01453   // in the MOAB instance. Otherwise checking to see if an edge exists before
01454   // creating a new one if very slow. This is partly the reason that Skinner is
01455   // very slow.
01456   ErrorCode find_skin( Range tris, const int dim,
01457  //                         std::vector<std::vector<EntityHandle> > &skin_edges,
01458                          Range &skin_edges,
01459                          const bool temp_bool ) {    
01460 
01461     const bool local_debug = false;
01462     //Skinner tool(MBI());
01463     //Range skin_verts;
01464     //return tool.find_skin( tris, dim, skin_edges, temp_bool );
01465     //return tool.find_skin_vertices( tris, skin_verts, &skin_edges, true );
01466   
01467     if(1 != dim) return MB_NOT_IMPLEMENTED;
01468     if(MBTRI != MBI()->type_from_handle(tris.front())) return MB_NOT_IMPLEMENTED;
01469 
01470     skin_edges.clear();
01471     if(tris.empty()) return MB_ENTITY_NOT_FOUND;
01472 
01473     // This implementation gets some of its speed due to not checking for edges
01474     ErrorCode result;
01475     int n_edges;
01476     result = MBI()->get_number_entities_by_type( 0, MBEDGE, n_edges );
01477     assert(MB_SUCCESS == result);
01478     if(0 != n_edges) {
01479       Range temp_edges;
01480       result = MBI()->get_entities_by_type( 0, MBEDGE, temp_edges);
01481       assert(MB_SUCCESS == result);
01482       result = MBI()->list_entities( temp_edges );
01483       assert(MB_SUCCESS == result);
01484     }      
01485     assert(0 == n_edges);
01486 
01487     // Get connectivity. Do not create edges.
01488     edge *edges = new edge[3*tris.size()];
01489     int n_verts;
01490     int ii = 0;
01491     for(Range::iterator i=tris.begin(); i!=tris.end(); i++) {
01492       const EntityHandle *conn;
01493       result = MBI()->get_connectivity( *i, conn, n_verts );
01494       assert(MB_SUCCESS == result);
01495       assert(3 == n_verts);
01496       // shouldn't be degenerate
01497       assert(conn[0] != conn[1]);
01498       assert(conn[1] != conn[2]);
01499       assert(conn[2] != conn[0]);
01500       // make edges
01501       edges[3*ii+0].v0 = conn[0];
01502       edges[3*ii+0].v1 = conn[1];
01503       edges[3*ii+1].v0 = conn[1];
01504       edges[3*ii+1].v1 = conn[2];
01505       edges[3*ii+2].v0 = conn[2];
01506       edges[3*ii+2].v1 = conn[0];
01507       ii++;
01508     }
01509 
01510     // Change the first handle to be lowest
01511     for(unsigned int i=0; i<3*tris.size(); ++i) {
01512       if(edges[i].v0 > edges[i].v1) {
01513         EntityHandle temp = edges[i].v0;
01514         edges[i].v0 = edges[i].v1;
01515         edges[i].v1 = temp;
01516       }
01517     }
01518 
01519     // Sort by first handle, then second handle. Do not sort the extra edge on the
01520     // back.
01521     qsort(edges, 3*tris.size(), sizeof(struct edge), compare_edge);    
01522 
01523     // Go through array, saving edges that are not paired.
01524     for(unsigned int i=0; i<3*tris.size(); i++) {
01525       // If the last edge has not been paired, create it. This avoids overrunning
01526       // the edges array with i+1.
01527       if(3*tris.size()-1 == i) {
01528         const EntityHandle conn[2] = {edges[i].v0, edges[i].v1};
01529         EntityHandle edge;
01530         result = MBI()->create_element( MBEDGE, conn, 2, edge );
01531         assert(MB_SUCCESS == result);
01532         skin_edges.insert(edge);
01533     
01534       // If a match exists, skip ahead
01535       } else if(edges[i].v0==edges[i+1].v0 && edges[i].v1==edges[i+1].v1) {
01536         i++;
01537         // test to make sure surface is manifold
01538         if(3*tris.size() != i+1) { // avoid overrunning array
01539           while( edges[i].v0==edges[i+1].v0 && edges[i].v1==edges[i+1].v1 ) {
01540             if(local_debug) {
01541               std::cout << "find_skin WARNING: non-manifold edge" << std::endl;
01542               MBI()->list_entity( edges[i].v0 );
01543               MBI()->list_entity( edges[i].v1 );
01544             }
01545             ++i;
01546           }
01547         }
01548         // otherwise a skin edge has been found
01549       } else {
01550         const EntityHandle conn[2] = {edges[i].v0, edges[i].v1};
01551         EntityHandle edge;
01552         result = MBI()->create_element( MBEDGE, conn, 2, edge );
01553         if(gen::error(MB_SUCCESS!=result, "could not create edge element")) return result;
01554         skin_edges.insert( edge );
01555       } 
01556     }
01557     delete[] edges;
01558     return MB_SUCCESS;
01559   }
01560   /*  ErrorCode find_skin( Range tris, const int dim, Range &skin_edges, const bool temp ) {
01561     std::vector<std::vector<EntityHandle> > skin_edges_vctr;
01562     ErrorCode result = find_skin( tris, dim, skin_edges_vctr, temp );
01563     assert(MB_SUCCESS == result);
01564     for(std::vector<std::vector<EntityHandle> >::const_iterator i=skin_edges_vctr.begin();
01565         i!=skin_edges_vctr.end(); i++) {
01566       EntityHandle edge;
01567       result = MBI()->create_element( MBEDGE, &(*i)[0], 2, edge );
01568       if(MB_SUCCESS != result) std::cout << "result=" << result << std::endl;
01569       assert(MB_SUCCESS == result);
01570       skin_edges.insert( edge );
01571     }
01572     return MB_SUCCESS;
01573     }*/
01574 
01575 // calculate volume of polyhedron
01576 // Copied from DagMC, without index_by_handle. The dagmc function will
01577 // segfault if build_indices is not first called. For sealing there is
01578 // no need to build_indices.
01579   ErrorCode measure_volume( const EntityHandle volume, double& result, bool debug, bool verbose )
01580 {
01581   ErrorCode rval;
01582   std::vector<EntityHandle> surfaces, surf_volumes;
01583   result = 0.0;
01584   
01585  
01586 
01587   // get surfaces from volume
01588   rval = MBI()->get_child_meshsets( volume, surfaces );
01589   if (MB_SUCCESS != rval)
01590     {
01591       return rval;
01592     }
01593   if(debug) std::cout << "in measure_volume 1" << std::endl;
01594 
01595     // get surface senses
01596   std::vector<int> senses( surfaces.size() );
01597 
01598 
01599   if (rval != MB_SUCCESS)
01600     {
01601       std::cout << "Could not measure volume" << std::endl;
01602       std::cout << "This can happen for 2 reasons, there are no volumes" << std::endl;
01603       std::cout << "or the pointer to the Moab instance is poor" << std::endl;
01604       exit(rval);
01605     }
01606   if(debug) std::cout << surfaces.size() << " " << result << std::endl;
01607 
01608   rval = surface_sense( volume, surfaces.size(), &surfaces[0], &senses[0] );
01609   
01610 
01611   if(debug) std::cout << "in measure_volume 2" << std::endl;
01612 
01613   if (MB_SUCCESS != rval)
01614     {
01615       std::cerr << "ERROR: Surface-Volume relative sense not available. "
01616                 << "Cannot calculate volume." << std::endl;
01617       return rval;
01618     }
01619 
01620 
01621   for (unsigned i = 0; i < surfaces.size(); ++i) 
01622     {
01623       // skip non-manifold surfaces
01624       if (!senses[i])
01625         continue;
01626     
01627       // get triangles in surface
01628       Range triangles;
01629       rval = MBI()->get_entities_by_dimension( surfaces[i], 2, triangles );
01630       if (MB_SUCCESS != rval)
01631         {
01632           return rval;
01633         }
01634     if (!triangles.all_of_type(MBTRI)) 
01635       {
01636         std::cout << "WARNING: Surface " << geom_id_by_handle(surfaces[i])
01637                   << " contains non-triangle elements. Volume calculation may be incorrect."
01638                   << std::endl;
01639         triangles.clear();
01640         rval = MBI()->get_entities_by_type( surfaces[i], MBTRI, triangles );
01641         if (MB_SUCCESS != rval)
01642           {
01643             return rval;
01644           }
01645     }
01646     
01647       // calculate signed volume beneath surface (x 6.0)
01648     double surf_sum = 0.0;
01649     const EntityHandle *conn;
01650     int len;
01651     CartVect coords[3];
01652     for (Range::iterator j = triangles.begin(); j != triangles.end(); ++j) {
01653       rval = MBI()->get_connectivity( *j, conn, len, true );
01654       if (MB_SUCCESS != rval) return rval;
01655       assert(3 == len);
01656       rval = MBI()->get_coords( conn, 3, coords[0].array() );
01657       if (MB_SUCCESS != rval) return rval;
01658     
01659       coords[1] -= coords[0];
01660       coords[2] -= coords[0];
01661       surf_sum += (coords[0] % (coords[1] * coords[2]));
01662     }
01663     result += senses[i] * surf_sum;
01664   }
01665   
01666   result /= 6.0;
01667   if(debug)  std::cout << result << std::endl;
01668   return MB_SUCCESS;
01669 }
01670 
01674   ErrorCode get_signed_volume( const EntityHandle surf_set, double &signed_volume) {
01675     ErrorCode rval;
01676     Range tris;
01677     rval = MBI()->get_entities_by_type( surf_set, MBTRI, tris );
01678     if(MB_SUCCESS != rval) return rval;
01679     signed_volume = 0.0;
01680     const EntityHandle *conn;
01681     int len;
01682     CartVect coords[3];
01683     for (Range::iterator j = tris.begin(); j != tris.end(); ++j) {
01684       rval = MBI()->get_connectivity( *j, conn, len, true );           
01685       if (MB_SUCCESS != rval) return rval;
01686       assert(3 == len);                    
01687       rval = MBI()->get_coords( conn, 3, coords[0].array() );              
01688       if (MB_SUCCESS != rval) return rval;
01689                             
01690       coords[1] -= coords[0];                              
01691       coords[2] -= coords[0];                     
01692       signed_volume += (coords[0] % (coords[1] * coords[2]));       
01693     }                                     
01694     return MB_SUCCESS;
01695   }                 
01696 
01697 ErrorCode measure( const EntityHandle set, const Tag geom_tag, double &size, bool debug,  bool verbose ) {
01698     ErrorCode result;
01699     int dim;
01700     result = MBI()->tag_get_data( geom_tag, &set, 1, &dim );                  
01701     assert(MB_SUCCESS == result);
01702     if(0 == dim) {
01703       std::cout << "measure: cannot measure vertex" << std::endl;
01704       return MB_FAILURE;
01705 
01706     } else if(1 == dim) {
01707       std::vector<EntityHandle> vctr;
01708       result = arc::get_meshset( set, vctr );
01709       assert(MB_SUCCESS == result);
01710       size = length( vctr );
01711 
01712     } else if(2 == dim) {
01713       Range tris;
01714       result = MBI()->get_entities_by_type( set, MBTRI, tris );
01715       assert(MB_SUCCESS == result);
01716       size = triangle_area( tris );
01717 
01718     } else if(3 == dim) {
01719       
01720       result = measure_volume( set, size, debug, verbose );
01721       
01722       if(MB_SUCCESS != result) {
01723         std::cout << "result=" << result << " vol_id=" 
01724                   << gen::geom_id_by_handle(set) << std::endl;
01725       }
01726     } else {
01727       std::cout << "measure: incorrect dimension" << std::endl;
01728       return MB_FAILURE;
01729     }
01730     return MB_SUCCESS;
01731   }
01732 
01733   // From CGMA/builds/dbg/include/CubitDefines
01735   ErrorCode get_curve_surf_sense( const EntityHandle surf_set, const EntityHandle curve_set,
01736                                     int &sense, bool debug ) {
01737     std::vector<EntityHandle> surfs;
01738     std::vector<int> senses;
01739     ErrorCode rval;
01740     GeomTopoTool gt( MBI(), false);
01741     rval = gt.get_senses( curve_set, surfs, senses );
01742     if(gen::error(MB_SUCCESS!=rval,"failed to get_senses")) return rval;
01743     
01744     unsigned counter = 0;
01745     int edim;
01746     for(unsigned i=0; i<surfs.size(); i++) {
01747       edim=gt.dimension(surfs[i]);
01748       if( edim == -1)
01749       { 
01750         surfs.erase(surfs.begin()+i);
01751         senses.erase(senses.begin()+i);
01752         i=0;
01753       }
01754      if(surf_set==surfs[i]) {
01755         sense = senses[i];
01756         
01757         ++counter;
01758       }
01759     }
01760 
01761    if(debug)
01762    {
01763     for(unsigned int index=0; index < surfs.size() ; index++)
01764     {
01765      std::cout << "surf = " << geom_id_by_handle(surfs[index]) << std::endl;
01766      std::cout << "sense = " << senses[index] << std::endl;
01767     }
01768    }
01769     // special case: parent surface does not exist
01770     if(gen::error(0==counter,"failed to find a surf in sense list")) return MB_FAILURE;
01771 
01772    
01773 
01774      // special case: ambiguous
01775     if(1<counter) sense = SENSE_UNKNOWN;
01776     return MB_SUCCESS;
01777   }
01778 
01779   ErrorCode surface_sense( EntityHandle volume,
01780                            int num_surfaces,
01781                            const EntityHandle* surfaces,
01782                            int* senses_out )
01783   {
01784     std::vector<EntityHandle> surf_volumes( 2*num_surfaces );
01785     Tag senseTag = get_tag( "GEOM_SENSE_2", 2, MB_TAG_SPARSE, MB_TYPE_HANDLE, NULL, false );
01786     ErrorCode rval = MBI()->tag_get_data( senseTag , surfaces, num_surfaces, &surf_volumes[0] );
01787     if (MB_SUCCESS != rval)
01788       {
01789         return rval;
01790       }
01791   
01792     const EntityHandle* end = surfaces + num_surfaces;
01793     std::vector<EntityHandle>::const_iterator surf_vols = surf_volumes.begin();
01794     while (surfaces != end) 
01795       {
01796         EntityHandle forward = *surf_vols; ++surf_vols;
01797         EntityHandle reverse = *surf_vols; ++surf_vols;
01798         if (volume == forward) 
01799           {
01800             *senses_out = (volume != reverse); // zero if both, otherwise 1
01801           }
01802         else if (volume == reverse)
01803           {
01804             *senses_out = SENSE_UNKNOWN;
01805           }
01806         else 
01807           {
01808             return MB_ENTITY_NOT_FOUND;
01809           }
01810     
01811         ++surfaces;
01812         ++senses_out;
01813       }
01814   
01815     return MB_SUCCESS;
01816   }
01817 
01819   ErrorCode surface_sense( EntityHandle volume,
01820                              EntityHandle surface,
01821                              int& sense_out )
01822   {
01823     // get sense of surfaces wrt volumes
01824     EntityHandle surf_volumes[2];
01825     Tag senseTag = get_tag( "GEOM_SENSE_2", 2, MB_TAG_SPARSE, MB_TYPE_HANDLE, NULL, false );
01826     ErrorCode rval = MBI()->tag_get_data( senseTag , &surface, 1, surf_volumes );
01827     if (MB_SUCCESS != rval)
01828       {
01829         return rval;
01830       }
01831   
01832     if (surf_volumes[0] == volume)
01833       {
01834         sense_out = (surf_volumes[1] != volume); // zero if both, otherwise 1
01835       }
01836     else if (surf_volumes[1] == volume)
01837       {
01838         sense_out = SENSE_UNKNOWN;
01839       }
01840     else
01841       {
01842         return MB_ENTITY_NOT_FOUND;
01843       }
01844   
01845   return MB_SUCCESS;
01846  }
01847 
01848  Tag get_tag( const char* name, int size, TagType store,
01849                      DataType type, const void* def_value,
01850                      bool create_if_missing)
01851  {
01852    Tag retval = 0;
01853    unsigned flags = store|MB_TAG_CREAT;
01854    if (!create_if_missing)
01855      {
01856        flags |= MB_TAG_EXCL;
01857      }
01858    ErrorCode result = MBI()->tag_get_handle(name, size, type, retval, flags, def_value);
01859    if (create_if_missing && MB_SUCCESS != result)
01860      {
01861        std::cerr << "Couldn't find nor create tag named " << name << std::endl;
01862      }
01863 
01864    return retval;
01865  }
01866 
01867 
01868 
01869 ErrorCode delete_surface( EntityHandle surf, Tag geom_tag, Range tris, int id, bool debug, bool verbose ) {
01870 
01871   ErrorCode result;
01872 
01873   //measure area of the surface  
01874         double area;
01875         result = gen::measure( surf, geom_tag, area, debug, verbose );
01876         if(gen::error(MB_SUCCESS!=result,"could not measure area")) return result;
01877         assert(MB_SUCCESS == result);
01878 
01879   //remove triagngles from the surface
01880         result = MBI()->remove_entities( surf, tris);                          
01881         if(gen::error(MB_SUCCESS!=result,"could not remove tris")) return result;
01882         assert(MB_SUCCESS == result);
01883   //print information about the deleted surface if requested by the user
01884         if(debug) std::cout << "  deleted surface " << id
01885                   << ", area=" << area << " cm^2, n_facets=" << tris.size() << std::endl;
01886 
01887   // delete triangles from mesh data
01888         result = MBI()->delete_entities( tris );                               
01889         if(gen::error(MB_SUCCESS!=result,"could not delete tris")) return result;
01890         assert(MB_SUCCESS == result);
01891 
01892   //remove the sense data for the surface from its child curves
01893         result = remove_surf_sense_data( surf, debug);
01894         if(gen::error(MB_SUCCESS!=result, "could not remove surface's sense data")) return result;
01895 
01896   //remove the surface set itself
01897         result = MBI()->delete_entities( &(surf), 1);
01898         if(gen::error(MB_SUCCESS!=result,"could not delete surface set")) return result;
01899         assert(MB_SUCCESS == result);
01900 
01901 
01902   return MB_SUCCESS;
01903  }
01904 
01906 
01907 ErrorCode remove_surf_sense_data(EntityHandle del_surf, bool debug) {
01908  
01909   ErrorCode result;
01910   GeomTopoTool gt(MBI(), false);
01911     int edim = gt.dimension(del_surf);
01912 
01913     if(gen::error(edim!=2,"could not remove sense data: entity is of the wrong dimension")) return MB_FAILURE;
01914 
01915  // get the curves of the surface
01916         Range del_surf_curves;
01917         result = MBI() -> get_child_meshsets( del_surf, del_surf_curves);
01918         if(gen::error(MB_SUCCESS!=result,"could not get the curves of the surface to delete")) return result;
01919         if (debug) std::cout << "got the curves" << std::endl;
01920        
01921   
01922 
01923         if (debug)
01924         {
01925           std::cout << "number of curves to the deleted surface = " << del_surf_curves.size() << std::endl;
01926           for(unsigned int index =0 ; index < del_surf_curves.size() ; index++)
01927           {
01928             std::cout << "deleted surface's child curve id " << index << " = " << gen::geom_id_by_handle(del_surf_curves[index]) << std::endl;
01929           }        
01930         }
01931         //get the sense data for each curve
01932 
01933         //get sense_tag handles from MOAB
01934         Tag senseEnts, senseSenses;
01935         unsigned flags = MB_TAG_SPARSE;
01936      
01937         //get tag for the entities with sense data associated with a given moab entity
01938         result = MBI()-> tag_get_handle(GEOM_SENSE_N_ENTS_TAG_NAME, 0, MB_TYPE_HANDLE, senseEnts, flags);
01939         if(gen::error(MB_SUCCESS!=result, "could not get senseEnts tag")) return result;
01940         
01941         //get tag for the sense data associated with the senseEnts entities for a given moab entity
01942         result = MBI()-> tag_get_handle(GEOM_SENSE_N_SENSES_TAG_NAME, 0, MB_TYPE_INTEGER, senseSenses, flags);
01943         if(gen::error(MB_SUCCESS!=result,"could not get senseSenses tag")) return result;
01944         
01945         //initialize vectors for entities and sense data
01946         std::vector<EntityHandle> surfaces;
01947         std::vector<int> senses;
01948         for(Range::iterator i=del_surf_curves.begin(); i!=del_surf_curves.end(); i++ )
01949         {
01950        
01951         result = gt.get_senses(*i, surfaces, senses);
01952         if(gen::error(MB_SUCCESS!=result, "could not get the senses for the del_surf_curve")) return result;
01953           // if the surface to be deleted (del_surf) exists in the sense data (which it should), then remove it
01954           for(unsigned int index = 0; index < senses.size() ; index++)
01955           {
01956             if(surfaces[index]==del_surf)
01957             {
01958              surfaces.erase(surfaces.begin() + index);
01959              senses.erase(senses.begin() +index);
01960              index=-1;
01961             }
01962           }
01963           //remove existing sense entity data for the curve
01964           result= MBI()-> tag_delete_data( senseEnts, &*i, 1);
01965           if(gen::error(MB_SUCCESS!=result, "could not delete sense entity data")) return result;
01966        
01967           //remove existing sense data for the curve
01968           result = MBI()-> tag_delete_data(senseSenses, &*i, 1);
01969           if(gen::error(MB_SUCCESS!=result, "could not delete sense data")) return result;
01970 
01971           //reset the sense data for each curve 
01972           result = gt.set_senses( *i, surfaces, senses);
01973           if(gen::error(MB_SUCCESS!=result, "could not update sense data for surface deletion")) return result;
01974 
01975         }
01976 
01977  return MB_SUCCESS;
01978  } 
01979 
01981  ErrorCode combine_merged_curve_senses( std::vector<EntityHandle> &curves, Tag merge_tag, bool debug) {
01982 
01983   ErrorCode result;
01984   
01985   for(std::vector<EntityHandle>::iterator j=curves.begin(); j!=curves.end(); j++) {
01986 
01987 
01988   EntityHandle merged_curve;
01989   result = MBI() -> tag_get_data( merge_tag, &(*j), 1, &merged_curve);
01990   if(gen::error(MB_SUCCESS!=result && MB_TAG_NOT_FOUND!=result, "could not get the merge_tag data of the curve")) return result;
01991 
01992   if(MB_SUCCESS==result) { // we have found a merged curve pairing
01993     // add the senses from the curve_to_delete to curve_to keep
01994     // create vectors for the senses and surfaces of each curve
01995     std::vector<EntityHandle> curve_to_keep_surfs, curve_to_delete_surfs, combined_surfs;
01996     std::vector<int> curve_to_keep_senses, curve_to_delete_senses, combined_senses;
01997 
01998     //initialize GeomTopoTool.cpp instance in MOAB
01999     GeomTopoTool gt(MBI(), false);
02000     // get senses of the iterator curve and place them in the curve_to_delete vectors
02001     result = gt.get_senses( *j, curve_to_delete_surfs, curve_to_delete_senses);
02002     if(gen::error(MB_SUCCESS!=result, "could not get the surfs/senses of the curve to delete")) return result;
02003     // get surfaces/senses of the merged_curve and place them in the curve_to_keep vectors
02004     result = gt.get_senses( merged_curve, curve_to_keep_surfs, curve_to_keep_senses);
02005     if(gen::error(MB_SUCCESS!=result, "could not get the surfs/senses of the curve to delete")) return result;
02006     
02007     if(debug){
02008           std::cout << "curve to keep id = " << gen::geom_id_by_handle(merged_curve) << std::endl;
02009           std::cout << "curve to delete id = " << gen::geom_id_by_handle(*j) << std::endl;
02010           for(unsigned int index=0; index < curve_to_keep_surfs.size(); index++)
02011           {
02012            std::cout << "curve_to_keep_surf " << index << " id = " << gen::geom_id_by_handle(curve_to_keep_surfs[index]) << std::endl;
02013            std::cout << "curve_to_keep_sense " << index << " = " << curve_to_keep_senses[index] << std::endl;
02014            
02015           }
02016           for(unsigned int index=0; index < curve_to_keep_surfs.size(); index++)
02017           {
02018           
02019           std::cout << "curve_to_delete_surf " << index << " id = " << gen::geom_id_by_handle(curve_to_delete_surfs[index]) << std::endl;
02020           std::cout << "curve_to_delete_sense " << index << " = " << curve_to_delete_senses[index] << std::endl;
02021           }
02022     } // end of debug st.
02023 
02024     // combine the surface and sense data for both curves into the same vector 
02025 
02026     combined_surfs.insert( combined_surfs.end(), curve_to_keep_surfs.begin(), 
02027                                                        curve_to_keep_surfs.end() );
02028     combined_surfs.insert( combined_surfs.end(), curve_to_delete_surfs.begin(), 
02029                                                        curve_to_delete_surfs.end() );
02030     combined_senses.insert(combined_senses.end(),curve_to_keep_senses.begin(),
02031                                                        curve_to_keep_senses.end() );
02032     combined_senses.insert(combined_senses.end(),curve_to_delete_senses.begin(),
02033                                                        curve_to_delete_senses.end() );
02034 
02035     if(debug){
02036       std::cout << combined_surfs.size() << std::endl;
02037       std::cout << combined_senses.size() << std::endl;
02038       for(unsigned int index=0; index < combined_senses.size(); index++)
02039           {
02040 
02041            std::cout << "combined_surfs{"<< index << "] = " << gen::geom_id_by_handle(combined_surfs[index]) << std::endl;
02042            std::cout << "combined_sense["<< index << "] = " << combined_senses[index] << std::endl;
02043           }
02044     } // end debug st. 
02045 
02046     result = gt.set_senses(merged_curve, combined_surfs, combined_senses);
02047     if(gen::error(MB_SUCCESS!=result && MB_MULTIPLE_ENTITIES_FOUND!=result,"failed to set senses: "));
02048 
02049 
02050     
02051     // add the duplicate curve_to_keep to the vector of curves
02052     *j = merged_curve;
02053 
02054    } //end merge_tag result if st.
02055 
02056 
02057   } //end curves loop
02058 
02059 
02060 
02061 
02062   return MB_SUCCESS;
02063   }
02064 
02065 ErrorCode get_sealing_mesh_tags( double &facet_tol,
02066                            double &sme_resabs_tol,
02067                            Tag &geom_tag,
02068                            Tag &id_tag,
02069                            Tag &normal_tag,
02070                            Tag &merge_tag,
02071                            Tag &faceting_tol_tag,
02072                            Tag &geometry_resabs_tag,
02073                            Tag &size_tag,
02074                            Tag &orig_curve_tag) {
02075 
02076   ErrorCode result;
02077 
02078     result = MBI()->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1,
02079                                 MB_TYPE_INTEGER, geom_tag, MB_TAG_DENSE|MB_TAG_CREAT );
02080     assert( MB_SUCCESS == result );
02081     if ( result != MB_SUCCESS )
02082       {
02083         moab_printer(result);
02084       }
02085     result = MBI()->tag_get_handle( GLOBAL_ID_TAG_NAME, 1,
02086                                 MB_TYPE_INTEGER, id_tag, MB_TAG_DENSE|MB_TAG_CREAT);
02087     assert( MB_SUCCESS == result );
02088     if ( result != MB_SUCCESS )
02089       {
02090         moab_printer(result);
02091       }
02092     result = MBI()->tag_get_handle( "NORMAL", sizeof(CartVect), MB_TYPE_OPAQUE,
02093         normal_tag, MB_TAG_DENSE|MB_TAG_CREAT);
02094     assert( MB_SUCCESS == result );
02095     if ( result != MB_SUCCESS )
02096       {
02097         moab_printer(result);
02098       }
02099     result = MBI()->tag_get_handle( "MERGE", 1, MB_TYPE_HANDLE,
02100         merge_tag, MB_TAG_SPARSE|MB_TAG_CREAT );
02101     assert( MB_SUCCESS == result );
02102     if ( result != MB_SUCCESS )
02103       {
02104         moab_printer(result);
02105       } 
02106     result = MBI()->tag_get_handle( "FACETING_TOL", 1, MB_TYPE_DOUBLE,
02107         faceting_tol_tag , MB_TAG_SPARSE|MB_TAG_CREAT );
02108     assert( MB_SUCCESS == result );
02109     if ( result != MB_SUCCESS )
02110       {
02111         moab_printer(result);
02112       }
02113     result = MBI()->tag_get_handle( "GEOMETRY_RESABS", 1,     MB_TYPE_DOUBLE,
02114                              geometry_resabs_tag, MB_TAG_SPARSE|MB_TAG_CREAT  );
02115     assert( MB_SUCCESS == result );
02116     if ( result != MB_SUCCESS )
02117       {
02118         moab_printer(result);
02119       }
02120     result = MBI()->tag_get_handle( "GEOM_SIZE", 1, MB_TYPE_DOUBLE,
02121                                     size_tag, MB_TAG_DENSE|MB_TAG_CREAT  );
02122     assert( (MB_SUCCESS == result) );
02123     if ( result != MB_SUCCESS )
02124       {
02125         moab_printer(result);
02126       }
02127     int true_int = 1;    
02128     result = MBI()->tag_get_handle( "ORIG_CURVE", 1,
02129                                 MB_TYPE_INTEGER, orig_curve_tag, MB_TAG_DENSE|MB_TAG_CREAT, &true_int );
02130     assert( MB_SUCCESS == result );
02131     if ( result != MB_SUCCESS )
02132       {
02133         moab_printer(result);
02134       }
02135     // PROBLEM: MOAB is not consistent with file_set behavior. The tag may not be
02136     // on the file_set.
02137     Range file_set;
02138     result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &faceting_tol_tag,
02139                                                   NULL, 1, file_set );
02140 
02141     if(gen::error(MB_SUCCESS!=result,"could not get faceting_tol_tag"))
02142       {
02143         return result;
02144       }
02145 
02146     gen::error(file_set.empty(),"file set not found");
02147 
02148     if(gen::error(1!=file_set.size(),"Refacet with newer version of ReadCGM.")) 
02149       {
02150         return MB_FAILURE;
02151       }
02152 
02153     result = MBI()->tag_get_data( faceting_tol_tag, &file_set.front(), 1,  
02154                                   &facet_tol );
02155     assert(MB_SUCCESS == result);
02156     result = MBI()->tag_get_data( geometry_resabs_tag, &file_set.front(), 1,  
02157                                   &sme_resabs_tol );
02158     if(MB_SUCCESS != result)
02159       {
02160         std::cout <<  "absolute tolerance could not be read from file" << std::endl;
02161       }
02162 
02163 
02164 
02165   return MB_SUCCESS;
02166 
02167   }
02168 
02169   ErrorCode get_geometry_meshsets( Range geometry_sets[], Tag geom_tag, bool verbose) {
02170 
02171     ErrorCode result;
02172 
02173     // get all geometry sets
02174     for(unsigned dim=0; dim<4; dim++) 
02175       {
02176         void *val[] = {&dim};
02177         result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag,
02178                                                     val, 1, geometry_sets[dim] );
02179         assert(MB_SUCCESS == result);
02180 
02181         // make sure that sets TRACK membership and curves are ordered
02182         // MESHSET_TRACK_OWNER=0x1, MESHSET_SET=0x2, MESHSET_ORDERED=0x4      
02183 
02184         for(Range::iterator i=geometry_sets[dim].begin(); i!=geometry_sets[dim].end(); i++)
02185           {
02186             unsigned int options;
02187             result = MBI()->get_meshset_options(*i, options );
02188             assert(MB_SUCCESS == result);
02189     
02190             // if options are wrong change them
02191             if(dim==1) 
02192               {
02193                 if( !(MESHSET_TRACK_OWNER&options) || !(MESHSET_ORDERED&options) ) 
02194                   {
02195                     result = MBI()->set_meshset_options(*i, MESHSET_TRACK_OWNER|MESHSET_ORDERED);
02196                     assert(MB_SUCCESS == result);
02197                   }
02198               } 
02199             else 
02200               {
02201                 if( !(MESHSET_TRACK_OWNER&options) ) 
02202                   {        
02203                     result = MBI()->set_meshset_options(*i, MESHSET_TRACK_OWNER);
02204                     assert(MB_SUCCESS == result);
02205                   }
02206               }
02207           }
02208       }
02209  
02210     return result;
02211 
02212     }
02213 
02214   ErrorCode check_for_geometry_sets(Tag geom_tag, bool verbose){
02215 
02216     ErrorCode result;
02217         // go get all geometry sets
02218         Range geometry_sets[4];
02219         result = get_geometry_meshsets( geometry_sets, geom_tag, false);
02220         if(gen::error(MB_SUCCESS!=result,"could not get the geometry meshsets")) return result;
02221 
02222         //make sure they're there
02223         for(unsigned dim=0; dim<4; dim++){
02224  
02225           if(geometry_sets[dim].size() == 0) return MB_FAILURE;
02226 
02227         }
02228  
02229     return MB_SUCCESS;
02230   }
02231 
02232 
02233 } //EOL
02234 
02235 Interface *MBI()
02236 {
02237     static Core instance;
02238     return &instance;
02239 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines