MeshKit  1.0
post_process.cpp
Go to the documentation of this file.
00001 // ********************************************************************
00002 // Brandon Smith
00003 // August, 2009
00004 
00005 /* _curve_to_be_tested_for_watertightness_
00006       vert1 X X vert1
00007             | |
00008       vert2 X |
00009   surf1     | |    surf2
00010             | |
00011       vert3 X X vert2
00012             | |
00013       vert4 X X vert3                   */
00014 
00015 // input:  h5m filename, tolerance
00016 // output: watertight h5m
00017 
00018 // make CXXFLAGS=-g for debug
00019 // make CXXFLAGS=-pg for profiling
00020 
00021 #include <iostream>
00022 #include <sstream>
00023 #include <iomanip> // for setprecision
00024 #include <limits> // for min/max values
00025 #include <cassert>
00026 #include <cmath>
00027 #include <ctime>
00028 #include <vector>
00029 #include "moab/Core.hpp"
00030 #include "MBTagConventions.hpp"
00031 #include "moab/Range.hpp"
00032 #include "moab/Skinner.hpp"
00033 
00034 #include "meshkit/gen.hpp"
00035 #include "meshkit/arc.hpp"
00036 #include "meshkit/zip.hpp"
00037 #include "meshkit/cleanup.hpp"
00038 
00039 using namespace moab;
00040 
00041 ErrorCode delete_all_edges() {
00042   // delete all of the edges. Never keep edges. They are too hard to track and use
00043   // due to orientation and multiple entities errors when merging.
00044   ErrorCode result;
00045   Range edges;
00046   result = MBI()->get_entities_by_type( 0, MBEDGE, edges );
00047   assert(MB_SUCCESS == result);
00048   result = MBI()->delete_entities( edges );
00049   assert(MB_SUCCESS == result);
00050   return MB_SUCCESS;
00051 }
00052 
00053 // Input: unordered sets of curves that do not track ownership
00054 // Output: ordered sets of verts that do track ownership. All edges are deleted.
00055 ErrorCode prepare_curves(Range &curve_sets,
00056                            Tag geom_tag, Tag id_tag, Tag merge_tag,
00057                            double const MERGE_TOL, double const FACET_TOL) {
00058   ErrorCode result;
00059 
00060   // process each curve
00061   for(MBRange::iterator i=curve_sets.begin(); i!=curve_sets.end(); i++ ) {
00062     // get the curve id of the curve meshset
00063     int id;
00064     result = MBI()->tag_get_data( id_tag, &(*i), 1, &id );
00065     assert(MB_SUCCESS == result);
00066     std::cout << "curve " << id << std::endl;
00067 
00068     // get the range of edges of the curve meshset
00069     std::vector<EntityHandle> curve_edges;
00070     result = MBI()->get_entities_by_type( *i, MBEDGE, curve_edges );
00071     assert( MB_SUCCESS == result );
00072 
00073     /* Merge the endpoints of the curve and remove its edges if it is too small.
00074     Use the MERGE_TOL because these edges will be merged with the MERGE_TOL 
00075     during zipping anyhow. Doing this now removes small curves from zipping and
00076     reduces ambiguity. */
00077     if(MERGE_TOL > gen::length(curve_edges)) {
00078       std::cout << "  removed curve with length=" << gen::length(curve_edges) 
00079                 << " n_verts=" << curve_edges.size()+1 << std::endl;
00080 
00081       // get the endpoints of the curve
00082       Range endpt_sets;
00083       result = MBI()->get_child_meshsets( *i, endpt_sets );
00084       assert(MB_SUCCESS==result);
00085       std::cout << "  endpt_sets.size()=" << endpt_sets.size() << std::endl;
00086       if(endpt_sets.empty()) {
00087         assert(false);
00088       } else if(1 == endpt_sets.size()) {
00089         // nothing
00090       } else if(2 == endpt_sets.size()) {
00091         Range front_endpt, back_endpt;
00092         result = MBI()->get_entities_by_type( endpt_sets.front(), MBVERTEX, front_endpt);
00093         assert(MB_SUCCESS == result);
00094         assert(1 == front_endpt.size());
00095         result = MBI()->get_entities_by_type( endpt_sets.back(), MBVERTEX, back_endpt);
00096         assert(MB_SUCCESS == result);
00097         assert(1 == back_endpt.size());
00098         // merge the endpoints-ALWAYS CHECK TO AVOID MERGING THE SAME ENTITY!!!
00099         if(front_endpt[0] != back_endpt[0]) {
00100           result = MBI()->merge_entities( front_endpt[0], back_endpt[0], false, true);
00101           assert(MB_SUCCESS == result);
00102           // check for and remove degenerate edges caused by the merge
00103           Range edges;
00104           EntityHandle temp = front_endpt[0];
00105           result = MBI()->get_adjacencies( &temp, 1, 1, false, edges);
00106           assert(MB_SUCCESS == result);
00107           for(MBRange::iterator j=edges.begin(); j!=edges.end(); j++) {
00108             const EntityHandle *conn;
00109             int n_verts;
00110             result = MBI()->get_connectivity( *j, conn, n_verts); 
00111             assert(MB_SUCCESS == result);
00112             if(conn[0] == conn[1]) {
00113               result = MBI()->delete_entities( &(*j), 1 );
00114               assert(MB_SUCCESS == result);
00115             }
00116           }
00117         }
00118       } else {
00119         assert(false);
00120       }
00121       // It is possible that the endpoints themselves are orphaned. Should these
00122       // be deleted?
00123       
00124       // Remove the curve set. This also removes parent-child relationships.
00125       result = MBI()->delete_entities( &(*i), 1);
00126       assert(MB_SUCCESS == result);
00127       i = curve_sets.erase(i) - 1;
00128       continue;
00129     }
00130 
00131     // convert the curve of edges into a curve of verts 
00132     std::vector<EntityHandle> ordered_verts;
00133     result = gen::ordered_verts_from_ordered_edges( curve_edges, ordered_verts);
00134     assert(MB_SUCCESS == result);
00135 
00136     // the edges are no longer needed
00137     result = MBI()->delete_entities( &curve_edges[0], curve_edges.size() );
00138     assert(MB_SUCCESS == result);
00139 
00140     // replace the unordered edges with the ordered verts
00141     result = arc::set_meshset( *i, ordered_verts );
00142     assert(MB_SUCCESS == result);
00143   }
00144 
00145   return MB_SUCCESS;
00146 }
00147 
00148 /* Isolate the failure by removing the curve and loop that failed. The zip_loop
00149    function will be called again on the remaining loops and curves. */
00150 ErrorCode remove_failed_loop_and_curve( std::vector<std::vector<EntityHandle> > &skin,
00151                                           std::vector<std::vector<EntityHandle> > &curves,
00152                                           std::vector<int> &curve_ids,
00153                                           std::vector<EntityHandle> &curve_sets,
00154                                           //Range &curve_sets,
00155                                           const unsigned int loop,
00156                                           const unsigned int curve ) {
00157   skin.erase( skin.begin()+loop );
00158   curves.erase( curves.begin()+curve );
00159   curve_ids.erase( curve_ids.begin()+curve ); 
00160   curve_sets.erase( curve_sets.begin()+curve );
00161   std::cout << "remove_failed_loop: removed loop " << loop << std::endl;
00162   return MB_SUCCESS;
00163 }
00164 
00165   // input: surface sets, ordered curve sets,
00166   // output: skin arcs corresponding to curves are added to parent surface sets
00167 ErrorCode prepare_surfaces(Range &surface_sets,
00168                              Tag geom_tag, Tag id_tag, Tag normal_tag, Tag merge_tag,
00169                              const double SME_RESABS_TOL, const double FACET_TOL, 
00170                                const double MERGE_TOL) {
00171     
00172     ErrorCode result;
00173 
00174     // loop over each surface meshset
00175     for(MBRange::iterator i=surface_sets.begin(); i!=surface_sets.end(); i++ ) {
00176 
00177       // get the surf id of the surface meshset
00178       int surf_id;
00179       result = MBI()->tag_get_data( id_tag, &(*i), 1, &surf_id );
00180       assert(MB_SUCCESS == result);
00181       std::cout << "  surf id=" << surf_id << std::endl;
00182 
00183       // get facets of the surface meshset
00184       Range tris;
00185       result = MBI()->get_entities_by_type( *i, MBTRI, tris );
00186       assert(MB_SUCCESS == result);
00187 
00188       // Get the curves sets
00189       std::vector<EntityHandle> curve_sets, unmerged_curve_sets;
00190       result = MBI()->get_child_meshsets( *i, curve_sets );
00191       assert(MB_SUCCESS==result);
00192 
00193       // Update the curve_sets with that contain entity_to_delete curves with their
00194       // entity_to_keep curves. Unmerged_curve_sets will end up holding the curves
00195       // of this surface that are not merged with another curve in this surface.
00196       for(std::vector<EntityHandle>::iterator j=curve_sets.begin();
00197           j!=curve_sets.end(); j++) {
00198         EntityHandle merged_curve, curve;
00199         result = MBI()->tag_get_data( merge_tag, &(*j), 1, &merged_curve );
00200         assert(MB_TAG_NOT_FOUND==result || MB_SUCCESS==result);
00201         if(MB_TAG_NOT_FOUND==result) {
00202           curve = *j;
00203         } else if(MB_SUCCESS == result) {
00204           std::cout << "  curve " << gen::geom_id_by_handle(*j) 
00205                     << " is entity_to_delete" << std::endl;
00206           curve = merged_curve;
00207           // should parent-childs be updated for the entity_to_keep?
00208         } else {
00209           std::cout << "prepare_surfaces: result=" << result << std::endl;
00210           return result;        
00211         }
00212       
00213         // If the curve is in unmerged_curve_sets, then remove it. Otherwise add it.
00214         std::vector<EntityHandle>::iterator k=find(unmerged_curve_sets.begin(),
00215           unmerged_curve_sets.end(), curve);
00216         if(unmerged_curve_sets.end() == k) {
00217           //std::cout << "  curve " << gen::geom_id_by_handle(*k) 
00218           //          << " is entity_to_keep" << std::endl;
00219           unmerged_curve_sets.push_back(curve);
00220         } else {
00221           unmerged_curve_sets.erase(k);
00222         }
00223       }
00224 
00225       // If all of the curves are merged, remove the surfaces facets.
00226       if(unmerged_curve_sets.empty()) {
00227         result = MBI()->remove_entities( *i, tris);                                           
00228         assert(MB_SUCCESS == result);
00229         std::cout << "  removed " << tris.size() << " facets and deleted surface" << std::endl;
00230         result = MBI()->delete_entities( tris );                                              
00231         assert(MB_SUCCESS == result);
00232         // remove the surface set itself
00233         result = MBI()->delete_entities( &(*i), 1);
00234         assert(MB_SUCCESS == result);
00235         i = surface_sets.erase(i) - 1;
00236         continue;
00237       }
00238 
00239       // Try zipping without curves that are merged with each other
00240       curve_sets.swap(unmerged_curve_sets);
00241 
00242       // Save the normals of the facets. These will later be used to determine if
00243       // the tri became inverted.
00244       result = gen::save_normals( tris, normal_tag );
00245       assert(MB_SUCCESS == result);
00246   
00247       // get the range of skin edges from the range of facets
00248       Skinner tool(MBI());
00249       Range skin_edges;
00250 
00251       // merge the vertices of the skin
00252       // BRANDON: For some reason cgm2moab does not do this? This was the 
00253       // problem with mod13 surf 881. Two skin verts were coincident. A tol=1e-10
00254       // found the verts, but tol=0 did not.
00255       Range skin_verts;
00256       result = MBI()->get_adjacencies( skin_edges, 0, false, skin_verts, 
00257                                        MBInterface::UNION );
00258       assert(MB_SUCCESS == result);
00259       result = gen::merge_vertices( skin_verts, SME_RESABS_TOL );         
00260       if (MB_SUCCESS != result) {
00261         std::cout << "result= " << result << std::endl;                  
00262         std::cout << "SURFACE_ZIPPING_FAILURE: could not merge vertices, surf_id="   
00263                   << surf_id << std::endl;  
00264         continue;
00265       }  
00266                                                       
00267       // Create loops with the skin edges.  
00268       std::vector< std::vector<EntityHandle> > skin_loops_of_edges;
00269       if(MB_SUCCESS != result) {
00270         std::cout << "SURFACE_ZIPPING_FAILURE: could not create loops for surf_id=" 
00271                   << surf_id << std::endl;
00272         continue;
00273       }
00274       std::cout << "    surf has " << skin_loops_of_edges.size() 
00275                 << " skin loop(s)." << std::endl;
00276     
00277       // Convert the loops of skin edges to loops of skin verts.
00278       std::vector< std::vector<EntityHandle> > skin(skin_loops_of_edges.size());
00279       for(unsigned int j=0; j<skin_loops_of_edges.size(); j++) {
00280         result = gen::ordered_verts_from_ordered_edges( skin_loops_of_edges[j], skin[j] );
00281         assert(MB_SUCCESS == result);
00282         // check to make sure that the loop is closed
00283         assert(skin[j].front() == skin[j].back());
00284       }
00285 
00286       // edges are no longer needed       
00287       result = delete_all_edges();
00288       assert(MB_SUCCESS == result);
00289 
00290       /* Get the curves that are part of the surface. Use vectors to store all curve
00291          stuff so that we can remove curves from the set as they are zipped. */
00292       //curve_sets.clear();
00293 
00294       //result = MBI()->get_child_meshsets( *i, curve_sets );
00295       //assert(MB_SUCCESS==result);
00296       std::vector<int> curve_ids;
00297       int curve_id;
00298       std::vector<std::vector<EntityHandle> > curves;
00299       //for(MBRange::iterator j=curve_sets.begin(); j!=curve_sets.end(); j++) {
00300       //for(unsigned int j=0; j<curve_sets.size(); j++) {
00301       for(std::vector<EntityHandle>::iterator j=curve_sets.begin();
00302         j!=curve_sets.end(); j++) {
00303 
00304         // If a delete_curve, replace it with the keep_curve. This approach allows
00305         // for duplicates because we are using vectors instead of ranges. Note that
00306         // parent-child links also cannot store duplicate handles.
00307         EntityHandle merged_curve;
00308         //EntityHandle temp = curve_sets[j];
00309         result = MBI()->tag_get_data( merge_tag, &(*j), 1, &merged_curve );     
00310         assert(MB_TAG_NOT_FOUND==result || MB_SUCCESS==result);
00311         if(MB_SUCCESS == result) *j = merged_curve;
00312 
00313         // do not add a curve if it contains nothing
00314         //temp = curve_sets[j];
00315         result = MBI()->tag_get_data( id_tag, &(*j), 1, &curve_id );
00316         assert(MB_SUCCESS == result);
00317         std::cout << "  curve_id=" << curve_id << " handle=" << *j << std::endl;
00318         curve_ids.push_back(curve_id);
00319         std::vector<EntityHandle> curve;
00320         result = arc::get_meshset( *j, curve );
00321         assert(MB_SUCCESS == result);
00322         curves.push_back( curve );
00323       }
00324 
00325       // Keep zipping loops until each is either zipped or failed. This function
00326       // returns only after all loops are zipped or a failure occurs.
00327       while(!skin.empty()) {
00328         //result = zip_loop( normal_tag, FACET_TOL, MERGE_TOL, 
00329         //                 curves, skin, curve_ids, *i, curve_sets );
00330         if(MB_SUCCESS != result) {
00331           std::cout << "SURFACE_ZIPPING_FAILURE: could not zip surf_id=" << surf_id << std::endl;
00332         }
00333       }
00334 
00335       // mod13surf2996, 3028 and 2997 are adjacent to the same bad geometry (figure 8 loop)
00336       //assert(MB_SUCCESS==result || 2996==surf_id || 2997==surf_id || 3028==surf_id);
00337     }
00338     return MB_SUCCESS;
00339   }
00340 
00341   ErrorCode test_edges() {
00342     ErrorCode result;
00343     Range edges;
00344     result = MBI()->get_entities_by_dimension( 0, 1, edges );
00345     assert(MB_SUCCESS == result);
00346     MBI()->list_entities( edges );
00347     return MB_SUCCESS;
00348   }
00349 
00350 
00351   int main(int argc, char **argv) {
00352 
00353     // ******************************************************************
00354     // Load the h5m file and create tags.
00355     // ******************************************************************
00356 
00357     clock_t start_time = clock(), prep_time, zip_time;
00358     if(2 > argc) {
00359       std::cout << "usage: do not use" << std::endl;
00360       std::cout << "./post_process <input_file>" << std::endl;
00361       return 1;
00362     }
00363     ErrorCode result;
00364     std::string input_name = argv[1];
00365 
00366     // The root name does not have an extension
00367     std::string root_name = argv[1];
00368     int len = root_name.length();
00369     root_name.erase(len - 4);
00370     const double MERGE_TOL = 1e-3; // should this depend on FACET_TOL? 
00371 
00372     // load the input file
00373     EntityHandle input_meshset;
00374     result = MBI()->create_meshset( MESHSET_SET, input_meshset );
00375     assert(MB_SUCCESS == result);
00376     if(std::string::npos != input_name.find("h5m")) {
00377       result = MBI()->load_file( input_name.c_str(), &input_meshset );
00378       assert( MB_SUCCESS == result );
00379     } else {
00380       std::cout << "invalid input file: must be h5m" << std::endl;
00381       return 1;
00382     }
00383 
00384     // create tags
00385     Tag geom_tag, id_tag, sense_tag, normal_tag, merge_tag;
00386     result = MBI()->tag_get_handle( GEOM_DIMENSION_TAG_NAME, sizeof(int), 
00387                                 MB_TYPE_INTEGER, geom_tag, MB_TAG_DENSE, 0, 0 );
00388     assert( MB_SUCCESS == result );
00389     result = MBI()->tag_get_handle( GLOBAL_ID_TAG_NAME, sizeof(int), 
00390                                 MB_TYPE_INTEGER, id_tag,MB_TAG_DENSE, 0, 0 );
00391     assert( MB_SUCCESS == result );
00392     result = MBI()->tag_get_handle( "GEOM_SENSE_2", 2*sizeof(EntityHandle),
00393                                 MB_TYPE_HANDLE, sense_tag, MB_TAG_DENSE, 0, 0 );
00394     assert( MB_SUCCESS == result );
00395     result = MBI()->tag_get_handle( "NORMAL", sizeof(MBCartVect), 
00396                                 MB_TYPE_OPAQUE, normal_tag, MB_TAG_DENSE, 0, 0 );
00397     assert( MB_SUCCESS == result );
00398     result = MBI()->tag_get_handle( "MERGE", sizeof(EntityHandle),
00399                                 MB_TYPE_HANDLE, merge_tag, MB_TAG_SPARSE, 0, 0 );
00400     assert( MB_SUCCESS == result );
00401 
00402     // get all geometry sets
00403     Range geom_sets[4];
00404     for(unsigned dim=0; dim<4; dim++) {
00405       void *val[] = {&dim};
00406       result = MBI()->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag,
00407                                                     val, 1, geom_sets[dim] );
00408       assert(MB_SUCCESS == result);
00409       // make sure that sets TRACK membership and curves are ordered
00410       // MESHSET_TRACK_OWNER=0x1, MESHSET_SET=0x2, MESHSET_ORDERED=0x4
00411       for(MBRange::iterator i=geom_sets[dim].begin(); i!=geom_sets[dim].end(); i++) {
00412         unsigned int options;
00413         result = MBI()->get_meshset_options(*i, options );
00414         assert(MB_SUCCESS == result);
00415     
00416         // if options are wrong change them
00417         if(dim==1) {
00418           if( !(MESHSET_TRACK_OWNER&options) || !(MESHSET_ORDERED&options) ) {
00419             result = MBI()->set_meshset_options(*i, MESHSET_TRACK_OWNER|MESHSET_ORDERED);
00420             assert(MB_SUCCESS == result);
00421           }
00422         } else {
00423           if( !(MESHSET_TRACK_OWNER&options) ) {        
00424             result = MBI()->set_meshset_options(*i, MESHSET_TRACK_OWNER);
00425             assert(MB_SUCCESS == result);
00426           }
00427         }
00428       }
00429     }
00430     std::cout << geom_sets[3].size() << " volumes, " 
00431               << geom_sets[2].size() << " surfaces, and "
00432               << geom_sets[1].size() << " curves" << std::endl;  
00433 
00434     result = cleanup::delete_small_edges(geom_sets[2], MERGE_TOL);
00435     assert(MB_SUCCESS == result);
00436   
00437     std::string output_filename = root_name + "_tri.h5m";
00438     // PROBLEM: If I write the input meshset the writer returns MB_FAILURE.
00439     // This happens only if I delete vertices when merging.
00440     // result = MBI()->write_mesh( filename_new.c_str(), &input_meshset, 1);
00441     result = MBI()->write_mesh( output_filename.c_str() );
00442     if (MB_SUCCESS != result) std::cout << "result= " << result << std::endl;
00443     assert(MB_SUCCESS == result);
00444 
00445     zip_time = clock();
00446     std::cout << "zipping took " << (double) (zip_time-prep_time)/CLOCKS_PER_SEC 
00447               << " sec." << std::endl;
00448   
00449     return 0;  
00450   }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines