 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 );
48 ErrorCode write_to_coords( Range& elems, Tag tagh );
49
50 const int SOLID_SETNO = 100, FLUID_SETNO = 200;
51
52 Interface* mb;
53
54 const bool debug = true;
55
56 class DeformMeshRemap
57 {
58 public:
59 //! Enumerator for solid/fluid, master/slave
60 enum
61 {
62 MASTER = 0,
63 SLAVE,
64 SOLID,
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
141 Interface* mbImpl;
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
155 EntityHandle masterSet, slaveSet;
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
170 Tag xNew;
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
266 ErrorCode DeformMeshRemap::execute()
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
417 DeformMeshRemap::DeformMeshRemap( Interface* impl, ParallelComm* master, ParallelComm* slave )
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
425 DeformMeshRemap::~DeformMeshRemap() {}
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 )
473 dfr->add_set_no( DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, *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 )
478 dfr->add_set_no( DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, *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 )
483 dfr->add_set_no( DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, *vit );
484
485 po.getOptAllArgs( "ss", set_nos );
486 for( vector< int >::iterator vit = set_nos.begin(); vit != set_nos.end(); ++vit )
487 dfr->add_set_no( DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, *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
504 ErrorCode DeformMeshRemap::write_and_save( Range& ents,
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
525 ErrorCode DeformMeshRemap::write_to_coords( Range& elems, Tag tagh, Tag tmp_tag )
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
743 ErrorCode DeformMeshRemap::find_other_sets( int m_or_s, EntityHandle file_set )
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 }
1.9.1.