cub2h5m.cpp

Go to the documentation of this file.

#include "moab/GeomQueryTool.hpp"
#include "moab/GeomTopoTool.hpp"
#include "InitCGMA.hpp"
#include "CGMApp.hpp"
#include "moab/CN.hpp"
#include "moab/Core.hpp"
#include "moab/CartVect.hpp"
#include "moab/FileOptions.hpp"
#include "moab/Skinner.hpp"
#include "quads_to_tris.hpp"
#include <limits>
#include <cstdlib>
#include <sstream>
#include <ctime>

#define GF_CUBIT_FILE_TYPE    "CUBIT"
#define GF_STEP_FILE_TYPE     "STEP"
#define GF_IGES_FILE_TYPE     "IGES"
#define GF_ACIS_TXT_FILE_TYPE "ACIS_SAT"
#define GF_ACIS_BIN_FILE_TYPE "ACIS_SAB"
#define GF_OCC_BREP_FILE_TYPE "OCC"

using namespace moab;

void tokenize( const std::string& str, std::vector< std::string >& tokens, const char* delimiters )
{
    std::string::size_type last = str.find_first_not_of( delimiters, 0 );
    std::string::size_type pos  = str.find_first_of( delimiters, last );
    if( std::string::npos == pos )
        tokens.push_back( str );
    else
        while( std::string::npos != pos && std::string::npos != last )
        {
            tokens.push_back( str.substr( last, pos - last ) );
            last = str.find_first_not_of( delimiters, pos );
            pos  = str.find_first_of( delimiters, last );
            if( std::string::npos == pos ) pos = str.size();
        }
}

ErrorCode get_group_names( Interface* MBI,
                           const EntityHandle group_set,
                           const Tag nameTag,
                           std::vector< std::string >& grp_names )
{
    // get names
    char name0[NAME_TAG_SIZE];
    std::fill( name0, name0 + NAME_TAG_SIZE, '\0' );
    ErrorCode result = MBI->tag_get_data( nameTag, &group_set, 1, &name0 );
    if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return MB_FAILURE;

    if( MB_TAG_NOT_FOUND != result ) grp_names.push_back( std::string( name0 ) );

    return MB_SUCCESS;
}

// For each material, sum the volume. If the coordinates were updated for
// deformation, summarize the volume change.
ErrorCode summarize_cell_volume_change( Interface* MBI,
                                        const EntityHandle cgm_file_set,
                                        const Tag categoryTag,
                                        const Tag dimTag,
                                        const Tag sizeTag,
                                        const Tag nameTag,
                                        const Tag idTag,
                                        const bool conserve_mass,
                                        const bool debug )
{
    // get groups
    ErrorCode rval;
    const char group_category[CATEGORY_TAG_SIZE] = { "Group\0" };
    const void* const group_val[]                = { &group_category };
    Range groups;
    rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &categoryTag, group_val, 1, groups );
    if( MB_SUCCESS != rval ) return rval;

    // Get the maximum group id. This is so that new groups do not have
    // duplicate ids.
    int max_grp_id = -1;
    for( Range::const_iterator i = groups.begin(); i != groups.end(); ++i )
    {
        int grp_id;
        rval = MBI->tag_get_data( idTag, &( *i ), 1, &grp_id );
        if( MB_SUCCESS != rval ) return rval;
        if( max_grp_id < grp_id ) max_grp_id = grp_id;
    }
    if( conserve_mass )
    {
        std::cout << "  Altering group densities to conserve mass for each volume." << std::endl;
        std::cout << "    maximum group id=" << max_grp_id << std::endl;
    }

    for( Range::const_iterator i = groups.begin(); i != groups.end(); ++i )
    {
        // get group names
        std::vector< std::string > grp_names;
        rval = get_group_names( MBI, *i, nameTag, grp_names );
        if( MB_SUCCESS != rval ) return MB_FAILURE;

        // determine if it is a material group
        bool material_grp = false;
        int mat_id        = -1;
        double rho        = 0;
        for( std::vector< std::string >::const_iterator j = grp_names.begin(); j != grp_names.end(); ++j )
        {
            if( std::string::npos != ( *j ).find( "mat" ) && std::string::npos != ( *j ).find( "rho" ) )
            {
                material_grp = true;
                std::cout << "  material group: " << *j << std::endl;

                // get the density and material id
                std::vector< std::string > tokens;
                tokenize( *j, tokens, "_" );
                mat_id = atoi( tokens[1].c_str() );
                rho    = atof( tokens[3].c_str() );
            }
        }
        if( !material_grp ) continue;

        // get the volume sets of the material group
        const int three               = 3;
        const void* const three_val[] = { &three };
        Range vols;
        rval = MBI->get_entities_by_type_and_tag( *i, MBENTITYSET, &dimTag, three_val, 1, vols );
        if( MB_SUCCESS != rval ) return rval;

        // for each volume, sum predeformed and deformed volume
        double orig_grp_volume = 0, defo_grp_volume = 0;

        moab::GeomTopoTool gtt  = moab::GeomTopoTool( MBI, false );
        moab::GeomQueryTool gqt = moab::GeomQueryTool( &gtt );
        for( Range::const_iterator j = vols.begin(); j != vols.end(); ++j )
        {
            double defo_size = 0, orig_size = 0;
            rval = gqt.measure_volume( *j, defo_size );
            if( MB_SUCCESS != rval ) return rval;
            defo_grp_volume += defo_size;
            rval = MBI->tag_get_data( sizeTag, &( *j ), 1, &orig_size );
            if( MB_SUCCESS != rval ) return rval;
            orig_grp_volume += orig_size;

            // calculate a new density to conserve mass through the deformation
            if( !conserve_mass ) continue;
            double new_rho = rho * orig_size / defo_size;

            // create a group for the volume with modified density
            EntityHandle new_grp;
            rval = MBI->create_meshset( MESHSET_SET, new_grp );
            if( MB_SUCCESS != rval ) return rval;
            std::stringstream new_name_ss;
            new_name_ss << "mat_" << mat_id << "_rho_" << new_rho << "\0";
            std::string new_name;
            new_name_ss >> new_name;
            rval = MBI->tag_set_data( nameTag, &new_grp, 1, new_name.c_str() );
            if( MB_SUCCESS != rval ) return rval;
            max_grp_id++;
            rval = MBI->tag_set_data( idTag, &new_grp, 1, &max_grp_id );
            if( MB_SUCCESS != rval ) return rval;
            const char group_category2[CATEGORY_TAG_SIZE] = "Group\0";
            rval = MBI->tag_set_data( categoryTag, &new_grp, 1, group_category2 );
            if( MB_SUCCESS != rval ) return rval;

            // add the volume to the new group
            rval = MBI->add_entities( new_grp, &( *j ), 1 );
            if( MB_SUCCESS != rval ) return rval;

            // add the new grp to the cgm_file_set
            rval = MBI->add_entities( cgm_file_set, &new_grp, 1 );
            if( MB_SUCCESS != rval ) return rval;

            // remove the volume from the old group
            rval = MBI->remove_entities( *i, &( *j ), 1 );
            if( MB_SUCCESS != rval ) return rval;
            if( debug ) std::cout << "    new group: " << new_name << " id=" << max_grp_id << std::endl;
        }

        std::cout << "    orig_volume=" << orig_grp_volume << " defo_volume=" << defo_grp_volume
                  << " defo/orig=" << defo_grp_volume / orig_grp_volume << std::endl;
    }

    return MB_SUCCESS;
}

