Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
cgm2moab.cpp
Go to the documentation of this file.
1 #include "cgm2moab.hpp"
2 #include "moab/ProgOptions.hpp"
3 
4 #include "moab/Core.hpp"
5 #include "moab/Range.hpp"
6 #include "moab/CartVect.hpp"
7 #include "MBTagConventions.hpp"
8 #include "moab/GeomQueryTool.hpp"
9 #include "moab/GeomTopoTool.hpp"
10 
11 #include <sstream>
12 #include <iomanip>
13 #include <cstdlib>
14 #include <algorithm>
15 #include <cstdio>
16 
17 using namespace moab;
18 
19 bool verbose = false;
20 
21 void chkerr( Interface* mbi, ErrorCode code, int line, const char* file )
22 {
23 
24  if( code != MB_SUCCESS )
25  {
26  std::cerr << "Error: " << mbi->get_error_string( code ) << " (" << code << ")" << std::endl;
27  std::string message;
28  if( MB_SUCCESS == mbi->get_last_error( message ) && !message.empty() )
29  std::cerr << "Error message: " << message << std::endl;
30  std::string fname = file;
31  size_t slash = fname.find_last_of( '/' );
32  if( slash != fname.npos )
33  {
34  fname = fname.substr( slash + 1 );
35  }
36  std::cerr << "At " << fname << " line: " << line << std::endl;
37  std::exit( EXIT_FAILURE );
38  }
39 }
40 
41 void chkerr( Interface& mbi, ErrorCode code, int line, const char* file )
42 {
43  chkerr( &mbi, code, line, file );
44 }
45 
46 void chkerr( GeomTopoTool& gtt, ErrorCode code, int line, const char* file )
47 {
48  chkerr( gtt.get_moab_instance(), code, line, file );
49 }
50 
51 /**
52  * Estimate the volume of the surface (actually multiplied by a factor of six).
53  * See DagMC::measure_volume, from which this code is borrowed, for algorithmic details.
54  * For our purpose, all that matters is the signedness.
55  *
56  * @param offset Offset to apply to surface to avoid a zero result.
57  */
59  const EntityHandle surf_set,
60  const CartVect& offset,
61  double& signed_volume )
62 {
63  ErrorCode rval;
64  Range tris;
65  rval = MBI->get_entities_by_type( surf_set, MBTRI, tris );
66  if( MB_SUCCESS != rval ) return rval;
67 
68  signed_volume = 0.0;
69  const EntityHandle* conn;
70  int len;
71  CartVect coords[3];
72  for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
73  {
74  rval = MBI->get_connectivity( *j, conn, len, true );
75  if( MB_SUCCESS != rval ) return rval;
76  if( 3 != len ) return MB_INVALID_SIZE;
77  rval = MBI->get_coords( conn, 3, coords[0].array() );
78  if( MB_SUCCESS != rval ) return rval;
79 
80  // Apply offset to avoid calculating 0 for cases when the origin is in the
81  // plane of the surface.
82  for( unsigned int k = 0; k < 3; ++k )
83  {
84  coords[k] += offset;
85  }
86 
87  coords[1] -= coords[0];
88  coords[2] -= coords[0];
89  signed_volume += ( coords[0] % ( coords[1] * coords[2] ) );
90  }
91  return MB_SUCCESS;
92 }
93 
94 /**
95  * Replace the triangles in an old surface with those in a new surface, ensuring
96  * that their surfaces senses match before the replacement
97  */
99  EntityHandle old_surf,
100  EntityHandle old_file_set,
101  EntityHandle new_surf,
102  const Tag& senseTag )
103 {
104 
105  ErrorCode rval;
106 
107  // Get the signed volume for each surface representation. If a volume comes
108  // back as zero, it's probably because a planar surface passes through the
109  // origin. In that case, try applying random offsets until reasonable
110  // values are returned.
111 
112  CartVect offset; // start with no offset
113  const double min_vol = 0.1; // try again if abs(vol) < this value
114 
115  double old_vol = 0, new_vol = 0;
116 
117  bool success = false;
118  int num_attempts = 100;
119 
120  while( num_attempts-- > 0 )
121  {
122 
123  rval = get_signed_volume( mbi, old_surf, offset, old_vol );
124  if( MB_SUCCESS != rval ) return rval;
125 
126  rval = get_signed_volume( mbi, new_surf, offset, new_vol );
127  if( MB_SUCCESS != rval ) return rval;
128 
129  if( std::fabs( old_vol ) >= min_vol && std::fabs( new_vol ) >= min_vol )
130  {
131  success = true;
132  break;
133  }
134 
135  // haven't succeeded yet: randomize the offset vector
136  const int max_random = 10;
137  for( int i = 0; i < 3; ++i )
138  {
139  offset[i] = std::rand() % max_random;
140  }
141  }
142 
143  if( !success )
144  {
145  std::cerr << "Error: could not calculate a surface volume" << std::endl;
146  return MB_FAILURE;
147  }
148 
149  // If the sign is different, reverse the old surf senses so that both
150  // representations have the same signed volume.
151  if( ( old_vol < 0 && new_vol > 0 ) || ( old_vol > 0 && new_vol < 0 ) )
152  {
153 
154  EntityHandle old_surf_volumes[2];
155  rval = mbi->tag_get_data( senseTag, &old_surf, 1, old_surf_volumes );
156  if( MB_SUCCESS != rval ) return rval;
157 
158  std::swap( old_surf_volumes[0], old_surf_volumes[1] );
159 
160  rval = mbi->tag_set_data( senseTag, &old_surf, 1, old_surf_volumes );
161  if( MB_SUCCESS != rval ) return rval;
162  }
163 
164  int num_old_tris, num_new_tris;
165 
166  // Remove the tris from the old surf. Also remove them from the
167  // old_file_set because it is not TRACKING.
168  Range old_tris;
169  rval = mbi->get_entities_by_type( old_surf, MBTRI, old_tris );
170  if( MB_SUCCESS != rval ) return rval;
171  num_old_tris = old_tris.size();
172  rval = mbi->remove_entities( old_surf, old_tris );
173  if( MB_SUCCESS != rval ) return rval;
174  rval = mbi->remove_entities( old_file_set, old_tris );
175  if( MB_SUCCESS != rval ) return rval;
176  rval = mbi->delete_entities( old_tris );
177  if( MB_SUCCESS != rval ) return rval;
178 
179  // Add the new_surf's triangles to the old_surf
180  Range new_tris;
181  rval = mbi->get_entities_by_type( new_surf, MBTRI, new_tris );
182  if( MB_SUCCESS != rval ) return rval;
183  num_new_tris = new_tris.size();
184  rval = mbi->add_entities( old_surf, new_tris );
185  if( MB_SUCCESS != rval ) return rval;
186 
187  if( verbose )
188  {
189  std::cout << num_old_tris << " tris -> " << num_new_tris << " tris" << std::endl;
190  }
191 
192  return MB_SUCCESS;
193 }
194 
195 /**
196  * Given an "old" file and a "new" file, replace the facets in any surface of the old
197  * file with facets from the new file.
198  */
200  const EntityHandle old_file_set,
201  const EntityHandle new_file_set,
202  const Tag& idTag,
203  const Tag& dimTag,
204  const Tag& senseTag )
205 {
206  ErrorCode rval;
207 
208  const int two = 2;
209  const Tag tags[2] = { dimTag, idTag };
210  const void* tag_vals[2] = { &two, NULL };
211 
212  Range old_surfs;
213  rval = mbi->get_entities_by_type_and_tag( old_file_set, MBENTITYSET, &dimTag, tag_vals, 1, old_surfs );
214  if( MB_SUCCESS != rval ) return rval;
215 
216  int count = 0;
217 
218  for( Range::iterator i = old_surfs.begin(); i != old_surfs.end(); ++i )
219  {
220  EntityHandle old_surf = *i;
221 
222  int surf_id;
223  rval = mbi->tag_get_data( idTag, &old_surf, 1, &surf_id );
224  if( MB_SUCCESS != rval ) return rval;
225 
226  Range new_surf_range;
227  tag_vals[1] = &surf_id;
228  rval = mbi->get_entities_by_type_and_tag( new_file_set, MBENTITYSET, tags, tag_vals, 2, new_surf_range );
229  if( MB_SUCCESS != rval ) return rval;
230 
231  if( new_surf_range.size() != 1 )
232  {
233  if( new_surf_range.size() > 1 )
234  {
235  std::cerr << "Warning: surface " << surf_id << " has more than one representation in new file"
236  << std::endl;
237  }
238  continue;
239  }
240 
241  // Now we have found a surf in new_file_set to replace an old surf
242  EntityHandle new_surf = new_surf_range[0];
243 
244  // TODO: check for quads and convert to triangles
245 
246  if( verbose )
247  {
248  std::cout << "Surface " << surf_id << ": " << std::flush;
249  }
250  rval = replace_surface( mbi, old_surf, old_file_set, new_surf, senseTag );
251  if( MB_SUCCESS != rval ) return rval;
252  count++;
253  }
254 
255  std::cout << "Replaced " << count << " surface" << ( count == 1 ? "." : "s." ) << std::endl;
256 
257  return MB_SUCCESS;
258 }
259 
260 int main( int argc, char* argv[] )
261 {
262 
263  ProgOptions po( "cgm2moab: a tool for preprocessing CAD and mesh files for analysis" );
264 
265  std::string input_file;
266  std::string output_file = "dagmc_preproc_out.h5m";
267  int grid = 50;
268 
269  po.addOpt< void >( ",v", "Verbose output", &verbose );
270  po.addOpt< std::string >( "outmesh,o", "Specify output file name (default " + output_file + ")", &output_file );
271  po.addOpt< void >( "no-outmesh,", "Do not write an output mesh" );
272  po.addOpt< std::string >( ",m", "Specify alternate input mesh to override surfaces in input_file" );
273  po.addOpt< std::string >( "obb-vis,O", "Specify obb visualization output file (default none)" );
274  po.addOpt< int >( "obb-vis-divs", "Resolution of obb visualization grid (default 50)", &grid );
275  po.addOpt< void >( "obb-stats,S", "Print obb statistics. With -v, print verbose statistics." );
276  po.addOpt< std::vector< int > >( "vols,V",
277  "Specify a set of volumes (applies to --obb_vis and --obb_stats, default all)" );
278  po.addOptionHelpHeading( "Options for loading CAD files" );
279  po.addOpt< double >( "ftol,f", "Faceting distance tolerance", po.add_cancel_opt );
280  po.addOpt< double >( "ltol,l", "Faceting edge length tolerance", po.add_cancel_opt );
281  po.addOpt< int >( "atol,a", "Faceting normal angle tolerance (degrees)", po.add_cancel_opt );
282  po.addOpt< void >( "all-warnings", "Verbose warnings about attributes and curve tolerances" );
283  po.addOpt< void >( "no-attribs", "Do not actuate CGM attributes" );
284  po.addOpt< void >( "fatal_curves", "Fatal Error when curves cannot be faceted" );
285 
286  po.addRequiredArg< std::string >( "input_file", "Path to input file for preprocessing", &input_file );
287 
288  po.parseCommandLine( argc, argv );
289 
290  /* Check that user has asked for at least one useful thing to be done */
291  bool obb_task = po.numOptSet( "obb-vis" ) || po.numOptSet( "obb-stats" );
292 
293  if( po.numOptSet( "no-outmesh" ) && !obb_task )
294  {
295  po.error( "Nothing to do. Please specify an OBB-related option, or remove --no_outmesh." );
296  }
297 
298  /* Load input file, with CAD processing options, if specified */
299  std::string options;
300 #define OPTION_APPEND( X ) \
301  { \
302  if( options.length() ) options += ";"; \
303  options += ( X ); \
304  }
305 
306  if( po.numOptSet( "no-attribs" ) )
307  {
308  OPTION_APPEND( "CGM_ATTRIBS=no" );
309  }
310 
311  if( po.numOptSet( "fatal_curves" ) )
312  {
313  OPTION_APPEND( "FATAL_ON_CURVES" );
314  }
315 
316  if( po.numOptSet( "all-warnings" ) )
317  {
318  OPTION_APPEND( "VERBOSE_CGM_WARNINGS" );
319  }
320 
321  // This is more roundabout than necessary, but we don't want *any* of the CGM-specific options
322  // to appear in the option string unless they were explicitly requested
323  double tol;
324  static const int tol_prec = 12;
325  if( po.getOpt( "ftol", &tol ) )
326  {
327  std::stringstream s;
328  s << "FACET_DISTANCE_TOLERANCE=" << std::setprecision( tol_prec ) << tol;
329  OPTION_APPEND( s.str() );
330  }
331 
332  if( po.getOpt( "ltol", &tol ) )
333  {
334  std::stringstream s;
335  s << "MAX_FACET_EDGE_LENGTH=" << std::setprecision( tol_prec ) << tol;
336  OPTION_APPEND( s.str() );
337  }
338 
339  int atol;
340  if( po.getOpt( "atol", &atol ) )
341  {
342  std::stringstream s;
343  s << "FACET_NORMAL_TOLERANCE=" << atol;
344  OPTION_APPEND( s.str() );
345  }
346 
347 #undef OPTION_APPEND
348 
349  /* Load main input file */
350  if( verbose )
351  {
352  std::cout << "Loading file " << input_file << std::endl;
353  if( options.length() ) std::cout << "Option string: " << options << std::endl;
354  }
355 
356  EntityHandle input_file_set;
357  ErrorCode ret;
358  Core mbi;
359 
360  ret = mbi.create_meshset( 0, input_file_set );
361  CHECKERR( mbi, ret );
362 
363  ret = mbi.load_file( input_file.c_str(), &input_file_set, options.c_str() );
364  if( ret == MB_UNHANDLED_OPTION )
365  {
366  std::cerr << "Warning: unhandled option while loading input_file, will proceed anyway" << std::endl;
367  }
368  else
369  {
370  CHECKERR( mbi, ret );
371  }
372 
373  /* Iterate through any -m alternate mesh files and replace surfaces */
374 
375  std::vector< std::string > m_list;
376  po.getOptAllArgs( ",m", m_list );
377 
378  if( m_list.size() > 0 )
379  {
380 
381  if( obb_task )
382  {
383  std::cerr << "Warning: using obb features in conjunction with -m may not work correctly!" << std::endl;
384  }
385 
386  // Create tags
387  Tag dimTag, idTag, senseTag;
389  CHECKERR( mbi, ret );
390 
391  idTag = mbi.globalId_tag();
392 
393  ret = mbi.tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, senseTag, MB_TAG_SPARSE | MB_TAG_CREAT );
394  CHECKERR( mbi, ret );
395 
396  for( std::vector< std::string >::iterator i = m_list.begin(); i != m_list.end(); ++i )
397  {
398  std::cout << "Loading alternate mesh file " << *i << std::endl;
399 
400  EntityHandle alt_file_set;
401  ret = mbi.create_meshset( 0, alt_file_set );
402  CHECKERR( mbi, ret );
403 
404  ret = mbi.load_file( ( *i ).c_str(), &alt_file_set );
405  CHECKERR( mbi, ret );
406 
407  if( verbose ) std::cout << "Merging input surfaces..." << std::flush;
408 
409  ret = merge_input_surfs( &mbi, input_file_set, alt_file_set, idTag, dimTag, senseTag );
410  CHECKERR( mbi, ret );
411 
412  if( verbose ) std::cout << "done." << std::endl;
413  }
414  }
415 
416  /* Write output file */
417 
418  if( !po.numOptSet( "no-outmesh" ) )
419  {
420  if( verbose )
421  {
422  std::cout << "Writing " << output_file << std::endl;
423  }
424  ret = mbi.write_file( output_file.c_str(), NULL, NULL, &input_file_set, 1 );
425  CHECKERR( mbi, ret );
426  }
427 
428  /* OBB statistics and visualization */
429  if( obb_task )
430  {
431 
432  if( verbose )
433  {
434  std::cout << "Loading data into GeomTopoTool" << std::endl;
435  }
436  GeomTopoTool* gtt = new GeomTopoTool( &mbi, false );
437  ret = gtt->find_geomsets();
438  CHECKERR( *gtt, ret );
439  ret = gtt->construct_obb_trees();
440  CHECKERR( *gtt, ret );
441  }
442  return 0;
443 }