MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "cgm2moab.hpp" 00002 #include "moab/ProgOptions.hpp" 00003 00004 #include "moab/Core.hpp" 00005 #include "moab/Range.hpp" 00006 #include "moab/CartVect.hpp" 00007 #include "MBTagConventions.hpp" 00008 #include "moab/GeomQueryTool.hpp" 00009 #include "moab/GeomTopoTool.hpp" 00010 00011 #include <sstream> 00012 #include <iomanip> 00013 #include <cstdlib> 00014 #include <algorithm> 00015 #include <cstdio> 00016 00017 using namespace moab; 00018 00019 bool verbose = false; 00020 00021 void chkerr( Interface* mbi, ErrorCode code, int line, const char* file ) 00022 { 00023 00024 if( code != MB_SUCCESS ) 00025 { 00026 std::cerr << "Error: " << mbi->get_error_string( code ) << " (" << code << ")" << std::endl; 00027 std::string message; 00028 if( MB_SUCCESS == mbi->get_last_error( message ) && !message.empty() ) 00029 std::cerr << "Error message: " << message << std::endl; 00030 std::string fname = file; 00031 size_t slash = fname.find_last_of( '/' ); 00032 if( slash != fname.npos ) { fname = fname.substr( slash + 1 ); } 00033 std::cerr << "At " << fname << " line: " << line << std::endl; 00034 std::exit( EXIT_FAILURE ); 00035 } 00036 } 00037 00038 void chkerr( Interface& mbi, ErrorCode code, int line, const char* file ) 00039 { 00040 chkerr( &mbi, code, line, file ); 00041 } 00042 00043 void chkerr( GeomTopoTool& gtt, ErrorCode code, int line, const char* file ) 00044 { 00045 chkerr( gtt.get_moab_instance(), code, line, file ); 00046 } 00047 00048 /** 00049 * Estimate the volume of the surface (actually multiplied by a factor of six). 00050 * See DagMC::measure_volume, from which this code is borrowed, for algorithmic details. 00051 * For our purpose, all that matters is the signedness. 00052 * 00053 * @param offset Offset to apply to surface to avoid a zero result. 00054 */ 00055 static ErrorCode get_signed_volume( Interface* MBI, const EntityHandle surf_set, const CartVect& offset, 00056 double& signed_volume ) 00057 { 00058 ErrorCode rval; 00059 Range tris; 00060 rval = MBI->get_entities_by_type( surf_set, MBTRI, tris ); 00061 if( MB_SUCCESS != rval ) return rval; 00062 00063 signed_volume = 0.0; 00064 const EntityHandle* conn; 00065 int len; 00066 CartVect coords[3]; 00067 for( Range::iterator j = tris.begin(); j != tris.end(); ++j ) 00068 { 00069 rval = MBI->get_connectivity( *j, conn, len, true ); 00070 if( MB_SUCCESS != rval ) return rval; 00071 if( 3 != len ) return MB_INVALID_SIZE; 00072 rval = MBI->get_coords( conn, 3, coords[0].array() ); 00073 if( MB_SUCCESS != rval ) return rval; 00074 00075 // Apply offset to avoid calculating 0 for cases when the origin is in the 00076 // plane of the surface. 00077 for( unsigned int k = 0; k < 3; ++k ) 00078 { 00079 coords[k] += offset; 00080 } 00081 00082 coords[1] -= coords[0]; 00083 coords[2] -= coords[0]; 00084 signed_volume += ( coords[0] % ( coords[1] * coords[2] ) ); 00085 } 00086 return MB_SUCCESS; 00087 } 00088 00089 /** 00090 * Replace the triangles in an old surface with those in a new surface, ensuring 00091 * that their surfaces senses match before the replacement 00092 */ 00093 static ErrorCode replace_surface( Interface* mbi, EntityHandle old_surf, EntityHandle old_file_set, 00094 EntityHandle new_surf, const Tag& senseTag ) 00095 { 00096 00097 ErrorCode rval; 00098 00099 // Get the signed volume for each surface representation. If a volume comes 00100 // back as zero, it's probably because a planar surface passes through the 00101 // origin. In that case, try applying random offsets until reasonable 00102 // values are returned. 00103 00104 CartVect offset; // start with no offset 00105 const double min_vol = 0.1; // try again if abs(vol) < this value 00106 00107 double old_vol = 0, new_vol = 0; 00108 00109 bool success = false; 00110 int num_attempts = 100; 00111 00112 while( num_attempts-- > 0 ) 00113 { 00114 00115 rval = get_signed_volume( mbi, old_surf, offset, old_vol ); 00116 if( MB_SUCCESS != rval ) return rval; 00117 00118 rval = get_signed_volume( mbi, new_surf, offset, new_vol ); 00119 if( MB_SUCCESS != rval ) return rval; 00120 00121 if( std::fabs( old_vol ) >= min_vol && std::fabs( new_vol ) >= min_vol ) 00122 { 00123 success = true; 00124 break; 00125 } 00126 00127 // haven't succeeded yet: randomize the offset vector 00128 const int max_random = 10; 00129 for( int i = 0; i < 3; ++i ) 00130 { 00131 offset[i] = std::rand() % max_random; 00132 } 00133 } 00134 00135 if( !success ) 00136 { 00137 std::cerr << "Error: could not calculate a surface volume" << std::endl; 00138 return MB_FAILURE; 00139 } 00140 00141 // If the sign is different, reverse the old surf senses so that both 00142 // representations have the same signed volume. 00143 if( ( old_vol < 0 && new_vol > 0 ) || ( old_vol > 0 && new_vol < 0 ) ) 00144 { 00145 00146 EntityHandle old_surf_volumes[2]; 00147 rval = mbi->tag_get_data( senseTag, &old_surf, 1, old_surf_volumes ); 00148 if( MB_SUCCESS != rval ) return rval; 00149 00150 std::swap( old_surf_volumes[0], old_surf_volumes[1] ); 00151 00152 rval = mbi->tag_set_data( senseTag, &old_surf, 1, old_surf_volumes ); 00153 if( MB_SUCCESS != rval ) return rval; 00154 } 00155 00156 int num_old_tris, num_new_tris; 00157 00158 // Remove the tris from the old surf. Also remove them from the 00159 // old_file_set because it is not TRACKING. 00160 Range old_tris; 00161 rval = mbi->get_entities_by_type( old_surf, MBTRI, old_tris ); 00162 if( MB_SUCCESS != rval ) return rval; 00163 num_old_tris = old_tris.size(); 00164 rval = mbi->remove_entities( old_surf, old_tris ); 00165 if( MB_SUCCESS != rval ) return rval; 00166 rval = mbi->remove_entities( old_file_set, old_tris ); 00167 if( MB_SUCCESS != rval ) return rval; 00168 rval = mbi->delete_entities( old_tris ); 00169 if( MB_SUCCESS != rval ) return rval; 00170 00171 // Add the new_surf's triangles to the old_surf 00172 Range new_tris; 00173 rval = mbi->get_entities_by_type( new_surf, MBTRI, new_tris ); 00174 if( MB_SUCCESS != rval ) return rval; 00175 num_new_tris = new_tris.size(); 00176 rval = mbi->add_entities( old_surf, new_tris ); 00177 if( MB_SUCCESS != rval ) return rval; 00178 00179 if( verbose ) { std::cout << num_old_tris << " tris -> " << num_new_tris << " tris" << std::endl; } 00180 00181 return MB_SUCCESS; 00182 } 00183 00184 /** 00185 * Given an "old" file and a "new" file, replace the facets in any surface of the old 00186 * file with facets from the new file. 00187 */ 00188 static ErrorCode merge_input_surfs( Interface* mbi, const EntityHandle old_file_set, const EntityHandle new_file_set, 00189 const Tag& idTag, const Tag& dimTag, const Tag& senseTag ) 00190 { 00191 ErrorCode rval; 00192 00193 const int two = 2; 00194 const Tag tags[2] = { dimTag, idTag }; 00195 const void* tag_vals[2] = { &two, NULL }; 00196 00197 Range old_surfs; 00198 rval = mbi->get_entities_by_type_and_tag( old_file_set, MBENTITYSET, &dimTag, tag_vals, 1, old_surfs ); 00199 if( MB_SUCCESS != rval ) return rval; 00200 00201 int count = 0; 00202 00203 for( Range::iterator i = old_surfs.begin(); i != old_surfs.end(); ++i ) 00204 { 00205 EntityHandle old_surf = *i; 00206 00207 int surf_id; 00208 rval = mbi->tag_get_data( idTag, &old_surf, 1, &surf_id ); 00209 if( MB_SUCCESS != rval ) return rval; 00210 00211 Range new_surf_range; 00212 tag_vals[1] = &surf_id; 00213 rval = mbi->get_entities_by_type_and_tag( new_file_set, MBENTITYSET, tags, tag_vals, 2, new_surf_range ); 00214 if( MB_SUCCESS != rval ) return rval; 00215 00216 if( new_surf_range.size() != 1 ) 00217 { 00218 if( new_surf_range.size() > 1 ) 00219 { 00220 std::cerr << "Warning: surface " << surf_id << " has more than one representation in new file" 00221 << std::endl; 00222 } 00223 continue; 00224 } 00225 00226 // Now we have found a surf in new_file_set to replace an old surf 00227 EntityHandle new_surf = new_surf_range[0]; 00228 00229 // TODO: check for quads and convert to triangles 00230 00231 if( verbose ) { std::cout << "Surface " << surf_id << ": " << std::flush; } 00232 rval = replace_surface( mbi, old_surf, old_file_set, new_surf, senseTag ); 00233 if( MB_SUCCESS != rval ) return rval; 00234 count++; 00235 } 00236 00237 std::cout << "Replaced " << count << " surface" << ( count == 1 ? "." : "s." ) << std::endl; 00238 00239 return MB_SUCCESS; 00240 } 00241 00242 int main( int argc, char* argv[] ) 00243 { 00244 00245 ProgOptions po( "cgm2moab: a tool for preprocessing CAD and mesh files for analysis" ); 00246 00247 std::string input_file; 00248 std::string output_file = "dagmc_preproc_out.h5m"; 00249 int grid = 50; 00250 00251 po.addOpt< void >( ",v", "Verbose output", &verbose ); 00252 po.addOpt< std::string >( "outmesh,o", "Specify output file name (default " + output_file + ")", &output_file ); 00253 po.addOpt< void >( "no-outmesh,", "Do not write an output mesh" ); 00254 po.addOpt< std::string >( ",m", "Specify alternate input mesh to override surfaces in input_file" ); 00255 po.addOpt< std::string >( "obb-vis,O", "Specify obb visualization output file (default none)" ); 00256 po.addOpt< int >( "obb-vis-divs", "Resolution of obb visualization grid (default 50)", &grid ); 00257 po.addOpt< void >( "obb-stats,S", "Print obb statistics. With -v, print verbose statistics." ); 00258 po.addOpt< std::vector< int > >( "vols,V", 00259 "Specify a set of volumes (applies to --obb_vis and --obb_stats, default all)" ); 00260 po.addOptionHelpHeading( "Options for loading CAD files" ); 00261 po.addOpt< double >( "ftol,f", "Faceting distance tolerance", po.add_cancel_opt ); 00262 po.addOpt< double >( "ltol,l", "Faceting edge length tolerance", po.add_cancel_opt ); 00263 po.addOpt< int >( "atol,a", "Faceting normal angle tolerance (degrees)", po.add_cancel_opt ); 00264 po.addOpt< void >( "all-warnings", "Verbose warnings about attributes and curve tolerances" ); 00265 po.addOpt< void >( "no-attribs", "Do not actuate CGM attributes" ); 00266 po.addOpt< void >( "fatal_curves", "Fatal Error when curves cannot be faceted" ); 00267 00268 po.addRequiredArg< std::string >( "input_file", "Path to input file for preprocessing", &input_file ); 00269 00270 po.parseCommandLine( argc, argv ); 00271 00272 /* Check that user has asked for at least one useful thing to be done */ 00273 bool obb_task = po.numOptSet( "obb-vis" ) || po.numOptSet( "obb-stats" ); 00274 00275 if( po.numOptSet( "no-outmesh" ) && !obb_task ) 00276 { po.error( "Nothing to do. Please specify an OBB-related option, or remove --no_outmesh." ); } 00277 00278 /* Load input file, with CAD processing options, if specified */ 00279 std::string options; 00280 #define OPTION_APPEND( X ) \ 00281 { \ 00282 if( options.length() ) options += ";"; \ 00283 options += ( X ); \ 00284 } 00285 00286 if( po.numOptSet( "no-attribs" ) ) { OPTION_APPEND( "CGM_ATTRIBS=no" ); } 00287 00288 if( po.numOptSet( "fatal_curves" ) ) { OPTION_APPEND( "FATAL_ON_CURVES" ); } 00289 00290 if( po.numOptSet( "all-warnings" ) ) { OPTION_APPEND( "VERBOSE_CGM_WARNINGS" ); } 00291 00292 // This is more roundabout than necessary, but we don't want *any* of the CGM-specific options 00293 // to appear in the option string unless they were explicitly requested 00294 double tol; 00295 static const int tol_prec = 12; 00296 if( po.getOpt( "ftol", &tol ) ) 00297 { 00298 std::stringstream s; 00299 s << "FACET_DISTANCE_TOLERANCE=" << std::setprecision( tol_prec ) << tol; 00300 OPTION_APPEND( s.str() ); 00301 } 00302 00303 if( po.getOpt( "ltol", &tol ) ) 00304 { 00305 std::stringstream s; 00306 s << "MAX_FACET_EDGE_LENGTH=" << std::setprecision( tol_prec ) << tol; 00307 OPTION_APPEND( s.str() ); 00308 } 00309 00310 int atol; 00311 if( po.getOpt( "atol", &atol ) ) 00312 { 00313 std::stringstream s; 00314 s << "FACET_NORMAL_TOLERANCE=" << atol; 00315 OPTION_APPEND( s.str() ); 00316 } 00317 00318 #undef OPTION_APPEND 00319 00320 /* Load main input file */ 00321 if( verbose ) 00322 { 00323 std::cout << "Loading file " << input_file << std::endl; 00324 if( options.length() ) std::cout << "Option string: " << options << std::endl; 00325 } 00326 00327 EntityHandle input_file_set; 00328 ErrorCode ret; 00329 Core mbi; 00330 00331 ret = mbi.create_meshset( 0, input_file_set ); 00332 CHECKERR( mbi, ret ); 00333 00334 ret = mbi.load_file( input_file.c_str(), &input_file_set, options.c_str() ); 00335 if( ret == MB_UNHANDLED_OPTION ) 00336 { std::cerr << "Warning: unhandled option while loading input_file, will proceed anyway" << std::endl; } 00337 else 00338 { 00339 CHECKERR( mbi, ret ); 00340 } 00341 00342 /* Iterate through any -m alternate mesh files and replace surfaces */ 00343 00344 std::vector< std::string > m_list; 00345 po.getOptAllArgs( ",m", m_list ); 00346 00347 if( m_list.size() > 0 ) 00348 { 00349 00350 if( obb_task ) 00351 { std::cerr << "Warning: using obb features in conjunction with -m may not work correctly!" << std::endl; } 00352 00353 // Create tags 00354 Tag dimTag, idTag, senseTag; 00355 ret = mbi.tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dimTag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00356 CHECKERR( mbi, ret ); 00357 00358 idTag = mbi.globalId_tag(); 00359 00360 ret = mbi.tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, senseTag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00361 CHECKERR( mbi, ret ); 00362 00363 for( std::vector< std::string >::iterator i = m_list.begin(); i != m_list.end(); ++i ) 00364 { 00365 std::cout << "Loading alternate mesh file " << *i << std::endl; 00366 00367 EntityHandle alt_file_set; 00368 ret = mbi.create_meshset( 0, alt_file_set ); 00369 CHECKERR( mbi, ret ); 00370 00371 ret = mbi.load_file( ( *i ).c_str(), &alt_file_set ); 00372 CHECKERR( mbi, ret ); 00373 00374 if( verbose ) std::cout << "Merging input surfaces..." << std::flush; 00375 00376 ret = merge_input_surfs( &mbi, input_file_set, alt_file_set, idTag, dimTag, senseTag ); 00377 CHECKERR( mbi, ret ); 00378 00379 if( verbose ) std::cout << "done." << std::endl; 00380 } 00381 } 00382 00383 /* Write output file */ 00384 00385 if( !po.numOptSet( "no-outmesh" ) ) 00386 { 00387 if( verbose ) { std::cout << "Writing " << output_file << std::endl; } 00388 ret = mbi.write_file( output_file.c_str(), NULL, NULL, &input_file_set, 1 ); 00389 CHECKERR( mbi, ret ); 00390 } 00391 00392 /* OBB statistics and visualization */ 00393 if( obb_task ) 00394 { 00395 00396 if( verbose ) { std::cout << "Loading data into GeomTopoTool" << std::endl; } 00397 GeomTopoTool* gtt = new GeomTopoTool( &mbi, false ); 00398 ret = gtt->find_geomsets(); 00399 CHECKERR( *gtt, ret ); 00400 ret = gtt->construct_obb_trees(); 00401 CHECKERR( *gtt, ret ); 00402 } 00403 return 0; 00404 }