Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
datacoupler_test.cpp
Go to the documentation of this file.
1 #include "moab/Core.hpp"
2 #include "moab/CpuTimer.hpp"
3 #include "DataCoupler.hpp"
4 #include "ElemUtil.hpp"
5 #include "TestUtil.hpp"
6 
7 #include <iostream>
8 #include <iomanip>
9 #include <sstream>
10 #include <cassert>
11 
12 #ifdef MOAB_HAVE_MPI
13 #include "moab/ParallelComm.hpp"
14 #include "MBParallelConventions.h"
15 #endif
16 
17 using namespace moab;
18 
19 #define PRINT_LAST_ERROR \
20  if( MB_SUCCESS != result ) \
21  { \
22  std::string tmp_str; \
23  std::cout << "Failure; message:" << std::endl; \
24  mbImpl->get_last_error( tmp_str ); \
25  std::cout << tmp_str << std::endl; \
26  MPI_Abort( MPI_COMM_WORLD, result ); \
27  return result; \
28  }
29 
30 // Print usage
31 void print_usage( char** argv )
32 {
33  std::cerr << "Usage: ";
34  std::cerr << argv[0]
35  << " -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] "
36  "[-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] [-outfile "
37  "<out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]"
38  << std::endl;
39  std::cerr << " -meshes" << std::endl;
40  std::cerr << " Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
41  std::cerr << " -itag" << std::endl;
42  std::cerr << " Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
43  std::cerr << " -gnorm" << std::endl;
44  std::cerr << " Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
45  std::cerr << " tag \"<gnorm_tag>_normf\" on the mesh set. Do this for all meshes." << std::endl;
46  std::cerr << " -ssnorm" << std::endl;
47  std::cerr << " Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
48  std::cerr << " tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset. Subsets "
49  "are selected"
50  << std::endl;
51  std::cerr << " using criteria in <ssnorm_selection>. Do this for all meshes." << std::endl;
52  std::cerr << " -ropts" << std::endl;
53  std::cerr << " Read in the mesh files using options in <roptions>." << std::endl;
54  std::cerr << " -outfile" << std::endl;
55  std::cerr << " Write out target mesh to <out_file>." << std::endl;
56  std::cerr << " -wopts" << std::endl;
57  std::cerr << " Write out mesh files using options in <woptions>." << std::endl;
58  std::cerr << " -dbgout" << std::endl;
59  std::cerr << " Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
60  std::cerr << " -eps" << std::endl;
61  std::cerr << " epsilon" << std::endl;
62  std::cerr << " -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl;
63 }
64 
65 #ifdef MOAB_HAVE_HDF5
66 
67 ErrorCode get_file_options( int argc,
68  char** argv,
69  std::vector< std::string >& meshFiles,
70  DataCoupler::Method& method,
71  std::string& interpTag,
72  std::string& gNormTag,
73  std::string& ssNormTag,
74  std::vector< const char* >& ssTagNames,
75  std::vector< const char* >& ssTagValues,
76  std::string& readOpts,
77  std::string& outFile,
78  std::string& writeOpts,
79  std::string& dbgFile,
80  bool& help,
81  double& epsilon );
82 
83 // ErrorCode get_file_options(int argc, char **argv,
84 // std::vector<const char*> &filenames,
85 // std::string &tag_name,
86 // std::string &out_fname,
87 // std::string &opts);
88 
89 #ifdef MOAB_HAVE_MPI
90 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, bool print_results );
91 #endif
92 
93 ErrorCode test_interpolation( Interface* mbImpl,
94  DataCoupler::Method method,
95  std::string& interpTag,
96  std::string& gNormTag,
97  std::string& ssNormTag,
98  std::vector< const char* >& ssTagNames,
99  std::vector< const char* >& ssTagValues,
100  std::vector< ParallelComm* >& pcs,
101  double& instant_time,
102  double& pointloc_time,
103  double& interp_time,
104  double& gnorm_time,
105  double& ssnorm_time,
106  double& toler );
107 
108 void reduceMax( double& v )
109 {
110  double buf;
111 
112  MPI_Barrier( MPI_COMM_WORLD );
113  MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
114 
115  v = buf;
116 }
117 
118 int main( int argc, char** argv )
119 {
120 #ifdef MOAB_HAVE_MPI
121  // Need to init MPI first, to tell how many procs and rank
122  int err = MPI_Init( &argc, &argv );
123  if( err != 0 )
124  {
125  std::cout << "MPI Initialization did not succeed.\n";
126  exit( 1 );
127  }
128 #endif
129 
130  std::vector< const char* > ssTagNames, ssTagValues;
131  std::vector< std::string > meshFiles;
132  std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
134 
135  ErrorCode result = MB_SUCCESS;
136  bool help = false;
137  double toler = 5.e-10;
138  result = get_file_options( argc, argv, meshFiles, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues,
139  readOpts, outFile, writeOpts, dbgFile, help, toler );
140 
141  if( result != MB_SUCCESS || help )
142  {
143  print_usage( argv );
144 #ifdef MOAB_HAVE_MPI
145  err = MPI_Finalize();
146  if( err != 0 )
147  {
148  std::cout << "MPI Initialization did not succeed.\n";
149  exit( 1 );
150  }
151 #endif
152  return 1;
153  }
154 
155  int nprocs = 1, rank = 0;
156 
157 #ifdef MOAB_HAVE_MPI
158  err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
159  err = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
160  std::vector< ParallelComm* > pcs( meshFiles.size() );
161 #endif
162 
163  // Redirect stdout and stderr if dbgFile is not null
164  if( !dbgFile.empty() )
165  {
166  std::stringstream dfname;
167  dfname << dbgFile << rank << ".txt";
168  if( !std::freopen( dfname.str().c_str(), "a", stdout ) ) return 2;
169  if( !std::freopen( dfname.str().c_str(), "a", stderr ) ) return 2;
170  }
171 
172  // Create MOAB instance based on that
173  Interface* mbImpl = new( std::nothrow ) Core();
174  if( NULL == mbImpl ) return 1;
175 
176  // Read in mesh(es)
177 
178  // Create root sets for each mesh using moab
179  std::vector< EntityHandle > roots( meshFiles.size() );
180 
181  for( unsigned int i = 0; i < meshFiles.size(); i++ )
182  {
183  std::string newReadopts;
184  std::ostringstream extraOpt;
185 #ifdef MOAB_HAVE_MPI
186  pcs[i] = new ParallelComm( mbImpl, MPI_COMM_WORLD );
187  int index = pcs[i]->get_id();
188  extraOpt << ";PARALLEL_COMM=" << index;
189  newReadopts = readOpts + extraOpt.str();
190 #endif
191 
192  result = mbImpl->create_meshset( MESHSET_SET, roots[i] );
194  result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );
196  }
197 
198 #ifdef MOAB_HAVE_MPI
199  result = report_iface_ents( mbImpl, pcs, true );
201 #endif
202 
203  double instant_time = 0.0, pointloc_time = 0.0, interp_time = 0.0, gnorm_time = 0.0, ssnorm_time = 0.0;
204  // Test interpolation and global normalization and subset normalization
205 
206  result = test_interpolation( mbImpl, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues, pcs,
207  instant_time, pointloc_time, interp_time, gnorm_time, ssnorm_time, toler );
209 
210  reduceMax( instant_time );
211  reduceMax( pointloc_time );
212  reduceMax( interp_time );
213 
214  if( 0 == rank )
215  printf( "\nMax time : %g %g %g (inst loc interp -- %d procs)\n", instant_time, pointloc_time, interp_time,
216  nprocs );
217 
218  // Output mesh
219  if( !outFile.empty() )
220  {
221  Range partSets;
222  // Only save the target mesh
223  partSets.insert( (EntityHandle)roots[1] );
224  std::string newwriteOpts = writeOpts;
225  std::ostringstream extraOpt;
226 #ifdef MOAB_HAVE_MPI
227  extraOpt << ";PARALLEL_COMM=" << 1;
228  newwriteOpts += extraOpt.str();
229 #endif
230  result = mbImpl->write_file( outFile.c_str(), NULL, newwriteOpts.c_str(), partSets );
232  std::cout << "Wrote " << outFile << std::endl;
233  std::cout << "mbcoupler_test complete." << std::endl;
234  }
235 
236 #ifdef MOAB_HAVE_MPI
237  for( unsigned int i = 0; i < meshFiles.size(); i++ )
238  delete pcs[i];
239 #endif
240 
241  delete mbImpl;
242  // May be leaking iMeshInst, don't care since it's end of program. Remove above deletes?
243 
244 #ifdef MOAB_HAVE_MPI
245  err = MPI_Finalize();
246 #endif
247 
248  return 0;
249 }
250 
251 #ifdef MOAB_HAVE_MPI
252 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, const bool print_results )
253 {
254  Range iface_ents[6];
255  ErrorCode result = MB_SUCCESS, tmp_result;
256 
257  // Now figure out which vertices are shared
258  for( unsigned int p = 0; p < pcs.size(); p++ )
259  {
260  for( int i = 0; i < 4; i++ )
261  {
262  tmp_result = pcs[p]->get_iface_entities( -1, i, iface_ents[i] );
263 
264  if( MB_SUCCESS != tmp_result )
265  {
266  std::cerr << "get_iface_entities returned error on proc " << pcs[p]->proc_config().proc_rank()
267  << "; message: " << std::endl;
268  std::string last_error;
269  result = mbImpl->get_last_error( last_error );
270  if( last_error.empty() )
271  std::cerr << "(none)" << std::endl;
272  else
273  std::cerr << last_error << std::endl;
274  result = tmp_result;
275  }
276  if( 0 != i ) iface_ents[4].merge( iface_ents[i] );
277  }
278  }
279 
280  // Report # iface entities
281  result = mbImpl->get_adjacencies( iface_ents[4], 0, false, iface_ents[5], Interface::UNION );
282 
283  int rank;
284  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
285 
286  if( print_results || iface_ents[0].size() != iface_ents[5].size() )
287  {
288  std::cerr << "Proc " << rank << " iface entities: " << std::endl;
289  for( int i = 0; i < 4; i++ )
290  std::cerr << " " << iface_ents[i].size() << " " << i << "d iface entities." << std::endl;
291  std::cerr << " (" << iface_ents[5].size() << " verts adj to other iface ents)" << std::endl;
292  }
293 
294  return result;
295 }
296 #endif
297 
298 // Check first character for a '-'.
299 // Return true if one is found. False otherwise.
300 bool check_for_flag( const char* str )
301 {
302  if( '-' == str[0] )
303  return true;
304  else
305  return false;
306 }
307 
308 // New get_file_options() function with added possibilities for mbcoupler_test.
309 ErrorCode get_file_options( int argc,
310  char** argv,
311  std::vector< std::string >& meshFiles,
312  DataCoupler::Method& method,
313  std::string& interpTag,
314  std::string& gNormTag,
315  std::string& ssNormTag,
316  std::vector< const char* >& ssTagNames,
317  std::vector< const char* >& ssTagValues,
318  std::string& readOpts,
319  std::string& outFile,
320  std::string& writeOpts,
321  std::string& dbgFile,
322  bool& help,
323  double& epsilon )
324 {
325  // Initialize some of the outputs to null values indicating not present
326  // in the argument list.
327  gNormTag = "";
328  ssNormTag = "";
329  readOpts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_"
330  "RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME";
331  outFile = "";
332  writeOpts = "PARALLEL=WRITE_PART;CPUTIME";
333  dbgFile = "";
334  std::string defaultDbgFile = argv[0]; // The executable name will be the default debug output file.
335 
336  // These will indicate if we've gotten our required parameters at the end of parsing.
337  bool haveMeshes = false;
338  bool haveInterpTag = false;
339 
340  // Loop over the values in argv pulling out an parsing each one
341  int npos = 1;
342 
343  if( argc > 1 && argv[1] == std::string( "-h" ) )
344  {
345  help = true;
346  return MB_SUCCESS;
347  }
348 
349  while( npos < argc )
350  {
351  if( argv[npos] == std::string( "-meshes" ) )
352  {
353  // Parse out the mesh filenames
354  npos++;
355  int numFiles = 2;
356  meshFiles.resize( numFiles );
357  for( int i = 0; i < numFiles; i++ )
358  {
359  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
360  meshFiles[i] = argv[npos++];
361  else
362  {
363  std::cerr << " ERROR - missing correct number of mesh filenames" << std::endl;
364  return MB_FAILURE;
365  }
366  }
367 
368  haveMeshes = true;
369  }
370  else if( argv[npos] == std::string( "-itag" ) )
371  {
372  // Parse out the interpolation tag
373  npos++;
374  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
375  interpTag = argv[npos++];
376  else
377  {
378  std::cerr << " ERROR - missing <interp_tag>" << std::endl;
379  return MB_FAILURE;
380  }
381 
382  haveInterpTag = true;
383  }
384  else if( argv[npos] == std::string( "-meth" ) )
385  {
386  // Parse out the interpolation tag
387  npos++;
388  if( argv[npos][0] == '0' )
389  method = DataCoupler::CONSTANT;
390  else if( argv[npos][0] == '1' )
391  method = DataCoupler::LINEAR_FE;
392  else if( argv[npos][0] == '2' )
393  method = DataCoupler::QUADRATIC_FE;
394  else if( argv[npos][0] == '3' )
395  method = DataCoupler::SPECTRAL;
396  else
397  {
398  std::cerr << " ERROR - unrecognized method number " << method << std::endl;
399  return MB_FAILURE;
400  }
401  npos++;
402  }
403  else if( argv[npos] == std::string( "-eps" ) )
404  {
405  // Parse out the tolerance
406  npos++;
407  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
408  epsilon = atof( argv[npos++] );
409  else
410  {
411  std::cerr << " ERROR - missing <epsilon>" << std::endl;
412  return MB_FAILURE;
413  }
414  }
415  else if( argv[npos] == std::string( "-gnorm" ) )
416  {
417  // Parse out the global normalization tag
418  npos++;
419  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
420  gNormTag = argv[npos++];
421  else
422  {
423  std::cerr << " ERROR - missing <gnorm_tag>" << std::endl;
424  return MB_FAILURE;
425  }
426  }
427  else if( argv[npos] == std::string( "-ssnorm" ) )
428  {
429  // Parse out the subset normalization tag and selection criteria
430  npos++;
431  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
432  ssNormTag = argv[npos++];
433  else
434  {
435  std::cerr << " ERROR - missing <ssnorm_tag>" << std::endl;
436  return MB_FAILURE;
437  }
438 
439  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
440  {
441  char* opts = argv[npos++];
442  char sep1[1] = { ';' };
443  char sep2[1] = { '=' };
444  bool end_vals_seen = false;
445  std::vector< char* > tmpTagOpts;
446 
447  // First get the options
448  for( char* i = strtok( opts, sep1 ); i; i = strtok( 0, sep1 ) )
449  tmpTagOpts.push_back( i );
450 
451  // Parse out the name and val or just name.
452  for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
453  {
454  char* e = strtok( tmpTagOpts[j], sep2 );
455  ssTagNames.push_back( e );
456  e = strtok( 0, sep2 );
457  if( e != NULL )
458  {
459  // We have a value
460  if( end_vals_seen )
461  {
462  // ERROR we should not have a value after none are seen
463  std::cerr << " ERROR - new value seen after end of values in "
464  "<ssnorm_selection>"
465  << std::endl;
466  return MB_FAILURE;
467  }
468  // Otherwise get the value string from e and convert it to an int
469  int* valp = new int;
470  *valp = atoi( e );
471  ssTagValues.push_back( (const char*)valp );
472  }
473  else
474  {
475  // Otherwise there is no '=' so push a null on the list
476  end_vals_seen = true;
477  ssTagValues.push_back( (const char*)0 );
478  }
479  }
480  }
481  else
482  {
483  std::cerr << " ERROR - missing <ssnorm_selection>" << std::endl;
484  return MB_FAILURE;
485  }
486  }
487  else if( argv[npos] == std::string( "-ropts" ) )
488  {
489  // Parse out the mesh file read options
490  npos++;
491  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
492  readOpts = argv[npos++];
493  else
494  {
495  std::cerr << " ERROR - missing <roptions>" << std::endl;
496  return MB_FAILURE;
497  }
498  }
499  else if( argv[npos] == std::string( "-outfile" ) )
500  {
501  // Parse out the output file name
502  npos++;
503  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
504  outFile = argv[npos++];
505  else
506  {
507  std::cerr << " ERROR - missing <out_file>" << std::endl;
508  return MB_FAILURE;
509  }
510  }
511  else if( argv[npos] == std::string( "-wopts" ) )
512  {
513  // Parse out the output file write options
514  npos++;
515  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
516  writeOpts = argv[npos++];
517  else
518  {
519  std::cerr << " ERROR - missing <woptions>" << std::endl;
520  return MB_FAILURE;
521  }
522  }
523  else if( argv[npos] == std::string( "-dbgout" ) )
524  {
525  // Parse out the debug output file name.
526  // If no name then use the default.
527  npos++;
528  if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
529  dbgFile = argv[npos++];
530  else
531  dbgFile = defaultDbgFile;
532  }
533  else
534  {
535  // Unrecognized parameter. Skip it and move along.
536  std::cerr << " ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
537  std::cerr << " Skipping..." << std::endl;
538  npos++;
539  }
540  }
541 
542  if( !haveMeshes )
543  {
544  meshFiles.resize( 2 );
545  meshFiles[0] = std::string( TestDir + "unittest/64bricks_1khex.h5m" );
546  meshFiles[1] = std::string( TestDir + "unittest/64bricks_12ktet.h5m" );
547  std::cout << "Mesh files not entered; using default files " << meshFiles[0] << " and " << meshFiles[1]
548  << std::endl;
549  }
550 
551  if( !haveInterpTag )
552  {
553  interpTag = "vertex_field";
554  std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
555  }
556 
557 #ifdef MOAB_HAVE_HDF5
558  if( 1 == argc )
559  {
560  std::cout << "No arguments given; using output file dum.h5m." << std::endl;
561  outFile = "dum.h5m";
562  }
563 #endif
564 
565  return MB_SUCCESS;
566 }
567 
568 // End new get_file_options()
569 
570 ErrorCode test_interpolation( Interface* mbImpl,
571  DataCoupler::Method method,
572  std::string& interpTag,
573  std::string& /* gNormTag */,
574  std::string& /* ssNormTag */,
575  std::vector< const char* >& /* ssTagNames */,
576  std::vector< const char* >& /* ssTagValues */,
577  std::vector< ParallelComm* >& pcs,
578  double& instant_time,
579  double& pointloc_time,
580  double& interp_time,
581  double& /* gnorm_time */,
582  double& /* ssnorm_time */,
583  double& toler )
584 {
585  assert( method >= DataCoupler::CONSTANT && method <= DataCoupler::SPECTRAL );
586 
587  // Source is 1st mesh, target is 2nd
588  Range src_elems, targ_elems, targ_verts;
589  ErrorCode result = pcs[0]->get_part_entities( src_elems, 3 );
591 
592  CpuTimer timer;
593 
594  // Instantiate a coupler, which also initializes the tree
595  DataCoupler dc( mbImpl, src_elems, 0, pcs[0] );
596 
597  // Initialize spectral elements, if they exist
598  // bool specSou = false, specTar = false;
599  // result = mbc.initialize_spectral_elements((EntityHandle)roots[0], (EntityHandle)roots[1],
600  // specSou, specTar);
601 
602  instant_time = timer.time_since_birth();
603 
604  // Get points from the target mesh to interpolate
605  // We have to treat differently the case when the target is a spectral mesh
606  // In that case, the points of interest are the GL points, not the vertex nodes
607  std::vector< double > vpos; // This will have the positions we are interested in
608  int numPointsOfInterest = 0;
609 #ifdef MOAB_HAVE_MPI
610  result = pcs[1]->get_part_entities( targ_elems, 3 );
611 #endif
613 
614  // First get all vertices adj to partition entities in target mesh
615  if( DataCoupler::CONSTANT == method )
616  targ_verts = targ_elems;
617  else
618  result = mbImpl->get_adjacencies( targ_elems, 0, false, targ_verts, Interface::UNION );
620 
621 #ifdef MOAB_HAVE_MPI
622  // Then get non-owned verts and subtract
623  Range tmp_verts;
624  result = pcs[1]->get_pstatus_entities( 0, PSTATUS_NOT_OWNED, tmp_verts );
626  targ_verts = subtract( targ_verts, tmp_verts );
627 #endif
628  // Get position of these entities; these are the target points
629  numPointsOfInterest = (int)targ_verts.size();
630  vpos.resize( 3 * targ_verts.size() );
631  result = mbImpl->get_coords( targ_verts, &vpos[0] );
633 
634  // Locate those points in the source mesh
635 #ifdef MOAB_HAVE_MPI
636  std::cout << "rank " << pcs[0]->proc_config().proc_rank();
637 #endif
638  std::cout << " points of interest: " << numPointsOfInterest << "\n";
639  result = dc.locate_points( &vpos[0], numPointsOfInterest, toler );
641 
642  pointloc_time = timer.time_elapsed();
643 
644  // Now interpolate tag onto target points
645  std::vector< double > field( numPointsOfInterest );
646 
647  result = dc.interpolate( method, interpTag, &field[0] );
649 
650  interp_time = timer.time_elapsed();
651 
652  // Set field values as tag on target vertices
653  // Use original tag
654  Tag tag;
655  result = mbImpl->tag_get_handle( interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag );
657  result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );
659 
660  // Done
661  return MB_SUCCESS;
662 }
663 
664 #else
665 
666 int main( int /*argc*/, char** argv )
667 {
668  print_usage( argv );
669  return 0;
670 }
671 
672 #endif