DeformMeshRemap.cpp

Go to the documentation of this file.

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