// We cannot build an OBB tree if all of a volume's surfaces have no facets.
// To prevent this, remove the cgm surface set if the cub surface set exists,
// but had its faced removed (due to dead elements). Remember that the cgm_file_set
// is not TRACKING.
ErrorCode remove_empty_cgm_surfs_and_vols( Interface* MBI,
                                           const EntityHandle cgm_file_set,
                                           const Tag idTag,
                                           const Tag dimTag,
                                           const bool /*debug */ )
{

    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int n_tris;
        result = MBI->get_number_entities_by_type( *i, MBTRI, n_tris );
        if( MB_SUCCESS != result ) return result;

        if( 0 == n_tris )
        {
            int surf_id;
            result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
            assert( MB_SUCCESS == result );

            Range parent_vols;
            result = MBI->get_parent_meshsets( *i, parent_vols );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = parent_vols.begin(); j != parent_vols.end(); ++j )
            {
                result = MBI->remove_parent_child( *j, *i );
                assert( MB_SUCCESS == result );
            }
            Range child_curves;
            result = MBI->get_child_meshsets( *i, child_curves );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = child_curves.begin(); j != child_curves.end(); ++j )
            {
                result = MBI->remove_parent_child( *i, *j );
                assert( MB_SUCCESS == result );
            }
            result = MBI->remove_entities( cgm_file_set, &( *i ), 1 );
            assert( MB_SUCCESS == result );

            // Is the set contained anywhere else? If the surface is in a CUBIT group,
            // such as "unmerged_surfs" it will cause write_mesh to fail. This should
            // be a MOAB bug.
            Range all_sets;
            result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = all_sets.begin(); j != all_sets.end(); ++j )
            {
                if( MBI->contains_entities( *j, &( *i ), 1 ) )
                {
                    result = MBI->remove_entities( *j, &( *i ), 1 );
                    assert( MB_SUCCESS == result );
                }
            }

            result = MBI->delete_entities( &( *i ), 1 );
            assert( MB_SUCCESS == result );
            std::cout << "  Surface " << surf_id << ": removed because all of its mesh faces have been removed"
                      << std::endl;
        }
    }

    // get all volumes
    const int three               = 3;
    const void* const three_val[] = { &three };
    Range cgm_vols;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, cgm_vols );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator i = cgm_vols.begin(); i != cgm_vols.end(); ++i )
    {
        // get the volume's number of surfaces
        int n_surfs;
        result = MBI->num_child_meshsets( *i, &n_surfs );
        assert( MB_SUCCESS == result );

        // if no surfaces remain, remove the volume
        if( 0 == n_surfs )
        {
            int vol_id;
            result = MBI->tag_get_data( idTag, &( *i ), 1, &vol_id );
            assert( MB_SUCCESS == result );
            // Is the set contained anywhere else? If the surface is in a CUBIT group,
            // such as "unmerged_surfs" it will cause write_mesh to fail. This should
            // be a MOAB bug.
            Range all_sets;
            result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = all_sets.begin(); j != all_sets.end(); ++j )
            {
                if( MBI->contains_entities( *j, &( *i ), 1 ) )
                {
                    result = MBI->remove_entities( *j, &( *i ), 1 );
                    assert( MB_SUCCESS == result );
                }
            }
            result = MBI->delete_entities( &( *i ), 1 );
            assert( MB_SUCCESS == result );
            std::cout << "  Volume " << vol_id << ": removed because all of its surfaces have been removed"
                      << std::endl;
        }
    }
    return MB_SUCCESS;
}

// Given parent volume senses, an id, and a set handle, this function creates a
// new surface set with dimension, geometry category, id, and sense tags.
ErrorCode build_new_surface( Interface* MBI,
                             EntityHandle& new_surf,
                             const EntityHandle forward_parent_vol,
                             const EntityHandle reverse_parent_vol,
                             const int new_surf_id,
                             const Tag dimTag,
                             const Tag idTag,
                             const Tag categoryTag,
                             const Tag senseTag )
{

    ErrorCode result;
    result = MBI->create_meshset( 0, new_surf );
    if( MB_SUCCESS != result ) return result;
    if( 0 != forward_parent_vol )
    {
        result = MBI->add_parent_child( forward_parent_vol, new_surf );
        if( MB_SUCCESS != result ) return result;
    }
    if( 0 != reverse_parent_vol )
    {
        result = MBI->add_parent_child( reverse_parent_vol, new_surf );
        if( MB_SUCCESS != result ) return result;
    }
    const int two = 2;
    result        = MBI->tag_set_data( dimTag, &new_surf, 1, &two );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_set_data( idTag, &new_surf, 1, &new_surf_id );
    if( MB_SUCCESS != result ) return result;
    const char geom_category[CATEGORY_TAG_SIZE] = { "Surface\0" };
    result                                      = MBI->tag_set_data( categoryTag, &new_surf, 1, &geom_category );
    if( MB_SUCCESS != result ) return result;
    EntityHandle vols[2] = { forward_parent_vol, reverse_parent_vol };
    result               = MBI->tag_set_data( senseTag, &new_surf, 1, vols );
    if( MB_SUCCESS != result ) return result;

    return MB_SUCCESS;
}

// Given a face, orient it outward wrt its adjacent mesh element.
// Each face must be adjacent to exactly one mesh element.
ErrorCode orient_faces_outward( Interface* MBI, const Range& faces, const bool /*debug*/ )
{

    ErrorCode result;
    for( Range::const_iterator i = faces.begin(); i != faces.end(); ++i )
    {
        Range adj_elem;
        result = MBI->get_adjacencies( &( *i ), 1, 3, false, adj_elem );
        if( MB_SUCCESS != result ) return result;
        if( 1 != adj_elem.size() ) return MB_INVALID_SIZE;

        // get connectivity for element and face
        const EntityHandle* elem_conn;
        int elem_n_nodes;
        result = MBI->get_connectivity( adj_elem.front(), elem_conn, elem_n_nodes );
        if( MB_SUCCESS != result ) return result;
        const EntityHandle* face_conn;
        int face_n_nodes;
        result = MBI->get_connectivity( *i, face_conn, face_n_nodes );
        if( MB_SUCCESS != result ) return result;

        // Get the sense of the face wrt the element
        EntityType elem_type = MBI->type_from_handle( adj_elem.front() );
        EntityType face_type = MBI->type_from_handle( *i );
        int face_num, offset;
        int sense          = 0;
        const int face_dim = CN::Dimension( face_type );
        int rval = CN::SideNumber( elem_type, elem_conn, face_conn, face_n_nodes, face_dim, face_num, sense, offset );
        if( 0 != rval ) return MB_FAILURE;

        // If the face is not oriented outward wrt the element, reverse it
        if( -1 == sense )
        {
            EntityHandle new_face_conn[4] = { face_conn[3], face_conn[2], face_conn[1], face_conn[0] };
            result                        = MBI->set_connectivity( *i, new_face_conn, 4 );
            if( MB_SUCCESS != result ) return result;
        }
    }
    return MB_SUCCESS;
}

