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