mbcoupler_test.cpp

Go to the documentation of this file.

// MOAB includes
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "moab/Core.hpp"
#include "Coupler.hpp"
#include "moab_mpi.h"
#include "ElemUtil.hpp"
#include "TestUtil.hpp"

// C++ includes
#include 
#include 
#include 
#include 

using namespace moab;

/*
  Sample usages:
    1) P_0 interpolation: ./mbcoupler_test -meshes   -itag  -meth
  0 -outfile  2) P_1 interpolation: ./mbcoupler_test -meshes   -itag
   -meth 1 -outfile  3) P_1 interpolation with epsilon control: ./mbcoupler_test
  -meshes   -itag  -meth 1 -eps  -outfile  3) P_0 interpolation with global normalization: ./mbcoupler_test
  -meshes   -itag  -meth 0 -gnorm  -outfile
   4) P_1 interpolation with subset normalization: ./mbcoupler_test -meshes 
   -itag  -meth 1 -ssnorm  
  -outfile  5) P_1 interpolation for meshes on sphere  ./mbcoupler_test -meshes  
    -itag    -meth 4 -outfile  
*/

// Print usage
void print_usage()
{
    std::cerr << "Usage: ";
    std::cerr << "mbcoupler_test -meshes   -itag  [-gnorm "
                 "] [-ssnorm  ] [-ropts ] "
                 "[-outfile  [-wopts ]] [-dbgout []]"
              << std::endl;
    std::cerr << "    -meshes" << std::endl;
    std::cerr << "        Read in mesh files  and ." << std::endl;
    std::cerr << "    -itag" << std::endl;
    std::cerr << "        Interpolate tag  from source mesh to target mesh." << std::endl;
    std::cerr << "    -gnorm" << std::endl;
    std::cerr << "        Normalize the value of tag  over then entire mesh and save to" << std::endl;
    std::cerr << "        tag \"_normf\" on the mesh set.  Do this for all meshes." << std::endl;
    std::cerr << "    -ssnorm" << std::endl;
    std::cerr << "        Normalize the value of tag  over subsets of a mesh and save to" << std::endl;
    std::cerr << "        tag \"_normf\" on the Entity Set for each subset.  Subsets "
                 "are selected"
              << std::endl;
    std::cerr << "        using criteria in .  Do this for all meshes." << std::endl;
    std::cerr << "    -ropts" << std::endl;
    std::cerr << "        Read in the mesh files using options in ." << std::endl;
    std::cerr << "    -outfile" << std::endl;
    std::cerr << "        Write out target mesh to ." << std::endl;
    std::cerr << "    -wopts" << std::endl;
    std::cerr << "        Write out mesh files using options in ." << std::endl;
    std::cerr << "    -dbgout" << std::endl;
    std::cerr << "        Write stdout and stderr streams to the file \'.txt\'." << std::endl;
    std::cerr << "    -eps" << std::endl;
    std::cerr << "        epsilon" << std::endl;
    std::cerr << "    -meth  (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL, 4=SPHERICAL)" << std::endl;
}

#ifdef MOAB_HAVE_HDF5

ErrorCode get_file_options( int argc,
                            char** argv,
                            int nprocs,
                            int rank,
                            std::vector< std::string >& meshFiles,
                            Coupler::Method& method,
                            std::string& interpTag,
                            std::string& gNormTag,
                            std::string& ssNormTag,
                            std::vector< const char* >& ssTagNames,
                            std::vector< const char* >& ssTagValues,
                            std::string& readOpts,
                            std::string& outFile,
                            std::string& writeOpts,
                            std::string& dbgFile,
                            bool& help,
                            double& epsilon );

ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, bool print_results );

ErrorCode test_interpolation( Interface* mbImpl,
                              Coupler::Method method,
                              std::string& interpTag,
                              std::string& gNormTag,
                              std::string& ssNormTag,
                              std::vector< const char* >& ssTagNames,
                              std::vector< const char* >& ssTagValues,
                              EntityHandle* roots,
                              std::vector< ParallelComm* >& pcs,
                              double& instant_time,
                              double& pointloc_time,
                              double& interp_time,
                              double& gnorm_time,
                              double& ssnorm_time,
                              double& toler );

void reduceMax( double& v )
{
    double buf;

    MPI_Barrier( MPI_COMM_WORLD );
    MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );

    v = buf;
}

int main( int argc, char** argv )
{
    // Need to init MPI first, to tell how many procs and rank
    int ierr = MPI_Init( &argc, &argv );
    assert( MPI_SUCCESS == ierr );

    std::vector< const char* > ssTagNames, ssTagValues;
    std::vector< std::string > meshFiles;
    std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
    Coupler::Method method = Coupler::CONSTANT;

    ErrorCode result;
    bool help    = false;
    double toler = 5.e-10;
    int nprocs, rank;
    ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    assert( MPI_SUCCESS == ierr );
    ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    assert( MPI_SUCCESS == ierr );

    result = get_file_options( argc, argv, nprocs, rank, meshFiles, method, interpTag, gNormTag, ssNormTag, ssTagNames,
                               ssTagValues, readOpts, outFile, writeOpts, dbgFile, help, toler );

    if( result != MB_SUCCESS || help )
    {
        print_usage();
        ierr = MPI_Finalize();
        return ierr;
    }

    // Redirect stdout and stderr if dbgFile is not null
    if( !dbgFile.empty() )
    {
        std::stringstream dfname;
        dfname << dbgFile << rank << ".txt";
        if( !std::freopen( dfname.str().c_str(), "a", stdout ) ) return -1;
        if( !std::freopen( dfname.str().c_str(), "a", stderr ) ) return -1;
    }

    // Create MOAB instance based on that
    Interface* mbImpl = new( std::nothrow ) Core();
    if( NULL == mbImpl ) return 1;

    // Read in mesh(es)
    std::vector< ParallelComm* > pcs( meshFiles.size() );

    // Create root sets for each mesh using the iMesh API. Then pass these
    // to the load_file functions to be populated.
    EntityHandle* roots = (EntityHandle*)malloc( sizeof( EntityHandle ) * meshFiles.size() );

    for( unsigned int i = 0; i < meshFiles.size(); i++ )
    {
        pcs[i]    = new ParallelComm( mbImpl, MPI_COMM_WORLD );
        int index = pcs[i]->get_id();
        std::string newReadopts;
        std::ostringstream extraOpt;
        extraOpt << ";PARALLEL_COMM=" << index;
        newReadopts = readOpts + extraOpt.str();

        result = mbImpl->create_meshset( MESHSET_SET, roots[i] );MB_CHK_ERR( result );

        result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );MB_CHK_ERR( result );
        // result = rps[i]->load_file(meshFiles[i].c_str(), &roots[i],
        // FileOptions(readOpts.c_str()));MB_CHK_ERR(result);
    }

    result = report_iface_ents( mbImpl, pcs, true );MB_CHK_ERR( result );

    double instant_time = 0.0, pointloc_time = 0.0, interp_time = 0.0, gnorm_time = 0.0, ssnorm_time = 0.0;
    // Test interpolation and global normalization and subset normalization

    result = test_interpolation( mbImpl, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues, roots, pcs,
                                 instant_time, pointloc_time, interp_time, gnorm_time, ssnorm_time, toler );MB_CHK_ERR( result );

    reduceMax( instant_time );
    reduceMax( pointloc_time );
    reduceMax( interp_time );

    if( 0 == rank )
        printf( "\nMax time : %g %g %g (inst loc interp -- %d procs)\n", instant_time, pointloc_time, interp_time,
                nprocs );

    // Output mesh
    if( !outFile.empty() )
    {
        Range partSets;
        // Only save the target mesh
        partSets.insert( (EntityHandle)roots[1] );
        std::string newwriteOpts;
        std::ostringstream extraOpt;
        if( nprocs > 1 ) extraOpt << ";PARALLEL_COMM=" << 1;
        newwriteOpts = writeOpts + extraOpt.str();
        result       = mbImpl->write_file( outFile.c_str(), NULL, newwriteOpts.c_str(), partSets );MB_CHK_ERR( result );
        if( 0 == rank )
        {
            std::cout << "Wrote " << outFile << std::endl;
            std::cout << "mbcoupler_test complete." << std::endl;
        }
    }

    free( roots );

    for( unsigned int i = 0; i < meshFiles.size(); i++ )
        delete pcs[i];

    delete mbImpl;

    ierr = MPI_Finalize();
    assert( MPI_SUCCESS == ierr );
    return 0;
}

ErrorCode report_iface_ents( Interface* mbImpl, std::vector< ParallelComm* >& pcs, const bool print_results )
{
    Range iface_ents[6];
    ErrorCode result = MB_SUCCESS, tmp_result;

    // Now figure out which vertices are shared
    for( unsigned int p = 0; p < pcs.size(); p++ )
    {
        for( int i = 0; i < 4; i++ )
        {
            tmp_result = pcs[p]->get_iface_entities( -1, i, iface_ents[i] );

            if( MB_SUCCESS != tmp_result )
            {
                std::cerr << "get_iface_entities returned error on proc " << pcs[p]->proc_config().proc_rank()
                          << "; message: " << std::endl;
                std::string last_error;
                result = mbImpl->get_last_error( last_error );
                if( last_error.empty() )
                    std::cerr << "(none)" << std::endl;
                else
                    std::cerr << last_error << std::endl;
                result = tmp_result;
            }
            if( 0 != i ) iface_ents[4].merge( iface_ents[i] );
        }
    }

    // Report # iface entities
    result = mbImpl->get_adjacencies( iface_ents[4], 0, false, iface_ents[5], Interface::UNION );MB_CHK_ERR( result );

    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    if( print_results || iface_ents[0].size() != iface_ents[5].size() )
    {
        std::cerr << "Proc " << rank << " iface entities: " << std::endl;
        for( int i = 0; i < 4; i++ )
            std::cerr << "    " << iface_ents[i].size() << " " << i << "d iface entities." << std::endl;
        std::cerr << "    (" << iface_ents[5].size() << " verts adj to other iface ents)" << std::endl;
    }

    return result;
}

// Check first character for a '-'.
// Return true if one is found. False otherwise.
bool check_for_flag( const char* str )
{
    if( '-' == str[0] )
        return true;
    else
        return false;
}