/* Isolate dead code in a preproc-killed block - sjackson 11/22/10 */
#if 0

/* qsort int comparison function */
int handle_compare(const void *a, const void *b)
{
  const EntityHandle *ia = (const EntityHandle *)a; // casting pointer types
  const EntityHandle *ib = (const EntityHandle *)b;
  return *ia - *ib;
  /* integer comparison: returns negative if b > a
   and positive if a > b */
}

// qsort face comparison function. assume each face has 4 nodes
int compare_face(const void *a, const void *b)
{
  EntityHandle *ia = (EntityHandle *)a;
  EntityHandle *ib = (EntityHandle *)b;
  if(*ia == *ib)
  {
    if(*(ia+1) == *(ib+1))
    {
      if(*(ia+2) == *(ib+2))
      {
        return (int)(*(ia+3) - *(ib+3));
      }
      else
      {
        return (int)(*(ia+2) - *(ib+2));
      }
    }
    else
    {
      return (int)(*(ia+1) - *(ib+1));
    }
  }
  else
  {
    return (int)(*ia - *ib);
  }
}

// Use this to get quad faces from hex elems.
ErrorCode skin_hex_elems(Interface *MBI, Range elems, const int dim,
    Range &faces )
{
  // get faces of each hex
  const int nodes_per_face = 4;
  const unsigned int faces_per_elem = 6;
  unsigned int n_faces = faces_per_elem*elems.size();
  EntityHandle f[n_faces][nodes_per_face];
  ErrorCode result;
  int counter = 0;
  for(Range::iterator i=elems.begin(); i!=elems.end(); ++i)
  {
    Range elem_faces;
    result = MBI->get_adjacencies( &(*i), 1, 2, true, elem_faces );
    if(MB_SUCCESS != result) return result;
    if(faces_per_elem != elem_faces.size()) return MB_INVALID_SIZE;
    for(Range::iterator j=elem_faces.begin(); j!=elem_faces.end(); ++j)
    {
      const EntityHandle *conn;
      int n_nodes;
      ErrorCode result = MBI->get_connectivity( *j, conn, n_nodes );
      if(MB_SUCCESS != result) return result;
      if(nodes_per_face != n_nodes) return MB_INVALID_SIZE;
      // Sort the node handles of the face
      for(int k=0; k<nodes_per_face; ++k) f[counter][k] = conn[k];
      qsort( &f[counter][0], nodes_per_face, sizeof(EntityHandle), handle_compare );
      ++counter;
    }
  }

  // Sort the faces by the first node handle, then second node, then third node...
  qsort( &f[0][0], n_faces, nodes_per_face*sizeof(EntityHandle), compare_face );

  // if a face has 1 or more duplicates, it is not on the skin
  faces.clear();
  for(unsigned int i=0; i<n_faces; ++i)
  {
    // if the last face is tested, it must be on the skin
    if(n_faces-1 == i)
    {
      Range face_handle;
      result = MBI->get_adjacencies( &(f[i][0]), nodes_per_face, 2, false, face_handle );
      if(MB_SUCCESS != result) return result;
      if(1 != face_handle.size()) return MB_INVALID_SIZE;
      faces.insert( face_handle.front() );
      // Due to sort, if a duplicate exists it must be next
    }
    else if( f[i][0]==f[i+1][0] && f[i][1]==f[i+1][1] &&
        f[i][2]==f[i+1][2] && f[i][3]==f[i+1][3] )
    {
      ++i;
      while( f[i][0]==f[i+1][0] && f[i][1]==f[i+1][1] &&
          f[i][2]==f[i+1][2] && f[i][3]==f[i+1][3] )
      {
        std::cout << "    skin WARNING: non-manifold face" << std::endl;
        ++i;
      }
      // otherwise it is on the skin
    }
    else
    {
      Range face_handle;
      result = MBI->get_adjacencies( &(f[i][0]), nodes_per_face, 2, false, face_handle );
      if(MB_SUCCESS != result) return result;
      if(1 != face_handle.size()) return MB_INVALID_SIZE;
      faces.insert( face_handle.front() );
    }
  }

  return MB_SUCCESS;
}

#endif  // dead code isolation
// Given a 1D array of data, axis labels, title, and number of bins, create a
// histogram.
void plot_histogram( const std::string& title,
                     const std::string& x_axis_label,
                     const std::string& y_axis_label,
                     const int n_bins,
                     const double data[],
                     const int n_data )
{
    // find max and min
    double min = std::numeric_limits< double >::max();
    double max = -std::numeric_limits< double >::max();
    for( int i = 0; i < n_data; ++i )
    {
        if( min > data[i] ) min = data[i];
        if( max < data[i] ) max = data[i];
    }

    // make bins for histogram
    double bin_width = ( max - min ) / n_bins;
    std::vector< int > bins( n_bins );
    for( int i = 0; i < n_bins; ++i )
        bins[i] = 0;

    // fill the bins
    for( int i = 0; i < n_data; ++i )
    {
        double diff = data[i] - min;
        int bin     = diff / bin_width;
        if( 9 < bin ) bin = 9;  // cheap fix for numerical precision error
        if( 0 > bin ) bin = 0;  // cheap fix for numerical precision error
        ++bins[bin];
    }

    // create bars
    int max_bin = 0;
    for( int i = 0; i < n_bins; ++i )
        if( max_bin < bins[i] ) max_bin = bins[i];
    int bar_height;
    int max_bar_chars = 72;
    std::vector< std::string > bars( n_bins );
    for( int i = 0; i < n_bins; ++i )
    {
        bar_height = ( max_bar_chars * bins[i] ) / max_bin;
        for( int j = 0; j < bar_height; ++j )
            bars[i] += "*";
    }

    // print histogram header
    std::cout << std::endl;
    std::cout << "                                 " << title << std::endl;

    // print results
    std::cout.width( 15 );
    std::cout << min << std::endl;
    for( int i = 0; i < n_bins; ++i )
    {
        std::cout.width( 15 );
        std::cout << min + ( ( i + 1 ) * bin_width );
        std::cout.width( max_bar_chars );
        std::cout << bars[i] << bins[i] << std::endl;
    }

    // print histogram footer
    std::cout.width( 15 );
    std::cout << y_axis_label;
    std::cout.width( max_bar_chars );
    std::cout << " " << x_axis_label << std::endl;
    std::cout << std::endl;
}

