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