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