umr.cpp

Go to the documentation of this file.

#include <iostream>
#include <vector>
#include <string>
#include <iomanip>
#include <fstream>
#include <ctime>
#include <cmath>
#include <cassert>
#include <limits>

#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/verdict/VerdictWrapper.hpp"
#include "moab/NestedRefine.hpp"

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

/* Exit values */
#define SUCCESS         0
#define USAGE_ERROR     1
#define NOT_IMPLEMENTED 2

using namespace moab;

static void print_usage( const char* name, std::ostream& stream )
{
    stream << "Usage: " << name << " <options> <input_file> [<input_file2> ...]" << std::endl
           << "Options: " << std::endl
           << "\t-h             - Print this help text and exit." << std::endl
           << "\t-d             - Dimension of the mesh." << std::endl
           << "\t-n             - Exact or a maximum number of levels for the hierarchy. Default 1." << std::endl
           << "\t-L             - Degree of refinement for each level. Pass an array or a number. It "
              "is mandatory to pass dimension and num_levels before to use this option. If this flag "
              "is not used, a default of deg 2 refinement is used. "
           << std::endl
           << "\t-V             - Pass a desired volume (absolute) . This will generate a hierarchy "
              "such that the maximum volume is reduced to the given amount approximately. The length "
              "of the hierarchy can be constrained further if a maximum number of levels is passed. "
              "It is mandatory to pass the dimension for this option. "
           << std::endl
           << "\t-q             - Prints out the maximum volume of the mesh and exits the program. "
              "This option can be used as a guide to volume constrainted mesh hierarchy later. "
           << std::endl
           << "\t-t             - Print out the time taken to generate hierarchy." << std::endl
           << "\t-s             - Print out the mesh sizes of each level of the generated hierarchy." << std::endl
           << "\t-o             - Specify true for output files for the mesh levels of the hierarchy." << std::endl
    //<< "\t-O option      - Specify read option." << std::endl
#ifdef MOAB_HAVE_MPI
           << "\t-p[0|1|2]      - Read in parallel[0], optionally also doing resolve_shared_ents (1) "
              "and exchange_ghosts (2)"
           << std::endl
#endif
        ;
    exit( USAGE_ERROR );
}

static void usage_error( const char* name )
{
    print_usage( name, std::cerr );
#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
    exit( USAGE_ERROR );
}

bool parse_id_list( const char* string, int dim, int nval, std::vector< int >& results );

bool make_opts_string( std::vector< std::string > options, std::string& opts );

ErrorCode get_degree_seq( Core& mb, EntityHandle fileset, int dim, double desired_vol, int& num_levels,
                          std::vector< int >& level_degs );

ErrorCode get_max_volume( Core& mb, EntityHandle fileset, int dim, double& vmax );

