Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
OptimizeMeshMesquite.cpp
Go to the documentation of this file.
1 #include "Mesquite.hpp"
2 #include "MsqIBase.hpp"
3 #include "MsqIGeom.hpp"
4 #include "MsqIMesh.hpp"
5 #include "MBiMesh.hpp"
6 #include "MeshImpl.hpp"
7 #include "moab/Core.hpp"
8 #include "moab/Skinner.hpp"
9 #include "moab/LloydSmoother.hpp"
10 #include "FacetModifyEngine.hpp"
11 #include "MsqError.hpp"
12 #include "InstructionQueue.hpp"
13 #include "PatchData.hpp"
14 #include "TerminationCriterion.hpp"
15 #include "QualityAssessor.hpp"
16 #include "SphericalDomain.hpp"
17 #include "PlanarDomain.hpp"
18 #include "MeshWriter.hpp"
19 
20 #include "moab/Core.hpp"
21 #include "moab/ProgOptions.hpp"
22 #ifdef MOAB_HAVE_MPI
23 #include "moab/ParallelComm.hpp"
24 #endif
25 
26 // algorithms
27 #include "IdealWeightInverseMeanRatio.hpp"
28 #include "TMPQualityMetric.hpp"
29 #include "AspectRatioGammaQualityMetric.hpp"
30 #include "ConditionNumberQualityMetric.hpp"
31 #include "VertexConditionNumberQualityMetric.hpp"
32 #include "MultiplyQualityMetric.hpp"
33 #include "LPtoPTemplate.hpp"
34 #include "LInfTemplate.hpp"
35 #include "PMeanPTemplate.hpp"
36 #include "SteepestDescent.hpp"
37 #include "FeasibleNewton.hpp"
38 #include "QuasiNewton.hpp"
39 #include "ConjugateGradient.hpp"
40 #include "SmartLaplacianSmoother.hpp"
41 #include "Randomize.hpp"
42 
43 #include "ReferenceMesh.hpp"
44 #include "RefMeshTargetCalculator.hpp"
45 #include "TShapeB1.hpp"
46 #include "TQualityMetric.hpp"
47 #include "IdealShapeTarget.hpp"
48 
49 #include <sys/stat.h>
50 #include <iostream>
51 using std::cerr;
52 using std::cout;
53 using std::endl;
54 
55 #include "iBase.h"
56 
57 using namespace MBMesquite;
58 
59 static int print_iGeom_error( const char* desc, int err, iGeom_Instance geom, const char* file, int line )
60 {
61  char buffer[1024];
62  iGeom_getDescription( geom, buffer, sizeof( buffer ) );
63  buffer[sizeof( buffer ) - 1] = '\0';
64 
65  std::cerr << "ERROR: " << desc << std::endl
66  << " Error code: " << err << std::endl
67  << " Error desc: " << buffer << std::endl
68  << " At : " << file << ':' << line << std::endl;
69 
70  return -1; // must always return false or CHECK macro will break
71 }
72 
73 static int print_iMesh_error( const char* desc, int err, iMesh_Instance mesh, const char* file, int line )
74 {
75  char buffer[1024];
76  iMesh_getDescription( mesh, buffer, sizeof( buffer ) );
77  buffer[sizeof( buffer ) - 1] = '\0';
78 
79  std::cerr << "ERROR: " << desc << std::endl
80  << " Error code: " << err << std::endl
81  << " Error desc: " << buffer << std::endl
82  << " At : " << file << ':' << line << std::endl;
83 
84  return -1; // must always return false or CHECK macro will break
85 }
86 
87 #define CHECK_IGEOM( STR ) \
88  if( err != iBase_SUCCESS ) return print_iGeom_error( STR, err, geom, __FILE__, __LINE__ )
89 
90 #define CHECK_IMESH( STR ) \
91  if( err != iBase_SUCCESS ) return print_iMesh_error( STR, err, instance, __FILE__, __LINE__ )
92 
93 const std::string default_file_name = std::string( MESH_DIR ) + std::string( "/surfrandomtris-4part.h5m" );
94 const std::string default_geometry_file_name = std::string( MESH_DIR ) + std::string( "/surfrandom.facet" );
95 
96 std::vector< double > solution_indicator;
97 
98 int write_vtk_mesh( Mesh* mesh, const char* filename );
99 
100 // Construct a MeshTSTT from the file
101 int get_imesh_mesh( MBMesquite::Mesh**, const char* file_name, int dimension );
102 
103 // Construct a MeshImpl from the file
104 int get_native_mesh( MBMesquite::Mesh**, const char* file_name, int dimension );
105 
106 int get_itaps_domain( MeshDomain**, const char* );
107 int get_mesquite_domain( MeshDomain**, const char* );
108 
109 // Run FeasibleNewton solver
110 int run_global_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value = 1e-4 );
111 
112 // Run SmoothLaplacian solver
113 int run_local_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value = 1e-3 );
114 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value = 1e-3 );
115 
116 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err );
117 
118 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err );
119 
120 bool file_exists( const std::string& name )
121 {
122  struct stat buffer;
123  return ( stat( name.c_str(), &buffer ) == 0 );
124 }
125 
126 int main( int argc, char* argv[] )
127 {
128  MBMesquite::MsqPrintError err( cout );
129  // command line arguments
130  std::string file_name, geometry_file_name;
131  bool use_native = false;
132  int dimension = 2;
133 
134 #ifdef MOAB_HAVE_MPI
135 // MPI_Init(&argc, &argv);
136 #endif
137  ProgOptions opts;
138 
139  opts.addOpt< void >( std::string( "native,N" ), std::string( "Use native representation (default=false)" ),
140  &use_native );
141  opts.addOpt< int >( std::string( "dim,d" ), std::string( "Topological dimension of the mesh (default=2)" ),
142  &dimension );
143  opts.addOpt< std::string >( std::string( "input_geo,i" ),
144  std::string( "Input file name for the geometry (pattern=*.stl, *.facet)" ),
145  &geometry_file_name );
146  opts.addOpt< std::string >( std::string( "input_mesh,m" ),
147  std::string( "Input file name for the mesh (pattern=*.vtk, *.h5m)" ), &file_name );
148 
149  opts.parseCommandLine( argc, argv );
150 
151  if( !geometry_file_name.length() )
152  {
153  file_name = default_file_name;
154  geometry_file_name = default_geometry_file_name;
155  cout << "No file specified: Using defaults.\n";
156  }
157  cout << "\t Mesh filename = " << file_name << endl;
158  cout << "\t Geometry filename = " << geometry_file_name << endl;
159 
160  // Vector3D pnt(0,0,0);
161  // Vector3D s_norm(0,0,1);
162  // PlanarDomain plane(s_norm, pnt);
163 
164  // PlanarDomain plane( PlanarDomain::XY );
165  MeshDomain* plane;
166  int ierr;
167  if( !file_exists( geometry_file_name ) ) geometry_file_name = "";
168  ierr = get_itaps_domain( &plane, geometry_file_name.c_str() ); // MB_CHK_ERR(ierr);
169 
170  // Try running a global smoother on the mesh
171  Mesh* mesh = 0;
172  if( use_native )
173  {
174  ierr = get_native_mesh( &mesh, file_name.c_str(), dimension ); // MB_CHK_ERR(ierr);
175  }
176  else
177  {
178  ierr = get_imesh_mesh( &mesh, file_name.c_str(), dimension ); // MB_CHK_ERR(ierr);
179  }
180 
181  if( !mesh )
182  {
183  std::cerr << "Failed to load input file. Aborting." << std::endl;
184  return 1;
185  }
186 
187  MeshDomainAssoc mesh_and_domain( mesh, plane );
188 
189  // ierr = write_vtk_mesh( mesh, "original.vtk");MB_CHK_ERR(ierr);
190  // cout << "Wrote \"original.vtk\"" << endl;
191 
192  // run_global_smoother( mesh_and_domain, err );
193  // if (err) return 1;
194 
195  // Try running a local smoother on the mesh
196  // Mesh* meshl=0;
197  // if(use_native)
198  // ierr = get_native_mesh(&meshl, file_name.c_str(), dimension);
199  // else
200  // ierr = get_imesh_mesh(&meshl, file_name.c_str(), dimension);
201  // if (!mesh || ierr) {
202  // std::cerr << "Failed to load input file. Aborting." << std::endl;
203  // return 1;
204  // }
205 
206  // MeshDomainAssoc mesh_and_domain_local(meshl, plane);
207 
208  // run_solution_mesh_optimizer( mesh_and_domain, err );
209  // if (err) return 1;
210 
211  run_local_smoother( mesh_and_domain, err, 1e-4 ); // MB_CHK_ERR(err);
212  if( err ) return 1;
213 
214  double reps = 5e-2;
215  for( int iter = 0; iter < 5; iter++ )
216  {
217 
218  if( !( iter % 2 ) )
219  {
220  run_local_smoother2( mesh_and_domain, err,
221  reps * 10 ); // CHECK_IMESH("local smoother2 failed");
222  if( err ) return 1;
223  }
224 
225  // run_global_smoother( mesh_and_domain, err, reps );MB_CHK_ERR(ierr);
226 
227  run_solution_mesh_optimizer( mesh_and_domain,
228  err ); // CHECK_IMESH("solution mesh optimizer failed");
229  if( err ) return 1;
230 
231  reps *= 0.01;
232  }
233 
234  run_local_smoother2( mesh_and_domain, err, 1e-4 ); // CHECK_IMESH("local smoother2 failed");
235  if( err ) return 1;
236 
237  // run_quality_optimizer( mesh_and_domain, err );MB_CHK_ERR(ierr);
238 
239  // run_local_smoother( mesh_and_domain, err );MB_CHK_ERR(ierr);
240 
241  // Delete MOAB instance
242  delete mesh;
243  delete plane;
244 
245 #ifdef MOAB_HAVE_MPI
246  // MPI_Finalize();
247 #endif
248 
249  return 0;
250 }
251 
252 int run_global_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
253 {
254  // double OF_value = 1e-6;
255 
256  Mesh* mesh = mesh_and_domain.get_mesh();
257  MeshDomain* domain = mesh_and_domain.get_domain();
258 
259  // creates an intruction queue
260  InstructionQueue queue1;
261 
262  // creates a mean ratio quality metric ...
263  IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio( err );
264  // ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
265  // TMPQualityMetric* mean_ratio = new TMPQualityMetric();
266 
267  // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
268  if( err ) return 1;
269  mean_ratio->set_averaging_method( QualityMetric::RMS, err );
270  if( err ) return 1;
271 
272  // ... and builds an objective function with it
273  LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
274  // LInfTemplate* obj_func = new LInfTemplate(mean_ratio);
275  if( err ) return 1;
276 
277  // creates the feas newt optimization procedures
278  // ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
279  // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
280  SteepestDescent* pass1 = new SteepestDescent( obj_func );
281  pass1->use_element_on_vertex_patch();
282  pass1->do_block_coordinate_descent_optimization();
283  pass1->use_global_patch();
284  if( err ) return 1;
285 
286  QualityAssessor stop_qa( mean_ratio );
287 
288  // **************Set stopping criterion****************
289  TerminationCriterion tc_inner;
290  tc_inner.add_absolute_vertex_movement( OF_value );
291  if( err ) return 1;
292  TerminationCriterion tc_outer;
293  tc_inner.add_iteration_limit( 10 );
294  tc_outer.add_iteration_limit( 5 );
295  tc_outer.set_debug_output_level( 3 );
296  tc_inner.set_debug_output_level( 3 );
297  pass1->set_inner_termination_criterion( &tc_inner );
298  pass1->set_outer_termination_criterion( &tc_outer );
299 
300  queue1.add_quality_assessor( &stop_qa, err );
301  if( err ) return 1;
302 
303  // adds 1 pass of pass1 to mesh_set1
304  queue1.set_master_quality_improver( pass1, err );
305  if( err ) return 1;
306 
307  queue1.add_quality_assessor( &stop_qa, err );
308  if( err ) return 1;
309 
310  // launches optimization on mesh_set
311  if( domain )
312  {
313  queue1.run_instructions( &mesh_and_domain, err );
314  }
315  else
316  {
317  queue1.run_instructions( mesh, err );
318  }
319  if( err ) return 1;
320 
321  // Construct a MeshTSTT from the file
322  int ierr = write_vtk_mesh( mesh, "feasible-newton-result.vtk" );
323  if( ierr ) return 1;
324  // MeshWriter::write_vtk(mesh, "feasible-newton-result.vtk", err);
325  // if (err) return 1;
326  cout << "Wrote \"feasible-newton-result.vtk\"" << endl;
327 
328  // print_timing_diagnostics( cout );
329  return 0;
330 }
331 
332 int write_vtk_mesh( Mesh* mesh, const char* filename )
333 {
334  moab::Interface* mbi =
335  reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl;
336 
337  mbi->write_file( filename );
338 
339  return 0;
340 }
341 
342 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value );
343 int run_local_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
344 {
345  Mesh* mesh = mesh_and_domain.get_mesh();
346  moab::Interface* mbi =
347  reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl;
348 
349  moab::Tag fixed;
350  moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
351  moab::Range cells;
352  rval = mbi->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" );
353 
354  moab::LloydSmoother lloyd( mbi, 0, cells, 0, 0 /*fixed*/ );
355 
356  lloyd.perform_smooth();
357 
358  // run_local_smoother2(mesh_and_domain, err, OF_value);
359 
360  // Construct a MeshTSTT from the file
361  int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" );
362  if( ierr ) return 1;
363  // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
364  // if (err) return 1;
365  cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
366  return 0;
367 }
368 
369 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
370 {
371  // double OF_value = 1e-5;
372 
373  Mesh* mesh = mesh_and_domain.get_mesh();
374  MeshDomain* domain = mesh_and_domain.get_domain();
375 
376  // creates an intruction queue
377  InstructionQueue queue1;
378 
379  // creates a mean ratio quality metric ...
380  // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
381  ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
382  // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
383  if( err ) return 1;
384  // mean_ratio->set_gradient_type(QualityMetric::NUMERICAL_GRADIENT);
385  // mean_ratio->set_hessian_type(QualityMetric::NUMERICAL_HESSIAN);
386  mean_ratio->set_averaging_method( QualityMetric::RMS, err );
387  if( err ) return 1;
388 
389  // ... and builds an objective function with it
390  LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
391  if( err ) return 1;
392 
393  if( false )
394  {
395  InstructionQueue qtmp;
396  Randomize rand( -0.005 );
397  TerminationCriterion sc_rand;
398  sc_rand.add_iteration_limit( 2 );
399  rand.set_outer_termination_criterion( &sc_rand );
400  qtmp.set_master_quality_improver( &rand, err );
401  if( err ) return 1;
402  if( domain )
403  {
404  qtmp.run_instructions( &mesh_and_domain, err );
405  }
406  else
407  {
408  qtmp.run_instructions( mesh, err );
409  }
410  if( err ) return 1;
411  }
412 
413  // creates the smart laplacian optimization procedures
414  SmartLaplacianSmoother* pass1 = new SmartLaplacianSmoother( obj_func );
415  // SteepestDescent* pass1 = new SteepestDescent( obj_func );
416 
417  QualityAssessor stop_qa( mean_ratio );
418 
419  // **************Set stopping criterion****************
420  TerminationCriterion tc_inner;
421  tc_inner.add_absolute_vertex_movement( OF_value );
422  TerminationCriterion tc_outer;
423  tc_outer.add_iteration_limit( 10 );
424  pass1->set_inner_termination_criterion( &tc_inner );
425  pass1->set_outer_termination_criterion( &tc_outer );
426 
427  queue1.add_quality_assessor( &stop_qa, err );
428  if( err ) return 1;
429 
430  // adds 1 pass of pass1 to mesh_set
431  queue1.set_master_quality_improver( pass1, err );
432  if( err ) return 1;
433 
434  queue1.add_quality_assessor( &stop_qa, err );
435  if( err ) return 1;
436 
437  // launches optimization on mesh_set
438  if( domain )
439  {
440  queue1.run_instructions( &mesh_and_domain, err );
441  }
442  else
443  {
444  queue1.run_instructions( mesh, err );
445  }
446  if( err ) return 1;
447 
448  // Construct a MeshTSTT from the file
449  int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" );
450  if( ierr ) return 1;
451  // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
452  // if (err) return 1;
453  cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
454 
455  // print_timing_diagnostics( cout );
456  return 0;
457 }
458 
459 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
460 {
461  Mesh* mesh = mesh_and_domain.get_mesh();
462  MeshDomain* domain = mesh_and_domain.get_domain();
463 
464  // creates an intruction queue
465  InstructionQueue q;
466 
467  // do optimization
468  const double eps = 0.01;
469  IdealShapeTarget w;
470  TShapeB1 mu;
471  TQualityMetric metric( &w, &mu );
472  PMeanPTemplate func( 1.0, &metric );
473 
474  SteepestDescent solver( &func );
475  solver.use_element_on_vertex_patch();
476  solver.do_jacobi_optimization();
477 
478  TerminationCriterion inner, outer;
479  inner.add_absolute_vertex_movement( 0.5 * eps );
480  outer.add_absolute_vertex_movement( 0.5 * eps );
481 
482  QualityAssessor qa( &metric );
483 
484  q.add_quality_assessor( &qa, err );
485  if( err ) return 1;
486  q.set_master_quality_improver( &solver, err );
487  if( err ) return 1;
488  q.add_quality_assessor( &qa, err );
489  if( err ) return 1;
490 
491  // launches optimization on mesh_set
492  if( domain )
493  {
494  q.run_instructions( &mesh_and_domain, err );
495  }
496  else
497  {
498  q.run_instructions( mesh, err );
499  }
500  if( err ) return 1;
501 
502  // Construct a MeshTSTT from the file
503  int ierr = write_vtk_mesh( mesh, "ideal-shape-result.vtk" );
504  if( ierr ) return 1;
505  // MeshWriter::write_vtk(mesh, "ideal-shape-result.vtk", err);
506  // if (err) return 1;
507  cout << "Wrote \"ideal-shape-result.vtk\"" << endl;
508 
509  print_timing_diagnostics( cout );
510  return 0;
511 }
512 
513 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
514 {
515  double OF_value = 0.01;
516 
517  Mesh* mesh = mesh_and_domain.get_mesh();
518  MeshDomain* domain = mesh_and_domain.get_domain();
519 
520  // creates an intruction queue
521  InstructionQueue queue1;
522 
523  // creates a mean ratio quality metric ...
524  // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
525  ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
526  // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
527  // AspectRatioGammaQualityMetric* mean_ratio = new AspectRatioGammaQualityMetric();
528 
529  // ElementSolIndQM* solindqm = new ElementSolIndQM(solution_indicator);
530  // MultiplyQualityMetric* mean_ratio = new MultiplyQualityMetric(new
531  // VertexConditionNumberQualityMetric(), solindqm, err);
532  // ElementSolIndQM* mean_ratio = solindqm;
533 
534  // LocalSizeQualityMetric* mean_ratio = new LocalSizeQualityMetric();
535 
536  mean_ratio->set_averaging_method( QualityMetric::SUM_OF_RATIOS_SQUARED, err );
537  if( err ) return 1;
538 
539  // ... and builds an objective function with it
540  LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
541  if( err ) return 1;
542 
543  // // creates the feas newt optimization procedures
544  ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
545  // QuasiNewton* pass1 = new QuasiNewton( obj_func );
546  // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
547  pass1->use_global_patch();
548 
549  QualityAssessor stop_qa( mean_ratio );
550 
551  // **************Set stopping criterion****************
552  TerminationCriterion tc_inner;
553  tc_inner.add_absolute_vertex_movement( OF_value );
554  if( err ) return 1;
555  TerminationCriterion tc_outer;
556  tc_inner.add_iteration_limit( 20 );
557  tc_outer.add_iteration_limit( 5 );
558  pass1->set_inner_termination_criterion( &tc_inner );
559  pass1->set_outer_termination_criterion( &tc_outer );
560  pass1->set_debugging_level( 3 );
561 
562  queue1.add_quality_assessor( &stop_qa, err );
563  if( err ) return 1;
564 
565  // adds 1 pass of pass1 to mesh_set1
566  queue1.set_master_quality_improver( pass1, err );
567  if( err ) return 1;
568 
569  queue1.add_quality_assessor( &stop_qa, err );
570  if( err ) return 1;
571 
572  // launches optimization on mesh_set
573  if( domain )
574  {
575  queue1.run_instructions( &mesh_and_domain, err );
576  }
577  else
578  {
579  queue1.run_instructions( mesh, err );
580  }
581  if( err ) return 1;
582 
583  // Construct a MeshTSTT from the file
584  int ierr = write_vtk_mesh( mesh, "solution-mesh-result.vtk" );
585  if( ierr ) return 1;
586  // MeshWriter::write_vtk(mesh, "solution-mesh-result.vtk", err);
587  // if (err) return 1;
588  cout << "Wrote \"solution-mesh-result.vtk\"" << endl;
589 
590  print_timing_diagnostics( cout );
591  return 0;
592 }
593 
594 int get_imesh_mesh( MBMesquite::Mesh** mesh, const char* file_name, int dimension )
595 {
596  int err;
597  iMesh_Instance instance = 0;
598  iMesh_newMesh( NULL, &instance, &err, 0 );
599  CHECK_IMESH( "Creation of mesh instance failed" );
600 
601  iBase_EntitySetHandle root_set;
602  iMesh_getRootSet( instance, &root_set, &err );
603  CHECK_IMESH( "Could not get root set" );
604 
605  iMesh_load( instance, root_set, file_name, 0, &err, strlen( file_name ), 0 );
606  CHECK_IMESH( "Could not load mesh from file" );
607 
608  iBase_TagHandle fixed_tag;
609  iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) );
610  // if (iBase_SUCCESS != err)
611  {
612  // get the skin vertices of those cells and mark them as fixed; we don't want to fix the
613  // vertices on a part boundary, but since we exchanged a layer of ghost cells, those
614  // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move
615  // those anyway use MOAB's skinner class to find the skin
616 
617  // get all vertices and cells
618  // make tag to specify fixed vertices, since it's input to the algorithm; use a default
619  // value of non-fixed so we only need to set the fixed tag for skin vertices
620  moab::Interface* mbi = reinterpret_cast< MBiMesh* >( instance )->mbImpl;
621  moab::EntityHandle currset = 0;
622  moab::Tag fixed;
623  int def_val = 0;
624  err = 0;
625  moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed,
626  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
627  moab::Range verts, cells, skin_verts;
628  rval = mbi->get_entities_by_type( currset, moab::MBVERTEX, verts );MB_CHK_SET_ERR( rval, "Querying vertices failed" );
629  rval = mbi->get_entities_by_dimension( currset, dimension, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" );
630  std::cout << "Found " << verts.size() << " vertices and " << cells.size() << " elements" << std::endl;
631 
632  moab::Skinner skinner( mbi );
633  rval = skinner.find_skin( currset, cells, true, skin_verts );MB_CHK_SET_ERR( rval,
634  "Finding the skin of the mesh failed" ); // 'true' param indicates we want
635  // vertices back, not cells
636 
637  std::vector< int > fix_tag( skin_verts.size(), 1 ); // initialized to 1 to indicate fixed
638  rval = mbi->tag_set_data( fixed, skin_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" );
639  std::cout << "Found " << skin_verts.size() << " vertices on the skin of the domain." << std::endl;
640 
641  // fix_tag.resize(verts.size(),0);
642  // rval = mbi->tag_get_data(fixed, verts, &fix_tag[0]); MB_CHK_SET_ERR(rval, "Getting tag
643  // data failed");
644 
645  iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) );
646  CHECK_IMESH( "Getting tag handle (fixed) failed" );
647 
648  // Set some arbitrary solution indicator
649  moab::Tag solindTag;
650  double def_val_dbl = 0.0;
651  rval = mbi->tag_get_handle( "solution_indicator", 1, moab::MB_TYPE_DOUBLE, solindTag,
652  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val_dbl );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
653  solution_indicator.resize( cells.size(), 0.01 );
654  for( unsigned i = 0; i < cells.size() / 4; i++ )
655  solution_indicator[i] = 0.1;
656  for( unsigned i = cells.size() / 4; i < 2 * cells.size() / 4; i++ )
657  solution_indicator[i] = 0.5;
658  for( unsigned i = 2 * cells.size() / 4; i < 3 * cells.size() / 4; i++ )
659  solution_indicator[i] = 0.5;
660  for( unsigned i = 3 * cells.size() / 4; i < cells.size(); i++ )
661  solution_indicator[i] = 0.5;
662 
663  rval = mbi->tag_set_data( solindTag, cells, &solution_indicator[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" );
664  }
665 
666  MsqError ierr;
667  MBMesquite::MsqIMesh* imesh =
668  new MBMesquite::MsqIMesh( instance, root_set, ( dimension == 3 ? iBase_REGION : iBase_FACE ), ierr,
669  &fixed_tag );
670  if( MSQ_CHKERR( ierr ) )
671  {
672  delete imesh;
673  cerr << err << endl;
674  err = iBase_FAILURE;
675  CHECK_IMESH( "Creation of MsqIMesh instance failed" );
676  return 0;
677  }
678 
679  *mesh = imesh;
680  return iBase_SUCCESS;
681 }
682 
683 int get_native_mesh( MBMesquite::Mesh** mesh, const char* file_name, int )
684 {
685  MsqError err;
686  MBMesquite::MeshImpl* imesh = new MBMesquite::MeshImpl();
687  imesh->read_vtk( file_name, err );
688  if( err )
689  {
690  cerr << err << endl;
691  exit( 3 );
692  }
693  *mesh = imesh;
694 
695  return iBase_SUCCESS;
696 }
697 
698 int get_itaps_domain( MeshDomain** odomain, const char* filename )
699 {
700 
701  if( filename == 0 || strlen( filename ) == 0 )
702  {
703  *odomain = new PlanarDomain( PlanarDomain::XY );
704  return 0;
705  }
706 
707  int err;
708  iGeom_Instance geom;
709  iGeom_newGeom( "", &geom, &err, 0 );
710  CHECK_IGEOM( "ERROR: iGeom creation failed" );
711 
712 #ifdef MOAB_HAVE_CGM_FACET
713  FacetModifyEngine::set_modify_enabled( CUBIT_TRUE );
714 #endif
715 
716  iGeom_load( geom, filename, 0, &err, strlen( filename ), 0 );
717  CHECK_IGEOM( "Cannot load the geometry" );
718 
719  iBase_EntitySetHandle root_set;
720  iGeom_getRootSet( geom, &root_set, &err );
721  CHECK_IGEOM( "getRootSet failed!" );
722 
723  // print out the number of entities
724  std::cout << "Model contents: " << std::endl;
725  const char* gtype[] = { "vertices: ", "edges: ", "faces: ", "regions: " };
726  int nents[4];
727  for( int i = 0; i <= 3; ++i )
728  {
729  iGeom_getNumOfType( geom, root_set, i, &nents[i], &err );
730  CHECK_IGEOM( "Error: problem getting entities after gLoad." );
731  std::cout << gtype[i] << nents[i] << std::endl;
732  }
733 
734  iBase_EntityHandle* hd_geom_ents;
735  int csize = 0, sizealloc = 0;
736  if( nents[3] > 0 )
737  {
738  hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[2] );
739  csize = nents[2];
740  iGeom_getEntities( geom, root_set, 2, &hd_geom_ents, &csize, &sizealloc, &err );
741  }
742  else
743  {
744  hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[1] );
745  csize = nents[1];
746  iGeom_getEntities( geom, root_set, 1, &hd_geom_ents, &csize, &sizealloc, &err );
747  }
748  CHECK_IGEOM( "ERROR: Could not get entities" );
749 
750  *odomain = new MsqIGeom( geom, hd_geom_ents[0] );
751  return iBase_SUCCESS;
752 }