// This is a helper function that creates data and labels for histograms.
void generate_plots( const double orig[], const double defo[], const int n_elems, const std::string& time_step )
{

    // find volume ratio then max and min
    double* ratio = new double[n_elems];
    for( int i = 0; i < n_elems; ++i )
        ratio[i] = ( defo[i] - orig[i] ) / orig[i];

    plot_histogram( "Predeformed Element Volume", "Num_Elems", "Volume [cc]", 10, orig, n_elems );
    std::string title = "Element Volume Change Ratio at Time Step " + time_step;
    plot_histogram( title, "Num_Elems", "Volume Ratio", 10, ratio, n_elems );
    delete[] ratio;
}

// Given four nodes, calculate the tet volume.
inline static double tet_volume( const CartVect& v0, const CartVect& v1, const CartVect& v2, const CartVect& v3 )
{
    return 1. / 6. * ( ( ( v1 - v0 ) * ( v2 - v0 ) ) % ( v3 - v0 ) );
}

// Measure and tet volume are taken from measure.cpp
double measure( Interface* MBI, const EntityHandle element )
{
    EntityType type = MBI->type_from_handle( element );

    const EntityHandle* conn;
    int num_vertices;
    ErrorCode result = MBI->get_connectivity( element, conn, num_vertices );
    if( MB_SUCCESS != result ) return result;

    std::vector< CartVect > coords( num_vertices );
    result = MBI->get_coords( conn, num_vertices, coords[0].array() );
    if( MB_SUCCESS != result ) return result;

    switch( type )
    {
        case MBEDGE:
            return ( coords[0] - coords[1] ).length();
        case MBTRI:
            return 0.5 * ( ( coords[1] - coords[0] ) * ( coords[2] - coords[0] ) ).length();
        case MBQUAD:
            num_vertices = 4;
        case MBPOLYGON: {
            CartVect mid( 0, 0, 0 );
            for( int i = 0; i < num_vertices; ++i )
                mid += coords[i];
            mid /= num_vertices;

            double sum = 0.0;
            for( int i = 0; i < num_vertices; ++i )
            {
                int j = ( i + 1 ) % num_vertices;
                sum += ( ( mid - coords[i] ) * ( mid - coords[j] ) ).length();
            }
            return 0.5 * sum;
        }
        case MBTET:
            return tet_volume( coords[0], coords[1], coords[2], coords[3] );
        case MBPYRAMID:
            return tet_volume( coords[0], coords[1], coords[2], coords[4] ) +
                   tet_volume( coords[0], coords[2], coords[3], coords[4] );
        case MBPRISM:
            return tet_volume( coords[0], coords[1], coords[2], coords[5] ) +
                   tet_volume( coords[3], coords[5], coords[4], coords[0] ) +
                   tet_volume( coords[1], coords[4], coords[5], coords[0] );
        case MBHEX:
            return tet_volume( coords[0], coords[1], coords[3], coords[4] ) +
                   tet_volume( coords[7], coords[3], coords[6], coords[4] ) +
                   tet_volume( coords[4], coords[5], coords[1], coords[6] ) +
                   tet_volume( coords[1], coords[6], coords[3], coords[4] ) +
                   tet_volume( coords[2], coords[6], coords[3], coords[1] );
        default:
            return 0.0;
    }
}

/* Calculate the signed volumes beneath the surface (x 6.0). Use the triangle's
 canonical sense. Do not take sense tags into account. Code taken from
 DagMC::measure_volume.

 Special Case: If the surface is planar, and the plane includes the origin,
 the signed volume will be ~0. If the signed volume is ~0 then offset everything
 by a random amount and try again. */
ErrorCode get_signed_volume( Interface* MBI,
                             const EntityHandle surf_set,
                             const CartVect& offset,
                             double& signed_volume )
{
    ErrorCode rval;
    Range tris;
    rval = MBI->get_entities_by_type( surf_set, MBTRI, tris );
    if( MB_SUCCESS != rval ) return rval;
    signed_volume = 0.0;
    const EntityHandle* conn;
    int len;
    CartVect coords[3];
    for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
    {
        rval = MBI->get_connectivity( *j, conn, len, true );
        if( MB_SUCCESS != rval ) return rval;
        if( 3 != len ) return MB_INVALID_SIZE;
        rval = MBI->get_coords( conn, 3, coords[0].array() );
        if( MB_SUCCESS != rval ) return rval;

        // Apply offset to avoid calculating 0 for cases when the origin is in the
        // plane of the surface.
        for( unsigned int k = 0; k < 3; ++k )
        {
            coords[k][0] += offset[0];
            coords[k][1] += offset[1];
            coords[k][2] += offset[2];
        }

        coords[1] -= coords[0];
        coords[2] -= coords[0];
        signed_volume += ( coords[0] % ( coords[1] * coords[2] ) );
    }
    return MB_SUCCESS;
}