int main( int argc, char* argv[] )
{
    int proc_id = 0, size = 1;
#ifdef MOAB_HAVE_MPI
    MPI_Init( &argc, &argv );
    MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
    MPI_Comm_size( MPI_COMM_WORLD, &size );
#endif

    int num_levels = 0, dim = 0;
    std::vector< int > level_degrees;
    bool optimize    = false;
    bool do_flag     = true;
    bool print_times = false, print_size = false, output = false;
    bool parallel = false, resolve_shared = false, exchange_ghosts = false;
    bool printusage = false, parselevels = true;
    bool qc_vol = false, only_quality = false;
    double cvol = 0;
    std::string infile;

    int i = 1;
    while( i < argc )
    {
        if( !argv[i][0] && 0 == proc_id )
        {
            usage_error( argv[0] );
#ifdef MOAB_HAVE_MPI
            MPI_Finalize();
#endif
        }

        if( do_flag && argv[i][0] == '-' )
        {
            switch( argv[i][1] )
            {
                    // do flag arguments:
                case '-':
                    do_flag = false;
                    break;
                case 't':
                    print_times = true;
                    break;
                case 's':
                    print_size = true;
                    break;
                case 'h':
                case 'H':
                    print_usage( argv[0], std::cerr );
                    printusage = true;
                    break;
                case 'd':
                    dim = atol( argv[i + 1] );
                    ++i;
                    break;
                case 'n':
                    num_levels = atol( argv[i + 1] );
                    ++i;
                    break;
                case 'L':
                    if( dim != 0 && num_levels != 0 )
                    {
                        parselevels = parse_id_list( argv[i + 1], dim, num_levels, level_degrees );
                        ++i;
                    }
                    else
                    {
                        print_usage( argv[0], std::cerr );
                        printusage = true;
                    }
                    break;
                case 'V':
                    qc_vol = true;
                    cvol   = strtod( argv[i + 1], NULL );
                    ++i;
                    break;
                case 'q':
                    only_quality = true;
                    break;
                case 'o':
                    output = true;
                    break;
                case 'O':
                    optimize = true;
                    break;
#ifdef MOAB_HAVE_MPI
                case 'p':
                    parallel       = true;
                    resolve_shared = true;
                    if( argv[i][1] == '1' ) exchange_ghosts = true;
                    break;
#endif
                default:
                    ++i;
                    switch( argv[i - 1][1] )
                    {
                            // case 'O':  read_opts.push_back(argv[i]); break;
                        default:
                            std::cerr << "Invalid option: " << argv[i] << std::endl;
                    }
            }
            ++i;
        }
        // do file names
        else
        {
            infile = argv[i];
            ++i;
        }
    }

    if( infile.empty() && !printusage ) print_usage( argv[0], std::cerr );

    if( !parselevels ) exit( USAGE_ERROR );

#ifdef MOAB_HAVE_MPI
    parallel       = true;
    resolve_shared = true;
#endif

    ErrorCode error;

    Core* moab    = new Core;
    Interface* mb = (Interface*)moab;
    EntityHandle fileset;

    // Create a fileset
    error = mb->create_meshset( MESHSET_SET, fileset );MB_CHK_ERR( error );

    // Set the read options for parallel file loading
    std::vector< std::string > read_opts, write_opts;
    std::string read_options, write_options;

    if( parallel && size > 1 )
    {
        read_opts.push_back( "PARALLEL=READ_PART" );
        read_opts.push_back( "PARTITION=PARALLEL_PARTITION" );
        if( resolve_shared ) read_opts.push_back( "PARALLEL_RESOLVE_SHARED_ENTS" );
        if( exchange_ghosts ) read_opts.push_back( "PARALLEL_GHOSTS=3.0.1" );

        /*  if (output)
            {
              write_opts.push_back(";;PARALLEL=WRITE_PART");
              std::cout<<"Here"<<std::endl;
            }*/
    }

    if( !make_opts_string( read_opts, read_options ) || !make_opts_string( write_opts, write_options ) )
    {
#ifdef MOAB_HAVE_MPI
        MPI_Finalize();
#endif
        return USAGE_ERROR;
    }

    // Load file
    std::cout << "READ OPTS=" << (char*)read_options.c_str() << std::endl;
    error = mb->load_file( infile.c_str(), &fileset, read_options.c_str() );MB_CHK_ERR( error );

    // Create the nestedrefine instance

#ifdef MOAB_HAVE_MPI
    ParallelComm* pc   = new ParallelComm( mb, MPI_COMM_WORLD );
    NestedRefine* uref = new NestedRefine( moab, pc, fileset );
#else
    NestedRefine* uref = new NestedRefine( moab );
#endif

    std::vector< EntityHandle > lsets;

    // std::cout<<"Read input file"<<std::endl;

    if( only_quality )
    {
        double vmax;
        error = get_max_volume( *moab, fileset, dim, vmax );MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
        int rank = 0;
        MPI_Comm_rank( MPI_COMM_WORLD, &rank );
        if( rank == 0 ) std::cout << "Maximum global volume = " << vmax << std::endl;
#else
        std::cout << "Maximum volume = " << vmax << std::endl;
#endif
        exit( SUCCESS );
    }

    // If a volume constraint is given, find an optimal degree sequence to reach the desired volume
    // constraint.
    if( qc_vol )
    {
        error = get_degree_seq( *moab, fileset, dim, cvol, num_levels, level_degrees );MB_CHK_ERR( error );

        if( dim == 0 ) print_usage( argv[0], std::cerr );
    }

    if( num_levels == 0 ) num_levels = 1;

    int* ldeg = new int[num_levels];
    // std::cout<<"level_degrees.size = "<<level_degrees.size()<<std::endl;
    if( level_degrees.empty() )
    {
        for( int l = 0; l < num_levels; l++ )
            ldeg[l] = 2;
    }
    else
    {
        for( int l = 0; l < num_levels; l++ )
            ldeg[l] = level_degrees[l];
    }

    std::cout << "Starting hierarchy generation" << std::endl;

    std::cout << "opt = " << optimize << std::endl;

    error = uref->generate_mesh_hierarchy( num_levels, ldeg, lsets, optimize );MB_CHK_ERR( error );

    if( print_times )
    {
        std::cout << "Finished hierarchy generation in " << uref->timeall.tm_total << "  secs" << std::endl;
        if( parallel )
        {
            std::cout << "Time taken for refinement " << uref->timeall.tm_refine << "  secs" << std::endl;
            std::cout << "Time taken for resolving shared interface " << uref->timeall.tm_resolve << "  secs"
                      << std::endl;
        }
    }
    else
        std::cout << "Finished hierarchy generation " << std::endl;

    if( print_size )
    {
        Range all_ents, ents[4];
        error = mb->get_entities_by_handle( fileset, all_ents );MB_CHK_ERR( error );

        for( int k = 0; k < 4; k++ )
            ents[k] = all_ents.subset_by_dimension( k );

        std::cout << std::endl;
        if( qc_vol )
        {
            double volume;
            error = get_max_volume( *moab, fileset, dim, volume );MB_CHK_ERR( error );
            std::cout << "Mesh size for level 0"
                      << "  :: nverts = " << ents[0].size() << ", nedges = " << ents[1].size()
                      << ", nfaces = " << ents[2].size() << ", ncells = " << ents[3].size() << " :: Vmax = " << volume
                      << std::endl;
        }
        else
            std::cout << "Mesh size for level 0"
                      << "  :: nverts = " << ents[0].size() << ", nedges = " << ents[1].size()
                      << ", nfaces = " << ents[2].size() << ", ncells = " << ents[3].size() << std::endl;

        for( int l = 0; l < num_levels; l++ )
        {
            all_ents.clear();
            ents[0].clear();
            ents[1].clear();
            ents[2].clear();
            ents[3].clear();
            error = mb->get_entities_by_handle( lsets[l + 1], all_ents );MB_CHK_ERR( error );

            for( int k = 0; k < 4; k++ )
                ents[k] = all_ents.subset_by_dimension( k );

            // std::cout<<std::endl;

            if( qc_vol )
            {
                double volume;
                error = get_max_volume( *moab, lsets[l + 1], dim, volume );MB_CHK_ERR( error );
                std::cout << "Mesh size for level " << l + 1 << "  :: nverts = " << ents[0].size()
                          << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size()
                          << ", ncells = " << ents[3].size() << " :: Vmax = " << volume << std::endl;
            }
            else
                std::cout << "Mesh size for level " << l + 1 << "  :: nverts = " << ents[0].size()
                          << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size()
                          << ", ncells = " << ents[3].size() << std::endl;
        }
    }

    if( output )
    {
        for( int l = 0; l < num_levels; l++ )
        {
            std::string::size_type idx1 = infile.find_last_of( "\\/" );
            std::string::size_type idx2 = infile.find_last_of( "." );
            std::string file            = infile.substr( idx1 + 1, idx2 - idx1 - 1 );
            std::stringstream out;
            if( parallel )
                out << "_ML_" << l + 1 << ".h5m";
            else
                out << "_ML_" << l + 1 << ".vtk";
            file                    = file + out.str();
            const char* output_file = file.c_str();
#ifdef MOAB_HAVE_MPI
            error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART", &lsets[l + 1], 1 );MB_CHK_ERR( error );
#else
            error = mb->write_file( output_file, 0, NULL, &lsets[l + 1], 1 );MB_CHK_ERR( error );
#endif
            //   const char* output_file = file.c_str();
            //   error =  mb->write_file(output_file, 0, write_options.c_str(), &lsets[l+1],
            //   1);MB_CHK_ERR(error);
            //    mb->list_entity(lsets[l+1]);
            //   mb->write_file(output_file, 0, "PARALLEL=WRITE_PART", &lsets[l+1], 1);
        }
    }

    delete uref;
#ifdef MOAB_HAVE_MPI
    delete pc;
#endif

    delete[] ldeg;
    delete moab;

#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif

    exit( SUCCESS );
}

