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