// get_file_options() function with added possibilities for mbcoupler_test.
ErrorCode get_file_options( int argc,
                            char** argv,
                            int nprocs,
                            int rank,
                            std::vector< std::string >& meshFiles,
                            Coupler::Method& method,
                            std::string& interpTag,
                            std::string& gNormTag,
                            std::string& ssNormTag,
                            std::vector< const char* >& ssTagNames,
                            std::vector< const char* >& ssTagValues,
                            std::string& readOpts,
                            std::string& outFile,
                            std::string& writeOpts,
                            std::string& dbgFile,
                            bool& help,
                            double& epsilon )
{
    // Initialize some of the outputs to null values indicating not present
    // in the argument list.
    gNormTag                   = "";
    ssNormTag                  = "";
    readOpts                   = ( nprocs > 1 ? "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;"
                                                "PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME"
                                              : "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;"
                                                "PARALLEL_RESOLVE_SHARED_ENTS;CPUTIME" );
    outFile                    = "";
    writeOpts                  = ( nprocs > 1 ? "PARALLEL=WRITE_PART;CPUTIME" : "" );
    dbgFile                    = "";
    std::string defaultDbgFile = argv[0];  // The executable name will be the default debug output file.

    // These will indicate if we've gotten our required parameters at the end of parsing.
    bool haveMeshes    = false;
    bool haveInterpTag = false;

    // Loop over the values in argv pulling out an parsing each one
    int npos = 1;

    if( argc > 1 && argv[1] == std::string( "-h" ) )
    {
        help = true;
        return MB_SUCCESS;
    }

    while( npos < argc )
    {
        if( argv[npos] == std::string( "-meshes" ) )
        {
            // Parse out the mesh filenames
            npos++;
            int numFiles = 2;
            meshFiles.resize( numFiles );
            for( int i = 0; i < numFiles; i++ )
            {
                if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                    meshFiles[i] = argv[npos++];
                else
                {
                    std::cerr << "    ERROR - missing correct number of mesh filenames" << std::endl;
                    return MB_FAILURE;
                }
            }

            haveMeshes = true;
        }
        else if( argv[npos] == std::string( "-itag" ) )
        {
            // Parse out the interpolation tag
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                interpTag = argv[npos++];
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }

            haveInterpTag = true;
        }
        else if( argv[npos] == std::string( "-meth" ) )
        {
            // Parse out the interpolation tag
            npos++;
            if( argv[npos][0] == '0' )
                method = Coupler::CONSTANT;
            else if( argv[npos][0] == '1' )
                method = Coupler::LINEAR_FE;
            else if( argv[npos][0] == '2' )
                method = Coupler::QUADRATIC_FE;
            else if( argv[npos][0] == '3' )
                method = Coupler::SPECTRAL;
            else if( argv[npos][0] == '4' )
                method = Coupler::SPHERICAL;
            else
            {
                std::cerr << "    ERROR - unrecognized method number " << method << std::endl;
                return MB_FAILURE;
            }
            npos++;
        }
        else if( argv[npos] == std::string( "-eps" ) )
        {
            // Parse out the tolerance
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                epsilon = atof( argv[npos++] );
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }
        }
        else if( argv[npos] == std::string( "-gnorm" ) )
        {
            // Parse out the global normalization tag
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                gNormTag = argv[npos++];
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }
        }
        else if( argv[npos] == std::string( "-ssnorm" ) )
        {
            // Parse out the subset normalization tag and selection criteria
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                ssNormTag = argv[npos++];
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }

            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
            {
                char* opts         = argv[npos++];
                char sep1[1]       = { ';' };
                char sep2[1]       = { '=' };
                bool end_vals_seen = false;
                std::vector< char* > tmpTagOpts;

                // First get the options
                for( char* i = strtok( opts, sep1 ); i; i = strtok( 0, sep1 ) )
                    tmpTagOpts.push_back( i );

                // Parse out the name and val or just name.
                for( unsigned int j = 0; j < tmpTagOpts.size(); j++ )
                {
                    char* e = strtok( tmpTagOpts[j], sep2 );
                    ssTagNames.push_back( e );
                    e = strtok( 0, sep2 );
                    if( e != NULL )
                    {
                        // We have a value
                        if( end_vals_seen )
                        {
                            // ERROR we should not have a value after none are seen
                            std::cerr << "    ERROR - new value seen after end of values in "
                                         ""
                                      << std::endl;
                            return MB_FAILURE;
                        }
                        // Otherwise get the value string from e and convert it to an int
                        int* valp = new int;
                        *valp     = atoi( e );
                        ssTagValues.push_back( (const char*)valp );
                    }
                    else
                    {
                        // Otherwise there is no '=' so push a null on the list
                        end_vals_seen = true;
                        ssTagValues.push_back( (const char*)0 );
                    }
                }
            }
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }
        }
        else if( argv[npos] == std::string( "-ropts" ) )
        {
            // Parse out the mesh file read options
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                readOpts = argv[npos++];
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }
        }
        else if( argv[npos] == std::string( "-outfile" ) )
        {
            // Parse out the output file name
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                outFile = argv[npos++];
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }
        }
        else if( argv[npos] == std::string( "-wopts" ) )
        {
            // Parse out the output file write options
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                writeOpts = argv[npos++];
            else
            {
                std::cerr << "    ERROR - missing " << std::endl;
                return MB_FAILURE;
            }
        }
        else if( argv[npos] == std::string( "-dbgout" ) )
        {
            // Parse out the debug output file name.
            // If no name then use the default.
            npos++;
            if( ( npos < argc ) && ( !check_for_flag( argv[npos] ) ) )
                dbgFile = argv[npos++];
            else
                dbgFile = defaultDbgFile;
        }
        else
        {
            // Unrecognized parameter.  Skip it and move along.
            std::cerr << "    ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
            std::cerr << "            Skipping..." << std::endl;
            npos++;
        }
    }

    if( !haveMeshes )
    {
        meshFiles.resize( 2 );
        meshFiles[0] = std::string( TestDir + "unittest/64bricks_1khex.h5m" );
        meshFiles[1] = std::string( TestDir + "unittest/64bricks_12ktet.h5m" );
        if( 0 == rank )
            std::cout << "Mesh files not entered; using default files " << meshFiles[0] << " and " << meshFiles[1]
                      << std::endl;
    }

    if( !haveInterpTag )
    {
        interpTag = "vertex_field";
        std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
    }