ErrorCode get_degree_seq( Core& mb, EntityHandle fileset, int dim, double desired_vol, int& num_levels,
                          std::vector< int >& level_degs )
{
    // Find max volume
    double vmax_global;
    ErrorCode error = get_max_volume( mb, fileset, dim, vmax_global );MB_CHK_ERR( error );

    int init_nl = num_levels;
    num_levels  = 0;
    level_degs.clear();

    // Find a sequence that reduces the maximum volume to desired.
    // desired_vol should be a fraction or absolute value ?

    // double remV = vmax_global*desired_vol/dim;
    double Vdesired = desired_vol;
    double try_x;
    double remV    = vmax_global;
    int degs[3][3] = { { 5, 3, 2 }, { 25, 9, 4 }, { 0, 27, 8 } };

    if( dim == 1 || dim == 2 )
    {
        while( remV - Vdesired >= 0 )
        {
            try_x = degs[dim - 1][0];
            if( ( remV / try_x - Vdesired ) >= 0 )
            {
                level_degs.push_back( 5 );
                num_levels += 1;
                remV /= try_x;
            }
            else
            {
                try_x = degs[dim - 1][1];
                if( ( remV / try_x - Vdesired ) >= 0 )
                {
                    level_degs.push_back( 3 );
                    num_levels += 1;
                    remV /= try_x;
                }
                else
                {
                    try_x = degs[dim - 1][2];
                    if( ( remV / try_x - Vdesired ) >= 0 )
                    {
                        level_degs.push_back( 2 );
                        num_levels += 1;
                        remV /= try_x;
                    }
                    else
                        break;
                }
            }
        }
    }
    else
    {
        while( remV - Vdesired >= 0 )
        {
            try_x = degs[dim - 1][1];
            if( ( remV / try_x - Vdesired ) >= 0 )
            {
                level_degs.push_back( 3 );
                num_levels += 1;
                remV /= try_x;
            }
            else
            {
                try_x = degs[dim - 1][2];
                if( ( remV / try_x - Vdesired ) >= 0 )
                {
                    level_degs.push_back( 2 );
                    num_levels += 1;
                    remV /= try_x;
                }
                else
                    break;
            }
        }
    }

    if( init_nl != 0 && init_nl < num_levels )
    {
        for( int i = level_degs.size(); i >= init_nl; i-- )
            level_degs.pop_back();
        num_levels = init_nl;
    }

    return MB_SUCCESS;
}

