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