OptimizeMeshMesquite.cpp

Go to the documentation of this file.

#include "Mesquite.hpp"
#include "MsqIBase.hpp"
#include "MsqIGeom.hpp"
#include "MsqIMesh.hpp"
#include "MBiMesh.hpp"
#include "MeshImpl.hpp"
#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/LloydSmoother.hpp"
#include "FacetModifyEngine.hpp"
#include "MsqError.hpp"
#include "InstructionQueue.hpp"
#include "PatchData.hpp"
#include "TerminationCriterion.hpp"
#include "QualityAssessor.hpp"
#include "SphericalDomain.hpp"
#include "PlanarDomain.hpp"
#include "MeshWriter.hpp"
 
#include "moab/Core.hpp"
#include "moab/ProgOptions.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#endif
 
// algorithms
#include "IdealWeightInverseMeanRatio.hpp"
#include "TMPQualityMetric.hpp"
#include "AspectRatioGammaQualityMetric.hpp"
#include "ConditionNumberQualityMetric.hpp"
#include "VertexConditionNumberQualityMetric.hpp"
#include "MultiplyQualityMetric.hpp"
#include "LPtoPTemplate.hpp"
#include "LInfTemplate.hpp"
#include "PMeanPTemplate.hpp"
#include "SteepestDescent.hpp"
#include "FeasibleNewton.hpp"
#include "QuasiNewton.hpp"
#include "ConjugateGradient.hpp"
#include "SmartLaplacianSmoother.hpp"
#include "Randomize.hpp"
 
#include "ReferenceMesh.hpp"
#include "RefMeshTargetCalculator.hpp"
#include "TShapeB1.hpp"
#include "TQualityMetric.hpp"
#include "IdealShapeTarget.hpp"
 
#include <sys/stat.h>
#include <iostream>
using std::cerr;
using std::cout;
using std::endl;
 
#include "iBase.h"
 
using namespace MBMesquite;
 
static int print_iGeom_error( const char* desc, int err, iGeom_Instance geom, const char* file, int line )
{
 char buffer[1024];
 iGeom_getDescription( geom, buffer, sizeof( buffer ) );
 buffer[sizeof( buffer ) - 1] = '\0';
 
 std::cerr << "ERROR: " << desc << std::endl
 << " Error code: " << err << std::endl
 << " Error desc: " << buffer << std::endl
 << " At : " << file << ':' << line << std::endl;
 
 return -1; // must always return false or CHECK macro will break
}
 
static int print_iMesh_error( const char* desc, int err, iMesh_Instance mesh, const char* file, int line )
{
 char buffer[1024];
 iMesh_getDescription( mesh, buffer, sizeof( buffer ) );
 buffer[sizeof( buffer ) - 1] = '\0';
 
 std::cerr << "ERROR: " << desc << std::endl
 << " Error code: " << err << std::endl
 << " Error desc: " << buffer << std::endl
 << " At : " << file << ':' << line << std::endl;
 
 return -1; // must always return false or CHECK macro will break
}
 
#define CHECK_IGEOM( STR ) \
 if( err != iBase_SUCCESS ) return print_iGeom_error( STR, err, geom, __FILE__, __LINE__ )
 
#define CHECK_IMESH( STR ) \
 if( err != iBase_SUCCESS ) return print_iMesh_error( STR, err, instance, __FILE__, __LINE__ )
 
const std::string default_file_name = std::string( MESH_DIR ) + std::string( "/surfrandomtris-4part.h5m" );
const std::string default_geometry_file_name = std::string( MESH_DIR ) + std::string( "/surfrandom.facet" );
 
std::vector< double > solution_indicator;
 
int write_vtk_mesh( Mesh* mesh, const char* filename );
 
// Construct a MeshTSTT from the file
int get_imesh_mesh( MBMesquite::Mesh**, const char* file_name, int dimension );
 
// Construct a MeshImpl from the file
int get_native_mesh( MBMesquite::Mesh**, const char* file_name, int dimension );
 
int get_itaps_domain( MeshDomain**, const char* );
int get_mesquite_domain( MeshDomain**, const char* );
 
// Run FeasibleNewton solver
int run_global_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value = 1e-4 );
 
// Run SmoothLaplacian solver
int run_local_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value = 1e-3 );
int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value = 1e-3 );
 
int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err );
 
int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err );
 
bool file_exists( const std::string& name )
{
 struct stat buffer;
 return ( stat( name.c_str(), &buffer ) == 0 );
}
 