ErrorCode get_max_volume( Core& mb, EntityHandle fileset, int dim, double& vmax )
{
    ErrorCode error;
    VerdictWrapper vw( &mb );
    QualityType q;

    switch( dim )
    {
        case 1:
            q = MB_LENGTH;
            break;
        case 2:
            q = MB_AREA;
            break;
        case 3:
            q = MB_VOLUME;
            break;
        default:
            return MB_FAILURE;
            break;
    }

    // Get all entities of the highest dimension which is passed as a command line argument.
    Range allents, owned;
    error = mb.get_entities_by_handle( fileset, allents );MB_CHK_ERR( error );
    owned = allents.subset_by_dimension( dim );MB_CHK_ERR( error );

    // Get all owned entities
#ifdef MOAB_HAVE_MPI
    int size = 1;
    MPI_Comm_size( MPI_COMM_WORLD, &size );
    int mpi_err;
    if( size > 1 )
    {
        // filter the entities not owned, so we do not process them more than once
        ParallelComm* pcomm = moab::ParallelComm::get_pcomm( &mb, 0 );
        Range current       = owned;
        owned.clear();
        error = pcomm->filter_pstatus( current, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned );
        if( error != MB_SUCCESS )
        {
            MPI_Finalize();
            return MB_FAILURE;
        }
    }
#endif

    double vmax_local = 0;
    // Find the maximum volume of an entity in the owned mesh
    for( Range::iterator it = owned.begin(); it != owned.end(); it++ )
    {
        double volume;
        error = vw.quality_measure( *it, q, volume );MB_CHK_ERR( error );
        if( volume > vmax_local ) vmax_local = volume;
    }

    // Get the global maximum
    double vmax_global = vmax_local;
#ifdef MOAB_HAVE_MPI
    mpi_err = MPI_Reduce( &vmax_local, &vmax_global, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
    if( mpi_err )
    {
        MPI_Finalize();
        return MB_FAILURE;
    }
#endif

    vmax = vmax_global;

    return MB_SUCCESS;
}

bool parse_id_list( const char* string, int dim, int nval, std::vector< int >& results )
{
    bool okay   = true;
    char* mystr = strdup( string );
    for( const char* ptr = strtok( mystr, "," ); ptr; ptr = strtok( 0, "," ) )
    {
        char* endptr;
        int val = strtol( ptr, &endptr, 0 );

        if( dim == 1 || dim == 2 )
        {
            if( val != 2 && val != 3 && val != 5 )
            {
                std::cerr << "Not a valid degree for the passed dimension" << std::endl;
                okay = false;
                break;
            }
        }
        else if( dim == 3 )
        {
            if( val != 2 && val != 3 )
            {
                std::cerr << "Not a valid degree for the passed dimension" << std::endl;
                okay = false;
                break;
            }
        }

        if( endptr == ptr || val <= 0 )
        {
            std::cerr << "Not a valid id: " << ptr << std::endl;
            okay = false;
            break;
        }

        results.push_back( val );
    }

    if( (int)results.size() < nval )
    {
        for( int i = results.size(); i <= nval - 1; i++ )
            results.push_back( results[0] );
    }
    else if( (int)results.size() > nval )
    {
        for( int i = results.size(); i > nval; i-- )
            results.pop_back();
    }

    free( mystr );
    return okay;
}

bool make_opts_string( std::vector< std::string > options, std::string& opts )
{
    opts.clear();
    if( options.empty() ) return true;

    // choose a separator character
    std::vector< std::string >::const_iterator i;
    char separator             = '\0';
    const char* alt_separators = ";+,:\t\n";
    for( const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr )
    {
        bool seen = false;
        for( i = options.begin(); i != options.end(); ++i )
            if( i->find( *sep_ptr, 0 ) != std::string::npos )
            {
                seen = true;
                break;
            }
        if( !seen )
        {
            separator = *sep_ptr;
            break;
        }
    }
    if( !separator )
    {
        std::cerr << "Error: cannot find separator character for options string" << std::endl;
        return false;
    }
    if( separator != ';' )
    {
        opts = ";";
        opts += separator;
    }

    // concatenate options
    i = options.begin();
    opts += *i;
    for( ++i; i != options.end(); ++i )
    {
        opts += separator;
        opts += *i;
    }

    return true;
}