// The cgm and cub surfaces may not have the same sense. Create triangles that
// represent the quads in the cub surface. Calculate the signed volume of both
// the cgm and cub surface. If they are different, change the cgm sense so that
// it matches the sense of the cub surface.
ErrorCode fix_surface_senses( Interface* MBI,
                              const EntityHandle cgm_file_set,
                              const EntityHandle cub_file_set,
                              const Tag idTag,
                              const Tag dimTag,
                              const Tag senseTag,
                              const bool debug )
{
    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;
    for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int surf_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
        if( MB_SUCCESS != result ) return result;

        // Find the meshed surface set with the same id
        Range cub_surf;
        const Tag tags[]             = { idTag, dimTag };
        const void* const tag_vals[] = { &surf_id, &two };
        result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_surf );
        if( MB_SUCCESS != result ) return result;
        if( 1 != cub_surf.size() )
        {
            std::cout << "  Surface " << surf_id << ": no meshed representation found, using CAD representation instead"
                      << std::endl;
            continue;
        }

        // Get tris that represent the quads of the cub surf
        Range quads;
        result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
        if( MB_SUCCESS != result ) return result;
        Range cub_tris;
        result = make_tris_from_quads( MBI, quads, cub_tris );
        if( MB_SUCCESS != result ) return result;

        // Add the tris to the same surface meshset as the quads are inside.
        result = MBI->add_entities( cub_surf.front(), cub_tris );
        if( MB_SUCCESS != result ) return result;

        // get the signed volume for each surface representation. Keep trying until
        // The signed volumes are much greater than numerical precision. Planar
        // surfaces will have a signed volume of zero if the plane goes through the
        // origin, unless we apply an offset.
        const int n_attempts        = 100;
        const int max_random        = 10;
        const double min_signed_vol = 0.1;
        double cgm_signed_vol, cub_signed_vol;
        for( int j = 0; j < n_attempts; ++j )
        {
            cgm_signed_vol = 0;
            cub_signed_vol = 0;
            CartVect offset( std::rand() % max_random, std::rand() % max_random, std::rand() % max_random );
            result = get_signed_volume( MBI, *i, offset, cgm_signed_vol );
            if( MB_SUCCESS != result ) return result;
            result = get_signed_volume( MBI, cub_surf.front(), offset, cub_signed_vol );
            if( MB_SUCCESS != result ) return result;
            if( debug )
                std::cout << "  surf_id=" << surf_id << " cgm_signed_vol=" << cgm_signed_vol
                          << " cub_signed_vol=" << cub_signed_vol << std::endl;
            if( fabs( cgm_signed_vol ) > min_signed_vol && fabs( cub_signed_vol ) > min_signed_vol ) break;
            if( n_attempts == j + 1 )
            {
                std::cout << "error: signed volume could not be calculated unambiguously" << std::endl;
                return MB_FAILURE;
            }
        }

        // If the sign is different, reverse the cgm senses so that both
        // representations have the same signed volume.
        if( ( cgm_signed_vol < 0 && cub_signed_vol > 0 ) || ( cgm_signed_vol > 0 && cub_signed_vol < 0 ) )
        {
            EntityHandle cgm_surf_volumes[2], reversed_cgm_surf_volumes[2];
            result = MBI->tag_get_data( senseTag, &( *i ), 1, cgm_surf_volumes );
            if( MB_SUCCESS != result ) return result;
            if( MB_SUCCESS != result ) return result;

            reversed_cgm_surf_volumes[0] = cgm_surf_volumes[1];
            reversed_cgm_surf_volumes[1] = cgm_surf_volumes[0];

            result = MBI->tag_set_data( senseTag, &( *i ), 1, reversed_cgm_surf_volumes );
            if( MB_SUCCESS != result ) return result;
        }
    }

    return MB_SUCCESS;
}

// The quads in the cub_file_set have been updated for dead elements. For each
// cgm_surf, if there exists a cub_surf with the same id, replace the cgm tris
// with cub_tris (created from the quads). Note the a surface that is not
// meshed (in cub file) will not be effected.
ErrorCode replace_faceted_cgm_surfs( Interface* MBI,
                                     const EntityHandle cgm_file_set,
                                     const EntityHandle cub_file_set,
                                     const Tag idTag,
                                     const Tag dimTag,
                                     const bool debug )
{
    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int surf_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
        if( MB_SUCCESS != result ) return result;
        if( debug ) std::cout << "surf_id=" << surf_id << std::endl;

        // Find the meshed surface set with the same id
        Range cub_surf;
        const Tag tags[]             = { idTag, dimTag };
        const void* const tag_vals[] = { &surf_id, &two };
        result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_surf );
        if( MB_SUCCESS != result ) return result;
        if( 1 != cub_surf.size() )
        {
            std::cout << "  Surface " << surf_id << ": no meshed representation found, using CAD representation instead"
                      << std::endl;
            continue;
        }

        // Get tris that represent the quads of the cub surf
        Range quads;
        result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
        if( MB_SUCCESS != result ) return result;

        Range cub_tris;
        result = make_tris_from_quads( MBI, quads, cub_tris );
        if( MB_SUCCESS != result ) return result;

        // Remove the tris from the cgm surf. Don't forget to remove them from the
        // cgm_file_set because it is not TRACKING.
        Range cgm_tris;
        result = MBI->get_entities_by_type( *i, MBTRI, cgm_tris );
        if( MB_SUCCESS != result ) return result;
        result = MBI->remove_entities( *i, cgm_tris );
        if( MB_SUCCESS != result ) return result;
        result = MBI->remove_entities( cgm_file_set, cgm_tris );
        if( MB_SUCCESS != result ) return result;
        result = MBI->delete_entities( cgm_tris );
        if( MB_SUCCESS != result ) return result;

        // Add the cub_tris to the cgm_surf
        result = MBI->add_entities( *i, cub_tris );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

