Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
skin.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <set>
3 #include <limits>
4 #include <ctime>
5 #include <vector>
6 #include <cstdlib>
7 #include <cstdio>
8 #include <cassert>
9 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
10 #include <unistd.h>
11 #include <sys/types.h>
12 #include <sys/stat.h>
13 #endif
14 #include <fcntl.h>
15 #include "moab/Interface.hpp"
16 #include "MBTagConventions.hpp"
17 #include "moab/Core.hpp"
18 #include "moab/Range.hpp"
19 #include "moab/Skinner.hpp"
20 #include "moab/AdaptiveKDTree.hpp"
21 #include "moab/CN.hpp"
22 
23 using namespace moab;
24 
25 static void get_time_mem( double& tot_time, double& tot_mem );
26 
27 // Different platforms follow different conventions for usage
28 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
29 #include <sys/resource.h>
30 #endif
31 #ifdef SOLARIS
32 extern "C" int getrusage( int, struct rusage* );
33 #ifndef RUSAGE_SELF
34 #include </usr/ucbinclude/sys/rusage.h>
35 #endif
36 #endif
37 
38 const char DEFAULT_FIXED_TAG[] = "fixed";
39 const int MIN_EDGE_LEN_DENOM = 4;
40 
41 #define CHKERROR( A ) \
42  do \
43  { \
44  if( MB_SUCCESS != ( A ) ) \
45  { \
46  std::cerr << "Internal error at line " << __LINE__ << std::endl; \
47  return 3; \
48  } \
49  } while( false )
50 
51 static ErrorCode merge_duplicate_vertices( Interface&, double epsilon );
52 static ErrorCode min_edge_length( Interface&, double& result );
53 
54 static void usage( const char* argv0, bool help = false )
55 {
56  std::ostream& str = help ? std::cout : std::cerr;
57 
58  str << "Usage: " << argv0
59  << " [-b <block_num> [-b ...] ] [-l] [-m] [-M <n>] [-p] [-s <sideset_num>] [-S] [-t|-T "
60  "<name>] [-w] [-v|-V <n>]"
61  << " <input_file> [<output_file>]" << std::endl;
62  str << "Help : " << argv0 << " -h" << std::endl;
63  if( !help ) exit( 1 );
64 
65  str << "Options: " << std::endl;
66  str << "-a : Compute skin using vert-elem adjacencies (more memory, less time)." << std::endl;
67  str << "-b <block_num> : Compute skin only for material set/block <block_num>." << std::endl;
68  str << "-p : Print cpu & memory performance." << std::endl;
69  str << "-s <sideset_num> : Put skin in neumann set/sideset <sideset_num>." << std::endl;
70  str << "-S : Look for and use structured mesh information to speed up skinning." << std::endl;
71  str << "-t : Set '" << DEFAULT_FIXED_TAG << "' tag on skin vertices." << std::endl;
72  str << "-T <name> : Create tag with specified name and set to 1 on skin vertices." << std::endl;
73  str << "-w : Write out whole mesh (otherwise just writes skin)." << std::endl;
74  str << "-m : consolidate duplicate vertices" << std::endl;
75  str << "-M <n> : consolidate duplicate vertices with specified tolerance. "
76  "(Default: min_edge_length/"
77  << MIN_EDGE_LEN_DENOM << ")" << std::endl;
78  str << "-l : List total numbers of entities and vertices in skin." << std::endl;
79  exit( 0 );
80 }
81 
82 int main( int argc, char* argv[] )
83 {
84  int i = 1;
85  std::vector< int > matsets;
86  int neuset_num = -1;
87  bool write_tag = false, write_whole_mesh = false;
88  bool print_perf = false;
89  bool use_vert_elem_adjs = false;
90  bool merge_vertices = false;
91  double merge_epsilon = -1;
92  bool list_skin = false;
93  bool use_scd = false;
94  const char* fixed_tag = DEFAULT_FIXED_TAG;
95  const char *input_file = 0, *output_file = 0;
96 
97  bool no_more_flags = false;
98  char* endptr = 0;
99  long block = 0; // initialize to eliminate compiler warning
100  while( i < argc )
101  {
102  if( !no_more_flags && argv[i][0] == '-' )
103  {
104  const int f = i++;
105  for( int j = 1; argv[f][j]; ++j )
106  {
107  switch( argv[f][j] )
108  {
109  case 'a':
110  use_vert_elem_adjs = true;
111  break;
112  case 'p':
113  print_perf = true;
114  break;
115  case 't':
116  write_tag = true;
117  break;
118  case 'w':
119  write_whole_mesh = true;
120  break;
121  case 'm':
122  merge_vertices = true;
123  break;
124  case '-':
125  no_more_flags = true;
126  break;
127  case 'h':
128  usage( argv[0], true );
129  break;
130  case 'l':
131  list_skin = true;
132  break;
133  case 'S':
134  use_scd = true;
135  break;
136  case 'b':
137  if( i == argc || 0 >= ( block = strtol( argv[i], &endptr, 0 ) ) || *endptr )
138  {
139  std::cerr << "Expected positive integer following '-b' flag" << std::endl;
140  usage( argv[0] );
141  }
142  matsets.push_back( (int)block );
143  ++i;
144  break;
145  case 's':
146  if( i == argc || 0 >= ( neuset_num = strtol( argv[i], &endptr, 0 ) ) || *endptr )
147  {
148  std::cerr << "Expected positive integer following '-s' flag" << std::endl;
149  usage( argv[0] );
150  }
151  ++i;
152  break;
153  case 'T':
154  if( i == argc || argv[i][0] == '-' )
155  {
156  std::cerr << "Expected tag name following '-T' flag" << std::endl;
157  usage( argv[0] );
158  }
159  fixed_tag = argv[i++];
160  break;
161  case 'M':
162  if( i == argc || 0.0 > ( merge_epsilon = strtod( argv[i], &endptr ) ) || *endptr )
163  {
164  std::cerr << "Expected positive numeric value following '-M' flag" << std::endl;
165  usage( argv[0] );
166  }
167  merge_vertices = true;
168  ++i;
169  break;
170  default:
171  std::cerr << "Unrecognized flag: '" << argv[f][j] << "'" << std::endl;
172  usage( argv[0] );
173  break;
174  }
175  }
176  }
177  else if( input_file && output_file )
178  {
179  std::cerr << "Extra argument: " << argv[i] << std::endl;
180  usage( argv[0] );
181  }
182  else if( input_file )
183  {
184  output_file = argv[i++];
185  }
186  else
187  {
188  input_file = argv[i++];
189  }
190  }
191 
192  if( !input_file )
193  {
194  std::cerr << "No input file specified" << std::endl;
195  usage( argv[0] );
196  }
197 
198  ErrorCode result;
199  Core mbimpl;
200  Interface* iface = &mbimpl;
201 
202  if( print_perf )
203  {
204  double tmp_time1, tmp_mem1;
205  get_time_mem( tmp_time1, tmp_mem1 );
206  std::cout << "Before reading: cpu time = " << tmp_time1 << ", memory = " << tmp_mem1 / 1.0e6 << "MB."
207  << std::endl;
208  }
209 
210  // read input file
211  result = iface->load_mesh( input_file );
212  if( MB_SUCCESS != result )
213  {
214  std::cerr << "Failed to load \"" << input_file << "\"." << std::endl;
215  return 2;
216  }
217  std::cerr << "Read \"" << input_file << "\"" << std::endl;
218  if( print_perf )
219  {
220  double tmp_time2, tmp_mem2;
221  get_time_mem( tmp_time2, tmp_mem2 );
222  std::cout << "After reading: cpu time = " << tmp_time2 << ", memory = " << tmp_mem2 / 1.0e6 << "MB."
223  << std::endl;
224  }
225 
226  if( merge_vertices )
227  {
228  if( merge_epsilon < 0.0 )
229  {
230  if( MB_SUCCESS != min_edge_length( *iface, merge_epsilon ) )
231  {
232  std::cerr << "Error determining minimum edge length" << std::endl;
233  return 1;
234  }
235  merge_epsilon /= MIN_EDGE_LEN_DENOM;
236  }
237  if( MB_SUCCESS != merge_duplicate_vertices( *iface, merge_epsilon ) )
238  {
239  std::cerr << "Error merging duplicate vertices" << std::endl;
240  return 1;
241  }
242  }
243 
244  // get entities of largest dimension
245  int dim = 4;
246  Range entities;
247  while( entities.empty() && dim > 1 )
248  {
249  dim--;
250  result = iface->get_entities_by_dimension( 0, dim, entities );CHKERROR( result );
251  }
252 
253  Range skin_ents;
254  Tag matset_tag = 0, neuset_tag = 0;
255  result = iface->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, matset_tag );
256  if( MB_SUCCESS != result ) return 1;
257  result = iface->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, neuset_tag );
258  if( MB_SUCCESS != result ) return 1;
259 
260  if( matsets.empty() )
261  skin_ents = entities;
262  else
263  {
264  // get all entities in the specified blocks
265  if( 0 == matset_tag )
266  {
267  std::cerr << "Couldn't find any material sets in this mesh." << std::endl;
268  return 1;
269  }
270 
271  for( std::vector< int >::iterator vit = matsets.begin(); vit != matsets.end(); ++vit )
272  {
273  int this_matset = *vit;
274  const void* this_matset_ptr = &this_matset;
275  Range this_range, ent_range;
276  result =
277  iface->get_entities_by_type_and_tag( 0, MBENTITYSET, &matset_tag, &this_matset_ptr, 1, this_range );
278  if( MB_SUCCESS != result )
279  {
280  std::cerr << "Trouble getting material set #" << *vit << std::endl;
281  return 1;
282  }
283  else if( this_range.empty() )
284  {
285  std::cerr << "Warning: couldn't find material set " << *vit << std::endl;
286  continue;
287  }
288 
289  result = iface->get_entities_by_dimension( *this_range.begin(), dim, ent_range, true );
290  if( MB_SUCCESS != result ) continue;
291  skin_ents.merge( ent_range );
292  }
293  }
294 
295  if( skin_ents.empty() )
296  {
297  std::cerr << "No entities for which to compute skin; exiting." << std::endl;
298  return 1;
299  }
300 
301  if( use_vert_elem_adjs )
302  {
303  // make a call which we know will generate vert-elem adjs
304  Range dum_range;
305  result = iface->get_adjacencies( &( *skin_ents.begin() ), 1, 1, false, dum_range );
306  if( MB_SUCCESS != result ) return 1;
307  }
308 
309  double tmp_time = 0.0, tmp_mem = 0.0;
310  if( print_perf )
311  {
312  get_time_mem( tmp_time, tmp_mem );
313  std::cout << "Before skinning: cpu time = " << tmp_time << ", memory = " << tmp_mem / 1.0e6 << "MB."
314  << std::endl;
315  }
316 
317  // skin the mesh
318  Range forward_lower, reverse_lower;
319  Skinner tool( iface );
320  if( use_scd )
321  result = tool.find_skin( 0, skin_ents, false, forward_lower, NULL, false, true, true );
322  else
323  result = tool.find_skin( 0, skin_ents, false, forward_lower, &reverse_lower );
324  Range boundary;
325  boundary.merge( forward_lower );
326  boundary.merge( reverse_lower );
327  if( MB_SUCCESS != result || boundary.empty() )
328  {
329  std::cerr << "Mesh skinning failed." << std::endl;
330  return 3;
331  }
332 
333  if( list_skin )
334  {
335  Range skin_verts;
336  result = iface->get_adjacencies( boundary, 0, true, skin_verts, Interface::UNION );
337  std::cout << "Skin has ";
338  if( skin_ents.num_of_dimension( 3 ) )
339  std::cout << boundary.num_of_dimension( 2 ) << " faces and ";
340  else if( skin_ents.num_of_dimension( 2 ) )
341  std::cout << boundary.num_of_dimension( 1 ) << " edges and ";
342  std::cout << skin_verts.size() << " vertices." << std::endl;
343  }
344  if( write_tag )
345  {
346  // get tag handle
347  Tag tag;
348  int zero = 0;
349  result = iface->tag_get_handle( fixed_tag, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );CHKERROR( result );
350 
351  // Set tags
352  std::vector< int > ones;
353  Range bverts;
354  result = iface->get_adjacencies( boundary, 0, false, bverts, Interface::UNION );
355  if( MB_SUCCESS != result )
356  {
357  std::cerr << "Trouble getting vertices on boundary." << std::endl;
358  return 1;
359  }
360  ones.resize( bverts.size(), 1 );
361  result = iface->tag_set_data( tag, bverts, &ones[0] );CHKERROR( result );
362  }
363 
364  if( -1 != neuset_num )
365  {
366  // create a neumann set with these entities
367  if( 0 == neuset_tag )
368  {
369  result = iface->tag_get_handle( "NEUMANN_SET_TAG_NAME", 1, MB_TYPE_INTEGER, neuset_tag,
370  MB_TAG_SPARSE | MB_TAG_CREAT );
371  if( MB_SUCCESS != result || 0 == neuset_tag ) return 1;
372  }
373 
374  // always create a forward neumann set, assuming we have something in the set
375  EntityHandle forward_neuset = 0;
376  result = iface->create_meshset( MESHSET_SET, forward_neuset );
377  if( MB_SUCCESS != result || 0 == forward_neuset ) return 1;
378  result = iface->tag_set_data( neuset_tag, &forward_neuset, 1, &neuset_num );
379  if( MB_SUCCESS != result ) return 1;
380 
381  if( !forward_lower.empty() )
382  {
383  result = iface->add_entities( forward_neuset, forward_lower );
384  if( MB_SUCCESS != result ) return 1;
385  }
386  if( !reverse_lower.empty() )
387  {
388  EntityHandle reverse_neuset = 1;
389  result = iface->create_meshset( MESHSET_SET, reverse_neuset );
390  if( MB_SUCCESS != result || 0 == forward_neuset ) return 1;
391 
392  result = iface->add_entities( reverse_neuset, reverse_lower );
393  if( MB_SUCCESS != result ) return 1;
394  Tag sense_tag;
395  int dum_sense = 0;
396  result = iface->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
397  &dum_sense );
398  if( result != MB_SUCCESS ) return 1;
399  int sense_val = -1;
400  result = iface->tag_set_data( neuset_tag, &reverse_neuset, 1, &sense_val );
401  if( MB_SUCCESS != result ) return 0;
402  result = iface->add_entities( forward_neuset, &reverse_neuset, 1 );
403  if( MB_SUCCESS != result ) return 0;
404  }
405  }
406 
407  if( NULL != output_file && write_whole_mesh )
408  {
409 
410  // write output file
411  result = iface->write_mesh( output_file );
412  if( MB_SUCCESS != result )
413  {
414  std::cerr << "Failed to write \"" << output_file << "\"." << std::endl;
415  return 2;
416  }
417  std::cerr << "Wrote \"" << output_file << "\"" << std::endl;
418  }
419  else if( NULL != output_file )
420  {
421  // write only skin; write them as one set
422  EntityHandle skin_set;
423  result = iface->create_meshset( MESHSET_SET, skin_set );
424  if( MB_SUCCESS != result ) return 1;
425  result = iface->add_entities( skin_set, forward_lower );
426  if( MB_SUCCESS != result ) return 1;
427  result = iface->add_entities( skin_set, reverse_lower );
428  if( MB_SUCCESS != result ) return 1;
429 
430  int dum = 10000;
431  result = iface->tag_set_data( matset_tag, &skin_set, 1, &dum );
432  if( MB_SUCCESS != result ) return 1;
433 
434  result = iface->write_mesh( output_file, &skin_set, 1 );
435  if( MB_SUCCESS != result )
436  {
437  std::cerr << "Failed to write \"" << output_file << "\"." << std::endl;
438  return 2;
439  }
440  std::cerr << "Wrote \"" << output_file << "\"" << std::endl;
441  }
442 
443  if( print_perf )
444  {
445  double tot_time, tot_mem;
446  get_time_mem( tot_time, tot_mem );
447  std::cout << "Total cpu time = " << tot_time << " seconds." << std::endl;
448  std::cout << "Total skin cpu time = " << tot_time - tmp_time << " seconds." << std::endl;
449  std::cout << "Total memory = " << tot_mem / 1024 << " MB." << std::endl;
450  std::cout << "Total skin memory = " << ( tot_mem - tmp_mem ) / 1024 << " MB." << std::endl;
451  std::cout << "Entities: " << std::endl;
452  iface->list_entities( 0, 0 );
453  }
454 
455  return 0;
456 }
457 
458 #if defined( _MSC_VER ) || defined( __MINGW32__ )
459 void get_time_mem( double& tot_time, double& tot_mem )
460 {
461  tot_time = (double)clock() / CLOCKS_PER_SEC;
462  tot_mem = 0;
463 }
464 #else
465 void get_time_mem( double& tot_time, double& tot_mem )
466 {
467  struct rusage r_usage;
468  getrusage( RUSAGE_SELF, &r_usage );
469  double utime = (double)r_usage.ru_utime.tv_sec + ( (double)r_usage.ru_utime.tv_usec / 1.e6 );
470  double stime = (double)r_usage.ru_stime.tv_sec + ( (double)r_usage.ru_stime.tv_usec / 1.e6 );
471  tot_time = utime + stime;
472  tot_mem = 0;
473  if( 0 != r_usage.ru_maxrss )
474  {
475  tot_mem = (double)r_usage.ru_maxrss;
476  }
477  else
478  {
479  // this machine doesn't return rss - try going to /proc
480  // print the file name to open
481  char file_str[4096], dum_str[4096];
482  int file_ptr = open( "/proc/self/stat", O_RDONLY );
483  int file_len = read( file_ptr, file_str, sizeof( file_str ) - 1 );
484  if( file_len <= 0 )
485  {
486  close( file_ptr );
487  return;
488  }
489 
490  close( file_ptr );
491  file_str[file_len] = '\0';
492  // read the preceding fields and the ones we really want...
493  int dum_int;
494  unsigned int dum_uint, vm_size, rss;
495  int num_fields = sscanf( file_str,
496  "%d " // pid
497  "%s " // comm
498  "%c " // state
499  "%d %d %d %d %d " // ppid, pgrp, session, tty, tpgid
500  "%u %u %u %u %u " // flags, minflt, cminflt, majflt, cmajflt
501  "%d %d %d %d %d %d " // utime, stime, cutime, cstime, counter, priority
502  "%u %u " // timeout, itrealvalue
503  "%d " // starttime
504  "%u %u", // vsize, rss
505  &dum_int, dum_str, dum_str, &dum_int, &dum_int, &dum_int, &dum_int, &dum_int,
506  &dum_uint, &dum_uint, &dum_uint, &dum_uint, &dum_uint, &dum_int, &dum_int, &dum_int,
507  &dum_int, &dum_int, &dum_int, &dum_uint, &dum_uint, &dum_int, &vm_size, &rss );
508  if( num_fields == 24 ) tot_mem = ( (double)vm_size );
509  }
510 }
511 #endif
512 
513 ErrorCode min_edge_length( Interface& moab, double& result )
514 {
515  double sqr_result = std::numeric_limits< double >::max();
516 
517  ErrorCode rval;
518  Range entities;
519  rval = moab.get_entities_by_handle( 0, entities );
520  if( MB_SUCCESS != rval ) return rval;
521  Range::iterator i = entities.upper_bound( MBVERTEX );
522  entities.erase( entities.begin(), i );
523  i = entities.lower_bound( MBENTITYSET );
524  entities.erase( i, entities.end() );
525 
526  std::vector< EntityHandle > storage;
527  for( i = entities.begin(); i != entities.end(); ++i )
528  {
529  EntityType t = moab.type_from_handle( *i );
530  const EntityHandle* conn;
531  int conn_len, indices[2];
532  rval = moab.get_connectivity( *i, conn, conn_len, true, &storage );
533  if( MB_SUCCESS != rval ) return rval;
534 
535  int num_edges = CN::NumSubEntities( t, 1 );
536  for( int j = 0; j < num_edges; ++j )
537  {
538  CN::SubEntityVertexIndices( t, 1, j, indices );
539  EntityHandle v[2] = { conn[indices[0]], conn[indices[1]] };
540  if( v[0] == v[1] ) continue;
541 
542  double c[6];
543  rval = moab.get_coords( v, 2, c );
544  if( MB_SUCCESS != rval ) return rval;
545 
546  c[0] -= c[3];
547  c[1] -= c[4];
548  c[2] -= c[5];
549  double len_sqr = c[0] * c[0] + c[1] * c[1] + c[2] * c[2];
550  if( len_sqr < sqr_result ) sqr_result = len_sqr;
551  }
552  }
553 
554  result = sqrt( sqr_result );
555  return MB_SUCCESS;
556 }
557 
558 ErrorCode merge_duplicate_vertices( Interface& moab, const double epsilon )
559 {
560  ErrorCode rval;
561  Range verts;
562  rval = moab.get_entities_by_type( 0, MBVERTEX, verts );
563  if( MB_SUCCESS != rval ) return rval;
564 
565  AdaptiveKDTree tree( &moab );
566  EntityHandle root;
567  rval = tree.build_tree( verts, &root );
568  if( MB_SUCCESS != rval )
569  {
570  fprintf( stderr, "Failed to build kD-tree.\n" );
571  return rval;
572  }
573 
574  std::set< EntityHandle > dead_verts;
575  std::vector< EntityHandle > leaves;
576  for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
577  {
578  double coords[3];
579  rval = moab.get_coords( &*i, 1, coords );
580  if( MB_SUCCESS != rval ) return rval;
581 
582  leaves.clear();
583  ;
584  rval = tree.distance_search( coords, epsilon, leaves, epsilon, epsilon );
585  if( MB_SUCCESS != rval ) return rval;
586 
587  Range near;
588  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
589  {
590  Range tmp;
591  rval = moab.get_entities_by_type( *j, MBVERTEX, tmp );
592  if( MB_SUCCESS != rval ) return rval;
593  near.merge( tmp.begin(), tmp.end() );
594  }
595 
596  Range::iterator v = near.find( *i );
597  assert( v != near.end() );
598  near.erase( v );
599 
600  EntityHandle merge = 0;
601  for( Range::iterator j = near.begin(); j != near.end(); ++j )
602  {
603  if( *j < *i && dead_verts.find( *j ) != dead_verts.end() ) continue;
604 
605  double coords2[3];
606  rval = moab.get_coords( &*j, 1, coords2 );
607  if( MB_SUCCESS != rval ) return rval;
608 
609  coords2[0] -= coords[0];
610  coords2[1] -= coords[1];
611  coords2[2] -= coords[2];
612  double dsqr = coords2[0] * coords2[0] + coords2[1] * coords2[1] + coords2[2] * coords2[2];
613  if( dsqr <= epsilon * epsilon )
614  {
615  merge = *j;
616  break;
617  }
618  }
619 
620  if( merge )
621  {
622  dead_verts.insert( *i );
623  rval = moab.merge_entities( merge, *i, false, true );
624  if( MB_SUCCESS != rval ) return rval;
625  }
626  }
627 
628  if( dead_verts.empty() )
629  std::cout << "No duplicate/coincident vertices." << std::endl;
630  else
631  std::cout << "Merged and deleted " << dead_verts.size() << " vertices "
632  << "coincident within " << epsilon << std::endl;
633 
634  return MB_SUCCESS;
635 }