MOAB: Mesh Oriented datABase  (version 5.4.0)
DeformMeshRemap.cpp
Go to the documentation of this file.
00001 /** @example DeformMeshRemap.cpp
00002  * Description: Account for mesh deformation of a solid due to structural mechanics\n
00003  * In this example there are two meshes, a "master" and "slave" mesh.  In the master mesh,
00004  * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
00005  * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
00006  * remapped to those new positions.  Then mesh positions and state variables are transferred
00007  * to the slave mesh, mimicing another mesh used by some other physics.
00008  *
00009  * To run: ./DeformMeshRemap [<master_meshfile> <slave_meshfile>]\n
00010  * (default values can run if users don't specify the mesh files)
00011  */
00012 
00013 #include "moab/Core.hpp"
00014 #include "moab/Range.hpp"
00015 #include "moab/Skinner.hpp"
00016 #include "moab/LloydSmoother.hpp"
00017 #include "moab/ProgOptions.hpp"
00018 #include "moab/BoundBox.hpp"
00019 #include "moab/SpatialLocator.hpp"
00020 #include "MBTagConventions.hpp"
00021 #include "DataCoupler.hpp"
00022 
00023 #ifdef USE_MPI
00024 #include "moab/ParallelComm.hpp"
00025 #endif
00026 
00027 #include <iostream>
00028 #include <set>
00029 #include <sstream>
00030 #include <cassert>
00031 
00032 using namespace moab;
00033 using namespace std;
00034 
00035 #ifndef MESH_DIR
00036 #define MESH_DIR "."
00037 #endif
00038 
00039 ErrorCode read_file( string& fname,
00040                      EntityHandle& seth,
00041                      Range& solids,
00042                      Range& solid_elems,
00043                      Range& fluids,
00044                      Range& fluid_elems );
00045 void deform_func( const BoundBox& bbox, double* xold, double* xnew );
00046 ErrorCode deform_master( Range& fluid_elems, Range& solid_elems, Tag& xnew );
00047 ErrorCode smooth_master( int dim, Tag xnew, EntityHandle& master, Range& fluids );
00048 ErrorCode write_to_coords( Range& elems, Tag tagh );
00049 
00050 const int SOLID_SETNO = 100, FLUID_SETNO = 200;
00051 
00052 Interface* mb;
00053 
00054 const bool debug = true;
00055 
00056 class DeformMeshRemap
00057 {
00058   public:
00059     //! Enumerator for solid/fluid, master/slave
00060     enum
00061     {
00062         MASTER = 0,
00063         SLAVE,
00064         SOLID,
00065         FLUID
00066     };
00067 
00068     //! Constructor
00069     //! If master is NULL, the MOAB part is run in serial;
00070     //! If slave is NULL but the master isn't, the slave is copied from the master
00071     //! Create communicators using moab::ParallelComm::get_pcomm
00072     DeformMeshRemap( Interface* impl, ParallelComm* master = NULL, ParallelComm* slave = NULL );
00073 
00074     //! Destructor
00075     ~DeformMeshRemap();
00076 
00077     //! Execute the deformed mesh process
00078     ErrorCode execute();
00079 
00080     //! Add a set number
00081     ErrorCode add_set_no( int m_or_s, int fluid_or_solid, int set_no );
00082 
00083     //! Remove a set number
00084     ErrorCode remove_set_no( int m_or_s, int fluid_or_solid, int set_no );
00085 
00086     //! Get the set numbers
00087     ErrorCode get_set_nos( int m_or_s, int fluid_or_solid, set< int >& set_nos ) const;
00088 
00089     //! Get the xNew tag handle
00090     inline Tag x_new() const
00091     {
00092         return xNew;
00093     }
00094 
00095     //! Get the tag name
00096     string x_new_name() const
00097     {
00098         return xNewName;
00099     }
00100 
00101     //! Set the tag name
00102     void x_new_name( const string& name )
00103     {
00104         xNewName = name;
00105     }
00106 
00107     //! Get/set the file name
00108     string get_file_name( int m_or_s ) const;
00109 
00110     //! Get/set the file name
00111     void set_file_name( int m_or_s, const string& name );
00112 
00113     //! Get/set the x displacement tag names
00114     string xdisp_name( int idx = 0 );
00115     void xdisp_name( const string& nm, int idx = 0 );
00116 
00117   private:
00118     //! Apply a known deformation to the solid elements, putting the results in the xNew tag; also
00119     //! write current coordinates to the xNew tag for fluid elements
00120     ErrorCode deform_master( Range& fluid_elems, Range& solid_elems, const char* tag_name = NULL );
00121 
00122     //! Read a file and establish proper ranges
00123     ErrorCode read_file( int m_or_s, string& fname, EntityHandle& seth );
00124 
00125     //! Write the input tag to the coordinates for the vertices in the input elems
00126     //! If non-zero tmp_tag is input, save coords to tmp_tag before over-writing with tag value
00127     ErrorCode write_to_coords( Range& elems, Tag tagh, Tag tmp_tag = 0 );
00128 
00129     //! Write the tag to the vertices, then save to the specified file
00130     //! If restore_coords is true, coords are restored to their initial state after file is written
00131     ErrorCode write_and_save( Range& ents,
00132                               EntityHandle seth,
00133                               Tag tagh,
00134                               const char* filename,
00135                               bool restore_coords = false );
00136 
00137     //! Find fluid/solid sets from complement of solid/fluid sets
00138     ErrorCode find_other_sets( int m_or_s, EntityHandle file_set );
00139 
00140     //! moab interface
00141     Interface* mbImpl;
00142 
00143 #ifdef USE_MPI
00144     //! ParallelComm for master, slave meshes
00145     ParallelComm *pcMaster, *pcSlave;
00146 #endif
00147 
00148     //! Material set numbers for fluid materials, for master/slave
00149     set< int > fluidSetNos[2];
00150 
00151     //! Material set numbers for solid materials, for master/slave
00152     set< int > solidSetNos[2];
00153 
00154     //! Sets defining master/slave meshes
00155     EntityHandle masterSet, slaveSet;
00156 
00157     //! Sets in master/slave meshes
00158     Range fluidSets[2], solidSets[2];
00159 
00160     //! Elements in master/slave meshes
00161     Range fluidElems[2], solidElems[2];
00162 
00163     //! Filenames for master/slave meshes
00164     string masterFileName, slaveFileName;
00165 
00166     //! Tag from file, might be 3
00167     Tag xDisp[3];
00168 
00169     //! Tag used for new positions
00170     Tag xNew;
00171 
00172     //! Tag name used to read disps from file
00173     string xDispNames[3];
00174 
00175     //! Tag name used for new positions
00176     string xNewName;
00177 };
00178 
00179 //! Add a set number
00180 inline ErrorCode DeformMeshRemap::add_set_no( int m_or_s, int f_or_s, int set_no )
00181 {
00182     set< int >* this_set;
00183     assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
00184     if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
00185 
00186     switch( f_or_s )
00187     {
00188         case FLUID:
00189             this_set = &fluidSetNos[m_or_s];
00190             break;
00191         case SOLID:
00192             this_set = &solidSetNos[m_or_s];
00193             break;
00194         default:
00195             assert( false && "f_or_s should be FLUID or SOLID." );
00196             return MB_FAILURE;
00197     }
00198 
00199     this_set->insert( set_no );
00200 
00201     return MB_SUCCESS;
00202 }
00203 
00204 //! Remove a set number
00205 inline ErrorCode DeformMeshRemap::remove_set_no( int m_or_s, int f_or_s, int set_no )
00206 {
00207     set< int >* this_set;
00208     assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
00209     if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
00210     switch( f_or_s )
00211     {
00212         case FLUID:
00213             this_set = &fluidSetNos[m_or_s];
00214             break;
00215         case SOLID:
00216             this_set = &solidSetNos[m_or_s];
00217             break;
00218         default:
00219             assert( false && "f_or_s should be FLUID or SOLID." );
00220             return MB_FAILURE;
00221     }
00222     set< int >::iterator sit = this_set->find( set_no );
00223     if( sit != this_set->end() )
00224     {
00225         this_set->erase( *sit );
00226         return MB_SUCCESS;
00227     }
00228 
00229     return MB_FAILURE;
00230 }
00231 
00232 //! Get the set numbers
00233 inline ErrorCode DeformMeshRemap::get_set_nos( int m_or_s, int f_or_s, set< int >& set_nos ) const
00234 {
00235     const set< int >* this_set;
00236     assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
00237     if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
00238     switch( f_or_s )
00239     {
00240         case FLUID:
00241             this_set = &fluidSetNos[m_or_s];
00242             break;
00243         case SOLID:
00244             this_set = &solidSetNos[m_or_s];
00245             break;
00246         default:
00247             assert( false && "f_or_s should be FLUID or SOLID." );
00248             return MB_FAILURE;
00249     }
00250 
00251     set_nos = *this_set;
00252 
00253     return MB_SUCCESS;
00254 }
00255 
00256 inline string DeformMeshRemap::xdisp_name( int idx )
00257 {
00258     return xDispNames[idx];
00259 }
00260 
00261 void DeformMeshRemap::xdisp_name( const string& nm, int idx )
00262 {
00263     xDispNames[idx] = nm;
00264 }
00265 
00266 ErrorCode DeformMeshRemap::execute()
00267 {
00268     // Read master/slave files and get fluid/solid material sets
00269     ErrorCode rval = read_file( MASTER, masterFileName, masterSet );MB_CHK_ERR( rval );
00270 
00271     if( solidSetNos[MASTER].empty() || fluidSetNos[MASTER].empty() )
00272     {
00273         rval = find_other_sets( MASTER, masterSet );MB_CHK_SET_ERR( rval, "Failed to find other sets in master mesh" );
00274     }
00275 
00276     bool have_slave = !( slaveFileName == "none" );
00277     if( have_slave )
00278     {
00279         rval = read_file( SLAVE, slaveFileName, slaveSet );MB_CHK_ERR( rval );
00280 
00281         if( solidSetNos[SLAVE].empty() || fluidSetNos[SLAVE].empty() )
00282         {
00283             rval = find_other_sets( SLAVE, slaveSet );MB_CHK_SET_ERR( rval, "Failed to find other sets in slave mesh" );
00284         }
00285     }
00286 
00287     if( debug ) cout << "Constructing data coupler/search tree on master mesh..." << endl;
00288 
00289     Range src_elems = solidElems[MASTER];
00290     src_elems.merge( fluidElems[MASTER] );
00291 
00292     // Initialize data coupler on source elements
00293     DataCoupler dc_master( mbImpl, src_elems, 0, NULL );
00294 
00295     Range tgt_verts;
00296     if( have_slave )
00297     {
00298         // Locate slave vertices in master, orig coords; do this with a data coupler, so you can
00299         // later interpolate
00300         Range tmp_range = solidElems[SLAVE];
00301         tmp_range.merge( fluidElems[SLAVE] );
00302         rval = mbImpl->get_adjacencies( tmp_range, 0, false, tgt_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get target verts" );
00303 
00304         // Locate slave vertices, caching results in dc
00305         if( debug ) cout << "Locating slave vertices in master mesh..." << endl;
00306         rval = dc_master.locate_points( tgt_verts );MB_CHK_SET_ERR( rval, "Point location of tgt verts failed" );
00307         int num_located = dc_master.spatial_locator()->local_num_located();
00308         if( num_located != (int)tgt_verts.size() )
00309         {
00310             rval = MB_FAILURE;
00311             cout << "Only " << num_located << " out of " << tgt_verts.size() << " target points successfully located."
00312                  << endl;
00313             return rval;
00314         }
00315     }
00316 
00317     // Deform the master's solid mesh, put results in a new tag
00318     if( debug ) cout << "Deforming fluid elements in master mesh..." << endl;
00319     rval = deform_master( fluidElems[MASTER], solidElems[MASTER], "xnew" );MB_CHK_ERR( rval );
00320 
00321     {  // To isolate the lloyd smoother & delete when done
00322         if( debug )
00323         {
00324             // Output the skin of smoothed elems, as a check
00325             // Get the skin; get facets, because we might need to filter on shared entities
00326             Skinner skinner( mbImpl );
00327             Range skin;
00328             rval = skinner.find_skin( 0, fluidElems[MASTER], false, skin );MB_CHK_SET_ERR( rval, "Unable to find skin" );
00329             EntityHandle skin_set;
00330             cout << "Writing skin_mesh.g and fluid_mesh.g." << endl;
00331             rval = mbImpl->create_meshset( MESHSET_SET, skin_set );MB_CHK_SET_ERR( rval, "Failed to create skin set" );
00332             rval = mbImpl->add_entities( skin_set, skin );MB_CHK_SET_ERR( rval, "Failed to add skin entities to set" );
00333             rval = mbImpl->write_file( "skin_mesh.vtk", NULL, NULL, &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failure to write skin set" );
00334             rval = mbImpl->remove_entities( skin_set, skin );MB_CHK_SET_ERR( rval, "Failed to remove skin entities from set" );
00335             rval = mbImpl->add_entities( skin_set, fluidElems[MASTER] );MB_CHK_SET_ERR( rval, "Failed to add fluid entities to set" );
00336             rval = mbImpl->write_file( "fluid_mesh.vtk", NULL, NULL, &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failure to write fluid set" );
00337             rval = mbImpl->delete_entities( &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failed to delete skin set" );
00338         }
00339 
00340         // Smooth the master mesh
00341         if( debug ) cout << "Smoothing fluid elements in master mesh..." << endl;
00342         LloydSmoother ll( mbImpl, NULL, fluidElems[MASTER], xNew );
00343         rval = ll.perform_smooth();MB_CHK_SET_ERR( rval, "Failed in lloyd smoothing" );
00344         cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
00345     }
00346 
00347     // Transfer xNew to coords, for master
00348     if( debug ) cout << "Transferring coords tag to vertex coordinates in master mesh..." << endl;
00349     rval = write_to_coords( solidElems[MASTER], xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to master fluid verts" );
00350     rval = write_to_coords( fluidElems[MASTER], xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to master fluid verts" );
00351 
00352     if( have_slave )
00353     {
00354         // Map new locations to slave
00355         // Interpolate xNew to slave points
00356         if( debug ) cout << "Interpolating new coordinates to slave vertices..." << endl;
00357         rval = dc_master.interpolate( (int)DataCoupler::VOLUME, "xnew" );MB_CHK_SET_ERR( rval, "Failed to interpolate target solution" );
00358         // Transfer xNew to coords, for slave
00359         if( debug ) cout << "Transferring coords tag to vertex coordinates in slave mesh..." << endl;
00360         rval = write_to_coords( tgt_verts, xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to slave verts" );
00361     }
00362 
00363     if( debug )
00364     {
00365         string str;
00366 #ifdef USE_MPI
00367         if( pcMaster && pcMaster->size() > 1 ) str = "PARALLEL=WRITE_PART";
00368 #endif
00369         if( debug ) cout << "Writing smoothed_master.h5m..." << endl;
00370         rval = mbImpl->write_file( "smoothed_master.h5m", NULL, str.c_str(), &masterSet, 1 );
00371 
00372         if( have_slave )
00373         {
00374 #ifdef USE_MPI
00375             str.clear();
00376             if( pcSlave && pcSlave->size() > 1 ) str = "PARALLEL=WRITE_PART";
00377 #endif
00378             if( debug ) cout << "Writing slave_interp.h5m..." << endl;
00379             rval = mbImpl->write_file( "slave_interp.h5m", NULL, str.c_str(), &slaveSet, 1 );
00380         }  // if have_slave
00381     }      // if debug
00382 
00383     if( debug ) dc_master.spatial_locator()->get_tree()->tree_stats().print();
00384 
00385     return MB_SUCCESS;
00386 }
00387 
00388 string DeformMeshRemap::get_file_name( int m_or_s ) const
00389 {
00390     switch( m_or_s )
00391     {
00392         case MASTER:
00393             return masterFileName;
00394         case SLAVE:
00395             return slaveFileName;
00396         default:
00397             assert( false && "m_or_s should be MASTER or SLAVE." );
00398             return string();
00399     }
00400 }
00401 
00402 void DeformMeshRemap::set_file_name( int m_or_s, const string& name )
00403 {
00404     switch( m_or_s )
00405     {
00406         case MASTER:
00407             masterFileName = name;
00408             break;
00409         case SLAVE:
00410             slaveFileName = name;
00411             break;
00412         default:
00413             assert( false && "m_or_s should be MASTER or SLAVE." );
00414     }
00415 }
00416 
00417 DeformMeshRemap::DeformMeshRemap( Interface* impl, ParallelComm* master, ParallelComm* slave )
00418     : mbImpl( impl ), pcMaster( master ), pcSlave( slave ), masterSet( 0 ), slaveSet( 0 ), xNew( 0 ), xNewName( "xnew" )
00419 {
00420     xDisp[0] = xDisp[1] = xDisp[2] = 0;
00421 
00422     if( !pcSlave && pcMaster ) pcSlave = pcMaster;
00423 }
00424 
00425 DeformMeshRemap::~DeformMeshRemap() {}
00426 
00427 int main( int argc, char** argv )
00428 {
00429     ErrorCode rval;
00430 
00431     ProgOptions po( "Deformed mesh options" );
00432     po.addOpt< string >( "master,m", "Specify the master meshfile name" );
00433     po.addOpt< string >( "worker,w", "Specify the slave/worker meshfile name, or 'none' (no quotes) if master only" );
00434     po.addOpt< string >( "d1,", "Tag name for displacement x or xyz" );
00435     po.addOpt< string >( "d2,", "Tag name for displacement y" );
00436     po.addOpt< string >( "d3,", "Tag name for displacement z" );
00437     po.addOpt< int >( "fm,", "Specify master fluid material set number(s). If none specified, "
00438                              "fluid sets derived from complement of solid sets." );
00439     po.addOpt< int >( "fs,", "Specify master solid material set number(s). If none specified, "
00440                              "solid sets derived from complement of fluid sets." );
00441     po.addOpt< int >( "sm,", "Specify slave fluid material set number(s). If none specified, fluid "
00442                              "sets derived from complement of solid sets." );
00443     po.addOpt< int >( "ss,", "Specify slave solid material set number(s). If none specified, solid "
00444                              "sets derived from complement of fluid sets." );
00445 
00446     po.parseCommandLine( argc, argv );
00447 
00448     mb = new Core();
00449 
00450     DeformMeshRemap* dfr;
00451 #ifdef USE_MPI
00452     ParallelComm* pc = new ParallelComm( mb, MPI_COMM_WORLD );
00453     dfr              = new DeformMeshRemap( mb, pc );
00454 #else
00455     dfr = new DeformMeshRemap( mb );
00456 #endif
00457 
00458     string masterf, slavef;
00459     if( !po.getOpt( "master", &masterf ) ) masterf = string( MESH_DIR ) + string( "/rodquad.g" );
00460     dfr->set_file_name( DeformMeshRemap::MASTER, masterf );
00461 
00462     if( !po.getOpt( "worker", &slavef ) ) slavef = string( MESH_DIR ) + string( "/rodtri.g" );
00463     dfr->set_file_name( DeformMeshRemap::SLAVE, slavef );
00464     if( slavef.empty() )
00465     {
00466         cerr << "Empty slave file name; if no slave, use filename 'none' (no quotes)." << endl;
00467         return 1;
00468     }
00469 
00470     vector< int > set_nos;
00471     po.getOptAllArgs( "fm", set_nos );
00472     for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
00473         dfr->add_set_no( DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, *vit );
00474     set_nos.clear();
00475 
00476     po.getOptAllArgs( "fs", set_nos );
00477     for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
00478         dfr->add_set_no( DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, *vit );
00479     set_nos.clear();
00480 
00481     po.getOptAllArgs( "sm", set_nos );
00482     for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
00483         dfr->add_set_no( DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, *vit );
00484 
00485     po.getOptAllArgs( "ss", set_nos );
00486     for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
00487         dfr->add_set_no( DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, *vit );
00488 
00489     string tnames[3];
00490     po.getOpt( "d1", &tnames[0] );
00491     po.getOpt( "d2", &tnames[1] );
00492     po.getOpt( "d3", &tnames[2] );
00493     for( int i = 0; i < 3; i++ )
00494         if( !tnames[i].empty() ) dfr->xdisp_name( tnames[i], i );
00495 
00496     rval = dfr->execute();
00497 
00498     delete dfr;
00499     delete mb;
00500 
00501     return rval;
00502 }
00503 
00504 ErrorCode DeformMeshRemap::write_and_save( Range& ents,
00505                                            EntityHandle seth,
00506                                            Tag tagh,
00507                                            const char* filename,
00508                                            bool restore_coords )
00509 {
00510     Tag tmp_tag = 0;
00511     ErrorCode rval;
00512     if( restore_coords ) rval = mbImpl->tag_get_handle( "", 3, MB_TYPE_DOUBLE, tmp_tag, MB_TAG_CREAT | MB_TAG_DENSE );
00513 
00514     rval = write_to_coords( ents, tagh, tmp_tag );MB_CHK_ERR( rval );
00515     rval = mbImpl->write_file( filename, NULL, NULL, &seth, 1 );MB_CHK_ERR( rval );
00516     if( restore_coords )
00517     {
00518         rval = write_to_coords( ents, tmp_tag );MB_CHK_ERR( rval );
00519         rval = mbImpl->tag_delete( tmp_tag );MB_CHK_ERR( rval );
00520     }
00521 
00522     return rval;
00523 }
00524 
00525 ErrorCode DeformMeshRemap::write_to_coords( Range& elems, Tag tagh, Tag tmp_tag )
00526 {
00527     // Write the tag to coordinates
00528     Range verts;
00529     ErrorCode rval = mbImpl->get_adjacencies( elems, 0, false, verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get adj vertices" );
00530     vector< double > coords( 3 * verts.size() );
00531 
00532     if( tmp_tag )
00533     {
00534         // Save the coords to tmp_tag first
00535         rval = mbImpl->get_coords( verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get tmp copy of coords" );
00536         rval = mbImpl->tag_set_data( tmp_tag, verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to save tmp copy of coords" );
00537     }
00538 
00539     rval = mbImpl->tag_get_data( tagh, verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get tag data" );
00540     rval = mbImpl->set_coords( verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to set coordinates" );
00541     return MB_SUCCESS;
00542 }
00543 
00544 void deform_func( const BoundBox& bbox, double* xold, double* xnew )
00545 {
00546     /*  Deformation function based on max delx and dely at top of rod
00547         const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
00548         // function: origin is at middle base of rod, and is .5 high
00549         // top of rod is (0,.55) on left and (.2,.6) on right
00550       double delx = 0.5*RODWIDTH;
00551 
00552       double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
00553       xnew[0] = xold[0] + yfrac * delx;
00554       xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
00555     */
00556     /* Deformation function based on fraction of bounding box dimension in each direction */
00557     double frac  = 0.01;  // taken from approximate relative deformation from LLNL Diablo of XX09 assys
00558     CartVect *xo = reinterpret_cast< CartVect* >( xold ), *xn = reinterpret_cast< CartVect* >( xnew );
00559     CartVect disp = frac * ( *xo - bbox.bMin );
00560     *xn           = *xo + disp;
00561 }
00562 
00563 ErrorCode DeformMeshRemap::deform_master( Range& fluid_elems, Range& solid_elems, const char* tag_name )
00564 {
00565     // Deform elements with an analytic function
00566     ErrorCode rval;
00567 
00568     // Get all the vertices and coords in the solid
00569     Range solid_verts, fluid_verts;
00570     rval = mbImpl->get_adjacencies( solid_elems, 0, false, solid_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get vertices" );
00571     vector< double > coords( 3 * solid_verts.size() ), new_coords( 3 * solid_verts.size() );
00572     rval = mbImpl->get_coords( solid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get vertex coords" );
00573     unsigned int num_verts = solid_verts.size();
00574 
00575     // Get or create the tag
00576     if( !xDispNames[0].empty() && !xDispNames[1].empty() && !xDispNames[2].empty() )
00577     {
00578         // 3 tags, specifying xyz individual data, integrate into one tag
00579         rval = mbImpl->tag_get_handle( ( tag_name ? tag_name : "" ), 3, MB_TYPE_DOUBLE, xNew,
00580                                        MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Failed to create xnew tag" );
00581         vector< double > disps( num_verts );
00582         for( int i = 0; i < 3; i++ )
00583         {
00584             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" );
00585             rval = mbImpl->tag_get_data( xDisp[i], solid_verts, &disps[0] );MB_CHK_SET_ERR( rval, "Failed to get xDisp tag values" );
00586             for( unsigned int j = 0; j < num_verts; j++ )
00587                 new_coords[3 * j + i] = coords[3 * j + i] + disps[j];
00588         }
00589     }
00590     else if( !xDispNames[0].empty() )
00591     {
00592         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" );
00593         xNew = xDisp[0];
00594         vector< double > disps( 3 * num_verts );
00595         rval = mbImpl->tag_get_data( xDisp[0], solid_verts, &disps[0] );
00596         for( unsigned int j = 0; j < 3 * num_verts; j++ )
00597             new_coords[j] = coords[j] + disps[j];
00598     }
00599     else
00600     {
00601         // Get the bounding box of the solid mesh
00602         BoundBox bbox;
00603         bbox.update( *mbImpl, solid_elems );
00604 
00605         for( unsigned int j = 0; j < num_verts; j++ )
00606             deform_func( bbox, &coords[3 * j], &new_coords[3 * j] );
00607     }
00608 
00609     if( debug )
00610     {
00611         double len = 0.0;
00612         for( unsigned int i = 0; i < num_verts; i++ )
00613         {
00614             CartVect dx    = CartVect( &new_coords[3 * i] ) - CartVect( &coords[3 * i] );
00615             double tmp_len = dx.length_squared();
00616             if( tmp_len > len ) len = tmp_len;
00617         }
00618         Range tmp_elems( fluid_elems );
00619         tmp_elems.merge( solid_elems );
00620         BoundBox box;
00621         box.update( *mbImpl, tmp_elems );
00622         double max_len =
00623             std::max( box.bMax[2] - box.bMin[2], std::max( box.bMax[1] - box.bMin[1], box.bMax[0] - box.bMin[0] ) );
00624 
00625         cout << "Max displacement = " << len << " (" << 100.0 * len / max_len << "% of max box length)" << endl;
00626     }
00627 
00628     if( !xNew )
00629     {
00630         rval = mbImpl->tag_get_handle( ( tag_name ? tag_name : "" ), 3, MB_TYPE_DOUBLE, xDisp[0],
00631                                        MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Failed to get xNew tag" );
00632         xNew = xDisp[0];
00633     }
00634 
00635     // Set the new tag to those coords
00636     rval = mbImpl->tag_set_data( xNew, solid_verts, &new_coords[0] );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
00637 
00638     // Get all the vertices and coords in the fluid, set xnew to them
00639     rval = mbImpl->get_adjacencies( fluid_elems, 0, false, fluid_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get vertices" );
00640     fluid_verts = subtract( fluid_verts, solid_verts );
00641 
00642     if( coords.size() < 3 * fluid_verts.size() ) coords.resize( 3 * fluid_verts.size() );
00643     rval = mbImpl->get_coords( fluid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get vertex coords" );
00644     rval = mbImpl->tag_set_data( xNew, fluid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to set xnew tag on fluid verts" );
00645 
00646     if( debug )
00647     {
00648         // Save deformed mesh coords to new file for visualizing
00649         Range tmp_range( fluidElems[MASTER] );
00650         tmp_range.merge( solidElems[MASTER] );
00651         rval = write_and_save( tmp_range, masterSet, xNew, "deformed_master.h5m", true );MB_CHK_ERR( rval );
00652     }
00653 
00654     return MB_SUCCESS;
00655 }
00656 
00657 ErrorCode DeformMeshRemap::read_file( int m_or_s, string& fname, EntityHandle& seth )
00658 {
00659     // Create meshset
00660     ErrorCode rval = mbImpl->create_meshset( 0, seth );MB_CHK_SET_ERR( rval, "Couldn't create master/slave set" );
00661     ostringstream options;
00662 #ifdef USE_MPI
00663     ParallelComm* pc = ( m_or_s == MASTER ? pcMaster : pcSlave );
00664     if( pc && pc->size() > 1 )
00665     {
00666         if( debug ) options << "DEBUG_IO=1;CPUTIME;";
00667         options << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
00668                 << "PARALLEL_GHOSTS=2.0.1;PARALLEL_COMM=" << pc->get_id();
00669     }
00670 #endif
00671     rval = mbImpl->load_file( fname.c_str(), &seth, options.str().c_str() );MB_CHK_SET_ERR( rval, "Couldn't load master/slave mesh" );
00672 
00673     if( *solidSetNos[m_or_s].begin() == -1 || *fluidSetNos[m_or_s].begin() == -1 ) return MB_SUCCESS;
00674 
00675     // Get material sets for solid/fluid
00676     Tag tagh;
00677     rval = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, tagh );MB_CHK_SET_ERR( rval, "Couldn't get material set tag name" );
00678     for( set< int >::iterator sit = solidSetNos[m_or_s].begin(); sit != solidSetNos[m_or_s].end(); ++sit )
00679     {
00680         Range sets;
00681         int set_no            = *sit;
00682         const void* setno_ptr = &set_no;
00683         rval                  = mbImpl->get_entities_by_type_and_tag( seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets );
00684         if( MB_SUCCESS != rval || sets.empty() )
00685         {
00686             MB_SET_ERR( MB_FAILURE, "Couldn't find solid set #" << *sit );
00687         }
00688         else
00689             solidSets[m_or_s].merge( sets );
00690     }
00691 
00692     // Get solid entities, and dimension
00693     Range tmp_range;
00694     for( Range::iterator rit = solidSets[m_or_s].begin(); rit != solidSets[m_or_s].end(); ++rit )
00695     {
00696         rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in solid" );
00697     }
00698     if( !tmp_range.empty() )
00699     {
00700         int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
00701         assert( dim > 0 && dim < 4 );
00702         solidElems[m_or_s] = tmp_range.subset_by_dimension( dim );
00703     }
00704 
00705     if( debug )
00706         cout << "Read " << solidElems[m_or_s].size() << " solid elements from " << solidSets[m_or_s].size()
00707              << " sets in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh." << endl;
00708 
00709     for( set< int >::iterator sit = fluidSetNos[m_or_s].begin(); sit != fluidSetNos[m_or_s].end(); ++sit )
00710     {
00711         Range sets;
00712         int set_no            = *sit;
00713         const void* setno_ptr = &set_no;
00714         rval                  = mbImpl->get_entities_by_type_and_tag( seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets );
00715         if( MB_SUCCESS != rval || sets.empty() )
00716         {
00717             MB_SET_ERR( MB_FAILURE, "Couldn't find fluid set #" << *sit );
00718         }
00719         else
00720             fluidSets[m_or_s].merge( sets );
00721     }
00722 
00723     // Get fluid entities, and dimension
00724     tmp_range.clear();
00725     for( Range::iterator rit = fluidSets[m_or_s].begin(); rit != fluidSets[m_or_s].end(); ++rit )
00726     {
00727         rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in fluid" );
00728     }
00729     if( !tmp_range.empty() )
00730     {
00731         int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
00732         assert( dim > 0 && dim < 4 );
00733         fluidElems[m_or_s] = tmp_range.subset_by_dimension( dim );
00734     }
00735 
00736     if( debug )
00737         cout << "Read " << fluidElems[m_or_s].size() << " fluid elements from " << fluidSets[m_or_s].size()
00738              << " sets in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh." << endl;
00739 
00740     return rval;
00741 }
00742 
00743 ErrorCode DeformMeshRemap::find_other_sets( int m_or_s, EntityHandle file_set )
00744 {
00745     // Solid or fluid sets are missing; find the other
00746     Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
00747 
00748     if( fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty() )
00749     {
00750         unfilled_sets  = &fluidSets[m_or_s];
00751         filled_sets    = &solidSets[m_or_s];
00752         unfilled_elems = &fluidElems[m_or_s];
00753         if( debug )
00754             cout << "Finding unspecified fluid elements in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh...";
00755     }
00756     else if( !fluidSets[m_or_s].empty() && solidSets[m_or_s].empty() )
00757     {
00758         filled_sets    = &fluidSets[m_or_s];
00759         unfilled_sets  = &solidSets[m_or_s];
00760         unfilled_elems = &solidElems[m_or_s];
00761         if( debug )
00762             cout << "Finding unspecified solid elements in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh...";
00763     }
00764 
00765     // Ok, we know the filled sets, now fill the unfilled sets, and the elems from those
00766     Tag tagh;
00767     ErrorCode rval = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, tagh );MB_CHK_SET_ERR( rval, "Couldn't get material set tag name" );
00768     Range matsets;
00769     rval = mbImpl->get_entities_by_type_and_tag( file_set, MBENTITYSET, &tagh, NULL, 1, matsets );
00770     if( MB_SUCCESS != rval || matsets.empty() )
00771     {
00772         MB_SET_ERR( MB_FAILURE, "Couldn't get any material sets" );
00773     }
00774     *unfilled_sets = subtract( matsets, *filled_sets );
00775     if( unfilled_sets->empty() )
00776     {
00777         MB_SET_ERR( MB_FAILURE, "Failed to find any unfilled material sets" );
00778     }
00779     Range tmp_range;
00780     for( Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); ++rit )
00781     {
00782         rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in unfilled set" );
00783     }
00784     int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
00785     assert( dim > 0 && dim < 4 );
00786     *unfilled_elems = tmp_range.subset_by_dimension( dim );
00787     if( unfilled_elems->empty() )
00788     {
00789         MB_SET_ERR( MB_FAILURE, "Failed to find any unfilled set entities" );
00790     }
00791 
00792     if( debug )
00793         cout << "found " << unfilled_sets->size() << " sets and " << unfilled_elems->size() << " elements." << endl;
00794 
00795     return MB_SUCCESS;
00796 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines