Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
DeformMeshRemap.cpp
Go to the documentation of this file.
1 /** @example DeformMeshRemap.cpp
2  * Description: Account for mesh deformation of a solid due to structural mechanics\n
3  * In this example there are two meshes, a "master" and "slave" mesh. In the master mesh,
4  * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
5  * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
6  * remapped to those new positions. Then mesh positions and state variables are transferred
7  * to the slave mesh, mimicing another mesh used by some other physics.
8  *
9  * To run: ./DeformMeshRemap [<master_meshfile> <slave_meshfile>]\n
10  * (default values can run if users don't specify the mesh files)
11  */
12 
13 #include "moab/Core.hpp"
14 #include "moab/Range.hpp"
15 #include "moab/Skinner.hpp"
16 #include "moab/LloydSmoother.hpp"
17 #include "moab/ProgOptions.hpp"
18 #include "moab/BoundBox.hpp"
19 #include "moab/SpatialLocator.hpp"
20 #include "MBTagConventions.hpp"
21 #include "DataCoupler.hpp"
22 
23 #ifdef USE_MPI
24 #include "moab/ParallelComm.hpp"
25 #endif
26 
27 #include <iostream>
28 #include <set>
29 #include <sstream>
30 #include <cassert>
31 
32 using namespace moab;
33 using namespace std;
34 
35 #ifndef MESH_DIR
36 #define MESH_DIR "."
37 #endif
38 
39 ErrorCode read_file( string& fname,
40  EntityHandle& seth,
41  Range& solids,
42  Range& solid_elems,
43  Range& fluids,
44  Range& fluid_elems );
45 void deform_func( const BoundBox& bbox, double* xold, double* xnew );
46 ErrorCode deform_master( Range& fluid_elems, Range& solid_elems, Tag& xnew );
47 ErrorCode smooth_master( int dim, Tag xnew, EntityHandle& master, Range& fluids );
49 
50 const int SOLID_SETNO = 100, FLUID_SETNO = 200;
51 
53 
54 const bool debug = true;
55 
57 {
58  public:
59  //! Enumerator for solid/fluid, master/slave
60  enum
61  {
62  MASTER = 0,
65  FLUID
66  };
67 
68  //! Constructor
69  //! If master is NULL, the MOAB part is run in serial;
70  //! If slave is NULL but the master isn't, the slave is copied from the master
71  //! Create communicators using moab::ParallelComm::get_pcomm
72  DeformMeshRemap( Interface* impl, ParallelComm* master = NULL, ParallelComm* slave = NULL );
73 
74  //! Destructor
75  ~DeformMeshRemap();
76 
77  //! Execute the deformed mesh process
78  ErrorCode execute();
79 
80  //! Add a set number
81  ErrorCode add_set_no( int m_or_s, int fluid_or_solid, int set_no );
82 
83  //! Remove a set number
84  ErrorCode remove_set_no( int m_or_s, int fluid_or_solid, int set_no );
85 
86  //! Get the set numbers
87  ErrorCode get_set_nos( int m_or_s, int fluid_or_solid, set< int >& set_nos ) const;
88 
89  //! Get the xNew tag handle
90  inline Tag x_new() const
91  {
92  return xNew;
93  }
94 
95  //! Get the tag name
96  string x_new_name() const
97  {
98  return xNewName;
99  }
100 
101  //! Set the tag name
102  void x_new_name( const string& name )
103  {
104  xNewName = name;
105  }
106 
107  //! Get/set the file name
108  string get_file_name( int m_or_s ) const;
109 
110  //! Get/set the file name
111  void set_file_name( int m_or_s, const string& name );
112 
113  //! Get/set the x displacement tag names
114  string xdisp_name( int idx = 0 );
115  void xdisp_name( const string& nm, int idx = 0 );
116 
117  private:
118  //! Apply a known deformation to the solid elements, putting the results in the xNew tag; also
119  //! write current coordinates to the xNew tag for fluid elements
120  ErrorCode deform_master( Range& fluid_elems, Range& solid_elems, const char* tag_name = NULL );
121 
122  //! Read a file and establish proper ranges
123  ErrorCode read_file( int m_or_s, string& fname, EntityHandle& seth );
124 
125  //! Write the input tag to the coordinates for the vertices in the input elems
126  //! If non-zero tmp_tag is input, save coords to tmp_tag before over-writing with tag value
127  ErrorCode write_to_coords( Range& elems, Tag tagh, Tag tmp_tag = 0 );
128 
129  //! Write the tag to the vertices, then save to the specified file
130  //! If restore_coords is true, coords are restored to their initial state after file is written
131  ErrorCode write_and_save( Range& ents,
132  EntityHandle seth,
133  Tag tagh,
134  const char* filename,
135  bool restore_coords = false );
136 
137  //! Find fluid/solid sets from complement of solid/fluid sets
138  ErrorCode find_other_sets( int m_or_s, EntityHandle file_set );
139 
140  //! moab interface
142 
143 #ifdef USE_MPI
144  //! ParallelComm for master, slave meshes
145  ParallelComm *pcMaster, *pcSlave;
146 #endif
147 
148  //! Material set numbers for fluid materials, for master/slave
149  set< int > fluidSetNos[2];
150 
151  //! Material set numbers for solid materials, for master/slave
152  set< int > solidSetNos[2];
153 
154  //! Sets defining master/slave meshes
156 
157  //! Sets in master/slave meshes
158  Range fluidSets[2], solidSets[2];
159 
160  //! Elements in master/slave meshes
161  Range fluidElems[2], solidElems[2];
162 
163  //! Filenames for master/slave meshes
164  string masterFileName, slaveFileName;
165 
166  //! Tag from file, might be 3
167  Tag xDisp[3];
168 
169  //! Tag used for new positions
171 
172  //! Tag name used to read disps from file
173  string xDispNames[3];
174 
175  //! Tag name used for new positions
176  string xNewName;
177 };
178 
179 //! Add a set number
180 inline ErrorCode DeformMeshRemap::add_set_no( int m_or_s, int f_or_s, int set_no )
181 {
182  set< int >* this_set;
183  assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
184  if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
185 
186  switch( f_or_s )
187  {
188  case FLUID:
189  this_set = &fluidSetNos[m_or_s];
190  break;
191  case SOLID:
192  this_set = &solidSetNos[m_or_s];
193  break;
194  default:
195  assert( false && "f_or_s should be FLUID or SOLID." );
196  return MB_FAILURE;
197  }
198 
199  this_set->insert( set_no );
200 
201  return MB_SUCCESS;
202 }
203 
204 //! Remove a set number
205 inline ErrorCode DeformMeshRemap::remove_set_no( int m_or_s, int f_or_s, int set_no )
206 {
207  set< int >* this_set;
208  assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
209  if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
210  switch( f_or_s )
211  {
212  case FLUID:
213  this_set = &fluidSetNos[m_or_s];
214  break;
215  case SOLID:
216  this_set = &solidSetNos[m_or_s];
217  break;
218  default:
219  assert( false && "f_or_s should be FLUID or SOLID." );
220  return MB_FAILURE;
221  }
222  set< int >::iterator sit = this_set->find( set_no );
223  if( sit != this_set->end() )
224  {
225  this_set->erase( *sit );
226  return MB_SUCCESS;
227  }
228 
229  return MB_FAILURE;
230 }
231 
232 //! Get the set numbers
233 inline ErrorCode DeformMeshRemap::get_set_nos( int m_or_s, int f_or_s, set< int >& set_nos ) const
234 {
235  const set< int >* this_set;
236  assert( ( m_or_s == MASTER || m_or_s == SLAVE ) && "m_or_s should be MASTER or SLAVE." );
237  if( m_or_s != MASTER && m_or_s != SLAVE ) return MB_INDEX_OUT_OF_RANGE;
238  switch( f_or_s )
239  {
240  case FLUID:
241  this_set = &fluidSetNos[m_or_s];
242  break;
243  case SOLID:
244  this_set = &solidSetNos[m_or_s];
245  break;
246  default:
247  assert( false && "f_or_s should be FLUID or SOLID." );
248  return MB_FAILURE;
249  }
250 
251  set_nos = *this_set;
252 
253  return MB_SUCCESS;
254 }
255 
256 inline string DeformMeshRemap::xdisp_name( int idx )
257 {
258  return xDispNames[idx];
259 }
260 
261 void DeformMeshRemap::xdisp_name( const string& nm, int idx )
262 {
263  xDispNames[idx] = nm;
264 }
265 
267 {
268  // Read master/slave files and get fluid/solid material sets
269  ErrorCode rval = read_file( MASTER, masterFileName, masterSet );MB_CHK_ERR( rval );
270 
271  if( solidSetNos[MASTER].empty() || fluidSetNos[MASTER].empty() )
272  {
273  rval = find_other_sets( MASTER, masterSet );MB_CHK_SET_ERR( rval, "Failed to find other sets in master mesh" );
274  }
275 
276  bool have_slave = !( slaveFileName == "none" );
277  if( have_slave )
278  {
279  rval = read_file( SLAVE, slaveFileName, slaveSet );MB_CHK_ERR( rval );
280 
281  if( solidSetNos[SLAVE].empty() || fluidSetNos[SLAVE].empty() )
282  {
283  rval = find_other_sets( SLAVE, slaveSet );MB_CHK_SET_ERR( rval, "Failed to find other sets in slave mesh" );
284  }
285  }
286 
287  if( debug ) cout << "Constructing data coupler/search tree on master mesh..." << endl;
288 
289  Range src_elems = solidElems[MASTER];
290  src_elems.merge( fluidElems[MASTER] );
291 
292  // Initialize data coupler on source elements
293  DataCoupler dc_master( mbImpl, src_elems, 0, NULL );
294 
295  Range tgt_verts;
296  if( have_slave )
297  {
298  // Locate slave vertices in master, orig coords; do this with a data coupler, so you can
299  // later interpolate
300  Range tmp_range = solidElems[SLAVE];
301  tmp_range.merge( fluidElems[SLAVE] );
302  rval = mbImpl->get_adjacencies( tmp_range, 0, false, tgt_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get target verts" );
303 
304  // Locate slave vertices, caching results in dc
305  if( debug ) cout << "Locating slave vertices in master mesh..." << endl;
306  rval = dc_master.locate_points( tgt_verts );MB_CHK_SET_ERR( rval, "Point location of tgt verts failed" );
307  int num_located = dc_master.spatial_locator()->local_num_located();
308  if( num_located != (int)tgt_verts.size() )
309  {
310  rval = MB_FAILURE;
311  cout << "Only " << num_located << " out of " << tgt_verts.size() << " target points successfully located."
312  << endl;
313  return rval;
314  }
315  }
316 
317  // Deform the master's solid mesh, put results in a new tag
318  if( debug ) cout << "Deforming fluid elements in master mesh..." << endl;
319  rval = deform_master( fluidElems[MASTER], solidElems[MASTER], "xnew" );MB_CHK_ERR( rval );
320 
321  { // To isolate the lloyd smoother & delete when done
322  if( debug )
323  {
324  // Output the skin of smoothed elems, as a check
325  // Get the skin; get facets, because we might need to filter on shared entities
326  Skinner skinner( mbImpl );
327  Range skin;
328  rval = skinner.find_skin( 0, fluidElems[MASTER], false, skin );MB_CHK_SET_ERR( rval, "Unable to find skin" );
329  EntityHandle skin_set;
330  cout << "Writing skin_mesh.g and fluid_mesh.g." << endl;
331  rval = mbImpl->create_meshset( MESHSET_SET, skin_set );MB_CHK_SET_ERR( rval, "Failed to create skin set" );
332  rval = mbImpl->add_entities( skin_set, skin );MB_CHK_SET_ERR( rval, "Failed to add skin entities to set" );
333  rval = mbImpl->write_file( "skin_mesh.vtk", NULL, NULL, &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failure to write skin set" );
334  rval = mbImpl->remove_entities( skin_set, skin );MB_CHK_SET_ERR( rval, "Failed to remove skin entities from set" );
335  rval = mbImpl->add_entities( skin_set, fluidElems[MASTER] );MB_CHK_SET_ERR( rval, "Failed to add fluid entities to set" );
336  rval = mbImpl->write_file( "fluid_mesh.vtk", NULL, NULL, &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failure to write fluid set" );
337  rval = mbImpl->delete_entities( &skin_set, 1 );MB_CHK_SET_ERR( rval, "Failed to delete skin set" );
338  }
339 
340  // Smooth the master mesh
341  if( debug ) cout << "Smoothing fluid elements in master mesh..." << endl;
342  LloydSmoother ll( mbImpl, NULL, fluidElems[MASTER], xNew );
343  rval = ll.perform_smooth();MB_CHK_SET_ERR( rval, "Failed in lloyd smoothing" );
344  cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
345  }
346 
347  // Transfer xNew to coords, for master
348  if( debug ) cout << "Transferring coords tag to vertex coordinates in master mesh..." << endl;
349  rval = write_to_coords( solidElems[MASTER], xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to master fluid verts" );
350  rval = write_to_coords( fluidElems[MASTER], xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to master fluid verts" );
351 
352  if( have_slave )
353  {
354  // Map new locations to slave
355  // Interpolate xNew to slave points
356  if( debug ) cout << "Interpolating new coordinates to slave vertices..." << endl;
357  rval = dc_master.interpolate( (int)DataCoupler::VOLUME, "xnew" );MB_CHK_SET_ERR( rval, "Failed to interpolate target solution" );
358  // Transfer xNew to coords, for slave
359  if( debug ) cout << "Transferring coords tag to vertex coordinates in slave mesh..." << endl;
360  rval = write_to_coords( tgt_verts, xNew );MB_CHK_SET_ERR( rval, "Failed writing tag to slave verts" );
361  }
362 
363  if( debug )
364  {
365  string str;
366 #ifdef USE_MPI
367  if( pcMaster && pcMaster->size() > 1 ) str = "PARALLEL=WRITE_PART";
368 #endif
369  if( debug ) cout << "Writing smoothed_master.h5m..." << endl;
370  rval = mbImpl->write_file( "smoothed_master.h5m", NULL, str.c_str(), &masterSet, 1 );
371 
372  if( have_slave )
373  {
374 #ifdef USE_MPI
375  str.clear();
376  if( pcSlave && pcSlave->size() > 1 ) str = "PARALLEL=WRITE_PART";
377 #endif
378  if( debug ) cout << "Writing slave_interp.h5m..." << endl;
379  rval = mbImpl->write_file( "slave_interp.h5m", NULL, str.c_str(), &slaveSet, 1 );
380  } // if have_slave
381  } // if debug
382 
383  if( debug ) dc_master.spatial_locator()->get_tree()->tree_stats().print();
384 
385  return MB_SUCCESS;
386 }
387 
388 string DeformMeshRemap::get_file_name( int m_or_s ) const
389 {
390  switch( m_or_s )
391  {
392  case MASTER:
393  return masterFileName;
394  case SLAVE:
395  return slaveFileName;
396  default:
397  assert( false && "m_or_s should be MASTER or SLAVE." );
398  return string();
399  }
400 }
401 
402 void DeformMeshRemap::set_file_name( int m_or_s, const string& name )
403 {
404  switch( m_or_s )
405  {
406  case MASTER:
407  masterFileName = name;
408  break;
409  case SLAVE:
410  slaveFileName = name;
411  break;
412  default:
413  assert( false && "m_or_s should be MASTER or SLAVE." );
414  }
415 }
416 
418  : mbImpl( impl ), pcMaster( master ), pcSlave( slave ), masterSet( 0 ), slaveSet( 0 ), xNew( 0 ), xNewName( "xnew" )
419 {
420  xDisp[0] = xDisp[1] = xDisp[2] = 0;
421 
422  if( !pcSlave && pcMaster ) pcSlave = pcMaster;
423 }
424 
426 
427 int main( int argc, char** argv )
428 {
429  ErrorCode rval;
430 
431  ProgOptions po( "Deformed mesh options" );
432  po.addOpt< string >( "master,m", "Specify the master meshfile name" );
433  po.addOpt< string >( "worker,w", "Specify the slave/worker meshfile name, or 'none' (no quotes) if master only" );
434  po.addOpt< string >( "d1,", "Tag name for displacement x or xyz" );
435  po.addOpt< string >( "d2,", "Tag name for displacement y" );
436  po.addOpt< string >( "d3,", "Tag name for displacement z" );
437  po.addOpt< int >( "fm,", "Specify master fluid material set number(s). If none specified, "
438  "fluid sets derived from complement of solid sets." );
439  po.addOpt< int >( "fs,", "Specify master solid material set number(s). If none specified, "
440  "solid sets derived from complement of fluid sets." );
441  po.addOpt< int >( "sm,", "Specify slave fluid material set number(s). If none specified, fluid "
442  "sets derived from complement of solid sets." );
443  po.addOpt< int >( "ss,", "Specify slave solid material set number(s). If none specified, solid "
444  "sets derived from complement of fluid sets." );
445 
446  po.parseCommandLine( argc, argv );
447 
448  mb = new Core();
449 
450  DeformMeshRemap* dfr;
451 #ifdef USE_MPI
452  ParallelComm* pc = new ParallelComm( mb, MPI_COMM_WORLD );
453  dfr = new DeformMeshRemap( mb, pc );
454 #else
455  dfr = new DeformMeshRemap( mb );
456 #endif
457 
458  string masterf, slavef;
459  if( !po.getOpt( "master", &masterf ) ) masterf = string( MESH_DIR ) + string( "/rodquad.g" );
460  dfr->set_file_name( DeformMeshRemap::MASTER, masterf );
461 
462  if( !po.getOpt( "worker", &slavef ) ) slavef = string( MESH_DIR ) + string( "/rodtri.g" );
463  dfr->set_file_name( DeformMeshRemap::SLAVE, slavef );
464  if( slavef.empty() )
465  {
466  cerr << "Empty slave file name; if no slave, use filename 'none' (no quotes)." << endl;
467  return 1;
468  }
469 
470  vector< int > set_nos;
471  po.getOptAllArgs( "fm", set_nos );
472  for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
474  set_nos.clear();
475 
476  po.getOptAllArgs( "fs", set_nos );
477  for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
479  set_nos.clear();
480 
481  po.getOptAllArgs( "sm", set_nos );
482  for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
484 
485  po.getOptAllArgs( "ss", set_nos );
486  for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
488 
489  string tnames[3];
490  po.getOpt( "d1", &tnames[0] );
491  po.getOpt( "d2", &tnames[1] );
492  po.getOpt( "d3", &tnames[2] );
493  for( int i = 0; i < 3; i++ )
494  if( !tnames[i].empty() ) dfr->xdisp_name( tnames[i], i );
495 
496  rval = dfr->execute();
497 
498  delete dfr;
499  delete mb;
500 
501  return rval;
502 }
503 
505  EntityHandle seth,
506  Tag tagh,
507  const char* filename,
508  bool restore_coords )
509 {
510  Tag tmp_tag = 0;
511  ErrorCode rval;
512  if( restore_coords ) rval = mbImpl->tag_get_handle( "", 3, MB_TYPE_DOUBLE, tmp_tag, MB_TAG_CREAT | MB_TAG_DENSE );
513 
514  rval = write_to_coords( ents, tagh, tmp_tag );MB_CHK_ERR( rval );
515  rval = mbImpl->write_file( filename, NULL, NULL, &seth, 1 );MB_CHK_ERR( rval );
516  if( restore_coords )
517  {
518  rval = write_to_coords( ents, tmp_tag );MB_CHK_ERR( rval );
519  rval = mbImpl->tag_delete( tmp_tag );MB_CHK_ERR( rval );
520  }
521 
522  return rval;
523 }
524 
526 {
527  // Write the tag to coordinates
528  Range verts;
529  ErrorCode rval = mbImpl->get_adjacencies( elems, 0, false, verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get adj vertices" );
530  vector< double > coords( 3 * verts.size() );
531 
532  if( tmp_tag )
533  {
534  // Save the coords to tmp_tag first
535  rval = mbImpl->get_coords( verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get tmp copy of coords" );
536  rval = mbImpl->tag_set_data( tmp_tag, verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to save tmp copy of coords" );
537  }
538 
539  rval = mbImpl->tag_get_data( tagh, verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get tag data" );
540  rval = mbImpl->set_coords( verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to set coordinates" );
541  return MB_SUCCESS;
542 }
543 
544 void deform_func( const BoundBox& bbox, double* xold, double* xnew )
545 {
546  /* Deformation function based on max delx and dely at top of rod
547  const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
548  // function: origin is at middle base of rod, and is .5 high
549  // top of rod is (0,.55) on left and (.2,.6) on right
550  double delx = 0.5*RODWIDTH;
551 
552  double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
553  xnew[0] = xold[0] + yfrac * delx;
554  xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
555  */
556  /* Deformation function based on fraction of bounding box dimension in each direction */
557  double frac = 0.01; // taken from approximate relative deformation from LLNL Diablo of XX09 assys
558  CartVect *xo = reinterpret_cast< CartVect* >( xold ), *xn = reinterpret_cast< CartVect* >( xnew );
559  CartVect disp = frac * ( *xo - bbox.bMin );
560  *xn = *xo + disp;
561 }
562 
563 ErrorCode DeformMeshRemap::deform_master( Range& fluid_elems, Range& solid_elems, const char* tag_name )
564 {
565  // Deform elements with an analytic function
566  ErrorCode rval;
567 
568  // Get all the vertices and coords in the solid
569  Range solid_verts, fluid_verts;
570  rval = mbImpl->get_adjacencies( solid_elems, 0, false, solid_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get vertices" );
571  vector< double > coords( 3 * solid_verts.size() ), new_coords( 3 * solid_verts.size() );
572  rval = mbImpl->get_coords( solid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get vertex coords" );
573  unsigned int num_verts = solid_verts.size();
574 
575  // Get or create the tag
576  if( !xDispNames[0].empty() && !xDispNames[1].empty() && !xDispNames[2].empty() )
577  {
578  // 3 tags, specifying xyz individual data, integrate into one tag
579  rval = mbImpl->tag_get_handle( ( tag_name ? tag_name : "" ), 3, MB_TYPE_DOUBLE, xNew,
580  MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Failed to create xnew tag" );
581  vector< double > disps( num_verts );
582  for( int i = 0; i < 3; i++ )
583  {
584  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" );
585  rval = mbImpl->tag_get_data( xDisp[i], solid_verts, &disps[0] );MB_CHK_SET_ERR( rval, "Failed to get xDisp tag values" );
586  for( unsigned int j = 0; j < num_verts; j++ )
587  new_coords[3 * j + i] = coords[3 * j + i] + disps[j];
588  }
589  }
590  else if( !xDispNames[0].empty() )
591  {
592  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" );
593  xNew = xDisp[0];
594  vector< double > disps( 3 * num_verts );
595  rval = mbImpl->tag_get_data( xDisp[0], solid_verts, &disps[0] );
596  for( unsigned int j = 0; j < 3 * num_verts; j++ )
597  new_coords[j] = coords[j] + disps[j];
598  }
599  else
600  {
601  // Get the bounding box of the solid mesh
602  BoundBox bbox;
603  bbox.update( *mbImpl, solid_elems );
604 
605  for( unsigned int j = 0; j < num_verts; j++ )
606  deform_func( bbox, &coords[3 * j], &new_coords[3 * j] );
607  }
608 
609  if( debug )
610  {
611  double len = 0.0;
612  for( unsigned int i = 0; i < num_verts; i++ )
613  {
614  CartVect dx = CartVect( &new_coords[3 * i] ) - CartVect( &coords[3 * i] );
615  double tmp_len = dx.length_squared();
616  if( tmp_len > len ) len = tmp_len;
617  }
618  Range tmp_elems( fluid_elems );
619  tmp_elems.merge( solid_elems );
620  BoundBox box;
621  box.update( *mbImpl, tmp_elems );
622  double max_len =
623  std::max( box.bMax[2] - box.bMin[2], std::max( box.bMax[1] - box.bMin[1], box.bMax[0] - box.bMin[0] ) );
624 
625  cout << "Max displacement = " << len << " (" << 100.0 * len / max_len << "% of max box length)" << endl;
626  }
627 
628  if( !xNew )
629  {
630  rval = mbImpl->tag_get_handle( ( tag_name ? tag_name : "" ), 3, MB_TYPE_DOUBLE, xDisp[0],
631  MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Failed to get xNew tag" );
632  xNew = xDisp[0];
633  }
634 
635  // Set the new tag to those coords
636  rval = mbImpl->tag_set_data( xNew, solid_verts, &new_coords[0] );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
637 
638  // Get all the vertices and coords in the fluid, set xnew to them
639  rval = mbImpl->get_adjacencies( fluid_elems, 0, false, fluid_verts, Interface::UNION );MB_CHK_SET_ERR( rval, "Failed to get vertices" );
640  fluid_verts = subtract( fluid_verts, solid_verts );
641 
642  if( coords.size() < 3 * fluid_verts.size() ) coords.resize( 3 * fluid_verts.size() );
643  rval = mbImpl->get_coords( fluid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to get vertex coords" );
644  rval = mbImpl->tag_set_data( xNew, fluid_verts, &coords[0] );MB_CHK_SET_ERR( rval, "Failed to set xnew tag on fluid verts" );
645 
646  if( debug )
647  {
648  // Save deformed mesh coords to new file for visualizing
649  Range tmp_range( fluidElems[MASTER] );
650  tmp_range.merge( solidElems[MASTER] );
651  rval = write_and_save( tmp_range, masterSet, xNew, "deformed_master.h5m", true );MB_CHK_ERR( rval );
652  }
653 
654  return MB_SUCCESS;
655 }
656 
657 ErrorCode DeformMeshRemap::read_file( int m_or_s, string& fname, EntityHandle& seth )
658 {
659  // Create meshset
660  ErrorCode rval = mbImpl->create_meshset( 0, seth );MB_CHK_SET_ERR( rval, "Couldn't create master/slave set" );
661  ostringstream options;
662 #ifdef USE_MPI
663  ParallelComm* pc = ( m_or_s == MASTER ? pcMaster : pcSlave );
664  if( pc && pc->size() > 1 )
665  {
666  if( debug ) options << "DEBUG_IO=1;CPUTIME;";
667  options << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
668  << "PARALLEL_GHOSTS=2.0.1;PARALLEL_COMM=" << pc->get_id();
669  }
670 #endif
671  rval = mbImpl->load_file( fname.c_str(), &seth, options.str().c_str() );MB_CHK_SET_ERR( rval, "Couldn't load master/slave mesh" );
672 
673  if( *solidSetNos[m_or_s].begin() == -1 || *fluidSetNos[m_or_s].begin() == -1 ) return MB_SUCCESS;
674 
675  // Get material sets for solid/fluid
676  Tag tagh;
677  rval = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, tagh );MB_CHK_SET_ERR( rval, "Couldn't get material set tag name" );
678  for( set< int >::iterator sit = solidSetNos[m_or_s].begin(); sit != solidSetNos[m_or_s].end(); ++sit )
679  {
680  Range sets;
681  int set_no = *sit;
682  const void* setno_ptr = &set_no;
683  rval = mbImpl->get_entities_by_type_and_tag( seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets );
684  if( MB_SUCCESS != rval || sets.empty() )
685  {
686  MB_SET_ERR( MB_FAILURE, "Couldn't find solid set #" << *sit );
687  }
688  else
689  solidSets[m_or_s].merge( sets );
690  }
691 
692  // Get solid entities, and dimension
693  Range tmp_range;
694  for( Range::iterator rit = solidSets[m_or_s].begin(); rit != solidSets[m_or_s].end(); ++rit )
695  {
696  rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in solid" );
697  }
698  if( !tmp_range.empty() )
699  {
700  int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
701  assert( dim > 0 && dim < 4 );
702  solidElems[m_or_s] = tmp_range.subset_by_dimension( dim );
703  }
704 
705  if( debug )
706  cout << "Read " << solidElems[m_or_s].size() << " solid elements from " << solidSets[m_or_s].size()
707  << " sets in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh." << endl;
708 
709  for( set< int >::iterator sit = fluidSetNos[m_or_s].begin(); sit != fluidSetNos[m_or_s].end(); ++sit )
710  {
711  Range sets;
712  int set_no = *sit;
713  const void* setno_ptr = &set_no;
714  rval = mbImpl->get_entities_by_type_and_tag( seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets );
715  if( MB_SUCCESS != rval || sets.empty() )
716  {
717  MB_SET_ERR( MB_FAILURE, "Couldn't find fluid set #" << *sit );
718  }
719  else
720  fluidSets[m_or_s].merge( sets );
721  }
722 
723  // Get fluid entities, and dimension
724  tmp_range.clear();
725  for( Range::iterator rit = fluidSets[m_or_s].begin(); rit != fluidSets[m_or_s].end(); ++rit )
726  {
727  rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in fluid" );
728  }
729  if( !tmp_range.empty() )
730  {
731  int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
732  assert( dim > 0 && dim < 4 );
733  fluidElems[m_or_s] = tmp_range.subset_by_dimension( dim );
734  }
735 
736  if( debug )
737  cout << "Read " << fluidElems[m_or_s].size() << " fluid elements from " << fluidSets[m_or_s].size()
738  << " sets in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh." << endl;
739 
740  return rval;
741 }
742 
744 {
745  // Solid or fluid sets are missing; find the other
746  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
747 
748  if( fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty() )
749  {
750  unfilled_sets = &fluidSets[m_or_s];
751  filled_sets = &solidSets[m_or_s];
752  unfilled_elems = &fluidElems[m_or_s];
753  if( debug )
754  cout << "Finding unspecified fluid elements in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh...";
755  }
756  else if( !fluidSets[m_or_s].empty() && solidSets[m_or_s].empty() )
757  {
758  filled_sets = &fluidSets[m_or_s];
759  unfilled_sets = &solidSets[m_or_s];
760  unfilled_elems = &solidElems[m_or_s];
761  if( debug )
762  cout << "Finding unspecified solid elements in " << ( m_or_s == MASTER ? "master" : "slave" ) << " mesh...";
763  }
764 
765  // Ok, we know the filled sets, now fill the unfilled sets, and the elems from those
766  Tag tagh;
767  ErrorCode rval = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, tagh );MB_CHK_SET_ERR( rval, "Couldn't get material set tag name" );
768  Range matsets;
769  rval = mbImpl->get_entities_by_type_and_tag( file_set, MBENTITYSET, &tagh, NULL, 1, matsets );
770  if( MB_SUCCESS != rval || matsets.empty() )
771  {
772  MB_SET_ERR( MB_FAILURE, "Couldn't get any material sets" );
773  }
774  *unfilled_sets = subtract( matsets, *filled_sets );
775  if( unfilled_sets->empty() )
776  {
777  MB_SET_ERR( MB_FAILURE, "Failed to find any unfilled material sets" );
778  }
779  Range tmp_range;
780  for( Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); ++rit )
781  {
782  rval = mbImpl->get_entities_by_handle( *rit, tmp_range, true );MB_CHK_SET_ERR( rval, "Failed to get entities in unfilled set" );
783  }
784  int dim = mbImpl->dimension_from_handle( *tmp_range.rbegin() );
785  assert( dim > 0 && dim < 4 );
786  *unfilled_elems = tmp_range.subset_by_dimension( dim );
787  if( unfilled_elems->empty() )
788  {
789  MB_SET_ERR( MB_FAILURE, "Failed to find any unfilled set entities" );
790  }
791 
792  if( debug )
793  cout << "found " << unfilled_sets->size() << " sets and " << unfilled_elems->size() << " elements." << endl;
794 
795  return MB_SUCCESS;
796 }