DeformMeshRemap.cpp

Description: Account for mesh deformation of a solid due to structural mechanics
In this example there are two meshes, a "master" and "slave" mesh. In the master mesh, the solid material is deformed, to mimic what happens when a solid heats up and deforms. The fluid mesh is smoothed to account for those deformations, and tags on the fluid are remapped to those new positions. Then mesh positions and state variables are transferred to the slave mesh, mimicing another mesh used by some other physics.

To run: ./DeformMeshRemap [<master_meshfile> <slave_meshfile>]
(default values can run if users don't specify the mesh files)

/** @example DeformMeshRemap.cpp
 * Description: Account for mesh deformation of a solid due to structural mechanics\n
 * In this example there are two meshes, a "master" and "slave" mesh.  In the master mesh,
 * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
 * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
 * remapped to those new positions.  Then mesh positions and state variables are transferred
 * to the slave mesh, mimicing another mesh used by some other physics.
 *
 * To run: ./DeformMeshRemap [<master_meshfile> <slave_meshfile>]\n
 * (default values can run if users don't specify the mesh files)
 */

#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/Skinner.hpp"
#include "moab/LloydSmoother.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/BoundBox.hpp"
#include "moab/SpatialLocator.hpp"
#include "MBTagConventions.hpp"
#include "DataCoupler.hpp"

#ifdef USE_MPI
#include "moab/ParallelComm.hpp"
#endif

#include <iostream>
#include <set>
#include <sstream>
#include <cassert>

using namespace moab;
using namespace std;

#ifndef MESH_DIR
#define MESH_DIR "."
#endif

ErrorCode read_file( string& fname,
                     EntityHandle& seth,
                     Range& solids,
                     Range& solid_elems,
                     Range& fluids,
                     Range& fluid_elems );
void deform_func( const BoundBox& bbox, double* xold, double* xnew );
ErrorCode deform_master( Range& fluid_elems, Range& solid_elems, Tag& xnew );
ErrorCode smooth_master( int dim, Tag xnew, EntityHandle& master, Range& fluids );
ErrorCode write_to_coords( Range& elems, Tag tagh );

const int SOLID_SETNO = 100, FLUID_SETNO = 200;

Interface* mb;

const bool debug = true;

class DeformMeshRemap
{
  public:
    //! Enumerator for solid/fluid, master/slave
    enum
    {
        MASTER = 0,
        SLAVE,
        SOLID,
        FLUID
    };

    //! Constructor
    //! If master is NULL, the MOAB part is run in serial;
    //! If slave is NULL but the master isn't, the slave is copied from the master
    //! Create communicators using moab::ParallelComm::get_pcomm
    DeformMeshRemap( Interface* impl, ParallelComm* master = NULL, ParallelComm* slave = NULL );

    //! Destructor
    ~DeformMeshRemap();

    //! Execute the deformed mesh process
    ErrorCode execute();

    //! Add a set number
    ErrorCode add_set_no( int m_or_s, int fluid_or_solid, int set_no );

    //! Remove a set number
    ErrorCode remove_set_no( int m_or_s, int fluid_or_solid, int set_no );

    //! Get the set numbers
    ErrorCode get_set_nos( int m_or_s, int fluid_or_solid, set< int >& set_nos ) const;

    //! Get the xNew tag handle
    inline Tag x_new() const
    {
        return xNew;
    }

    //! Get the tag name
    string x_new_name() const
    {
        return xNewName;
    }

    //! Set the tag name
    void x_new_name( const string& name )
    {
        xNewName = name;
    }

    //! Get/set the file name
    string get_file_name( int m_or_s ) const;

    //! Get/set the file name
    void set_file_name( int m_or_s, const string& name );

    //! Get/set the x displacement tag names
    string xdisp_name( int idx = 0 );
    void xdisp_name( const string& nm, int idx = 0 );

  private:
    //! Apply a known deformation to the solid elements, putting the results in the xNew tag; also
    //! write current coordinates to the xNew tag for fluid elements
    ErrorCode deform_master( Range& fluid_elems, Range& solid_elems, const char* tag_name = NULL );

    //! Read a file and establish proper ranges
    ErrorCode read_file( int m_or_s, string& fname, EntityHandle& seth );

    //! Write the input tag to the coordinates for the vertices in the input elems
    //! If non-zero tmp_tag is input, save coords to tmp_tag before over-writing with tag value
    ErrorCode write_to_coords( Range& elems, Tag tagh, Tag tmp_tag = 0 );

    //! Write the tag to the vertices, then save to the specified file
    //! If restore_coords is true, coords are restored to their initial state after file is written
    ErrorCode write_and_save( Range& ents,
                              EntityHandle seth,
                              Tag tagh,
                              const char* filename,
                              bool restore_coords = false );

    //! Find fluid/solid sets from complement of solid/fluid sets
    ErrorCode find_other_sets( int m_or_s, EntityHandle file_set );

    //! moab interface
    Interface* mbImpl;

#ifdef USE_MPI
    //! ParallelComm for master, slave meshes
    ParallelComm *pcMaster, *pcSlave;
#endif

    //! Material set numbers for fluid materials, for master/slave
    set< int > fluidSetNos[2];

    //! Material set numbers for solid materials, for master/slave
    set< int > solidSetNos[2];

    //! Sets defining master/slave meshes
    EntityHandle masterSet, slaveSet;

    //! Sets in master/slave meshes
    Range fluidSets[2], solidSets[2];

    //! Elements in master/slave meshes
    Range fluidElems[2], solidElems[2];

    //! Filenames for master/slave meshes
    string masterFileName, slaveFileName;

    //! Tag from file, might be 3
    Tag xDisp[3];

    //! Tag used for new positions
    Tag xNew;

    //! Tag name used to read disps from file
    string xDispNames[3];

    //! Tag name used for new positions
    string xNewName;
};

//! Add a set number
inline ErrorCode DeformMeshRemap::add_set_no( int m_or_s, int f_or_s, int set_no )
{
    set< int >* this_set;
    assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
    if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;

    switch( f_or_s )
    {
        case FLUID:
            this_set = &fluidSetNos[m_or_s];
            break;
        case SOLID:
            this_set = &solidSetNos[m_or_s];
            break;
        default:
            assert( false && "f_or_s should be FLUID or SOLID." );
            return MB_FAILURE;
    }

    this_set->insert( set_no );

    return MB_SUCCESS;
}

//! Remove a set number
inline ErrorCode DeformMeshRemap::remove_set_no( int m_or_s, int f_or_s, int set_no )
{
    set< int >* this_set;
    assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
    if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
    switch( f_or_s )
    {
        case FLUID:
            this_set = &fluidSetNos[m_or_s];
            break;
        case SOLID:
            this_set = &solidSetNos[m_or_s];
            break;
        default:
            assert( false && "f_or_s should be FLUID or SOLID." );
            return MB_FAILURE;
    }
    set< int >::iterator sit = this_set->find( set_no );
    if( sit != this_set->end() )
    {
        this_set->erase( *sit );
        return MB_SUCCESS;
    }

    return MB_FAILURE;
}

//! Get the set numbers
inline ErrorCode DeformMeshRemap::get_set_nos( int m_or_s, int f_or_s, set< int >& set_nos ) const
{
    const set< int >* this_set;
    assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
    if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
    switch( f_or_s )
    {
        case FLUID:
            this_set = &fluidSetNos[m_or_s];
            break;
        case SOLID:
            this_set = &solidSetNos[m_or_s];
            break;
        default:
            assert( false && "f_or_s should be FLUID or SOLID." );
            return MB_FAILURE;
    }

    set_nos = *this_set;

    return MB_SUCCESS;
}

inline string DeformMeshRemap::xdisp_name( int idx )
{
    return xDispNames[idx];
}

void DeformMeshRemap::xdisp_name( const string& nm, int idx )
{
    xDispNames[idx] = nm;
}

ErrorCode DeformMeshRemap::execute()
{
    // Read master/slave files and get fluid/solid material sets
    ErrorCode rval = read_file( MASTER, masterFileName, masterSet );MB_CHK_ERR( rval );

    if( solidSetNos[MASTER].empty() || fluidSetNos[MASTER].empty() )
    {
        rval = find_other_sets( MASTER, masterSet );MB_CHK_SET_ERR( rval, "Failed to find other sets in master mesh" );
    }

    bool have_slave = !( slaveFileName == "none" );
    if( have_slave )
    {
        rval = read_file( SLAVE, slaveFileName, slaveSet );MB_CHK_ERR( rval );

        if( solidSetNos[SLAVE].empty() || fluidSetNos[SLAVE].empty() )
        {
            rval = find_other_sets( SLAVE, slaveSet );MB_CHK_SET_ERR( rval, "Failed to find other sets in slave mesh" );
        }
    }

    if( debug ) cout << "Constructing data coupler/search tree on master mesh..." << endl;

    Range src_elems = solidElems[MASTER];
    src_elems.merge( fluidElems[MASTER] );

    // Initialize data coupler on source elements
    DataCoupler dc_master( mbImpl, src_elems, 0, NULL );

    Range tgt_verts;
    if( have_slave )
    {
        // Locate slave vertices in master, orig coords; do this with a data coupler, so you can
        // later interpolate
        Range tmp_range = solidElems[SLAVE];
        tmp_range.merge( fluidElems[SLAVE] );
        rval = mbImpl->get_adjacencies( tmp_range, 0, false, tgt_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get target verts" );

        // Locate slave vertices, caching results in dc
        if( debug ) cout << "Locating slave vertices in master mesh..." << endl;
        rval = dc_master.locate_points( tgt_verts );MB_CHK_SET_ERR( rval, "Point location of tgt verts failed" );
        int num_located = dc_master.spatial_locator()->local_num_located();
        if( num_located != (int)tgt_verts.size() )
        {
            rval = MB_FAILURE;
            cout << "Only " << num_located << " out of " << tgt_verts.size() << " target points successfully located."
                 << endl;
            return rval;
        }
    }

    // Deform the master's solid mesh, put results in a new tag
    if( debug ) cout << "Deforming fluid elements in master mesh..." << endl;
    rval = deform_master( fluidElems[MASTER], solidElems[MASTER], "xnew" );MB_CHK_ERR( rval );

    {  // To isolate the lloyd smoother & delete when done
        if( debug )
        {
            // Output the skin of smoothed elems, as a check
            // Get the skin; get facets, because we might need to filter on shared entities
            Skinner skinner( mbImpl );
            Range skin;
            rval = skinner.find_skin( 0, fluidElems[MASTER], false, skin );MB_CHK_SET_ERR( rval, "Unable to find skin" );
            EntityHandle skin_set;
            cout << "Writing skin_mesh.g and fluid_mesh.g." << endl;
            rval = mbImpl->create_meshset( MESHSET_SET, skin_set );MB_CHK_SET_ERR( rval, "Failed to create skin set" );
            rval = mbImpl->add_entities( skin_set, skin );MB_CHK_SET_ERR( rval, "Failed to add skin entities to set" );
            rval = mbImpl->write_file( "skin_mesh.vtk", NULL, NULL, &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failure to write skin set" );
            rval = mbImpl->remove_entities( skin_set, skin );MB_CHK_SET_ERR( rval, "Failed to remove skin entities from set" );
            rval = mbImpl->add_entities( skin_set, fluidElems[MASTER] );MB_CHK_SET_ERR( rval, "Failed to add fluid entities to set" );
            rval = mbImpl->write_file( "fluid_mesh.vtk", NULL, NULL, &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failure to write fluid set" );
            rval = mbImpl->delete_entities( &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failed to delete skin set" );
        }

        // Smooth the master mesh
        if( debug ) cout << "Smoothing fluid elements in master mesh..." << endl;
        LloydSmoother ll( mbImpl, NULL, fluidElems[MASTER], xNew );
        rval = ll.perform_smooth();MB_CHK_SET_ERR( rval, "Failed in lloyd smoothing" );
        cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
    }

    // Transfer xNew to coords, for master
    if( debug ) cout << "Transferring coords tag to vertex coordinates in master mesh..." << endl;
    rval = write_to_coords( solidElems[MASTER], xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to master fluid verts" );
    rval = write_to_coords( fluidElems[MASTER], xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to master fluid verts" );

    if( have_slave )
    {
        // Map new locations to slave
        // Interpolate xNew to slave points
        if( debug ) cout << "Interpolating new coordinates to slave vertices..." << endl;
        rval = dc_master.interpolate( (int)DataCoupler::VOLUME, "xnew" );MB_CHK_SET_ERR( rval, "Failed to interpolate target solution" );
        // Transfer xNew to coords, for slave
        if( debug ) cout << "Transferring coords tag to vertex coordinates in slave mesh..." << endl;
        rval = write_to_coords( tgt_verts, xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to slave verts" );
    }

    if( debug )
    {
        string str;
#ifdef USE_MPI
        if( pcMaster && pcMaster->size() > 1 ) str = "PARALLEL=WRITE_PART";
#endif
        if( debug ) cout << "Writing smoothed_master.h5m..." << endl;
        rval = mbImpl->write_file( "smoothed_master.h5m", NULL, str.c_str(), &masterSet, 1 );

        if( have_slave )
        {
#ifdef USE_MPI
            str.clear();
            if( pcSlave && pcSlave->size() > 1 ) str = "PARALLEL=WRITE_PART";
#endif
            if( debug ) cout << "Writing slave_interp.h5m..." << endl;
            rval = mbImpl->write_file( "slave_interp.h5m", NULL, str.c_str(), &slaveSet, 1 );
        }  // if have_slave
    }      // if debug

    if( debug ) dc_master.spatial_locator()->get_tree()->tree_stats().print();

    return MB_SUCCESS;
}

string DeformMeshRemap::get_file_name( int m_or_s ) const
{
    switch( m_or_s )
    {
        case MASTER:
            return masterFileName;
        case SLAVE:
            return slaveFileName;
        default:
            assert( false && "m_or_s should be MASTER or SLAVE." );
            return string();
    }
}

void DeformMeshRemap::set_file_name( int m_or_s, const string& name )
{
    switch( m_or_s )
    {
        case MASTER:
            masterFileName = name;
            break;
        case SLAVE:
            slaveFileName = name;
            break;
        default:
            assert( false && "m_or_s should be MASTER or SLAVE." );
    }
}

DeformMeshRemap::DeformMeshRemap( Interface* impl, ParallelComm* master, ParallelComm* slave )
    : mbImpl( impl ), pcMaster( master ), pcSlave( slave ), masterSet( 0 ), slaveSet( 0 ), xNew( 0 ), xNewName( "xnew" )
{
    xDisp[0] = xDisp[1] = xDisp[2] = 0;

    if( !pcSlave && pcMaster ) pcSlave = pcMaster;
}

DeformMeshRemap::~DeformMeshRemap() {}

int main( int argc, char** argv )
{
    ErrorCode rval;

    ProgOptions po( "Deformed mesh options" );
    po.addOpt< string >( "master,m", "Specify the master meshfile name" );
    po.addOpt< string >( "worker,w", "Specify the slave/worker meshfile name, or 'none' (no quotes) if master only" );
    po.addOpt< string >( "d1,", "Tag name for displacement x or xyz" );
    po.addOpt< string >( "d2,", "Tag name for displacement y" );
    po.addOpt< string >( "d3,", "Tag name for displacement z" );
    po.addOpt< int >( "fm,", "Specify master fluid material set number(s). If none specified, "
                             "fluid sets derived from complement of solid sets." );
    po.addOpt< int >( "fs,", "Specify master solid material set number(s). If none specified, "
                             "solid sets derived from complement of fluid sets." );
    po.addOpt< int >( "sm,", "Specify slave fluid material set number(s). If none specified, fluid "
                             "sets derived from complement of solid sets." );
    po.addOpt< int >( "ss,", "Specify slave solid material set number(s). If none specified, solid "
                             "sets derived from complement of fluid sets." );

    po.parseCommandLine( argc, argv );

    mb = new Core();

    DeformMeshRemap* dfr;
#ifdef USE_MPI
    ParallelComm* pc = new ParallelComm( mb, MPI_COMM_WORLD );
    dfr              = new DeformMeshRemap( mb, pc );
#else
    dfr = new DeformMeshRemap( mb );
#endif

    string masterf, slavef;
    if( !po.getOpt( "master", &masterf ) ) masterf = string( MESH_DIR ) + string( "/rodquad.g" );
    dfr->set_file_name( DeformMeshRemap::MASTER, masterf );

    if( !po.getOpt( "worker", &slavef ) ) slavef = string( MESH_DIR ) + string( "/rodtri.g" );
    dfr->set_file_name( DeformMeshRemap::SLAVE, slavef );
    if( slavef.empty() )
    {
        cerr << "Empty slave file name; if no slave, use filename 'none' (no quotes)." << endl;
        return 1;
    }

    vector< int > set_nos;
    po.getOptAllArgs( "fm", set_nos );
    for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
        dfr->add_set_no( DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, *vit );
    set_nos.clear();

    po.getOptAllArgs( "fs", set_nos );
    for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
        dfr->add_set_no( DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, *vit );
    set_nos.clear();

    po.getOptAllArgs( "sm", set_nos );
    for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
        dfr->add_set_no( DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, *vit );

    po.getOptAllArgs( "ss", set_nos );
    for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
        dfr->add_set_no( DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, *vit );

    string tnames[3];
    po.getOpt( "d1", &tnames[0] );
    po.getOpt( "d2", &tnames[1] );
    po.getOpt( "d3", &tnames[2] );
    for( int i = 0; i < 3; i++ )
        if( !tnames[i].empty() ) dfr->xdisp_name( tnames[i], i );

    rval = dfr->execute();

    delete dfr;
    delete mb;

    return rval;
}

ErrorCode DeformMeshRemap::write_and_save( Range& ents,
                                           EntityHandle seth,
                                           Tag tagh,
                                           const char* filename,
                                           bool restore_coords )
{
    Tag tmp_tag = 0;
    ErrorCode rval;
    if( restore_coords ) rval = mbImpl->tag_get_handle( "", 3, MB_TYPE_DOUBLE, tmp_tag, MB_TAG_CREAT | MB_TAG_DENSE );

    rval = write_to_coords( ents, tagh, tmp_tag );MB_CHK_ERR( rval );
    rval = mbImpl->write_file( filename, NULL, NULL, &seth, 1 );MB_CHK_ERR( rval );
    if( restore_coords )
    {
        rval = write_to_coords( ents, tmp_tag );MB_CHK_ERR( rval );
        rval = mbImpl->tag_delete( tmp_tag );MB_CHK_ERR( rval );
    }

    return rval;
}

ErrorCode DeformMeshRemap::write_to_coords( Range& elems, Tag tagh, Tag tmp_tag )
{
    // Write the tag to coordinates
    Range verts;
    ErrorCode rval = mbImpl->get_adjacencies( elems, 0, false, verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get adj vertices" );
    vector< double > coords( 3 * verts.size() );

    if( tmp_tag )
    {
        // Save the coords to tmp_tag first
        rval = mbImpl->get_coords( verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get tmp copy of coords" );
        rval = mbImpl->tag_set_data( tmp_tag, verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to save tmp copy of coords" );
    }

    rval = mbImpl->tag_get_data( tagh, verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get tag data" );
    rval = mbImpl->set_coords( verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to set coordinates" );
    return MB_SUCCESS;
}

void deform_func( const BoundBox& bbox, double* xold, double* xnew )
{
    /*  Deformation function based on max delx and dely at top of rod
        const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
        // function: origin is at middle base of rod, and is .5 high
        // top of rod is (0,.55) on left and (.2,.6) on right
      double delx = 0.5*RODWIDTH;

      double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
      xnew[0] = xold[0] + yfrac * delx;
      xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
    */
    /* Deformation function based on fraction of bounding box dimension in each direction */
    double frac  = 0.01;  // taken from approximate relative deformation from LLNL Diablo of XX09 assys
    CartVect *xo = reinterpret_cast< CartVect* >( xold ), *xn = reinterpret_cast< CartVect* >( xnew );
    CartVect disp = frac * ( *xo - bbox.bMin );
    *xn           = *xo + disp;
}

ErrorCode DeformMeshRemap::deform_master( Range& fluid_elems, Range& solid_elems, const char* tag_name )
{
    // Deform elements with an analytic function
    ErrorCode rval;

    // Get all the vertices and coords in the solid
    Range solid_verts, fluid_verts;
    rval = mbImpl->get_adjacencies( solid_elems, 0, false, solid_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get vertices" );
    vector< double > coords( 3 * solid_verts.size() ), new_coords( 3 * solid_verts.size() );
    rval = mbImpl->get_coords( solid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get vertex coords" );
    unsigned int num_verts = solid_verts.size();

    // Get or create the tag
    if( !xDispNames[0].empty() && !xDispNames[1].empty() && !xDispNames[2].empty() )
    {
        // 3 tags, specifying xyz individual data, integrate into one tag
        rval = mbImpl->tag_get_handle( ( tag_name ? tag_name : "" ), 3, MB_TYPE_DOUBLE, xNew,
                                       MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Failed to create xnew tag" );
        vector< double > disps( num_verts );
        for( int i = 0; i < 3; i++ )
        {
            rval = mbImpl->tag_get_handle( xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i] );MB_CHK_SET_ERR( rval, "Failed to get xDisp tag" );
            rval = mbImpl->tag_get_data( xDisp[i], solid_verts, &disps[0] );MB_CHK_SET_ERR( rval, "Failed to get xDisp tag values" );
            for( unsigned int j = 0; j < num_verts; j++ )
                new_coords[3 * j + i] = coords[3 * j + i] + disps[j];
        }
    }
    else if( !xDispNames[0].empty() )
    {
        rval = mbImpl->tag_get_handle( xDispNames[0].c_str(), 3, MB_TYPE_DOUBLE, xDisp[0] );MB_CHK_SET_ERR( rval, "Failed to get first xDisp tag" );
        xNew = xDisp[0];
        vector< double > disps( 3 * num_verts );
        rval = mbImpl->tag_get_data( xDisp[0], solid_verts, &disps[0] );
        for( unsigned int j = 0; j < 3 * num_verts; j++ )
            new_coords[j] = coords[j] + disps[j];
    }
    else
    {
        // Get the bounding box of the solid mesh
        BoundBox bbox;
        bbox.update( *mbImpl, solid_elems );

        for( unsigned int j = 0; j < num_verts; j++ )
            deform_func( bbox, &coords[3 * j], &new_coords[3 * j] );
    }

    if( debug )
    {
        double len = 0.0;
        for( unsigned int i = 0; i < num_verts; i++ )
        {
            CartVect dx    = CartVect( &new_coords[3 * i] ) - CartVect( &coords[3 * i] );
            double tmp_len = dx.length_squared();
            if( tmp_len > len ) len = tmp_len;
        }
        Range tmp_elems( fluid_elems );
        tmp_elems.merge( solid_elems );
        BoundBox box;
        box.update( *mbImpl, tmp_elems );
        double max_len =
            std::max( box.bMax[2] - box.bMin[2], std::max( box.bMax[1] - box.bMin[1], box.bMax[0] - box.bMin[0] ) );

        cout << "Max displacement = " << len << " (" << 100.0 * len / max_len << "% of max box length)" << endl;
    }

    if( !xNew )
    {
        rval = mbImpl->tag_get_handle( ( tag_name ? tag_name : "" ), 3, MB_TYPE_DOUBLE, xDisp[0],
                                       MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Failed to get xNew tag" );
        xNew = xDisp[0];
    }

    // Set the new tag to those coords
    rval = mbImpl->tag_set_data( xNew, solid_verts, &new_coords[0] );MB_CHK_SET_ERR( rval, "Failed to set tag data" );

    // Get all the vertices and coords in the fluid, set xnew to them
    rval = mbImpl->get_adjacencies( fluid_elems, 0, false, fluid_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get vertices" );
    fluid_verts = subtract( fluid_verts, solid_verts );

    if( coords.size() < 3 * fluid_verts.size() ) coords.resize( 3 * fluid_verts.size() );
    rval = mbImpl->get_coords( fluid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get vertex coords" );
    rval = mbImpl->tag_set_data( xNew, fluid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to set xnew tag on fluid verts" );

    if( debug )
    {
        // Save deformed mesh coords to new file for visualizing
        Range tmp_range( fluidElems[MASTER] );
        tmp_range.merge( solidElems[MASTER] );
        rval = write_and_save( tmp_range, masterSet, xNew, "deformed_master.h5m", true );MB_CHK_ERR( rval );
    }

    return MB_SUCCESS;
}

ErrorCode DeformMeshRemap::read_file( int m_or_s, string& fname, EntityHandle& seth )
{
    // Create meshset
    ErrorCode rval = mbImpl->create_meshset( 0, seth );MB_CHK_SET_ERR( rval, "Couldn't create master/slave set" );
    ostringstream options;
#ifdef USE_MPI
    ParallelComm* pc = ( m_or_s == MASTER ? pcMaster : pcSlave );
    if( pc && pc->size() > 1 )
    {
        if( debug ) options << "DEBUG_IO=1;CPUTIME;";
        options << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
                << "PARALLEL_GHOSTS=2.0.1;PARALLEL_COMM=" << pc->get_id();
    }
#endif
    rval = mbImpl->load_file( fname.c_str(), &seth, options.str().c_str() );MB_CHK_SET_ERR( rval, "Couldn't load master/slave mesh" );

    if( *solidSetNos[m_or_s].begin() == -1 || *fluidSetNos[m_or_s].begin() == -1 ) return MB_SUCCESS;

    // Get material sets for solid/fluid
    Tag tagh;
    rval = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, tagh );MB_CHK_SET_ERR( rval, "Couldn't get material set tag name" );
    for( set< int >::iterator sit = solidSetNos[m_or_s].begin(); sit != solidSetNos[m_or_s].end(); ++sit )
    {
        Range sets;
        int set_no            = *sit;
        const void* setno_ptr = &set_no;
        rval                  = mbImpl->get_entities_by_type_and_tag( seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets );
        if( MB_SUCCESS != rval || sets.empty() )
        {
            MB_SET_ERR( MB_FAILURE, "Couldn't find solid set #" << *sit );
        }
        else
            solidSets[m_or_s].merge( sets );
    }

    // Get solid entities, and dimension
    Range tmp_range;
    for( Range::iterator rit = solidSets[m_or_s].begin(); rit != solidSets[m_or_s].end(); ++rit )
    {
        rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in solid" );
    }
    if( !tmp_range.empty() )
    {
        int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
        assert( dim > 0 && dim < 4 );
        solidElems[m_or_s] = tmp_range.subset_by_dimension( dim );
    }

    if( debug )
        cout << "Read " << solidElems[m_or_s].size() << " solid elements from " << solidSets[m_or_s].size()
             << " sets in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh." << endl;

    for( set< int >::iterator sit = fluidSetNos[m_or_s].begin(); sit != fluidSetNos[m_or_s].end(); ++sit )
    {
        Range sets;
        int set_no            = *sit;
        const void* setno_ptr = &set_no;
        rval                  = mbImpl->get_entities_by_type_and_tag( seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets );
        if( MB_SUCCESS != rval || sets.empty() )
        {
            MB_SET_ERR( MB_FAILURE, "Couldn't find fluid set #" << *sit );
        }
        else
            fluidSets[m_or_s].merge( sets );
    }

    // Get fluid entities, and dimension
    tmp_range.clear();
    for( Range::iterator rit = fluidSets[m_or_s].begin(); rit != fluidSets[m_or_s].end(); ++rit )
    {
        rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in fluid" );
    }
    if( !tmp_range.empty() )
    {
        int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
        assert( dim > 0 && dim < 4 );
        fluidElems[m_or_s] = tmp_range.subset_by_dimension( dim );
    }

    if( debug )
        cout << "Read " << fluidElems[m_or_s].size() << " fluid elements from " << fluidSets[m_or_s].size()
             << " sets in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh." << endl;

    return rval;
}

ErrorCode DeformMeshRemap::find_other_sets( int m_or_s, EntityHandle file_set )
{
    // Solid or fluid sets are missing; find the other
    Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;

    if( fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty() )
    {
        unfilled_sets  = &fluidSets[m_or_s];
        filled_sets    = &solidSets[m_or_s];
        unfilled_elems = &fluidElems[m_or_s];
        if( debug )
            cout << "Finding unspecified fluid elements in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh...";
    }
    else if( !fluidSets[m_or_s].empty() && solidSets[m_or_s].empty() )
    {
        filled_sets    = &fluidSets[m_or_s];
        unfilled_sets  = &solidSets[m_or_s];
        unfilled_elems = &solidElems[m_or_s];
        if( debug )
            cout << "Finding unspecified solid elements in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh...";
    }

    // Ok, we know the filled sets, now fill the unfilled sets, and the elems from those
    Tag tagh;
    ErrorCode rval = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, tagh );MB_CHK_SET_ERR( rval, "Couldn't get material set tag name" );
    Range matsets;
    rval = mbImpl->get_entities_by_type_and_tag( file_set, MBENTITYSET, &tagh, NULL, 1, matsets );
    if( MB_SUCCESS != rval || matsets.empty() )
    {
        MB_SET_ERR( MB_FAILURE, "Couldn't get any material sets" );
    }
    *unfilled_sets = subtract( matsets, *filled_sets );
    if( unfilled_sets->empty() )
    {
        MB_SET_ERR( MB_FAILURE, "Failed to find any unfilled material sets" );
    }
    Range tmp_range;
    for( Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); ++rit )
    {
        rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in unfilled set" );
    }
    int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
    assert( dim > 0 && dim < 4 );
    *unfilled_elems = tmp_range.subset_by_dimension( dim );
    if( unfilled_elems->empty() )
    {
        MB_SET_ERR( MB_FAILURE, "Failed to find any unfilled set entities" );
    }

    if( debug )
        cout << "found " << unfilled_sets->size() << " sets and " << unfilled_elems->size() << " elements." << endl;

    return MB_SUCCESS;
}