// Dead elements need removed from the simulation. This is done by removing them
// from their volume set and adding them to the implicit complement. New surfaces
// must be created for this.
// IF MODIFYING THIS CODE, BE AWARE THAT DEAD ELEMENTS CAN BE ADJACENT TO MORE
// THAN ONE SURFACE, MAKING THE ASSOCIATION BETWEEN NEWLY EXPOSED AND EXISTING
// SURFACES AMBIGUOUS.
ErrorCode add_dead_elems_to_impl_compl( Interface* MBI,
                                        const EntityHandle cgm_file_set,
                                        const EntityHandle cub_file_set,
                                        const Tag idTag,
                                        const Tag dimTag,
                                        const Tag categoryTag,
                                        const Tag senseTag,
                                        const bool debug )
{

    // Get the cgm surfaces
    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;

    // Get the maximum surface id. This is so that new surfaces do not have
    // duplicate ids.
    int max_surf_id = -1;
    for( Range::const_iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int surf_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
        if( MB_SUCCESS != result ) return result;
        if( max_surf_id < surf_id ) max_surf_id = surf_id;
    }
    std::cout << "  Maximum surface id=" << max_surf_id << std::endl;

    // For each cgm volume, does a cub volume with the same id exist?
    const int three               = 3;
    const void* const three_val[] = { &three };
    Range cgm_vols;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, cgm_vols );
    if( MB_SUCCESS != result ) return result;

    // get the corresponding cub volume
    for( Range::iterator i = cgm_vols.begin(); i != cgm_vols.end(); ++i )
    {
        int vol_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &vol_id );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        std::cout << "  Volume " << vol_id;

        // Find the meshed vol set with the same id
        Range cub_vol;
        const Tag tags[]             = { idTag, dimTag };
        const void* const tag_vals[] = { &vol_id, &three };
        result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_vol );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        if( 1 != cub_vol.size() )
        {
            std::cout << ": no meshed representation found" << std::endl;
            continue;
        }
        else
        {
            std::cout << std::endl;
        }

        // get the mesh elements of the volume.
        Range elems;
        result = MBI->get_entities_by_type( cub_vol.front(), MBHEX, elems );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        if( debug ) std::cout << "    found " << elems.size() << " hex elems" << std::endl;

        // skin the volumes
        Skinner tool( MBI );
        Range skin_faces;
        result = tool.find_skin( 0, elems, 2, skin_faces, true );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;

        // Reconcile the difference between faces of the cub file surfaces and skin
        // faces of the 3D mesh. The difference exists because dead elements have been
        // removed. Faces are divided into:
        //   cub_faces    (in) - the faces in the cub file surface
        //   skin_faces   (in) - the faces of the 3D mesh elements
        //   common_faces (out)- the faces common to both the cub file surface and 3D mesh
        //                       they are still adjacent to this volume
        //   old_faces    (out)- the faces of the cub surface not on the 3D mesh skin
        //                       they are no longer adjacent to this vol
        //   new_faces    (out)- the faces of the 3D mesh skin not on the cub surface
        //                       they are now adjacent to this volume

        // get cub child surfaces.
        Range cub_surfs;
        result = MBI->get_child_meshsets( cub_vol.front(), cub_surfs );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        for( Range::iterator j = cub_surfs.begin(); j != cub_surfs.end(); ++j )
        {
            int surf_id;
            result = MBI->tag_get_data( idTag, &( *j ), 1, &surf_id );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // get the quads on each surface
            Range cub_faces;
            result = MBI->get_entities_by_type( *j, MBQUAD, cub_faces );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // Meshed volumes must have meshed surfaces
            if( cub_faces.empty() )
            {
                std::cout << "    Surface " << surf_id << ": contains no meshed faces" << std::endl;
                // return MB_ENTITY_NOT_FOUND;
            }
            // get the faces common to both the skin and this surface
            Range common_faces = intersect( cub_faces, skin_faces );
            // find the surface faces not on the skin - these are old and need removed
            Range old_faces = subtract( cub_faces, common_faces );
            result          = MBI->remove_entities( *j, old_faces );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;

            // remove the common faces from the skin faces
            skin_faces = subtract( skin_faces, common_faces );
            // If no old faces exist we are done
            if( old_faces.empty() ) continue;
            std::cout << "    Surface " << surf_id << ": " << old_faces.size() << " old faces removed" << std::endl;
            // Place the old faces in a new surface, because they may still be adjacent
            // to 3D mesh in another volume. Get the parent vols of the surface.
            Range cgm_surf;
            const Tag tags2[]             = { idTag, dimTag };
            const void* const tag_vals2[] = { &surf_id, &two };
            result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, tags2, tag_vals2, 2, cgm_surf );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            if( 1 != cgm_surf.size() )
            {
                std::cout << "invalid size" << std::endl;
                return MB_INVALID_SIZE;
            }
            EntityHandle cgm_vols2[2], cub_vols[2];
            result = MBI->tag_get_data( senseTag, cgm_surf, &cgm_vols2 );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            result = MBI->tag_get_data( senseTag, &( *j ), 1, &cub_vols );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // for the new surf, replace the current volume with the impl compl vol.
            // This is because the faces that no longer exist will become adjacent to
            // the impl compl
            if( *i == cgm_vols2[0] )
            {
                cgm_vols2[0] = 0;
                cub_vols[0]  = 0;
            }
            if( *i == cgm_vols2[1] )
            {
                cgm_vols2[1] = 0;
                cub_vols[1]  = 0;
            }
            // If both sides of the surface are the impl comp, do not create the surface.
            if( 0 == cgm_vols2[0] && 0 == cgm_vols2[1] )
            {
                std::cout << "    New surface was not created for old faces because both parents "
                             "are impl_compl volume "
                          << std::endl;
                continue;
            }
            // build the new surface.
            EntityHandle new_cgm_surf, new_cub_surf;
            ++max_surf_id;
            result = build_new_surface( MBI, new_cgm_surf, cgm_vols2[0], cgm_vols2[1], max_surf_id, dimTag, idTag,
                                        categoryTag, senseTag );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            result = build_new_surface( MBI, new_cub_surf, cub_vols[0], cub_vols[1], max_surf_id, dimTag, idTag,
                                        categoryTag, senseTag );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // add the new surface to the file set and populate it with the old faces
            result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            assert( MB_SUCCESS == result );
            result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
            if( MB_SUCCESS != result ) return result;
            assert( MB_SUCCESS == result );
            result = MBI->add_entities( new_cub_surf, old_faces );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;

            std::cout << "    Surface " << max_surf_id << ": created for " << old_faces.size() << " old faces"
                      << std::endl;
        }

        // the remaining skin faces are newly exposed faces
        Range new_faces = skin_faces;

        // new skin faces must be assigned to a surface
        if( new_faces.empty() ) continue;
        std::cout << "    Surface " << max_surf_id + 1 << ": created for " << new_faces.size() << " new faces"
                  << std::endl;

        // Ensure that faces are oriented outwards
        result = orient_faces_outward( MBI, new_faces, debug );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;

        // Create the new surface.
        EntityHandle new_cgm_surf, new_cub_surf;
        ++max_surf_id;
        result = build_new_surface( MBI, new_cgm_surf, *i, 0, max_surf_id, dimTag, idTag, categoryTag, senseTag );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        result = build_new_surface( MBI, new_cub_surf, cub_vol.front(), 0, max_surf_id, dimTag, idTag, categoryTag,
                                    senseTag );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        // Insert the new surf into file sets and populate it with faces.
        result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        result = MBI->add_entities( new_cub_surf, new_faces );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

/* Get the type of a file.
 Return value is one of the above constants
 */
const char* get_geom_file_type( const char* filename );
const char* get_geom_fptr_type( FILE* file );

int is_cubit_file( FILE* file );
int is_step_file( FILE* file );
int is_iges_file( FILE* file );
int is_acis_txt_file( FILE* file );
int is_acis_bin_file( FILE* file );
int is_occ_brep_file( FILE* file );

double DEFAULT_DISTANCE = 0.001;
double DEFAULT_LEN      = 0.0;
int DEFAULT_NORM        = 5;

// load cub file
// load cgm file
// for each surface
//   convert cub surf quads to tris
//   get signed volume from cgm and cub surf         MUST COME BEFORE COORD UPDATE, NEEDS MBTRIS
//   reverse cgm surface sense if needed
//   replace cgm surf tris with cub surf tris
// measure volume of predeformed cub elements
// convert cub volumes sets to tracking so that dead elems are removed from vol sets
// update coordinates and delete dead elems
// measure volume of deformed cub elems
// print histogram of volume change
// for each cub volume
//   skin volume elems to get faces
//   for each child cub surface
//     assign old skin faces to a new surface in case they are adjacent to another volume
//   orient each skin face outward
//   assign new skin faces to a new surface
// for each surface
//   remove existing tris (from before the update)
//   convert quads to tris
// remove empty surfaces and volumes due to dead elements
int main( int argc, char* argv[] )
{
    clock_t start_time = clock();

    const bool debug      = false;
    const char* file_type = NULL;

    const char* cub_name  = 0;
    const char* exo_name  = 0;
    const char* out_name  = 0;
    const char* time_step = 0;
    const char* sat_name  = 0;
    double dist_tol = 0.001, len_tol = 0.0;
    int norm_tol = 5;

    if( 6 != argc && 9 != argc )
    {
        std::cerr << "To read meshed geometry for MOAB:" << std::endl;
        std::cerr << "$> <cub_file.cub> <acis_file.sat> <facet_tol> <output_file.h5m> conserve_mass<bool>" << std::endl;
        std::cerr << "To read meshed geometry for MOAB and update node coordinates:" << std::endl;
        std::cerr << "$> <cub_file.cub> <acis_file.sat> <facet_tol> <output_file.h5m> "
                     "<deformed_exo_file.e> time_step<int> check_vol_change<bool> conserve_mass<bool>"
                  << std::endl;
        exit( 4 );
    }

    // check filenames for proper suffix
    std::string temp;
    cub_name = argv[1];
    temp.assign( cub_name );
    if( std::string::npos == temp.find( ".cub" ) )
    {
        std::cerr << "cub_file does not have correct suffix" << std::endl;
        return 1;
    }
    sat_name = argv[2];  // needed because the cub file's embedded sat file does not have groups
    temp.assign( sat_name );
    if( std::string::npos == temp.find( ".sat" ) )
    {
        std::cerr << "sat_file does not have correct suffix" << std::endl;
        return 1;
    }
    out_name = argv[4];
    temp.assign( out_name );
    if( std::string::npos == temp.find( ".h5m" ) )
    {
        std::cerr << "out_file does not have correct suffix" << std::endl;
        return 1;
    }

    // get the facet tolerance
    dist_tol = atof( argv[3] );
    if( 0 > dist_tol || 1 < dist_tol )
    {
        std::cout << "error: facet_tolerance=" << dist_tol << std::endl;
        return 1;
    }

    // Should the nodes be updated?
    bool update_coords = false;
    if( 9 == argc )
    {
        exo_name = argv[5];
        temp.assign( exo_name );
        if( std::string::npos == temp.find( ".e" ) )
        {
            std::cerr << "e_file does not have correct suffix" << std::endl;
            return 1;
        }
        time_step     = argv[6];
        update_coords = true;
    }

    // Should the volume change be determined?
    bool determine_volume_change = false;
    if( 9 == argc )
    {
        temp.assign( argv[7] );
        if( std::string::npos != temp.find( "true" ) ) determine_volume_change = true;
    }

    // Should densities be changed to conserve mass?
    bool conserve_mass = false;
    if( 9 == argc )
    {
        temp.assign( argv[8] );
        if( std::string::npos != temp.find( "true" ) ) conserve_mass = true;
    }
    else if( 6 == argc )
    {
        temp.assign( argv[5] );
        if( std::string::npos != temp.find( "true" ) ) conserve_mass = true;
    }

    // Get CGM file type
    if( !file_type )
    {
        file_type = get_geom_file_type( cub_name );
        if( !file_type )
        {
            std::cerr << cub_name << " : unknown file type, try '-t'" << std::endl;
            exit( 1 );
        }
    }

    // Read the mesh from the cub file with Tqcdfr
    Core* MBI = new Core();
    ErrorCode result;
    EntityHandle cub_file_set;
    result = MBI->create_meshset( 0, cub_file_set );
    if( MB_SUCCESS != result ) return result;
    // Do not ignore the Cubit file version. In testing, a cub file from Cubit12
    // did not work.
    // char cub_options[256] = "120";
    // char cub_options[256] = "IGNORE_VERSION";
    // result = MBI->load_file(cub_name, &cub_file_set, cub_options, NULL, 0, 0);
    result = MBI->load_file( cub_name, &cub_file_set, 0, NULL, 0, 0 );
    if( MB_SUCCESS != result )
    {
        std::cout << "error: problem reading cub file" << std::endl;
        return result;
    }
    std::cout << "Mesh file read." << std::endl;

    // Read the ACIS file with ReadCGM
    char cgm_options[256];
    std::cout << "  facet tolerance=" << dist_tol << std::endl;
    sprintf( cgm_options,
             "CGM_ATTRIBS=yes;FACET_DISTANCE_TOLERANCE=%g;FACET_NORMAL_TOLERANCE=%d;MAX_FACET_EDGE_"
             "LENGTH=%g;",
             dist_tol, norm_tol, len_tol );
    EntityHandle cgm_file_set;
    result = MBI->create_meshset( 0, cgm_file_set );
    if( MB_SUCCESS != result ) return result;
    result = MBI->load_file( sat_name, &cgm_file_set, cgm_options, NULL, 0, 0 );
    if( MB_SUCCESS != result )
    {
        std::cout << "error: problem reading sat file" << std::endl;
        return result;
    }
    std::cout << "CAD file read." << std::endl;

    // Create tags
    Tag dimTag, idTag, categoryTag, senseTag, sizeTag, nameTag;
    result = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dimTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    idTag  = MBI->globalId_tag();
    result = MBI->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
                                  MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, senseTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_get_handle( "GEOM_SIZE", 1, MB_TYPE_DOUBLE, sizeTag, MB_TAG_DENSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, nameTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;

    // Create triangles from the quads of the cub surface sets and add them to the
    // cub surface sets. Get the signed volume of each surface for both cgm and
    // cub representations. Change the sense of the cgm representation to match
    // the cub representation.
    result = fix_surface_senses( MBI, cgm_file_set, cub_file_set, idTag, dimTag, senseTag, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Fixed CAD surface senses to match meshed surface senses." << std::endl;

    // Get the 3D elements in the cub file and measure their volume.
    Range orig_elems;
    std::vector< double > orig_size;
    if( determine_volume_change )
    {
        result = MBI->get_entities_by_dimension( 0, 3, orig_elems );
        if( MB_SUCCESS != result ) return result;
        orig_size.resize( orig_elems.size() );
        for( unsigned int i = 0; i < orig_elems.size(); ++i )
        {
            orig_size[i] = measure( MBI, orig_elems[i] );
        }
    }

    // Before updating the nodes and removing dead elements, force the cub volume
    // sets to track ownership so that dead elements will be deleted from the sets.
    const int three               = 3;
    const void* const three_val[] = { &three };
    Range cub_vols;
    result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, &dimTag, three_val, 1, cub_vols );
    if( MB_SUCCESS != result ) return result;
    for( Range::const_iterator i = cub_vols.begin(); i != cub_vols.end(); ++i )
    {
        result = MBI->set_meshset_options( *i, MESHSET_TRACK_OWNER );
        if( MB_SUCCESS != result ) return result;
    }

    // Tag volume sets with the undeformed size of each volume.
    Range vol_sets;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, vol_sets );
    if( MB_SUCCESS != result ) return result;

    moab::GeomTopoTool gtt  = moab::GeomTopoTool( MBI, false );
    moab::GeomQueryTool gqt = moab::GeomQueryTool( &gtt );

    for( Range::const_iterator i = vol_sets.begin(); i != vol_sets.end(); ++i )
    {
        double size;
        result = gqt.measure_volume( *i, size );
        if( MB_SUCCESS != result ) return result;
        result = MBI->tag_set_data( sizeTag, &( *i ), 1, &size );
        if( MB_SUCCESS != result ) return result;
    }

    // Update the coordinates if needed. Do not do this before checking surface
    // sense, because the coordinate update could deform the surfaces too much
    // to make an accurate comparison.
    // The cub node ids are unique because cgm vertex ids are tagged on the vertex
    // meshset and not the vertex itself.
    // result = MBI->delete_entities( &cub_file_set, 1 );
    // if(MB_SUCCESS != result) return result;
    // Assume dead elements exist until I think of something better.
    bool dead_elements_exist = true;
    if( update_coords )
    {
        // ReadNCDF my_ex_reader(MBI);
        char exo_options[120] = "tdata=coord,";
        strcat( exo_options, time_step );
        strcat( exo_options, ",set" );
        // FileOptions exo_opts(exo_options)  ;
        // opts = "tdata=coord, 100, sum, temp.exo";
        // result =  my_ex_reader.load_file(exo_name, cgm_file_set, exo_opts, NULL, 0 , 0);
        // result =  my_ex_reader.load_file(exo_name, cub_file_set, exo_opts, NULL, 0 , 0);
        // result = my_ex_reader.load_file(exo_name, &cub_file_set, exo_opts, NULL, 0 , 0);
        MBI->load_file( exo_name, &cub_file_set, exo_options );
        if( MB_SUCCESS != result )
        {
            std::string last_error;
            MBI->get_last_error( last_error );
            std::cout << "coordinate update failed, " << last_error << std::endl;
            return result;
        }
        std::cout << "Updated mesh nodes with deformed coordinates from exodus file." << std::endl;
    }

    if( determine_volume_change )
    {
        // Dead elements have been removed by the deformation. Get the elements that
        // still exist.
        Range defo_elems;
        result = MBI->get_entities_by_dimension( 0, 3, defo_elems );
        if( MB_SUCCESS != result ) return result;

        // Determine the volume of the elements now that a deformation has been
        // applied. Condense the original array by removing dead elements.
        double* orig_size_condensed = new double[defo_elems.size()];
        double* defo_size_condensed = new double[defo_elems.size()];
        int j                       = 0;
        for( unsigned int i = 0; i < orig_elems.size(); ++i )
        {
            if( orig_elems[i] == defo_elems[j] )
            {
                orig_size_condensed[j] = orig_size[i];
                defo_size_condensed[j] = measure( MBI, defo_elems[j] );
                ++j;
            }
        }
        generate_plots( orig_size_condensed, defo_size_condensed, defo_elems.size(), std::string( time_step ) );
        delete[] orig_size_condensed;  // can't use the same delete[] for both
        delete[] defo_size_condensed;
    }

    // Deal with dead elements. For now, add them to the impl_compl volume.
    // Extra surfaces are created to do this.
    if( update_coords && dead_elements_exist )
    {
        result = add_dead_elems_to_impl_compl( MBI, cgm_file_set, cub_file_set, idTag, dimTag, categoryTag, senseTag,
                                               debug );
        if( MB_SUCCESS != result ) return result;
        std::cout << "Placed dead elements in implicit complement volume and added required surfaces." << std::endl;
    }

    // The quads in the cub_file_set have been updated for dead elements. For each
    // cgm_surf, if there exists a cub_surf with the same id, replace the cgm tris
    // with cub_tris (created from the quads). Note the a surface that is not
    // meshed (in cub file) will not be effected.
    result = replace_faceted_cgm_surfs( MBI, cgm_file_set, cub_file_set, idTag, dimTag, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Replaced faceted CAD surfaces with meshed surfaces of triangles." << std::endl;

    result = remove_empty_cgm_surfs_and_vols( MBI, cgm_file_set, idTag, dimTag, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Removed surfaces and volumes that no longer have any triangles." << std::endl;

    // For each material, sum the volume. If the coordinates were updated for
    // deformation, summarize the volume change.
    result = summarize_cell_volume_change( MBI, cgm_file_set, categoryTag, dimTag, sizeTag, nameTag, idTag,
                                           conserve_mass, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Summarized the volume change of each material, with respect to the solid model." << std::endl;

    result = MBI->write_mesh( out_name, &cgm_file_set, 1 );
    if( MB_SUCCESS != result )
    {
        std::cout << "write mesh failed" << std::endl;
        return result;
    }
    std::cout << "Saved output file for mesh-based analysis." << std::endl;

    clock_t end_time = clock();
    std::cout << "  " << (double)( end_time - start_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
    std::cout << std::endl;

    return 0;
}

const char* get_geom_file_type( const char* name )
{
    FILE* file;
    const char* result = 0;

    file = fopen( name, "r" );
    if( file )
    {
        result = get_geom_fptr_type( file );
        fclose( file );
    }

    return result;
}

const char* get_geom_fptr_type( FILE* file )
{
    static const char* CUBIT_NAME = GF_CUBIT_FILE_TYPE;
    static const char* STEP_NAME  = GF_STEP_FILE_TYPE;
    static const char* IGES_NAME  = GF_IGES_FILE_TYPE;
    static const char* SAT_NAME   = GF_ACIS_TXT_FILE_TYPE;
    static const char* SAB_NAME   = GF_ACIS_BIN_FILE_TYPE;
    static const char* BREP_NAME  = GF_OCC_BREP_FILE_TYPE;

    if( is_cubit_file( file ) )
        return CUBIT_NAME;
    else if( is_step_file( file ) )
        return STEP_NAME;
    else if( is_iges_file( file ) )
        return IGES_NAME;
    else if( is_acis_bin_file( file ) )
        return SAB_NAME;
    else if( is_acis_txt_file( file ) )
        return SAT_NAME;
    else if( is_occ_brep_file( file ) )
        return BREP_NAME;
    else
        return 0;
}

int is_cubit_file( FILE* file )
{
    unsigned char buffer[4];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 4, 1, file ) && !memcmp( buffer, "CUBE", 4 );
}

int is_step_file( FILE* file )
{
    unsigned char buffer[9];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 9, 1, file ) && !memcmp( buffer, "ISO-10303", 9 );
}

int is_iges_file( FILE* file )
{
    unsigned char buffer[10];
    return !fseek( file, 72, SEEK_SET ) && fread( buffer, 10, 1, file ) && !memcmp( buffer, "S      1\r\n", 10 );
}

int is_acis_bin_file( FILE* file )
{
    char buffer[15];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 15, 1, file ) && !memcmp( buffer, "ACIS BinaryFile", 9 );
}

int is_acis_txt_file( FILE* file )
{
    char buffer[5];
    int version, length;

    if( fseek( file, 0, SEEK_SET ) || 2 != fscanf( file, "%d %*d %*d %*d %d ", &version, &length ) ) return 0;

    if( version < 1 || version > 0xFFFF ) return 0;

    // Skip application name
    if( fseek( file, length, SEEK_CUR ) ) return 0;

    // Read length of version string followed by first 5 characters
    if( 2 != fscanf( file, "%d %4s", &length, buffer ) ) return 0;

    return !strcmp( buffer, "ACIS" );
}

int is_occ_brep_file( FILE* file )
{
    unsigned char buffer[6];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 6, 1, file ) && !memcmp( buffer, "DBRep_", 6 );
}