MeshKit  1.0
check_watertight.cpp
Go to the documentation of this file.
00001 // ********************************************************************
00002 // Brandon Smith
00003 // August 2009
00004 
00005 // This is a function to test DagMC-style mesh for watertightness. For
00006 // now this will be a stand-alone code that uses MOAB. For volumes to 
00007 // be watertight, the facet edges of each surface must be matched 
00008 // one-to-one. By default this checks for topological watertightness.
00009 // To instead use a geometrical tolerance between two vertices that are
00010 // considered the same, pass in a tolerance.
00011 //
00012 // input:  h5m file name, tolerance(optional)
00013 // output: list of unmatched facet edges and their parent surfaces, by volume.
00014 
00015 // make CXXFLAGS=-g for debug
00016 // make CXXFLAGS=-pg for profiling
00017 
00018 /* Assume no MBEDGEs exist for fastest skinning
00019    For each volume:
00020      For each child surface:
00021        skin surface tris
00022        enter data into coords_and_id
00023      }
00024      match edges
00025    }
00026    Each surface is skinned twice, but the logic is simple and the memory handling 
00027    is easy.
00028 */   
00029 
00030 #include <iostream>
00031 #include <cmath>
00032 #include <cstdlib>
00033 #include <ctime>
00034 #include <vector>
00035 #include <set>
00036 #include <algorithm>
00037 #include "moab/Core.hpp"
00038 #include "MBTagConventions.hpp"
00039 #include "moab/Range.hpp"
00040 #include "moab/Skinner.hpp"
00041 
00042 #include "meshkit/cw_func.hpp"
00043 #include "meshkit/gen.hpp"
00044 #include "meshkit/arc.hpp"
00045 #include "meshkit/zip.hpp"
00046 using namespace moab;
00047 
00048 // struct to hold coordinates of skin edge, it's surface id, and a matched flag
00049 struct coords_and_id {
00050   double x1;
00051   double y1;
00052   double z1;
00053   double x2;
00054   double y2;
00055   double z2;
00056   int  surf_id;
00057   bool matched;
00058   EntityHandle vert1;
00059   EntityHandle vert2;
00060 };
00061 
00062 /* qsort struct comparision function */
00063 int compare_by_handle(const void *a, const void *b)
00064 {
00065   struct coords_and_id *ia = (struct coords_and_id *)a;
00066   struct coords_and_id *ib = (struct coords_and_id *)b;
00067   if(ia->vert1 == ib->vert1) 
00068   {
00069     return (int)(ia->vert2 - ib->vert2);
00070   } 
00071   else 
00072   {
00073     return (int)(ia->vert1 - ib->vert1);
00074   }
00075   /* float comparison: returns negative if b > a 
00076      and positive if a > b. We multiplied result by 100.0
00077      to preserve decimal fraction */
00078 } 
00079 
00080 /* qsort struct comparision function */
00081 // This is tricky because doubles always get rounded down to ints.
00082 int compare_by_coords(const void *a, const void *b)
00083 {
00084   struct coords_and_id *ia = (struct coords_and_id *)a;
00085   struct coords_and_id *ib = (struct coords_and_id *)b;
00086   if(ia->x1 == ib->x1) {
00087     if(ia->y1 == ib->y1) {
00088       if(ia->z1 == ib->z1) {
00089         if(ia->x2 == ib->x2) {
00090           if(ia->y2 == ib->y2) {
00091             if(ia->z2 == ib->z2) {
00092               return ia->surf_id - ib->surf_id;
00093             } else {
00094               return (ia->z2 > ib->z2) - (ia->z2 < ib->z2);
00095             }
00096           } else {
00097             return (ia->y2 > ib->y2) - (ia->y2 < ib->y2);
00098           }
00099         } else {
00100           return (ia->x2 > ib->x2) - (ia->x2 < ib->x2);
00101         }
00102       } else {
00103         return (ia->z1 > ib->z1) - (ia->z1 < ib->z1);
00104       }
00105     } else {
00106       return (ia->y1 > ib->y1) - (ia->y1 < ib->y1);;
00107     }
00108   } else {
00109     return (ia->x1 > ib->x1) - (ia->x1 < ib->x1);
00110   }
00111   /* float comparison: returns negative if b > a 
00112      and positive if a > b. We multiplied result by 100.0
00113      to preserve decimal fraction */
00114 } 
00115  
00116 int main(int argc, char **argv) 
00117 {
00118 
00119   // ******************************************************************
00120   // Load the h5m file and create tags.
00121   // ******************************************************************
00122 
00123   clock_t start_time;
00124   start_time = clock();
00125   // check input args
00126   
00127   if( argc < 2 || argc > 5) 
00128     {
00129     std::cout << "To check using topology of facet points:              " << std::endl;
00130     std::cout << "./check_watertight <filename> <verbose(true or false)>" << std::endl;
00131     std::cout << "To check using geometry tolerance of facet points:    " << std::endl;
00132     std::cout << "./check_watertight <filename> <verbose(true or false)> <tolerance>" << std::endl;
00133     return 1;
00134     }
00135 
00136   // load file and get tolerance from input argument
00137   ErrorCode result;
00138   std::string filename = argv[1]; //set filename
00139   EntityHandle input_set;
00140   result = MBI()->create_meshset( MESHSET_SET, input_set ); //create handle to meshset
00141   if(MB_SUCCESS != result)
00142     {
00143       return result;
00144     }
00145 
00146   result = MBI()->load_file( filename.c_str(), &input_set ); //load the file into the meshset
00147   if(MB_SUCCESS != result)
00148     {
00149       // failed to load the file
00150       std::cout << "could not load file" << std::endl;
00151       return result;
00152     }
00153 
00154   double tol = 1e-04; // tolerance for two verts to be considered the same
00155                       // used in geometry checking a model
00156                       // set to this default value, unless user specifies otherwize, see below
00157   bool check_topology, verbose;
00158 
00159   if(2 == argc) // set topological check
00160     {
00161       std::cout << "topology check" << std::endl;
00162       check_topology = true;
00163       verbose = false;
00164     } 
00165   else if (3 == argc)  // set topological check with different tolerance
00166     {
00167       std::cout << "topology check" << std::endl;
00168       check_topology = true;
00169       const std::string verbose_string = argv[2];
00170       verbose = ( 0==verbose_string.compare("true") );
00171     } 
00172   else // otherwise do geometry check
00173     {
00174       std::cout << "geometry check";
00175       check_topology = false;
00176       if(!argv[3]){
00177         std::cout << "No tolerance specified for geometric check. Using default of " << tol << std::endl;
00178       }
00179       else tol = atof( argv[3] );
00180       std::cout<< " tolerance=" << tol << std::endl;
00181       const std::string verbose_string = argv[2];
00182       verbose = ( 0==verbose_string.compare("true") );
00183 
00184     }
00185 
00186   // replaced much of this code with a more modular version in check_watertight_func for testing purposes 
00187   std::set<int> leaky_surfs, leaky_vols;
00188   bool sealed, test;
00189   test=false;
00190   // is the order of the optional variables going to be a problem?
00191   // (i.e. we 'skipped' the variable test)
00192   result=cw_func::check_mesh_for_watertightness( input_set, tol, sealed, test, verbose, check_topology);
00193   if(gen::error(MB_SUCCESS!=result, "could not check model for watertightness")) return result;
00194 
00195   clock_t end_time = clock();
00196   std::cout << (double) (end_time-start_time)/CLOCKS_PER_SEC << " seconds" << std::endl;
00197  
00198 }
00199 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines