datacoupler_test.cpp

Go to the documentation of this file.

#include "moab/Core.hpp"
#include "moab/CpuTimer.hpp"
#include "DataCoupler.hpp"
#include "ElemUtil.hpp"
#include "TestUtil.hpp"

#include <iostream>
#include <iomanip>
#include <sstream>
#include <cassert>

#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#endif

using namespace moab;

#define PRINT_LAST_ERROR                               \
    if( MB_SUCCESS != result )                         \
    {                                                  \
        std::string tmp_str;                           \
        std::cout << "Failure; message:" << std::endl; \
        mbImpl->get_last_error( tmp_str );             \
        std::cout << tmp_str << std::endl;             \
        MPI_Abort( MPI_COMM_WORLD, result );           \
        return result;                                 \
    }

// Print usage
void print_usage( char** argv )
{
    std::cerr << "Usage: ";
    std::cerr << argv[0]
              << " -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] "
                 "[-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] [-outfile "
                 "<out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]"
              << std::endl;
    std::cerr << "    -meshes" << std::endl;
    std::cerr << "        Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
    std::cerr << "    -itag" << std::endl;
    std::cerr << "        Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
    std::cerr << "    -gnorm" << std::endl;
    std::cerr << "        Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
    std::cerr << "        tag \"<gnorm_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 <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
    std::cerr << "        tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset.  Subsets "
                 "are selected"
              << std::endl;
    std::cerr << "        using criteria in <ssnorm_selection>.  Do this for all meshes." << std::endl;
    std::cerr << "    -ropts" << std::endl;
    std::cerr << "        Read in the mesh files using options in <roptions>." << std::endl;
    std::cerr << "    -outfile" << std::endl;
    std::cerr << "        Write out target mesh to <out_file>." << std::endl;
    std::cerr << "    -wopts" << std::endl;
    std::cerr << "        Write out mesh files using options in <woptions>." << std::endl;
    std::cerr << "    -dbgout" << std::endl;
    std::cerr << "        Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
    std::cerr << "    -eps" << std::endl;
    std::cerr << "        epsilon" << std::endl;
    std::cerr << "    -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl;
}

#ifdef MOAB_HAVE_HDF5

ErrorCode get_file_options( int argc,
                            char** argv,
                            std::vector< std::string >& meshFiles,
                            DataCoupler::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 get_file_options(int argc, char **argv,
//                           std::vector<const char*> &filenames,
//                           std::string &tag_name,
//                           std::string &out_fname,
//                           std::string &opts);

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

ErrorCode test_interpolation( Interface* mbImpl,
                              DataCoupler::Method method,
                              std::string& interpTag,
                              std::string& gNormTag,
                              std::string& ssNormTag,
                              std::vector< const char* >& ssTagNames,
                              std::vector< const char* >& ssTagValues,
                              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 )
{
#ifdef MOAB_HAVE_MPI
    // Need to init MPI first, to tell how many procs and rank
    int err = MPI_Init( &argc, &argv );
    if( err != 0 )
    {
        std::cout << "MPI Initialization did not succeed.\n";
        exit( 1 );
    }
#endif

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

    ErrorCode result = MB_SUCCESS;
    bool help        = false;
    double toler     = 5.e-10;
    result = get_file_options( argc, argv, meshFiles, method, interpTag, gNormTag, ssNormTag, ssTagNames, ssTagValues,
                               readOpts, outFile, writeOpts, dbgFile, help, toler );

    if( result != MB_SUCCESS || help )
    {
        print_usage( argv );
#ifdef MOAB_HAVE_MPI
        err = MPI_Finalize();
        if( err != 0 )
        {
            std::cout << "MPI Initialization did not succeed.\n";
            exit( 1 );
        }
#endif
        return 1;
    }

    int nprocs = 1, rank = 0;

#ifdef MOAB_HAVE_MPI
    err = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    err = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    std::vector< ParallelComm* > pcs( meshFiles.size() );
#endif

    // 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 2;
        if( !std::freopen( dfname.str().c_str(), "a", stderr ) ) return 2;
    }

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

    // Read in mesh(es)

    // Create root sets for each mesh using moab
    std::vector< EntityHandle > roots( meshFiles.size() );

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

        result = mbImpl->create_meshset( MESHSET_SET, roots[i] );
        PRINT_LAST_ERROR;
        result = mbImpl->load_file( meshFiles[i].c_str(), &roots[i], newReadopts.c_str() );
        PRINT_LAST_ERROR;
    }

#ifdef MOAB_HAVE_MPI
    result = report_iface_ents( mbImpl, pcs, true );
    PRINT_LAST_ERROR;
#endif

    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, pcs,
                                 instant_time, pointloc_time, interp_time, gnorm_time, ssnorm_time, toler );
    PRINT_LAST_ERROR;

    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 = writeOpts;
        std::ostringstream extraOpt;
#ifdef MOAB_HAVE_MPI
        extraOpt << ";PARALLEL_COMM=" << 1;
        newwriteOpts += extraOpt.str();
#endif
        result = mbImpl->write_file( outFile.c_str(), NULL, newwriteOpts.c_str(), partSets );
        PRINT_LAST_ERROR;
        std::cout << "Wrote " << outFile << std::endl;
        std::cout << "mbcoupler_test complete." << std::endl;
    }

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

    delete mbImpl;
    // May be leaking iMeshInst, don't care since it's end of program. Remove above deletes?

#ifdef MOAB_HAVE_MPI
    err = MPI_Finalize();
#endif

    return 0;
}

#ifdef MOAB_HAVE_MPI
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 );

    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;
}
#endif

// 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;
}

// New get_file_options() function with added possibilities for mbcoupler_test.
ErrorCode get_file_options( int argc,
                            char** argv,
                            std::vector< std::string >& meshFiles,
                            DataCoupler::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                   = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_"
                                 "RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME";
    outFile                    = "";
    writeOpts                  = "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 <interp_tag>" << 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 = DataCoupler::CONSTANT;
            else if( argv[npos][0] == '1' )
                method = DataCoupler::LINEAR_FE;
            else if( argv[npos][0] == '2' )
                method = DataCoupler::QUADRATIC_FE;
            else if( argv[npos][0] == '3' )
                method = DataCoupler::SPECTRAL;
            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 <epsilon>" << 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 <gnorm_tag>" << 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 <ssnorm_tag>" << 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 "
                                         "<ssnorm_selection>"
                                      << 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 <ssnorm_selection>" << 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 <roptions>" << 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 <out_file>" << 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 <woptions>" << 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" );
        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::cout << "No arguments given; using output file dum.h5m." << std::endl;
        outFile = "dum.h5m";
    }
#endif

    return MB_SUCCESS;
}

// End new get_file_options()

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

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

    CpuTimer timer;

    // Instantiate a coupler, which also initializes the tree
    DataCoupler dc( mbImpl, src_elems, 0, pcs[0] );

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

    instant_time = timer.time_since_birth();

    // 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;
#ifdef MOAB_HAVE_MPI
    result = pcs[1]->get_part_entities( targ_elems, 3 );
#endif
    PRINT_LAST_ERROR;

    // First get all vertices adj to partition entities in target mesh
    if( DataCoupler::CONSTANT == method )
        targ_verts = targ_elems;
    else
        result = mbImpl->get_adjacencies( targ_elems, 0, false, targ_verts, Interface::UNION );
    PRINT_LAST_ERROR;

#ifdef MOAB_HAVE_MPI
    // Then get non-owned verts and subtract
    Range tmp_verts;
    result = pcs[1]->get_pstatus_entities( 0, PSTATUS_NOT_OWNED, tmp_verts );
    PRINT_LAST_ERROR;
    targ_verts = subtract( targ_verts, tmp_verts );
#endif
    // 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] );
    PRINT_LAST_ERROR;

    // Locate those points in the source mesh
#ifdef MOAB_HAVE_MPI
    std::cout << "rank " << pcs[0]->proc_config().proc_rank();
#endif
    std::cout << " points of interest: " << numPointsOfInterest << "\n";
    result = dc.locate_points( &vpos[0], numPointsOfInterest, toler );
    PRINT_LAST_ERROR;

    pointloc_time = timer.time_elapsed();

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

    result = dc.interpolate( method, interpTag, &field[0] );
    PRINT_LAST_ERROR;

    interp_time = timer.time_elapsed();

    // Set field values as tag on target vertices
    // Use original tag
    Tag tag;
    result = mbImpl->tag_get_handle( interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag );
    PRINT_LAST_ERROR;
    result = mbImpl->tag_set_data( tag, targ_verts, &field[0] );
    PRINT_LAST_ERROR;

    // Done
    return MB_SUCCESS;
}

#else

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

#endif