MOAB: Mesh Oriented datABase  (version 5.2.1)
datacoupler_test.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/CpuTimer.hpp"
00003 #include "DataCoupler.hpp"
00004 #include "ElemUtil.hpp"
00005 #include "TestUtil.hpp"
00006 
00007 #include <iostream>
00008 #include <iomanip>
00009 #include <sstream>
00010 #include <assert.h>
00011 
00012 #ifdef MOAB_HAVE_MPI
00013 #include "moab/ParallelComm.hpp"
00014 #include "MBParallelConventions.h"
00015 #endif
00016 
00017 using namespace moab;
00018 
00019 #define PRINT_LAST_ERROR                               \
00020     if( MB_SUCCESS != result )                         \
00021     {                                                  \
00022         std::string tmp_str;                           \
00023         std::cout << "Failure; message:" << std::endl; \
00024         mbImpl->get_last_error( tmp_str );             \
00025         std::cout << tmp_str << std::endl;             \
00026         MPI_Abort( MPI_COMM_WORLD, result );           \
00027         return result;                                 \
00028     }
00029 
00030 // Print usage
00031 void print_usage( char** argv )
00032 {
00033     std::cerr << "Usage: ";
00034     std::cerr << argv[0]
00035               << " -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] "
00036                  "[-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] [-outfile "
00037                  "<out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]"
00038               << std::endl;
00039     std::cerr << "    -meshes" << std::endl;
00040     std::cerr << "        Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
00041     std::cerr << "    -itag" << std::endl;
00042     std::cerr << "        Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
00043     std::cerr << "    -gnorm" << std::endl;
00044     std::cerr << "        Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
00045     std::cerr << "        tag \"<gnorm_tag>_normf\" on the mesh set.  Do this for all meshes." << std::endl;
00046     std::cerr << "    -ssnorm" << std::endl;
00047     std::cerr << "        Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
00048     std::cerr << "        tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset.  Subsets "
00049                  "are selected"
00050               << std::endl;
00051     std::cerr << "        using criteria in <ssnorm_selection>.  Do this for all meshes." << std::endl;
00052     std::cerr << "    -ropts" << std::endl;
00053     std::cerr << "        Read in the mesh files using options in <roptions>." << std::endl;
00054     std::cerr << "    -outfile" << std::endl;
00055     std::cerr << "        Write out target mesh to <out_file>." << std::endl;
00056     std::cerr << "    -wopts" << std::endl;
00057     std::cerr << "        Write out mesh files using options in <woptions>." << std::endl;
00058     std::cerr << "    -dbgout" << std::endl;
00059     std::cerr << "        Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
00060     std::cerr << "    -eps" << std::endl;
00061     std::cerr << "        epsilon" << std::endl;
00062     std::cerr << "    -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl;
00063 }
00064 
00065 #ifdef MOAB_HAVE_HDF5
00066 
00067 ErrorCode get_file_options( int argc, char** argv, std::vector< std::string >& meshFiles, DataCoupler::Method& method,
00068                             std::string& interpTag, std::string& gNormTag, std::string& ssNormTag,
00069                             std::vector< const char* >& ssTagNames, std::vector< const char* >& ssTagValues,
00070                             std::string& readOpts, std::string& outFile, std::string& writeOpts, std::string& dbgFile,
00071                             bool& help, double& epsilon );
00072 
00073 // ErrorCode get_file_options(int argc, char **argv,
00074 //                           std::vector<const char*> &filenames,
00075 //                           std::string &tag_name,
00076 //                           std::string &out_fname,
00077 //                           std::string &opts);
00078 
00079 #ifdef MOAB_HAVE_MPI
00080 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, bool print_results );
00081 #endif
00082 
00083 ErrorCode test_interpolation( Interface* mbImpl, DataCoupler::Method method, std::string& interpTag,
00084                               std::string& gNormTag, std::string& ssNormTag, std::vector< const char* >& ssTagNames,
00085                               std::vector< const char* >& ssTagValues, std::vector< ParallelComm* >& pcs,
00086                               double& instant_time, double& pointloc_time, double& interp_time, double& gnorm_time,
00087                               double& ssnorm_time, double& toler );
00088 
00089 void reduceMax( double& v )
00090 {
00091     double buf;
00092 
00093     MPI_Barrier( MPI_COMM_WORLD );
00094     MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
00095 
00096     v = buf;
00097 }
00098 
00099 int main( int argc, char** argv )
00100 {
00101 #ifdef MOAB_HAVE_MPI
00102     // Need to init MPI first, to tell how many procs and rank
00103     int err = MPI_Init( &argc, &argv );
00104     if( err != 0 )
00105     {
00106         std::cout << "MPI Initialization did not succeed.\n";
00107         exit( 1 );
00108     }
00109 #endif
00110 
00111     std::vector< const char* > ssTagNames, ssTagValues;
00112     std::vector< std::string > meshFiles;
00113     std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
00114     DataCoupler::Method method = DataCoupler::CONSTANT;
00115 
00116     ErrorCode result = MB_SUCCESS;
00117     bool help        = false;
00118     double toler     = 5.e-10;
00119     result = get_file_options( argc, argv, meshFiles, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues,
00120                                readOpts, outFile, writeOpts, dbgFile, help, toler );
00121 
00122     if( result != MB_SUCCESS || help )
00123     {
00124         print_usage( argv );
00125 #ifdef MOAB_HAVE_MPI
00126         err = MPI_Finalize();
00127         if( err != 0 )
00128         {
00129             std::cout << "MPI Initialization did not succeed.\n";
00130             exit( 1 );
00131         }
00132 #endif
00133         return 1;
00134     }
00135 
00136     int nprocs = 1, rank = 0;
00137 
00138 #ifdef MOAB_HAVE_MPI
00139     err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00140     err = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00141     std::vector< ParallelComm* > pcs( meshFiles.size() );
00142 #endif
00143 
00144     // Redirect stdout and stderr if dbgFile is not null
00145     if( !dbgFile.empty() )
00146     {
00147         std::stringstream dfname;
00148         dfname << dbgFile << rank << ".txt";
00149         if( !std::freopen( dfname.str().c_str(), "a", stdout ) ) return 2;
00150         if( !std::freopen( dfname.str().c_str(), "a", stderr ) ) return 2;
00151     }
00152 
00153     // Create MOAB instance based on that
00154     Interface* mbImpl = new( std::nothrow ) Core();
00155     if( NULL == mbImpl ) return 1;
00156 
00157     // Read in mesh(es)
00158 
00159     // Create root sets for each mesh using moab
00160     std::vector< EntityHandle > roots( meshFiles.size() );
00161 
00162     for( unsigned int i = 0; i < meshFiles.size(); i++ )
00163     {
00164         std::string newReadopts;
00165         std::ostringstream extraOpt;
00166 #ifdef MOAB_HAVE_MPI
00167         pcs[i]    = new ParallelComm( mbImpl, MPI_COMM_WORLD );
00168         int index = pcs[i]->get_id();
00169         extraOpt << ";PARALLEL_COMM=" << index;
00170         newReadopts = readOpts + extraOpt.str();
00171 #endif
00172 
00173         result = mbImpl->create_meshset( MESHSET_SET, roots[i] );
00174         PRINT_LAST_ERROR;
00175         result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );
00176         PRINT_LAST_ERROR;
00177     }
00178 
00179 #ifdef MOAB_HAVE_MPI
00180     result = report_iface_ents( mbImpl, pcs, true );
00181     PRINT_LAST_ERROR;
00182 #endif
00183 
00184     double instant_time = 0.0, pointloc_time = 0.0, interp_time = 0.0, gnorm_time = 0.0, ssnorm_time = 0.0;
00185     // Test interpolation and global normalization and subset normalization
00186 
00187     result = test_interpolation( mbImpl, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues, pcs,
00188                                  instant_time, pointloc_time, interp_time, gnorm_time, ssnorm_time, toler );
00189     PRINT_LAST_ERROR;
00190 
00191     reduceMax( instant_time );
00192     reduceMax( pointloc_time );
00193     reduceMax( interp_time );
00194 
00195     if( 0 == rank )
00196         printf( "\nMax time : %g %g %g (inst loc interp -- %d procs)\n", instant_time, pointloc_time, interp_time,
00197                 nprocs );
00198 
00199     // Output mesh
00200     if( !outFile.empty() )
00201     {
00202         Range partSets;
00203         // Only save the target mesh
00204         partSets.insert( (EntityHandle)roots[1] );
00205         std::string newwriteOpts = writeOpts;
00206         std::ostringstream extraOpt;
00207 #ifdef MOAB_HAVE_MPI
00208         extraOpt << ";PARALLEL_COMM=" << 1;
00209         newwriteOpts += extraOpt.str();
00210 #endif
00211         result = mbImpl->write_file( outFile.c_str(), NULL, newwriteOpts.c_str(), partSets );
00212         PRINT_LAST_ERROR;
00213         std::cout << "Wrote " << outFile << std::endl;
00214         std::cout << "mbcoupler_test complete." << std::endl;
00215     }
00216 
00217 #ifdef MOAB_HAVE_MPI
00218     for( unsigned int i = 0; i < meshFiles.size(); i++ )
00219         delete pcs[i];
00220 #endif
00221 
00222     delete mbImpl;
00223     // May be leaking iMeshInst, don't care since it's end of program. Remove above deletes?
00224 
00225 #ifdef MOAB_HAVE_MPI
00226     err = MPI_Finalize();
00227 #endif
00228 
00229     return 0;
00230 }
00231 
00232 #ifdef MOAB_HAVE_MPI
00233 ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, const bool print_results )
00234 {
00235     Range iface_ents[6];
00236     ErrorCode result = MB_SUCCESS, tmp_result;
00237 
00238     // Now figure out which vertices are shared
00239     for( unsigned int p = 0; p < pcs.size(); p++ )
00240     {
00241         for( int i = 0; i < 4; i++ )
00242         {
00243             tmp_result = pcs[p]->get_iface_entities( -1, i, iface_ents[i] );
00244 
00245             if( MB_SUCCESS != tmp_result )
00246             {
00247                 std::cerr << "get_iface_entities returned error on proc " << pcs[p]->proc_config().proc_rank()
00248                           << "; message: " << std::endl;
00249                 std::string last_error;
00250                 result = mbImpl->get_last_error( last_error );
00251                 if( last_error.empty() )
00252                     std::cerr << "(none)" << std::endl;
00253                 else
00254                     std::cerr << last_error << std::endl;
00255                 result = tmp_result;
00256             }
00257             if( 0 != i ) iface_ents[4].merge( iface_ents[i] );
00258         }
00259     }
00260 
00261     // Report # iface entities
00262     result = mbImpl->get_adjacencies( iface_ents[4], 0, false, iface_ents[5], Interface::UNION );
00263 
00264     int rank;
00265     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00266 
00267     if( print_results || iface_ents[0].size() != iface_ents[5].size() )
00268     {
00269         std::cerr << "Proc " << rank << " iface entities: " << std::endl;
00270         for( int i = 0; i < 4; i++ )
00271             std::cerr << "    " << iface_ents[i].size() << " " << i << "d iface entities." << std::endl;
00272         std::cerr << "    (" << iface_ents[5].size() << " verts adj to other iface ents)" << std::endl;
00273     }
00274 
00275     return result;
00276 }
00277 #endif
00278 
00279 // Check first character for a '-'.
00280 // Return true if one is found. False otherwise.
00281 bool check_for_flag( const char* str )
00282 {
00283     if( '-' == str[0] )
00284         return true;
00285     else
00286         return false;
00287 }
00288 
00289 // New get_file_options() function with added possibilities for mbcoupler_test.
00290 ErrorCode get_file_options( int argc, char** argv, std::vector< std::string >& meshFiles, DataCoupler::Method& method,
00291                             std::string& interpTag, std::string& gNormTag, std::string& ssNormTag,
00292                             std::vector< const char* >& ssTagNames, std::vector< const char* >& ssTagValues,
00293                             std::string& readOpts, std::string& outFile, std::string& writeOpts, std::string& dbgFile,
00294                             bool& help, double& epsilon )
00295 {
00296     // Initialize some of the outputs to null values indicating not present
00297     // in the argument list.
00298     gNormTag  = "";
00299     ssNormTag = "";
00300     readOpts  = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_"
00301                "RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME";
00302     outFile                    = "";
00303     writeOpts                  = "PARALLEL=WRITE_PART;CPUTIME";
00304     dbgFile                    = "";
00305     std::string defaultDbgFile = argv[0];  // The executable name will be the default debug output file.
00306 
00307     // These will indicate if we've gotten our required parameters at the end of parsing.
00308     bool haveMeshes    = false;
00309     bool haveInterpTag = false;
00310 
00311     // Loop over the values in argv pulling out an parsing each one
00312     int npos = 1;
00313 
00314     if( argc > 1 && argv[1] == std::string( "-h" ) )
00315     {
00316         help = true;
00317         return MB_SUCCESS;
00318     }
00319 
00320     while( npos < argc )
00321     {
00322         if( argv[npos] == std::string( "-meshes" ) )
00323         {
00324             // Parse out the mesh filenames
00325             npos++;
00326             int numFiles = 2;
00327             meshFiles.resize( numFiles );
00328             for( int i = 0; i < numFiles; i++ )
00329             {
00330                 if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00331                     meshFiles[i] = argv[npos++];
00332                 else
00333                 {
00334                     std::cerr << "    ERROR - missing correct number of mesh filenames" << std::endl;
00335                     return MB_FAILURE;
00336                 }
00337             }
00338 
00339             haveMeshes = true;
00340         }
00341         else if( argv[npos] == std::string( "-itag" ) )
00342         {
00343             // Parse out the interpolation tag
00344             npos++;
00345             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00346                 interpTag = argv[npos++];
00347             else
00348             {
00349                 std::cerr << "    ERROR - missing <interp_tag>" << std::endl;
00350                 return MB_FAILURE;
00351             }
00352 
00353             haveInterpTag = true;
00354         }
00355         else if( argv[npos] == std::string( "-meth" ) )
00356         {
00357             // Parse out the interpolation tag
00358             npos++;
00359             if( argv[npos][0] == '0' )
00360                 method = DataCoupler::CONSTANT;
00361             else if( argv[npos][0] == '1' )
00362                 method = DataCoupler::LINEAR_FE;
00363             else if( argv[npos][0] == '2' )
00364                 method = DataCoupler::QUADRATIC_FE;
00365             else if( argv[npos][0] == '3' )
00366                 method = DataCoupler::SPECTRAL;
00367             else
00368             {
00369                 std::cerr << "    ERROR - unrecognized method number " << method << std::endl;
00370                 return MB_FAILURE;
00371             }
00372             npos++;
00373         }
00374         else if( argv[npos] == std::string( "-eps" ) )
00375         {
00376             // Parse out the tolerance
00377             npos++;
00378             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00379                 epsilon = atof( argv[npos++] );
00380             else
00381             {
00382                 std::cerr << "    ERROR - missing <epsilon>" << std::endl;
00383                 return MB_FAILURE;
00384             }
00385         }
00386         else if( argv[npos] == std::string( "-gnorm" ) )
00387         {
00388             // Parse out the global normalization tag
00389             npos++;
00390             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00391                 gNormTag = argv[npos++];
00392             else
00393             {
00394                 std::cerr << "    ERROR - missing <gnorm_tag>" << std::endl;
00395                 return MB_FAILURE;
00396             }
00397         }
00398         else if( argv[npos] == std::string( "-ssnorm" ) )
00399         {
00400             // Parse out the subset normalization tag and selection criteria
00401             npos++;
00402             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00403                 ssNormTag = argv[npos++];
00404             else
00405             {
00406                 std::cerr << "    ERROR - missing <ssnorm_tag>" << std::endl;
00407                 return MB_FAILURE;
00408             }
00409 
00410             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00411             {
00412                 char* opts         = argv[npos++];
00413                 char sep1[1]       = { ';' };
00414                 char sep2[1]       = { '=' };
00415                 bool end_vals_seen = false;
00416                 std::vector< char* > tmpTagOpts;
00417 
00418                 // First get the options
00419                 for( char* i = strtok( opts, sep1 ); i; i = strtok( 0, sep1 ) )
00420                     tmpTagOpts.push_back( i );
00421 
00422                 // Parse out the name and val or just name.
00423                 for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
00424                 {
00425                     char* e = strtok( tmpTagOpts[j], sep2 );
00426                     ssTagNames.push_back( e );
00427                     e = strtok( 0, sep2 );
00428                     if( e != NULL )
00429                     {
00430                         // We have a value
00431                         if( end_vals_seen )
00432                         {
00433                             // ERROR we should not have a value after none are seen
00434                             std::cerr << "    ERROR - new value seen after end of values in "
00435                                          "<ssnorm_selection>"
00436                                       << std::endl;
00437                             return MB_FAILURE;
00438                         }
00439                         // Otherwise get the value string from e and convert it to an int
00440                         int* valp = new int;
00441                         *valp     = atoi( e );
00442                         ssTagValues.push_back( (const char*)valp );
00443                     }
00444                     else
00445                     {
00446                         // Otherwise there is no '=' so push a null on the list
00447                         end_vals_seen = true;
00448                         ssTagValues.push_back( (const char*)0 );
00449                     }
00450                 }
00451             }
00452             else
00453             {
00454                 std::cerr << "    ERROR - missing <ssnorm_selection>" << std::endl;
00455                 return MB_FAILURE;
00456             }
00457         }
00458         else if( argv[npos] == std::string( "-ropts" ) )
00459         {
00460             // Parse out the mesh file read options
00461             npos++;
00462             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00463                 readOpts = argv[npos++];
00464             else
00465             {
00466                 std::cerr << "    ERROR - missing <roptions>" << std::endl;
00467                 return MB_FAILURE;
00468             }
00469         }
00470         else if( argv[npos] == std::string( "-outfile" ) )
00471         {
00472             // Parse out the output file name
00473             npos++;
00474             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00475                 outFile = argv[npos++];
00476             else
00477             {
00478                 std::cerr << "    ERROR - missing <out_file>" << std::endl;
00479                 return MB_FAILURE;
00480             }
00481         }
00482         else if( argv[npos] == std::string( "-wopts" ) )
00483         {
00484             // Parse out the output file write options
00485             npos++;
00486             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00487                 writeOpts = argv[npos++];
00488             else
00489             {
00490                 std::cerr << "    ERROR - missing <woptions>" << std::endl;
00491                 return MB_FAILURE;
00492             }
00493         }
00494         else if( argv[npos] == std::string( "-dbgout" ) )
00495         {
00496             // Parse out the debug output file name.
00497             // If no name then use the default.
00498             npos++;
00499             if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
00500                 dbgFile = argv[npos++];
00501             else
00502                 dbgFile = defaultDbgFile;
00503         }
00504         else
00505         {
00506             // Unrecognized parameter.  Skip it and move along.
00507             std::cerr << "    ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
00508             std::cerr << "            Skipping..." << std::endl;
00509             npos++;
00510         }
00511     }
00512 
00513     if( !haveMeshes )
00514     {
00515         meshFiles.resize( 2 );
00516         meshFiles[0] = std::string( TestDir + "/64bricks_1khex.h5m" );
00517         meshFiles[1] = std::string( TestDir + "/64bricks_12ktet.h5m" );
00518         std::cout << "Mesh files not entered; using default files " << meshFiles[0] << " and " << meshFiles[1]
00519                   << std::endl;
00520     }
00521 
00522     if( !haveInterpTag )
00523     {
00524         interpTag = "vertex_field";
00525         std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
00526     }
00527 
00528 #ifdef MOAB_HAVE_HDF5
00529     if( 1 == argc )
00530     {
00531         std::cout << "No arguments given; using output file dum.h5m." << std::endl;
00532         outFile = "dum.h5m";
00533     }
00534 #endif
00535 
00536     return MB_SUCCESS;
00537 }
00538 
00539 // End new get_file_options()
00540 
00541 ErrorCode test_interpolation( Interface* mbImpl, DataCoupler::Method method, std::string& interpTag,
00542                               std::string& /* gNormTag */, std::string& /* ssNormTag */,
00543                               std::vector< const char* >& /* ssTagNames */,
00544                               std::vector< const char* >& /* ssTagValues */, std::vector< ParallelComm* >& pcs,
00545                               double& instant_time, double& pointloc_time, double& interp_time,
00546                               double& /* gnorm_time */, double& /* ssnorm_time */, double& toler )
00547 {
00548     assert( method >= DataCoupler::CONSTANT && method <= DataCoupler::SPECTRAL );
00549 
00550     // Source is 1st mesh, target is 2nd
00551     Range src_elems, targ_elems, targ_verts;
00552     ErrorCode result = pcs[0]->get_part_entities( src_elems, 3 );
00553     PRINT_LAST_ERROR;
00554 
00555     CpuTimer timer;
00556 
00557     // Instantiate a coupler, which also initializes the tree
00558     DataCoupler dc( mbImpl, src_elems, 0, pcs[0] );
00559 
00560     // Initialize spectral elements, if they exist
00561     // bool specSou = false, specTar = false;
00562     // result = mbc.initialize_spectral_elements((EntityHandle)roots[0], (EntityHandle)roots[1],
00563     // specSou, specTar);
00564 
00565     instant_time = timer.time_since_birth();
00566 
00567     // Get points from the target mesh to interpolate
00568     // We have to treat differently the case when the target is a spectral mesh
00569     // In that case, the points of interest are the GL points, not the vertex nodes
00570     std::vector< double > vpos;  // This will have the positions we are interested in
00571     int numPointsOfInterest = 0;
00572 #ifdef MOAB_HAVE_MPI
00573     result = pcs[1]->get_part_entities( targ_elems, 3 );
00574 #endif
00575     PRINT_LAST_ERROR;
00576 
00577     // First get all vertices adj to partition entities in target mesh
00578     if( DataCoupler::CONSTANT == method )
00579         targ_verts = targ_elems;
00580     else
00581         result = mbImpl->get_adjacencies( targ_elems, 0, false, targ_verts, Interface::UNION );
00582     PRINT_LAST_ERROR;
00583 
00584 #ifdef MOAB_HAVE_MPI
00585     // Then get non-owned verts and subtract
00586     Range tmp_verts;
00587     result = pcs[1]->get_pstatus_entities( 0, PSTATUS_NOT_OWNED, tmp_verts );
00588     PRINT_LAST_ERROR;
00589     targ_verts = subtract( targ_verts, tmp_verts );
00590 #endif
00591     // Get position of these entities; these are the target points
00592     numPointsOfInterest = (int)targ_verts.size();
00593     vpos.resize( 3 * targ_verts.size() );
00594     result = mbImpl->get_coords( targ_verts, &vpos[0] );
00595     PRINT_LAST_ERROR;
00596 
00597     // Locate those points in the source mesh
00598 #ifdef MOAB_HAVE_MPI
00599     std::cout << "rank " << pcs[0]->proc_config().proc_rank();
00600 #endif
00601     std::cout << " points of interest: " << numPointsOfInterest << "\n";
00602     result = dc.locate_points( &vpos[0], numPointsOfInterest, toler );
00603     PRINT_LAST_ERROR;
00604 
00605     pointloc_time = timer.time_elapsed();
00606 
00607     // Now interpolate tag onto target points
00608     std::vector< double > field( numPointsOfInterest );
00609 
00610     result = dc.interpolate( method, interpTag, &field[0] );
00611     PRINT_LAST_ERROR;
00612 
00613     interp_time = timer.time_elapsed();
00614 
00615     // Set field values as tag on target vertices
00616     // Use original tag
00617     Tag tag;
00618     result = mbImpl->tag_get_handle( interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag );
00619     PRINT_LAST_ERROR;
00620     result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );
00621     PRINT_LAST_ERROR;
00622 
00623     // Done
00624     return MB_SUCCESS;
00625 }
00626 
00627 #else
00628 
00629 int main( int /*argc*/, char** argv )
00630 {
00631     print_usage( argv );
00632     return 0;
00633 }
00634 
00635 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines