Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
umr.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 #include <string>
4 #include <iomanip>
5 #include <fstream>
6 #if defined( __MINGW32__ )
7 #include <sys/time.h>
8 #else
9 #include <ctime>
10 #endif
11 
12 #include <ctime>
13 #include <cmath>
14 #include <cassert>
15 #include <cfloat>
16 
17 #include "moab/Core.hpp"
18 #include "moab/Interface.hpp"
20 #include "moab/NestedRefine.hpp"
21 
22 #ifdef MOAB_HAVE_MPI
23 #include "moab_mpi.h"
24 #include "moab/ParallelComm.hpp"
25 #include "MBParallelConventions.h"
26 #endif
27 
28 /* Exit values */
29 #define SUCCESS 0
30 #define USAGE_ERROR 1
31 #define NOT_IMPLEMENTED 2
32 
33 using namespace moab;
34 
35 static void print_usage( const char* name, std::ostream& stream )
36 {
37  stream << "Usage: " << name << " <options> <input_file> [<input_file2> ...]" << std::endl
38  << "Options: " << std::endl
39  << "\t-h - Print this help text and exit." << std::endl
40  << "\t-d - Dimension of the mesh." << std::endl
41  << "\t-n - Exact or a maximum number of levels for the hierarchy. Default 1." << std::endl
42  << "\t-L - Degree of refinement for each level. Pass an array or a number. It "
43  "is mandatory to pass dimension and num_levels before to use this option. If this flag "
44  "is not used, a default of deg 2 refinement is used. "
45  << std::endl
46  << "\t-V - Pass a desired volume (absolute) . This will generate a hierarchy "
47  "such that the maximum volume is reduced to the given amount approximately. The length "
48  "of the hierarchy can be constrained further if a maximum number of levels is passed. "
49  "It is mandatory to pass the dimension for this option. "
50  << std::endl
51  << "\t-q - Prints out the maximum volume of the mesh and exits the program. "
52  "This option can be used as a guide to volume constrainted mesh hierarchy later. "
53  << std::endl
54  << "\t-t - Print out the time taken to generate hierarchy." << std::endl
55  << "\t-s - Print out the mesh sizes of each level of the generated hierarchy." << std::endl
56  << "\t-o - Specify true for output files for the mesh levels of the hierarchy." << std::endl
57  //<< "\t-O option - Specify read option." << std::endl
58 #ifdef MOAB_HAVE_MPI
59  << "\t-p[0|1|2] - Read in parallel[0], optionally also doing resolve_shared_ents (1) "
60  "and exchange_ghosts (2)"
61  << std::endl
62 #endif
63  ;
64  exit( USAGE_ERROR );
65 }
66 
67 static void usage_error( const char* name )
68 {
69  print_usage( name, std::cerr );
70 #ifdef MOAB_HAVE_MPI
71  MPI_Finalize();
72 #endif
73  exit( USAGE_ERROR );
74 }
75 
76 bool parse_id_list( const char* string, int dim, int nval, std::vector< int >& results );
77 
78 bool make_opts_string( std::vector< std::string > options, std::string& opts );
79 
81  EntityHandle fileset,
82  int dim,
83  double desired_vol,
84  int& num_levels,
85  std::vector< int >& level_degs );
86 
87 ErrorCode get_max_volume( Core& mb, EntityHandle fileset, int dim, double& vmax );
88 
89 int main( int argc, char* argv[] )
90 {
91  int proc_id = 0, size = 1;
92 #ifdef MOAB_HAVE_MPI
93  MPI_Init( &argc, &argv );
94  MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
95  MPI_Comm_size( MPI_COMM_WORLD, &size );
96 #endif
97 
98  int num_levels = 0, dim = 0;
99  std::vector< int > level_degrees;
100  bool optimize = false;
101  bool do_flag = true;
102  bool print_times = false, print_size = false, output = false;
103  bool parallel = false, resolve_shared = false, exchange_ghosts = false;
104  bool printusage = false, parselevels = true;
105  bool qc_vol = false, only_quality = false;
106  double cvol = 0;
107  std::string infile;
108 
109  int i = 1;
110  while( i < argc )
111  {
112  if( !argv[i][0] && 0 == proc_id )
113  {
114  usage_error( argv[0] );
115 #ifdef MOAB_HAVE_MPI
116  MPI_Finalize();
117 #endif
118  }
119 
120  if( do_flag && argv[i][0] == '-' )
121  {
122  switch( argv[i][1] )
123  {
124  // do flag arguments:
125  case '-':
126  do_flag = false;
127  break;
128  case 't':
129  print_times = true;
130  break;
131  case 's':
132  print_size = true;
133  break;
134  case 'h':
135  case 'H':
136  print_usage( argv[0], std::cerr );
137  printusage = true;
138  break;
139  case 'd':
140  dim = atol( argv[i + 1] );
141  ++i;
142  break;
143  case 'n':
144  num_levels = atol( argv[i + 1] );
145  ++i;
146  break;
147  case 'L':
148  if( dim != 0 && num_levels != 0 )
149  {
150  parselevels = parse_id_list( argv[i + 1], dim, num_levels, level_degrees );
151  ++i;
152  }
153  else
154  {
155  print_usage( argv[0], std::cerr );
156  printusage = true;
157  }
158  break;
159  case 'V':
160  qc_vol = true;
161  cvol = strtod( argv[i + 1], NULL );
162  ++i;
163  break;
164  case 'q':
165  only_quality = true;
166  break;
167  case 'o':
168  output = true;
169  break;
170  case 'O':
171  optimize = true;
172  break;
173 #ifdef MOAB_HAVE_MPI
174  case 'p':
175  parallel = true;
176  resolve_shared = true;
177  if( argv[i][1] == '1' ) exchange_ghosts = true;
178  break;
179 #endif
180  default:
181  ++i;
182  switch( argv[i - 1][1] )
183  {
184  // case 'O': read_opts.push_back(argv[i]); break;
185  default:
186  std::cerr << "Invalid option: " << argv[i] << std::endl;
187  }
188  }
189  ++i;
190  }
191  // do file names
192  else
193  {
194  infile = argv[i];
195  ++i;
196  }
197  }
198 
199  if( infile.empty() && !printusage ) print_usage( argv[0], std::cerr );
200 
201  if( !parselevels ) exit( USAGE_ERROR );
202 
203 #ifdef MOAB_HAVE_MPI
204  parallel = true;
205  resolve_shared = true;
206 #endif
207 
209 
210  Core* moab = new Core;
211  Interface* mb = (Interface*)moab;
212  EntityHandle fileset;
213 
214  // Create a fileset
216 
217  // Set the read options for parallel file loading
218  std::vector< std::string > read_opts, write_opts;
219  std::string read_options, write_options;
220 
221  if( parallel && size > 1 )
222  {
223  read_opts.push_back( "PARALLEL=READ_PART" );
224  read_opts.push_back( "PARTITION=PARALLEL_PARTITION" );
225  if( resolve_shared ) read_opts.push_back( "PARALLEL_RESOLVE_SHARED_ENTS" );
226  if( exchange_ghosts ) read_opts.push_back( "PARALLEL_GHOSTS=3.0.1" );
227 
228  /* if (output)
229  {
230  write_opts.push_back(";;PARALLEL=WRITE_PART");
231  std::cout<<"Here"<<std::endl;
232  }*/
233  }
234 
235  if( !make_opts_string( read_opts, read_options ) || !make_opts_string( write_opts, write_options ) )
236  {
237 #ifdef MOAB_HAVE_MPI
238  MPI_Finalize();
239 #endif
240  return USAGE_ERROR;
241  }
242 
243  // Load file
244  std::cout << "READ OPTS=" << (char*)read_options.c_str() << std::endl;
245  error = mb->load_file( infile.c_str(), &fileset, read_options.c_str() );MB_CHK_ERR( error );
246 
247  // Create the nestedrefine instance
248 
249 #ifdef MOAB_HAVE_MPI
250  ParallelComm* pc = new ParallelComm( mb, MPI_COMM_WORLD );
251  NestedRefine* uref = new NestedRefine( moab, pc, fileset );
252 #else
253  NestedRefine* uref = new NestedRefine( moab );
254 #endif
255 
256  std::vector< EntityHandle > lsets;
257 
258  // std::cout<<"Read input file"<<std::endl;
259 
260  if( only_quality )
261  {
262  double vmax;
263  error = get_max_volume( *moab, fileset, dim, vmax );MB_CHK_ERR( error );
264 #ifdef MOAB_HAVE_MPI
265  int rank = 0;
266  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
267  if( rank == 0 ) std::cout << "Maximum global volume = " << vmax << std::endl;
268 #else
269  std::cout << "Maximum volume = " << vmax << std::endl;
270 #endif
271  exit( SUCCESS );
272  }
273 
274  // If a volume constraint is given, find an optimal degree sequence to reach the desired volume
275  // constraint.
276  if( qc_vol )
277  {
278  error = get_degree_seq( *moab, fileset, dim, cvol, num_levels, level_degrees );MB_CHK_ERR( error );
279 
280  if( dim == 0 ) print_usage( argv[0], std::cerr );
281  }
282 
283  if( num_levels == 0 ) num_levels = 1;
284 
285  int* ldeg = new int[num_levels];
286  // std::cout<<"level_degrees.size = "<<level_degrees.size()<<std::endl;
287  if( level_degrees.empty() )
288  {
289  for( int l = 0; l < num_levels; l++ )
290  ldeg[l] = 2;
291  }
292  else
293  {
294  for( int l = 0; l < num_levels; l++ )
295  ldeg[l] = level_degrees[l];
296  }
297 
298  std::cout << "Starting hierarchy generation" << std::endl;
299 
300  std::cout << "opt = " << optimize << std::endl;
301 
302  error = uref->generate_mesh_hierarchy( num_levels, ldeg, lsets, optimize );MB_CHK_ERR( error );
303 
304  if( print_times )
305  {
306  std::cout << "Finished hierarchy generation in " << uref->timeall.tm_total << " secs" << std::endl;
307  if( parallel )
308  {
309  std::cout << "Time taken for refinement " << uref->timeall.tm_refine << " secs" << std::endl;
310  std::cout << "Time taken for resolving shared interface " << uref->timeall.tm_resolve << " secs"
311  << std::endl;
312  }
313  }
314  else
315  std::cout << "Finished hierarchy generation " << std::endl;
316 
317  if( print_size )
318  {
319  Range all_ents, ents[4];
320  error = mb->get_entities_by_handle( fileset, all_ents );MB_CHK_ERR( error );
321 
322  for( int k = 0; k < 4; k++ )
323  ents[k] = all_ents.subset_by_dimension( k );
324 
325  std::cout << std::endl;
326  if( qc_vol )
327  {
328  double volume;
329  error = get_max_volume( *moab, fileset, dim, volume );MB_CHK_ERR( error );
330  std::cout << "Mesh size for level 0"
331  << " :: nverts = " << ents[0].size() << ", nedges = " << ents[1].size()
332  << ", nfaces = " << ents[2].size() << ", ncells = " << ents[3].size() << " :: Vmax = " << volume
333  << std::endl;
334  }
335  else
336  std::cout << "Mesh size for level 0"
337  << " :: nverts = " << ents[0].size() << ", nedges = " << ents[1].size()
338  << ", nfaces = " << ents[2].size() << ", ncells = " << ents[3].size() << std::endl;
339 
340  for( int l = 0; l < num_levels; l++ )
341  {
342  all_ents.clear();
343  ents[0].clear();
344  ents[1].clear();
345  ents[2].clear();
346  ents[3].clear();
347  error = mb->get_entities_by_handle( lsets[l + 1], all_ents );MB_CHK_ERR( error );
348 
349  for( int k = 0; k < 4; k++ )
350  ents[k] = all_ents.subset_by_dimension( k );
351 
352  // std::cout<<std::endl;
353 
354  if( qc_vol )
355  {
356  double volume;
357  error = get_max_volume( *moab, lsets[l + 1], dim, volume );MB_CHK_ERR( error );
358  std::cout << "Mesh size for level " << l + 1 << " :: nverts = " << ents[0].size()
359  << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size()
360  << ", ncells = " << ents[3].size() << " :: Vmax = " << volume << std::endl;
361  }
362  else
363  std::cout << "Mesh size for level " << l + 1 << " :: nverts = " << ents[0].size()
364  << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size()
365  << ", ncells = " << ents[3].size() << std::endl;
366  }
367  }
368 
369  if( output )
370  {
371  for( int l = 0; l < num_levels; l++ )
372  {
373  std::string::size_type idx1 = infile.find_last_of( "\\/" );
374  std::string::size_type idx2 = infile.find_last_of( "." );
375  std::string file = infile.substr( idx1 + 1, idx2 - idx1 - 1 );
376  std::stringstream out;
377  if( parallel )
378  out << "_ML_" << l + 1 << ".h5m";
379  else
380  out << "_ML_" << l + 1 << ".vtk";
381  file = file + out.str();
382  const char* output_file = file.c_str();
383 #ifdef MOAB_HAVE_MPI
384  error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART", &lsets[l + 1], 1 );MB_CHK_ERR( error );
385 #else
386  error = mb->write_file( output_file, 0, NULL, &lsets[l + 1], 1 );MB_CHK_ERR( error );
387 #endif
388  // const char* output_file = file.c_str();
389  // error = mb->write_file(output_file, 0, write_options.c_str(), &lsets[l+1],
390  // 1);MB_CHK_ERR(error);
391  // mb->list_entity(lsets[l+1]);
392  // mb->write_file(output_file, 0, "PARALLEL=WRITE_PART", &lsets[l+1], 1);
393  }
394  }
395 
396  delete uref;
397 #ifdef MOAB_HAVE_MPI
398  delete pc;
399 #endif
400 
401  delete[] ldeg;
402  delete moab;
403 
404 #ifdef MOAB_HAVE_MPI
405  MPI_Finalize();
406 #endif
407 
408  exit( SUCCESS );
409 }
410 
412  EntityHandle fileset,
413  int dim,
414  double desired_vol,
415  int& num_levels,
416  std::vector< int >& level_degs )
417 {
418  // Find max volume
419  double vmax_global;
420  ErrorCode error = get_max_volume( mb, fileset, dim, vmax_global );MB_CHK_ERR( error );
421 
422  int init_nl = num_levels;
423  num_levels = 0;
424  level_degs.clear();
425 
426  // Find a sequence that reduces the maximum volume to desired.
427  // desired_vol should be a fraction or absolute value ?
428 
429  // double remV = vmax_global*desired_vol/dim;
430  double Vdesired = desired_vol;
431  double try_x;
432  double remV = vmax_global;
433  int degs[3][3] = { { 5, 3, 2 }, { 25, 9, 4 }, { 0, 27, 8 } };
434 
435  if( dim == 1 || dim == 2 )
436  {
437  while( remV - Vdesired >= 0 )
438  {
439  try_x = degs[dim - 1][0];
440  if( ( remV / try_x - Vdesired ) >= 0 )
441  {
442  level_degs.push_back( 5 );
443  num_levels += 1;
444  remV /= try_x;
445  }
446  else
447  {
448  try_x = degs[dim - 1][1];
449  if( ( remV / try_x - Vdesired ) >= 0 )
450  {
451  level_degs.push_back( 3 );
452  num_levels += 1;
453  remV /= try_x;
454  }
455  else
456  {
457  try_x = degs[dim - 1][2];
458  if( ( remV / try_x - Vdesired ) >= 0 )
459  {
460  level_degs.push_back( 2 );
461  num_levels += 1;
462  remV /= try_x;
463  }
464  else
465  break;
466  }
467  }
468  }
469  }
470  else
471  {
472  while( remV - Vdesired >= 0 )
473  {
474  try_x = degs[dim - 1][1];
475  if( ( remV / try_x - Vdesired ) >= 0 )
476  {
477  level_degs.push_back( 3 );
478  num_levels += 1;
479  remV /= try_x;
480  }
481  else
482  {
483  try_x = degs[dim - 1][2];
484  if( ( remV / try_x - Vdesired ) >= 0 )
485  {
486  level_degs.push_back( 2 );
487  num_levels += 1;
488  remV /= try_x;
489  }
490  else
491  break;
492  }
493  }
494  }
495 
496  if( init_nl != 0 && init_nl < num_levels )
497  {
498  for( int i = level_degs.size(); i >= init_nl; i-- )
499  level_degs.pop_back();
500  num_levels = init_nl;
501  }
502 
503  return MB_SUCCESS;
504 }
505 
506 ErrorCode get_max_volume( Core& mb, EntityHandle fileset, int dim, double& vmax )
507 {
509  VerdictWrapper vw( &mb );
510  QualityType q;
511 
512  switch( dim )
513  {
514  case 1:
515  q = MB_LENGTH;
516  break;
517  case 2:
518  q = MB_AREA;
519  break;
520  case 3:
521  q = MB_VOLUME;
522  break;
523  default:
524  return MB_FAILURE;
525  break;
526  }
527 
528  // Get all entities of the highest dimension which is passed as a command line argument.
529  Range allents, owned;
530  error = mb.get_entities_by_handle( fileset, allents );MB_CHK_ERR( error );
531  owned = allents.subset_by_dimension( dim );MB_CHK_ERR( error );
532 
533  // Get all owned entities
534 #ifdef MOAB_HAVE_MPI
535  int size = 1;
536  MPI_Comm_size( MPI_COMM_WORLD, &size );
537  int mpi_err;
538  if( size > 1 )
539  {
540  // filter the entities not owned, so we do not process them more than once
542  Range current = owned;
543  owned.clear();
544  error = pcomm->filter_pstatus( current, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned );
545  if( error != MB_SUCCESS )
546  {
547  MPI_Finalize();
548  return MB_FAILURE;
549  }
550  }
551 #endif
552 
553  double vmax_local = 0;
554  // Find the maximum volume of an entity in the owned mesh
555  for( Range::iterator it = owned.begin(); it != owned.end(); it++ )
556  {
557  double volume;
558  error = vw.quality_measure( *it, q, volume );MB_CHK_ERR( error );
559  if( volume > vmax_local ) vmax_local = volume;
560  }
561 
562  // Get the global maximum
563  double vmax_global = vmax_local;
564 #ifdef MOAB_HAVE_MPI
565  mpi_err = MPI_Reduce( &vmax_local, &vmax_global, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
566  if( mpi_err )
567  {
568  MPI_Finalize();
569  return MB_FAILURE;
570  }
571 #endif
572 
573  vmax = vmax_global;
574 
575  return MB_SUCCESS;
576 }
577 
578 bool parse_id_list( const char* string, int dim, int nval, std::vector< int >& results )
579 {
580  bool okay = true;
581  char* mystr = strdup( string );
582  for( const char* ptr = strtok( mystr, "," ); ptr; ptr = strtok( 0, "," ) )
583  {
584  char* endptr;
585  int val = strtol( ptr, &endptr, 0 );
586 
587  if( dim == 1 || dim == 2 )
588  {
589  if( val != 2 && val != 3 && val != 5 )
590  {
591  std::cerr << "Not a valid degree for the passed dimension" << std::endl;
592  okay = false;
593  break;
594  }
595  }
596  else if( dim == 3 )
597  {
598  if( val != 2 && val != 3 )
599  {
600  std::cerr << "Not a valid degree for the passed dimension" << std::endl;
601  okay = false;
602  break;
603  }
604  }
605 
606  if( endptr == ptr || val <= 0 )
607  {
608  std::cerr << "Not a valid id: " << ptr << std::endl;
609  okay = false;
610  break;
611  }
612 
613  results.push_back( val );
614  }
615 
616  if( (int)results.size() < nval )
617  {
618  for( int i = results.size(); i <= nval - 1; i++ )
619  results.push_back( results[0] );
620  }
621  else if( (int)results.size() > nval )
622  {
623  for( int i = results.size(); i > nval; i-- )
624  results.pop_back();
625  }
626 
627  free( mystr );
628  return okay;
629 }
630 
631 bool make_opts_string( std::vector< std::string > options, std::string& opts )
632 {
633  opts.clear();
634  if( options.empty() ) return true;
635 
636  // choose a separator character
637  std::vector< std::string >::const_iterator i;
638  char separator = '\0';
639  const char* alt_separators = ";+,:\t\n";
640  for( const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr )
641  {
642  bool seen = false;
643  for( i = options.begin(); i != options.end(); ++i )
644  if( i->find( *sep_ptr, 0 ) != std::string::npos )
645  {
646  seen = true;
647  break;
648  }
649  if( !seen )
650  {
651  separator = *sep_ptr;
652  break;
653  }
654  }
655  if( !separator )
656  {
657  std::cerr << "Error: cannot find separator character for options string" << std::endl;
658  return false;
659  }
660  if( separator != ';' )
661  {
662  opts = ";";
663  opts += separator;
664  }
665 
666  // concatenate options
667  i = options.begin();
668  opts += *i;
669  for( ++i; i != options.end(); ++i )
670  {
671  opts += separator;
672  opts += *i;
673  }
674 
675  return true;
676 }