#ifdef MOAB_HAVE_HDF5
    if( 1 == argc )
    {
        std::stringstream dfname;
        dfname << "dum" << nprocs << ".h5m";
        outFile = dfname.str();
        if( 0 == rank ) std::cout << "No arguments given; using output file " << outFile << std::endl;
    }
#endif

    return MB_SUCCESS;
}

// End get_file_options()

ErrorCode test_interpolation( Interface* mbImpl,
                              Coupler::Method method,
                              std::string& interpTag,
                              std::string& gNormTag,
                              std::string& ssNormTag,
                              std::vector< const char* >& ssTagNames,
                              std::vector< const char* >& ssTagValues,
                              EntityHandle* roots,
                              std::vector< ParallelComm* >& pcs,
                              double& instant_time,
                              double& pointloc_time,
                              double& interp_time,
                              double& gnorm_time,
                              double& ssnorm_time,
                              double& toler )
{
    assert( method >= Coupler::CONSTANT && method <= Coupler::SPHERICAL );

    // Source is 1st mesh, target is 2nd
    Range src_elems, targ_elems, targ_verts;
    ErrorCode result = pcs[0]->get_part_entities( src_elems, 3 );MB_CHK_ERR( result );

    double start_time = MPI_Wtime();

    // Instantiate a coupler, which does not initialize the tree yet
    Coupler mbc( mbImpl, pcs[0], src_elems, 0, false );  // do not initialize tree yet
    if( Coupler::SPHERICAL == method ) mbc.set_spherical();
    mbc.initialize_tree();  // it is expensive, but do something different for spherical

    // Initialize spectral elements, if they exist
    bool specSou = false, specTar = false;
    result = mbc.initialize_spectral_elements( roots[0], roots[1], specSou, specTar );

    instant_time = MPI_Wtime();

    // Get points from the target mesh to interpolate
    // We have to treat differently the case when the target is a spectral mesh
    // In that case, the points of interest are the GL points, not the vertex nodes
    std::vector< double > vpos;  // This will have the positions we are interested in
    int numPointsOfInterest = 0;
    if( !specTar )
    {  // Usual case
        Range tmp_verts;

        // First get all vertices adj to partition entities in target mesh
        result = pcs[1]->get_part_entities( targ_elems, 3 );MB_CHK_ERR( result );
        if( Coupler::SPHERICAL == method ) result = pcs[1]->get_part_entities( targ_elems, 2 );MB_CHK_ERR( result );  // get the polygons/quads on a sphere.
        if( Coupler::CONSTANT == method )
            targ_verts = targ_elems;
        else
            result = mbImpl->get_adjacencies( targ_elems, 0, false, targ_verts, Interface::UNION );MB_CHK_ERR( result );

        // Then get non-owned verts and subtract
        result = pcs[1]->get_pstatus_entities( 0, PSTATUS_NOT_OWNED, tmp_verts );MB_CHK_ERR( result );
        targ_verts = subtract( targ_verts, tmp_verts );
        // get position of these entities; these are the target points
        numPointsOfInterest = (int)targ_verts.size();
        vpos.resize( 3 * targ_verts.size() );
        result = mbImpl->get_coords( targ_verts, &vpos[0] );MB_CHK_ERR( result );
        // Locate those points in the source mesh
        std::cout << "rank " << pcs[0]->proc_config().proc_rank() << " points of interest: " << numPointsOfInterest
                  << "\n";
        result = mbc.locate_points( &vpos[0], numPointsOfInterest, 0, toler );MB_CHK_ERR( result );
    }
    else
    {
        // In this case, the target mesh is spectral, we want values
        // interpolated on the GL positions; for each element, get the GL points, and construct
        // CartVect!!!
        result = pcs[1]->get_part_entities( targ_elems, 3 );MB_CHK_ERR( result );
        result = mbc.get_gl_points_on_elements( targ_elems, vpos, numPointsOfInterest );MB_CHK_ERR( result );
        std::cout << "rank " << pcs[0]->proc_config().proc_rank() << " points of interest: " << numPointsOfInterest
                  << "\n";
    }

    pointloc_time = MPI_Wtime();

    // Now interpolate tag onto target points
    std::vector< double > field( numPointsOfInterest );

    result = mbc.interpolate( method, interpTag, &field[0] );MB_CHK_ERR( result );

    interp_time = MPI_Wtime();

    // Do global normalization if specified
    if( !gNormTag.empty() )
    {
        // Normalize the source mesh
        result = mbc.normalize_mesh( roots[0], gNormTag.c_str(), Coupler::VOLUME, 4 );MB_CHK_ERR( result );

        // Normalize the target mesh
        result = mbc.normalize_mesh( roots[1], gNormTag.c_str(), Coupler::VOLUME, 4 );MB_CHK_ERR( result );
    }

    gnorm_time = MPI_Wtime();

    // Do subset normalization if specified

    if( !ssNormTag.empty() )
    {

        result = mbc.normalize_subset( roots[0], ssNormTag.c_str(), &ssTagNames[0], ssTagNames.size(), &ssTagValues[0],
                                       Coupler::VOLUME, 4 );MB_CHK_ERR( result );

        result = mbc.normalize_subset( roots[1], ssNormTag.c_str(), &ssTagNames[0], ssTagNames.size(), &ssTagValues[0],
                                       Coupler::VOLUME, 4 );MB_CHK_ERR( result );
    }

    ssnorm_time = MPI_Wtime();

    ssnorm_time -= gnorm_time;
    gnorm_time -= interp_time;
    interp_time -= pointloc_time;
    pointloc_time -= instant_time;
    instant_time -= start_time;

    // Set field values as tag on target vertices
    if( specSou )
    {
        // Create a new tag for the values on the target
        Tag tag;
        std::string newtag = interpTag + "_TAR";
        result = mbImpl->tag_get_handle( newtag.c_str(), 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( result );
        result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );MB_CHK_ERR( result );
    }
    else
    {
        if( !specTar )
        {
            // Use original tag
            Tag tag;
            result = mbImpl->tag_get_handle( interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag );MB_CHK_ERR( result );
            result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );MB_CHK_ERR( result );
        }
        else
        {
            // We have the field values computed at each GL points, on each element
            // in the target mesh
            // We need to create a new tag, on elements, of size _ntot, to hold
            // all those values.
            // So it turns out we need ntot. Maybe we can compute it from the
            // number of values computed, divided by number of elements
            int ntot = numPointsOfInterest / targ_elems.size();
            Tag tag;
            std::string newtag = interpTag + "_TAR";
            result = mbImpl->tag_get_handle( newtag.c_str(), ntot, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( result );
            result = mbImpl->tag_set_data( tag, targ_elems, &field[0] );MB_CHK_ERR( result );
        }
    }

    // Done
    return MB_SUCCESS;
}

#else

int main( int /*argc*/, char** /*argv*/ )
{
    print_usage();
    return 0;
}

#endif