OptimizeMeshMesquite.cpp

This example demonstrates mesh optimization using the Mesquite library. It shows how to load meshes using both native MOAB and iMesh interfaces, apply various mesh optimization algorithms from Mesquite, use different quality metrics (aspect ratio, condition number, etc.), perform global and local mesh smoothing, handle geometry constraints using domain definitions, and measure mesh quality before and after optimization. The example demonstrates integration between MOAB and Mesquite for advanced mesh optimization and quality improvement.

/** @example OptimizeMeshMesquite.cpp
 * This example demonstrates mesh optimization using the Mesquite library.
 * It shows how to load meshes using both native MOAB and iMesh interfaces,
 * apply various mesh optimization algorithms from Mesquite,
 * use different quality metrics (aspect ratio, condition number, etc.),
 * perform global and local mesh smoothing,
 * handle geometry constraints using domain definitions,
 * and measure mesh quality before and after optimization.
 * The example demonstrates integration between MOAB and Mesquite
 * for advanced mesh optimization and quality improvement.
 */
 
#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;
    MB_CHK_SET_ERR( mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed ), "Getting tag handle failed" );
    moab::Range cells;
    MB_CHK_SET_ERR( mbi->get_entities_by_dimension( 0, 2, cells ), "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::MB_CHK_SET_ERR( mbi->tag_get_handle( "fixed", 1, moab::MB_TYPE_INTEGER, fixed,
                                                   moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val ),
                              "Getting tag handle failed" );
        moab::Range verts, cells, skin_verts;
        MB_CHK_SET_ERR( mbi->get_entities_by_type( currset, moab::MBVERTEX, verts ), "Querying vertices failed" );
        MB_CHK_SET_ERR( mbi->get_entities_by_dimension( currset, dimension, cells ), "Querying elements failed" );
        std::cout << "Found " << verts.size() << " vertices and " << cells.size() << " elements" << std::endl;
 
        moab::Skinner skinner( mbi );
        MB_CHK_SET_ERR( skinner.find_skin( currset, cells, true, skin_verts ),
                        "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
        MB_CHK_SET_ERR( mbi->tag_set_data( fixed, skin_verts, &fix_tag[0] ), "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);
        // MB_CHK_SET_ERR( mbi->tag_get_data(fixed, verts, &fix_tag[0]), "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;
        MB_CHK_SET_ERR( mbi->tag_get_handle( "solution_indicator", 1, moab::MB_TYPE_DOUBLE, solindTag,
                                             moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val_dbl ),
                        "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;
 
        MB_CHK_SET_ERR( mbi->tag_set_data( solindTag, cells, &solution_indicator[0] ), "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;
}