int main( int argc, char* argv[] )
{
 MBMesquite::MsqPrintError err( cout );
 // command line arguments
 std::string file_name, geometry_file_name;
 bool use_native = false;
 int dimension = 2;
 
#ifdef MOAB_HAVE_MPI
// MPI_Init(&argc, &argv);
#endif
 ProgOptions opts;
 
 opts.addOpt< void >( std::string( "native,N" ), std::string( "Use native representation (default=false)" ),
 &use_native );
 opts.addOpt< int >( std::string( "dim,d" ), std::string( "Topological dimension of the mesh (default=2)" ),
 &dimension );
 opts.addOpt< std::string >( std::string( "input_geo,i" ),
 std::string( "Input file name for the geometry (pattern=*.stl, *.facet)" ),
 &geometry_file_name );
 opts.addOpt< std::string >( std::string( "input_mesh,m" ),
 std::string( "Input file name for the mesh (pattern=*.vtk, *.h5m)" ), &file_name );
 
 opts.parseCommandLine( argc, argv );
 
 if( !geometry_file_name.length() )
 {
 file_name = default_file_name;
 geometry_file_name = default_geometry_file_name;
 cout << "No file specified: Using defaults.\n";
 }
 cout << "\t Mesh filename = " << file_name << endl;
 cout << "\t Geometry filename = " << geometry_file_name << endl;
 
 // Vector3D pnt(0,0,0);
 // Vector3D s_norm(0,0,1);
 // PlanarDomain plane(s_norm, pnt);
 
 // PlanarDomain plane( PlanarDomain::XY );
 MeshDomain* plane;
 int ierr;
 if( !file_exists( geometry_file_name ) ) geometry_file_name = "";
 ierr = get_itaps_domain( &plane, geometry_file_name.c_str() ); // MB_CHK_ERR(ierr);
 
 // Try running a global smoother on the mesh
 Mesh* mesh = 0;
 if( use_native )
 {
 ierr = get_native_mesh( &mesh, file_name.c_str(), dimension ); // MB_CHK_ERR(ierr);
 }
 else
 {
 ierr = get_imesh_mesh( &mesh, file_name.c_str(), dimension ); // MB_CHK_ERR(ierr);
 }
 
 if( !mesh )
 {
 std::cerr << "Failed to load input file. Aborting." << std::endl;
 return 1;
 }
 
 MeshDomainAssoc mesh_and_domain( mesh, plane );
 
 // ierr = write_vtk_mesh( mesh, "original.vtk");MB_CHK_ERR(ierr);
 // cout << "Wrote \"original.vtk\"" << endl;
 
 // run_global_smoother( mesh_and_domain, err );
 // if (err) return 1;
 
 // Try running a local smoother on the mesh
 // Mesh* meshl=0;
 // if(use_native)
 // ierr = get_native_mesh(&meshl, file_name.c_str(), dimension);
 // else
 // ierr = get_imesh_mesh(&meshl, file_name.c_str(), dimension);
 // if (!mesh || ierr) {
 // std::cerr << "Failed to load input file. Aborting." << std::endl;
 // return 1;
 // }
 
 // MeshDomainAssoc mesh_and_domain_local(meshl, plane);
 
 // run_solution_mesh_optimizer( mesh_and_domain, err );
 // if (err) return 1;
 
 run_local_smoother( mesh_and_domain, err, 1e-4 ); // MB_CHK_ERR(err);
 if( err ) return 1;
 
 double reps = 5e-2;
 for( int iter = 0; iter < 5; iter++ )
 {
 
 if( !( iter % 2 ) )
 {
 run_local_smoother2( mesh_and_domain, err,
 reps * 10 ); // CHECK_IMESH("local smoother2 failed");
 if( err ) return 1;
 }
 
 // run_global_smoother( mesh_and_domain, err, reps );MB_CHK_ERR(ierr);
 
 run_solution_mesh_optimizer( mesh_and_domain,
 err ); // CHECK_IMESH("solution mesh optimizer failed");
 if( err ) return 1;
 
 reps *= 0.01;
 }
 
 run_local_smoother2( mesh_and_domain, err, 1e-4 ); // CHECK_IMESH("local smoother2 failed");
 if( err ) return 1;
 
 // run_quality_optimizer( mesh_and_domain, err );MB_CHK_ERR(ierr);
 
 // run_local_smoother( mesh_and_domain, err );MB_CHK_ERR(ierr);
 
 // Delete MOAB instance
 delete mesh;
 delete plane;
 
#ifdef MOAB_HAVE_MPI
 // MPI_Finalize();
#endif
 
 return 0;
}
 
int run_global_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
{
 // double OF_value = 1e-6;
 
 Mesh* mesh = mesh_and_domain.get_mesh();
 MeshDomain* domain = mesh_and_domain.get_domain();
 
 // creates an intruction queue
 InstructionQueue queue1;
 
 // creates a mean ratio quality metric ...
 IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio( err );
 // ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
 // TMPQualityMetric* mean_ratio = new TMPQualityMetric();
 
 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
 if( err ) return 1;
 mean_ratio->set_averaging_method( QualityMetric::RMS, err );
 if( err ) return 1;
 
 // ... and builds an objective function with it
 LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
 // LInfTemplate* obj_func = new LInfTemplate(mean_ratio);
 if( err ) return 1;
 
 // creates the feas newt optimization procedures
 // ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
 // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
 SteepestDescent* pass1 = new SteepestDescent( obj_func );
 pass1->use_element_on_vertex_patch();
 pass1->do_block_coordinate_descent_optimization();
 pass1->use_global_patch();
 if( err ) return 1;
 
 QualityAssessor stop_qa( mean_ratio );
 
 // **************Set stopping criterion****************
 TerminationCriterion tc_inner;
 tc_inner.add_absolute_vertex_movement( OF_value );
 if( err ) return 1;
 TerminationCriterion tc_outer;
 tc_inner.add_iteration_limit( 10 );
 tc_outer.add_iteration_limit( 5 );
 tc_outer.set_debug_output_level( 3 );
 tc_inner.set_debug_output_level( 3 );
 pass1->set_inner_termination_criterion( &tc_inner );
 pass1->set_outer_termination_criterion( &tc_outer );
 
 queue1.add_quality_assessor( &stop_qa, err );
 if( err ) return 1;
 
 // adds 1 pass of pass1 to mesh_set1
 queue1.set_master_quality_improver( pass1, err );
 if( err ) return 1;
 
 queue1.add_quality_assessor( &stop_qa, err );
 if( err ) return 1;
 
 // launches optimization on mesh_set
 if( domain )
 {
 queue1.run_instructions( &mesh_and_domain, err );
 }
 else
 {
 queue1.run_instructions( mesh, err );
 }
 if( err ) return 1;
 
 // Construct a MeshTSTT from the file
 int ierr = write_vtk_mesh( mesh, "feasible-newton-result.vtk" );
 if( ierr ) return 1;
 // MeshWriter::write_vtk(mesh, "feasible-newton-result.vtk", err);
 // if (err) return 1;
 cout << "Wrote \"feasible-newton-result.vtk\"" << endl;
 
 // print_timing_diagnostics( cout );
 return 0;
}
 
int write_vtk_mesh( Mesh* mesh, const char* filename )
{
 moab::Interface* mbi =
 reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl;
 
 mbi->write_file( filename );
 
 return 0;
}
 
int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value );
int run_local_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
{
 Mesh* mesh = mesh_and_domain.get_mesh();
 moab::Interface* mbi =
 reinterpret_cast< MBiMesh* >( dynamic_cast< MsqIMesh* >( mesh )->get_imesh_instance() )->mbImpl;
 
 moab::Tag fixed;
 moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
 moab::Range cells;
 rval = mbi->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" );
 
 moab::LloydSmoother lloyd( mbi, 0, cells, 0, 0 /*fixed*/ );
 
 lloyd.perform_smooth();
 
 // run_local_smoother2(mesh_and_domain, err, OF_value);
 
 // Construct a MeshTSTT from the file
 int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" );
 if( ierr ) return 1;
 // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
 // if (err) return 1;
 cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
 return 0;
}
 
