Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 "DataCoupler.hpp"
#ifdef USE_MPI
#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;
{
public:
//! Enumerator for solid/fluid, master/slave
enum
{
MASTER = 0,
SLAVE,
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
//! 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,
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;
}
{
// 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;
}
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();
#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" );
if( !po.getOpt( "worker", &slavef ) ) slavef = string( MESH_DIR ) + string( "/rodtri.g" );
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 )
set_nos.clear();
po.getOptAllArgs( "fs", set_nos );
for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
set_nos.clear();
po.getOptAllArgs( "sm", set_nos );
for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
po.getOptAllArgs( "ss", set_nos );
for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++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;
}
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;
}
{
// 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;
}