Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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; 133  DataCoupler::Method method = DataCoupler::CONSTANT; 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] ); 193  PRINT_LAST_ERROR; 194  result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() ); 195  PRINT_LAST_ERROR; 196  } 197  198 #ifdef MOAB_HAVE_MPI 199  result = report_iface_ents( mbImpl, pcs, true ); 200  PRINT_LAST_ERROR; 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 ); 208  PRINT_LAST_ERROR; 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 ); 231  PRINT_LAST_ERROR; 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 ); 590  PRINT_LAST_ERROR; 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 612  PRINT_LAST_ERROR; 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 ); 619  PRINT_LAST_ERROR; 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 ); 625  PRINT_LAST_ERROR; 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] ); 632  PRINT_LAST_ERROR; 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 ); 640  PRINT_LAST_ERROR; 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] ); 648  PRINT_LAST_ERROR; 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 ); 656  PRINT_LAST_ERROR; 657  result = mbImpl->tag_set_data( tag, targ_verts, &field[0] ); 658  PRINT_LAST_ERROR; 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