int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
{
 // double OF_value = 1e-5;
 
 Mesh* mesh = mesh_and_domain.get_mesh();
 MeshDomain* domain = mesh_and_domain.get_domain();
 
 // creates an intruction queue
 InstructionQueue queue1;
 
 // creates a mean ratio quality metric ...
 // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
 ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
 if( err ) return 1;
 // mean_ratio->set_gradient_type(QualityMetric::NUMERICAL_GRADIENT);
 // mean_ratio->set_hessian_type(QualityMetric::NUMERICAL_HESSIAN);
 mean_ratio->set_averaging_method( QualityMetric::RMS, err );
 if( err ) return 1;
 
 // ... and builds an objective function with it
 LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
 if( err ) return 1;
 
 if( false )
 {
 InstructionQueue qtmp;
 Randomize rand( -0.005 );
 TerminationCriterion sc_rand;
 sc_rand.add_iteration_limit( 2 );
 rand.set_outer_termination_criterion( &sc_rand );
 qtmp.set_master_quality_improver( &rand, err );
 if( err ) return 1;
 if( domain )
 {
 qtmp.run_instructions( &mesh_and_domain, err );
 }
 else
 {
 qtmp.run_instructions( mesh, err );
 }
 if( err ) return 1;
 }
 
 // creates the smart laplacian optimization procedures
 SmartLaplacianSmoother* pass1 = new SmartLaplacianSmoother( obj_func );
 // SteepestDescent* pass1 = new SteepestDescent( obj_func );
 
 QualityAssessor stop_qa( mean_ratio );
 
 // **************Set stopping criterion****************
 TerminationCriterion tc_inner;
 tc_inner.add_absolute_vertex_movement( OF_value );
 TerminationCriterion tc_outer;
 tc_outer.add_iteration_limit( 10 );
 pass1->set_inner_termination_criterion( &tc_inner );
 pass1->set_outer_termination_criterion( &tc_outer );
 
 queue1.add_quality_assessor( &stop_qa, err );
 if( err ) return 1;
 
 // adds 1 pass of pass1 to mesh_set
 queue1.set_master_quality_improver( pass1, err );
 if( err ) return 1;
 
 queue1.add_quality_assessor( &stop_qa, err );
 if( err ) return 1;
 
 // launches optimization on mesh_set
 if( domain )
 {
 queue1.run_instructions( &mesh_and_domain, err );
 }
 else
 {
 queue1.run_instructions( mesh, err );
 }
 if( err ) return 1;
 
 // Construct a MeshTSTT from the file
 int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk" );
 if( ierr ) return 1;
 // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
 // if (err) return 1;
 cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
 
 // print_timing_diagnostics( cout );
 return 0;
}
 
int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
{
 Mesh* mesh = mesh_and_domain.get_mesh();
 MeshDomain* domain = mesh_and_domain.get_domain();
 
 // creates an intruction queue
 InstructionQueue q;
 
 // do optimization
 const double eps = 0.01;
 IdealShapeTarget w;
 TShapeB1 mu;
 TQualityMetric metric( &w, &mu );
 PMeanPTemplate func( 1.0, &metric );
 
 SteepestDescent solver( &func );
 solver.use_element_on_vertex_patch();
 solver.do_jacobi_optimization();
 
 TerminationCriterion inner, outer;
 inner.add_absolute_vertex_movement( 0.5 * eps );
 outer.add_absolute_vertex_movement( 0.5 * eps );
 
 QualityAssessor qa( &metric );
 
 q.add_quality_assessor( &qa, err );
 if( err ) return 1;
 q.set_master_quality_improver( &solver, err );
 if( err ) return 1;
 q.add_quality_assessor( &qa, err );
 if( err ) return 1;
 
 // launches optimization on mesh_set
 if( domain )
 {
 q.run_instructions( &mesh_and_domain, err );
 }
 else
 {
 q.run_instructions( mesh, err );
 }
 if( err ) return 1;
 
 // Construct a MeshTSTT from the file
 int ierr = write_vtk_mesh( mesh, "ideal-shape-result.vtk" );
 if( ierr ) return 1;
 // MeshWriter::write_vtk(mesh, "ideal-shape-result.vtk", err);
 // if (err) return 1;
 cout << "Wrote \"ideal-shape-result.vtk\"" << endl;
 
 print_timing_diagnostics( cout );
 return 0;
}
 
int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
{
 double OF_value = 0.01;
 
 Mesh* mesh = mesh_and_domain.get_mesh();
 MeshDomain* domain = mesh_and_domain.get_domain();
 
 // creates an intruction queue
 InstructionQueue queue1;
 
 // creates a mean ratio quality metric ...
 // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
 ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
 // AspectRatioGammaQualityMetric* mean_ratio = new AspectRatioGammaQualityMetric();
 
 // ElementSolIndQM* solindqm = new ElementSolIndQM(solution_indicator);
 // MultiplyQualityMetric* mean_ratio = new MultiplyQualityMetric(new
 // VertexConditionNumberQualityMetric(), solindqm, err);
 // ElementSolIndQM* mean_ratio = solindqm;
 
 // LocalSizeQualityMetric* mean_ratio = new LocalSizeQualityMetric();
 
 mean_ratio->set_averaging_method( QualityMetric::SUM_OF_RATIOS_SQUARED, err );
 if( err ) return 1;
 
 // ... and builds an objective function with it
 LPtoPTemplate* obj_func = new LPtoPTemplate( mean_ratio, 1, err );
 if( err ) return 1;
 
 // // creates the feas newt optimization procedures
 ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
 // QuasiNewton* pass1 = new QuasiNewton( obj_func );
 // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
 pass1->use_global_patch();
 
 QualityAssessor stop_qa( mean_ratio );
 
 // **************Set stopping criterion****************
 TerminationCriterion tc_inner;
 tc_inner.add_absolute_vertex_movement( OF_value );
 if( err ) return 1;
 TerminationCriterion tc_outer;
 tc_inner.add_iteration_limit( 20 );
 tc_outer.add_iteration_limit( 5 );
 pass1->set_inner_termination_criterion( &tc_inner );
 pass1->set_outer_termination_criterion( &tc_outer );
 pass1->set_debugging_level( 3 );
 
 queue1.add_quality_assessor( &stop_qa, err );
 if( err ) return 1;
 
 // adds 1 pass of pass1 to mesh_set1
 queue1.set_master_quality_improver( pass1, err );
 if( err ) return 1;
 
 queue1.add_quality_assessor( &stop_qa, err );
 if( err ) return 1;
 
 // launches optimization on mesh_set
 if( domain )
 {
 queue1.run_instructions( &mesh_and_domain, err );
 }
 else
 {
 queue1.run_instructions( mesh, err );
 }
 if( err ) return 1;
 
 // Construct a MeshTSTT from the file
 int ierr = write_vtk_mesh( mesh, "solution-mesh-result.vtk" );
 if( ierr ) return 1;
 // MeshWriter::write_vtk(mesh, "solution-mesh-result.vtk", err);
 // if (err) return 1;
 cout << "Wrote \"solution-mesh-result.vtk\"" << endl;
 
 print_timing_diagnostics( cout );
 return 0;
}
 
int get_imesh_mesh( MBMesquite::Mesh** mesh, const char* file_name, int dimension )
{
 int err;
 iMesh_Instance instance = 0;
 iMesh_newMesh( NULL, &instance, &err, 0 );
 CHECK_IMESH( "Creation of mesh instance failed" );
 
 iBase_EntitySetHandle root_set;
 iMesh_getRootSet( instance, &root_set, &err );
 CHECK_IMESH( "Could not get root set" );
 
 iMesh_load( instance, root_set, file_name, 0, &err, strlen( file_name ), 0 );
 CHECK_IMESH( "Could not load mesh from file" );
 
 iBase_TagHandle fixed_tag;
 iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) );
 // if (iBase_SUCCESS != err)
 {
 // get the skin vertices of those cells and mark them as fixed; we don't want to fix the
 // vertices on a part boundary, but since we exchanged a layer of ghost cells, those
 // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move
 // those anyway use MOAB's skinner class to find the skin
 
 // get all vertices and cells
 // make tag to specify fixed vertices, since it's input to the algorithm; use a default
 // value of non-fixed so we only need to set the fixed tag for skin vertices
 moab::Interface* mbi = reinterpret_cast< MBiMesh* >( instance )->mbImpl;
 moab::EntityHandle currset = 0;
 moab::Tag fixed;
 int def_val = 0;
 err = 0;
 moab::ErrorCode rval = mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed,
 moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
 moab::Range verts, cells, skin_verts;
 rval = mbi->get_entities_by_type( currset, moab::MBVERTEX, verts );MB_CHK_SET_ERR( rval, "Querying vertices failed" );
 rval = mbi->get_entities_by_dimension( currset, dimension, cells );MB_CHK_SET_ERR( rval, "Querying elements failed" );
 std::cout << "Found " << verts.size() << " vertices and " << cells.size() << " elements" << std::endl;
 
 moab::Skinner skinner( mbi );
 rval = skinner.find_skin( currset, cells, true, skin_verts );MB_CHK_SET_ERR( rval,
 "Finding the skin of the mesh failed" ); // 'true' param indicates we want
 // vertices back, not cells
 
 std::vector< int > fix_tag( skin_verts.size(), 1 ); // initialized to 1 to indicate fixed
 rval = mbi->tag_set_data( fixed, skin_verts, &fix_tag[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" );
 std::cout << "Found " << skin_verts.size() << " vertices on the skin of the domain." << std::endl;
 
 // fix_tag.resize(verts.size(),0);
 // rval = mbi->tag_get_data(fixed, verts, &fix_tag[0]); MB_CHK_SET_ERR(rval, "Getting tag
 // data failed");
 
 iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen( "fixed" ) );
 CHECK_IMESH( "Getting tag handle (fixed) failed" );
 
 // Set some arbitrary solution indicator
 moab::Tag solindTag;
 double def_val_dbl = 0.0;
 rval = mbi->tag_get_handle( "solution_indicator", 1, moab::MB_TYPE_DOUBLE, solindTag,
 moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val_dbl );MB_CHK_SET_ERR( rval, "Getting tag handle failed" );
 solution_indicator.resize( cells.size(), 0.01 );
 for( unsigned i = 0; i < cells.size() / 4; i++ )
 solution_indicator[i] = 0.1;
 for( unsigned i = cells.size() / 4; i < 2 * cells.size() / 4; i++ )
 solution_indicator[i] = 0.5;
 for( unsigned i = 2 * cells.size() / 4; i < 3 * cells.size() / 4; i++ )
 solution_indicator[i] = 0.5;
 for( unsigned i = 3 * cells.size() / 4; i < cells.size(); i++ )
 solution_indicator[i] = 0.5;
 
 rval = mbi->tag_set_data( solindTag, cells, &solution_indicator[0] );MB_CHK_SET_ERR( rval, "Setting tag data failed" );
 }
 
 MsqError ierr;
 MBMesquite::MsqIMesh* imesh =
 new MBMesquite::MsqIMesh( instance, root_set, ( dimension == 3 ? iBase_REGION : iBase_FACE ), ierr,
 &fixed_tag );
 if( MSQ_CHKERR( ierr ) )
 {
 delete imesh;
 cerr << err << endl;
 err = iBase_FAILURE;
 CHECK_IMESH( "Creation of MsqIMesh instance failed" );
 return 0;
 }
 
 *mesh = imesh;
 return iBase_SUCCESS;
}
 
