MeshKit
1.0
|
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