int get_native_mesh( MBMesquite::Mesh** mesh, const char* file_name, int )
{
 MsqError err;
 MBMesquite::MeshImpl* imesh = new MBMesquite::MeshImpl();
 imesh->read_vtk( file_name, err );
 if( err )
 {
 cerr << err << endl;
 exit( 3 );
 }
 *mesh = imesh;
 
 return iBase_SUCCESS;
}
 
int get_itaps_domain( MeshDomain** odomain, const char* filename )
{
 
 if( filename == 0 || strlen( filename ) == 0 )
 {
 *odomain = new PlanarDomain( PlanarDomain::XY );
 return 0;
 }
 
 int err;
 iGeom_Instance geom;
 iGeom_newGeom( "", &geom, &err, 0 );
 CHECK_IGEOM( "ERROR: iGeom creation failed" );
 
#ifdef MOAB_HAVE_CGM_FACET
 FacetModifyEngine::set_modify_enabled( CUBIT_TRUE );
#endif
 
 iGeom_load( geom, filename, 0, &err, strlen( filename ), 0 );
 CHECK_IGEOM( "Cannot load the geometry" );
 
 iBase_EntitySetHandle root_set;
 iGeom_getRootSet( geom, &root_set, &err );
 CHECK_IGEOM( "getRootSet failed!" );
 
 // print out the number of entities
 std::cout << "Model contents: " << std::endl;
 const char* gtype[] = { "vertices: ", "edges: ", "faces: ", "regions: " };
 int nents[4];
 for( int i = 0; i <= 3; ++i )
 {
 iGeom_getNumOfType( geom, root_set, i, &nents[i], &err );
 CHECK_IGEOM( "Error: problem getting entities after gLoad." );
 std::cout << gtype[i] << nents[i] << std::endl;
 }
 
 iBase_EntityHandle* hd_geom_ents;
 int csize = 0, sizealloc = 0;
 if( nents[3] > 0 )
 {
 hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[2] );
 csize = nents[2];
 iGeom_getEntities( geom, root_set, 2, &hd_geom_ents, &csize, &sizealloc, &err );
 }
 else
 {
 hd_geom_ents = (iBase_EntityHandle*)malloc( sizeof( iBase_EntityHandle ) * nents[1] );
 csize = nents[1];
 iGeom_getEntities( geom, root_set, 1, &hd_geom_ents, &csize, &sizealloc, &err );
 }
 CHECK_IGEOM( "ERROR: Could not get entities" );
 
 *odomain = new MsqIGeom( geom, hd_geom_ents[0] );
 return iBase_SUCCESS;
}